3D Triangle with 2 Points, All Lengths, and Z of the 3rd Point Given

Question

I have a 3D triangle ABC with known Lengths AB, BC, and AC and known points A(x,y,z), B(x,y,z) and C(z). I need to be able to find the x and y coordinates of C using a non-iterative method.

I've been going through another post and was following along until the aa bb and c.z calulations in the answer. Unfortunately, I don't have enough reputation to comment on that post...

3d geometry: triangle 2 points known, find 3rd point

Any help is greatly appreciated!

I was just trying to show I only know the z coordinate of point C. I probably should have said I need to find point C(?,?,z) where z is given. — Lance Miller, Jun 28 '17 at 20:19
The side lengths and known vertices give you two spheres on which $C$ might lie. These spheres intersect the plane on which $C$ lies (if at all) in a pair of circles. Finding equations for these circles is a matter of plugging $C$’s known $z$-coordinate into the equations of these spheres, and then you’re reduced to computing the intersection of two circles in 2-D. — amd, Jun 28 '17 at 22:23

amd · Accepted Answer · 2017-06-30T21:03:37.050

The point $C$ lies somewhere on the intersection of a pair of spheres with centers $A$ and $B$. Since we also know that $C$ lies somewhere on the plane $z=z_C$, this reduces to finding the intersection of a pair of circles.

Let $r_A=|AC|$ and $r_B=|BC|$. The equations of the two spheres are $$(x-x_A)^2+(y-y_A)^2+(z-z_A)^2=r_A^2 \\ (x-x_B)^2+(y-y_B)^2+(z-z_B)^2=r_B^2.$$ If $|z_C-z_A|\gt r_A$ or $|z_C-z_B|\gt r_B$, then there’s no solution. Otherwise substitute $z=z_C$ into the above equations to get the equations of two circles: $$(x-x_A)^2+(y-y_A)^2=r_A^2-(z_C-z_A)^2 \\ (x-x_B)^2+(y-y_B)^2=r_B^2-(z_C-z_B)^2.$$ If you subtract one of these equations from the other you get a linear equation in $x$ and $y$, which you can solve for one or the other and then back-substitute into either circle equation to get a solution. The most common case will be two solutions, but if the vertices are lined up just right you might get one or an infinite number of them (the two circles coincide).

Alternatively, a direct calculation with vectors produces the points of intersection of the two circles. Let $r_A'^2=r_A^2-(z_C-z_A)^2$ and $r_B'^2=r_B^2-(z_C-z_B)^2$ be the radii of the two circles, and $d=\sqrt{(x_A-x_B)^2+(y_A-y_B)^2}$ the distance between their centers. If $r_A'+r_B'>d$, the circles clearly don’t intersect. If $d=0$, the circles are concentric, for which see below. Otherwise, the intersections of the circles are equally spaced on a line perpendicular to the line through the two centers.

Referring to the illustration above, by the Pythagorean theorem the distance of this line from the center of circle $A$ satisfies $$r_A'^2-w^2=r_B'^2-(d-w)^2$$ from which $$w={d^2+r_A'^2-r_B'^2\over2d}.$$ Again by the Pythagorean theorem, the distance $h$ along the perpendicular line to the intersection points is $h=\sqrt{r_A'^2-w^2}=\sqrt{r_B'^2-(d-w)^2}$. Taking unit vectors in the directions of the two lines, this gives for the $x$- and $y$- coordinates of the intersection points $$\frac wd(x_B-x_A,y_B-y_A)\pm\frac hd(y_A-y_B,x_B-x_A).\tag{*}$$ If $d=0$, then the $x$- and $y$- coordinates of points $A$ and $B$ are equal, so there’s one solution—$(x_A,y_A,z_C)$—when $r_A'=r_B'=0$ (the two spheres are tangent) or an infinite number when $r_A'=r_B'\gt0$. When $r_A'\ne r_B'$, there is, of course, no solution.

As an implementation note, observe that $\frac wd={d^2+r_A'^2-r_B'^2\over d^2}$ and similarly $h/d$ only involves $d^2$ after factoring in the contribution of $w^2$. All of the radii are squared in the actual computations as well, so in practice the only point at which you really need to take a square root is when computing $h$. Range checks can be delayed to avoid taking square roots: If $z_C$ is out of range, one or both of $r_A'^2$ or $r_B'^2$ will end up being negative; if $r_A'+r_B'\gt d$, then $h$ will end up being imaginary.

I originally considered this approach as well, but the problem is, the unit vector along $h$ is unknown, except for being perpendicular to $\overline{AB}$. This means the actual expressions to find the $x$ and $y$ coordinates for the point become as difficult to implement as before. Unless, of course, you have some nice, clever way to solve the actual coordinates for the point, perhaps? (Note: I am in no way disputing your answer, in either current or original form. You are quite correct. I am only saying that us programmers need more advice/help for implementing this in practice.) — Nominal Animal, Jun 30 '17 at 20:37
@NominalAnimal Errr... have another look at equation (*), which has that perpendicular vector in its second term. That doesn’t look very difficult to implement to me. Everything happens in the $z=z_C$ plane, making the problem two-dimensional. We don’t need a perpendicular to $\overline{AB}$ per se, but do want a perpendicular to the line through their projections onto that plane, i.e., a vector with $z=0$. Once you’ve rejected the no-solution cases and computed $x$ and $y$ as above, slap $z_C$ on the end to bring it back into 3-D. — amd, Jun 30 '17 at 21:05
!!! I think I now understand. Your second pair of equations can be treated as the 2D solution in the $z = z_C$ plane, the right sides corresponding to the squared radii of the two circles corresponding to the 3D spheres. Wow; why didn't I see that at all? I think this is a brilliant approach; let me go through the steps to see how practical code would look like. — Nominal Animal, Jun 30 '17 at 21:12
Oh yes; this is much simpler, code easier to implement and understand, and faster too. (I still want to show other programmers how to implement this, though, so I think I shall rewrite my answer.) @lance-miller, could you please accept this answer instead? — Nominal Animal, Jun 30 '17 at 21:43
@amd Thank you so much for your help! It's been so long since I've done anything like this that I didn't know how to work through your first response. — Lance Miller, Jul 01 '17 at 16:19

Nominal Animal · Answer 2 · 2017-06-30T22:05:52.703

amd's answer describes a much better method than what I suggested in this answer.

Simply put, amd pointed out that the problem simplifies to a 2D one. In 2D, the distance between $A$ and $C$ is $\sqrt{AC^2 - (z_C - z_A)^2}$, and the distance between $B$ and $C$ is $\sqrt{BC^2 - (z_C - z_B)^2}$.

amd's answer is very simple to implement. In pseudocode, it is

Function Find_third_point(A_x, A_y, A_z,
                          B_x, B_y, B_z,
                          AC, BC, C_z, epsilon):
    # 2D coordinates
    Bx = B_x - A_x
    By = B_y - A_y
    Ad2 = AC*AC - (C_z - A_z)*(C_z - A_z)
    Bd2 = BC*BC - (C_z - B_z)*(C_z - B_z)
    d2 = Bx*Bx + By*By
    If (d2 <= epsilon):
        Fail: "No solutions, or an infinite number of solutions"
    End if

    w_per_d = (d2 + Ar2 - Br2) / (d2 + d2)
    w2_per_d2 = w_per_d*w_per_d

    If (Ar2 < d2*w2_per_d2):
        Fail: "No solutions"
    End if
    h_per_d = sqrt(Ar2/d2 - w2_per_d2)

    Return Point(A_x + w_per_d*Bx - h_per_d*By,
                 A_y + w_per_d*By + h_per_d*Bx,
                 C_z),
           Point(A_x + w_per_d*Bx + h_per_d*By,
                 A_y + w_per_d*By - h_per_d*Bx,
                 C_z)

which only needs one square root, sixteen multiplications, two divisions, and 21 additions or subtractions; a brilliant solution.

Here is an awk test script that tests a million random points (or NUMBER if you add -v n=NUMBER to the command line):

#!/usr/bin/awk -f
function distance(xp1, yp1, zp1, xp2, yp2, zp2) { return sqrt((xp2-xp1)*(xp2-xp1) + (yp2-yp1)*(yp2-yp1) + (zp2-zp1)*(zp2-zp1)) }
BEGIN {
    srand();

    eps = 0.0005
    xmin = -9.999 ; ymin = -9.999 ; zmin = -9.999
    xmax = +9.999 ; ymax = +9.999 ; zmax = +9.999

    n = int(n)
    if (n < 1)
        n = 1000000;

    for (i = 0; i < n; i++) {

        if (!(i % 10000))
            printf "\r%.1f%% tests completed\033[K\r", 100.0*i/n

        while (1) {
            x1 = xmin + (xmax - xmin) * rand()
            y1 = ymin + (ymax - ymin) * rand()
            z1 = zmin + (zmax - zmin) * rand()
            x2 = xmin + (xmax - xmin) * rand()
            y2 = ymin + (ymax - ymin) * rand()
            z2 = zmin + (zmax - zmin) * rand()
            x3 = xmin + (xmax - xmin) * rand()
            y3 = ymin + (ymax - ymin) * rand()
            z3 = zmin + (zmax - zmin) * rand()
            dd12 = (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1)
            dd13 = (x3-x1)*(x3-x1) + (y3-y1)*(y3-y1) + (z3-z1)*(z3-z1)
            dd23 = (x3-x2)*(x3-x2) + (y3-y2)*(y3-y2) + (z3-z2)*(z3-z2)
            if (dd12 > eps*eps &&
                dd13 > eps*eps &&
                dd23 > eps*eps) break
        }

        # Preparation

        Bx = x2 - x1
        By = y2 - y1

        Cz = z3 - z1

        Ar2 = dd13 - (z3 - z1)*(z3 - z1)
        Br2 = dd23 - (z3 - z2)*(z3 - z2)

        d2 = Bx*Bx + By*By
        if (d2 <= eps*eps) {
            printf "A:   %7.3f %7.3f %7.3f\n", x1, y1, z1
            printf "B:   %7.3f %7.3f %7.3f\n", x2, y2, z2
            printf "C:   %7.3f %7.3f %7.3f\n", x3, y3, z3
            printf "AB: %11.6f\n", sqrt(dd12)
            printf "AC: %11.6f\n", sqrt(dd13)
            printf "BC: %11.6f\n", sqrt(dd23)
            printf "\n"
            printf "Cannot solve, A and B differ only by their z coordinates.\n"
            exit(1)
        }

        w_per_d = (d2 + Ar2 - Br2) / (d2 + d2)
        w2_per_d2 = w_per_d * w_per_d
        if (Ar2 < d2*w2_per_d2) {
            printf "A:   %7.3f %7.3f %7.3f\n", x1, y1, z1
            printf "B:   %7.3f %7.3f %7.3f\n", x2, y2, z2
            printf "C:   %7.3f %7.3f %7.3f\n", x3, y3, z3
            printf "AB: %11.6f\n", sqrt(dd12)
            printf "AC: %11.6f\n", sqrt(dd13)
            printf "BC: %11.6f\n", sqrt(dd23)
            printf "\n"
            printf "No solutions.\n"
            exit(1)
        }
        h_per_d = sqrt(Ar2/d2 - w2_per_d2)

        Cx1 = x1 + w_per_d * Bx - h_per_d * By
        Cy1 = y1 + w_per_d * By + h_per_d * Bx
        Cz1 = z3

        Cx2 = x1 + w_per_d * Bx + h_per_d * By
        Cy2 = y1 + w_per_d * By - h_per_d * Bx
        Cz2 = z3

        err11 = distance(Cx1,Cy1,Cz1,x1,y1,z1) - sqrt(dd13)
        err12 = distance(Cx1,Cy1,Cz1,x2,y2,z2) - sqrt(dd23)
        err21 = distance(Cx2,Cy2,Cz2,x1,y1,z1) - sqrt(dd13)
        err22 = distance(Cx2,Cy2,Cz2,x2,y2,z2) - sqrt(dd23)

        if (err11 > eps || err11 < -eps ||
            err12 > eps || err12 < -eps ||
            err21 > eps || err21 < -eps ||
            err22 > eps || err22 < -eps) {
            printf "A:   %7.3f %7.3f %7.3f\n", x1, y1, z1
            printf "B:   %7.3f %7.3f %7.3f\n", x2, y2, z2
            printf "C:   %7.3f %7.3f %7.3f\n", x3, y3, z3
            printf "AB: %11.6f\n", sqrt(dd12)
            printf "AC: %11.6f\n", sqrt(dd13)
            printf "BC: %11.6f\n", sqrt(dd23)
            printf "\n"
            printf "Found: %7.3f %7.3f %7.3f errors: %+.6f %+.6f %+.6f\n", Cx1, Cy1, Cz1, distance(Cx1,Cy1,Cz1,x3,y3,z3), distance(Cx1,Cy1,Cz1,x1,y1,z1) - sqrt(dd13), distance(Cx1,Cy1,Cz1,x2,y2,z2) - sqrt(dd23)
            printf "Found: %7.3f %7.3f %7.3f errors: %+.6f %+.6f %+.6f\n", Cx2, Cy2, Cz2, distance(Cx2,Cy2,Cz2,x3,y3,z3), distance(Cx2,Cy2,Cz2,x1,y1,z1) - sqrt(dd13), distance(Cx2,Cy2,Cz2,x2,y2,z2) - sqrt(dd23)
            exit(1)
        }
    }

    printf "All %d tests successful.\n", n
}

I am honestly ashamed that I didn't see the 2D simplification at all; it took quite a bit of prodding from amd for me to get it.

In the interests of completeness, below is my earlier, comparatively stupid, suggestion, that originally got accepted by the asker. (This is to show that if somebody with good math-fu shows you another way, it always pays to listen and understand: algorithmic enhancements always beat code optimizations.)

If point $A$ is at the origin, then we have $$\begin{cases} C_x^2 + C_y^2 + C_z^2 = AC^2 \\ (C_x - B_x)^2 + (C_y - B_y)^2 + (C_z - B_z)^2 = BC^2 \end{cases} \tag{1a}\label{1a}$$ i.e. $$\begin{cases} C_x^2 + C_y^2 = AC^2 - C_z^2 \\ C_x^2 - 2 B_x C_x + C_y^2 - 2 B_y C_y = BC^2 - (C_z - B_z)^2 - B_x^2 - B_y^2\end{cases}\tag{1b}\label{1b}$$ that we wish to solve for $C_x$ and $C_y$.

Solving $\eqref{1a}$ or $\eqref{1b}$ for $C_x$ and $C_y$ yields two solutions. I used Maple to look at what they are:

C_x = A_x+((B_x-A_x)*(AC^2-BC^2+(B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2-2*(C_z-A_z)*(B_z-A_z))-signum(-B_x+A_x)*signum(-B_y+A_y)*sqrt((B_y-A_y)^2*((B_x-A_x)^4+2*((B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2)*(AC^2+BC^2-(B_x-A_x)^2)-(AC^2-BC^2)^2-((B_y-A_y)^2+(B_z-A_z)^2)^2+4*(B_z-A_z)*(C_z-A_z)*(AC^2-BC^2+(B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2)-4*AC^2*(B_z-A_z)^2-4*((B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2)*(C_z-A_z)^2)))/(2*(B_x-A_x)^2+2*(B_y-A_y)^2)
C_y = A_y+((B_y-A_y)*(AC^2-BC^2+(B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2-2*(C_z-A_z)*(B_z-A_z))+sqrt((B_x-A_x)^2*((B_x-A_x)^4+2*((B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2)*(AC^2+BC^2-(B_x-A_x)^2)-(AC^2-BC^2)^2-((B_y-A_y)^2+(B_z-A_z)^2)^2+4*(B_z-A_z)*(C_z-A_z)*(AC^2-BC^2+(B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2)-4*AC^2*(B_z-A_z)^2-4*((B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2)*(C_z-A_z)^2)))/(2*(B_x-A_x)^2+2*(B_y-A_y)^2)

and

C_x = A_x+((B_x-A_x)*(AC^2-BC^2+(B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2-2*(C_z-A_z)*(B_z-A_z))+signum(-B_x+A_x)*signum(-B_y+A_y)*sqrt((B_y-A_y)^2*((B_x-A_x)^4+2*((B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2)*(AC^2+BC^2-(B_x-A_x)^2)-(AC^2-BC^2)^2-((B_y-A_y)^2+(B_z-A_z)^2)^2+4*(B_z-A_z)*(C_z-A_z)*(AC^2-BC^2+(B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2)-4*AC^2*(B_z-A_z)^2-4*((B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2)*(C_z-A_z)^2)))/(2*(B_x-A_x)^2+2*(B_y-A_y)^2)
C_y = A_y+((B_y-A_y)*(AC^2-BC^2+(B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2-2*(C_z-A_z)*(B_z-A_z))-sqrt((B_x-A_x)^2*((B_x-A_x)^4+2*((B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2)*(AC^2+BC^2-(B_x-A_x)^2)-(AC^2-BC^2)^2-((B_y-A_y)^2+(B_z-A_z)^2)^2+4*(B_z-A_z)*(C_z-A_z)*(AC^2-BC^2+(B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2)-4*AC^2*(B_z-A_z)^2-4*((B_x-A_x)^2+(B_y-A_y)^2+(B_z-A_z)^2)*(C_z-A_z)^2)))/(2*(B_x-A_x)^2+2*(B_y-A_y)^2)

where ^ denotes exponentiation, not binary exclusive-or operation; and signum() is zero if its argument is zero, -1 if its argument is negative, and +1 if its argument is positive.

Since the above forms are horrible, I fiddled with the expressions a bit to find a form a programmer might use. Mostly, I did this to show that even if the solution looks horrible, it may still be quite easy to use, if you first understand the solution and its requirements. (You do need to avoid things like division by zero, and verify your solutions numerically.)

Let's define

ABAB = AB * AB    # Squared distance from A to B
ACAC = AC * AC    # Squared distance from A to C
BCBC = BC * BC    # Squared distance from B to C
Bx = B_x - A_x    # Location of B relative to A
By = B_y - A_y
Bz = B_z - A_z
BxBx = Bx * Bx
ByBy = By * By
BzBz = Bz * Bz
Cz = C_z - A_z    # Z coordinate of C relative to A
#                   Note: ABAB = BxBx + ByBy + BzBz

so that Bx, By, Bz, and Cz are relative to the location of point $A$; our solutions for Cx and Cy will similarly be relative to $A$, so the actual solutions where $A$ is not at origin will be

C_x = A_x + Cx
C_y = A_y + Cy

Note that I use an _ only in the original coordinate variable names (where point $A$ is not at origin). Coordinate variables without a _ are relative to the location of point $A$.

If one of AB, AC, and AB are all zero, we have a degenerate triangle, a point. Then, C_x = A_x and C_y = A_y.

If AB = 0, but AC ≠ BC, there is no solution.

If AB = 0, and AC = BC > 0, any point $C$ in a circle of radius $AC$ around $A$ works.

If BC = 0, but AC ≠ AB or B_z ≠ C_z, there is no solution.

If BC = 0, AC = AB > 0, and B_z = C_z, then C_x = B_x and C_y = B_y.

If Bx = 0 and By = 0, but AC > 0 and BC > 0, there may be no solution (for example if C_z > AC + A_z). If there are any solutions, they form a circle around (A_x, A_y, C_z), perpendicular to the $z$ axis.

If Bx = 0, the solution simplifies to

    Cy = (ACAC + ABAB - BCBC - 2*Bz*Cz) / (2*By)
    Cx = +- sqrt(ACAC - Cy*Cy - Cz*Cz)

but note that if ACAC = Cy*Cy + Cz*Cz, there is just one solution, Cx = 0, and if ACAC < Cy*Cy + Cz*Cz, there are no (real) solutions.

If By = 0, the solution simplifies to

    Cx = (ACAC + ABAB - BCBC - 2*Bz*Cz) / (2*Bx)
    Cy = +- sqrt(ACAC - Cx*Cx - Cz*Cz)

but note that if ACAC = Cy*Cy + Cz*Cz, there is just one solution, Cy = 0, and if ACAC < Cx*Cx + Cz*Cz, there are no (real) solutions.

In the general case, where neither Bx or By is zero, calculate S:

S = BxBx * BxBx
  - (ACAC - BCBC) * (ACAC - BCBC)
  - (ByBy + BzBz) * (ByBy + BzBz)
  + 2 * ABAB * (ACAC + BCBC - BxBx)
  + 4*( Bz * Cz * (ABAB + ACAC - BCBC) - ACAC*BzBz - ABAB*Cz*Cz )

If S is negative, there are no real solutions. If S is zero, there is exactly one solution (the two solutions are the same). If S is positive, there are two solutions.

(Note that S*BxBx corresponds to the value inside the square root, sqrt(), in the horrible formulas above. I took the term outside, so that we do not try to calculate By/Bx outside the square root.)

Next, calculate temporary values

BB = 2 * (BxBx + ByBy)
BL = ACAC + ABAB - BCBC - 2 * Bz * Cz

If BB is zero, then points $A$ and $B$ have the same $x$ and $y$ coordinates, and there is no solution.

We also need a sign constant:

BS = +1  if Bx and By have the same sign
   = -1  if Bx and By have different signs

Note that if you use floating-point numbers, you do not need to inspect the actual signs; you can just use e.g.

If Bx*By > 0:
    BS = +1
Else:
    Bs = -1
End If

If we define BS so that it is only zero if both are zero, and otherwise only considers the sign of the nonzero value(s), this approach works for the (Bx=0 or By=0) cases as well. That is,

If (Bx > 0):
    If (By > 0):
        BS = +1
    Else if (By < 0):
        BS = -1
    Else:
        BS = +1
    End if
Else if (Bx < 0):
    If (By > 0):
        BS = -1
    Else if (By < 0):
        BS = +1
    Else:
        BS = -1
    End if
Else:
    If (By > 0):
        BS = +1
    Else if (By < 0):
        BS = -1
    Else:
        BS = 0
    End if
End if

The first solution, relative to point $A$, is

Cx1 = (Bx * BL - BS * sqrt( ByBy * S )) / BB
Cy1 = (By * BL + sqrt( BxBx * S )) / BB

and the second solution is

Cx2 = (Bx * BL + BS * sqrt( ByBy * S )) / BB
Cy2 = (By * BL - sqrt( BxBx * S )) / BB

To repeat, if $A$ is not at the origin, the solutions are

C_x = A_x + (Bx * BL - BS * sqrt( ByBy * S )) / BB
C_y = A_y + (By * BL + sqrt( BxBx * S )) / BB

and

C_x = A_x + (Bx * BL + BS * sqrt( ByBy * S )) / BB
C_y = A_y + (By * BL - sqrt( BxBx * S )) / BB

Here is an example Awk script, that generates three random points not too close to each other, and uses the general form above to solve the $x$ and $y$ coordinates for the third point, and verifies that the solution is correct. By default, the test is repeated for a million (different) points; you can specify the number of points to test using -v n=NUMBER.

#!/usr/bin/awk -f
function distance(xp1, yp1, zp1, xp2, yp2, zp2) { return sqrt((xp2-xp1)*(xp2-xp1) + (yp2-yp1)*(yp2-yp1) + (zp2-zp1)*(zp2-zp1)) }
BEGIN {
    srand();

    eps = 0.0005
    xmin = -9.999 ; ymin = -9.999 ; zmin = -9.999
    xmax = +9.999 ; ymax = +9.999 ; zmax = +9.999

    if (int(n) < 1)
        n = 1000000;
    else
        n = int(n);

    for (i = 0; i < n; i++) {

        if (!(i % 10000))
            printf "\r%.1f%% tests completed\033[K\r", 100.0*i/n

        while (1) {
            x1 = xmin + (xmax - xmin) * rand()
            y1 = ymin + (ymax - ymin) * rand()
            z1 = zmin + (zmax - zmin) * rand()
            x2 = xmin + (xmax - xmin) * rand()
            y2 = ymin + (ymax - ymin) * rand()
            z2 = zmin + (zmax - zmin) * rand()
            x3 = xmin + (xmax - xmin) * rand()
            y3 = ymin + (ymax - ymin) * rand()
            z3 = zmin + (zmax - zmin) * rand()
            dd12 = (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1)
            dd13 = (x3-x1)*(x3-x1) + (y3-y1)*(y3-y1) + (z3-z1)*(z3-z1)
            dd23 = (x3-x2)*(x3-x2) + (y3-y2)*(y3-y2) + (z3-z2)*(z3-z2)
            if (dd12 > eps*eps &&
                dd13 > eps*eps &&
                dd23 > eps*eps) break
        }

        # Preparation

        Bx = x2-x1 ; By = y2-y1 ; Bz = z2-z1 ; Cz = z3-z1
        BxBx = Bx*Bx ; ByBy = By*By ; BzBz = Bz*Bz
        ABAB = dd12 ; ACAC = dd13 ; BCBC = dd23

        # Calculations

        S = BxBx * BxBx - (ACAC - BCBC) * (ACAC - BCBC) - (ByBy + BzBz) * (ByBy + BzBz) + 2 * ABAB * (ACAC + BCBC - BxBx) + 4*( Bz * Cz * (ABAB + ACAC - BCBC) - ACAC*BzBz - ABAB*Cz*Cz )
        if (S < 0) {
            printf "A:   %7.3f %7.3f %7.3f\n", x1, y1, z1
            printf "B:   %7.3f %7.3f %7.3f\n", x2, y2, z2
            printf "C:   %7.3f %7.3f %7.3f\n", x3, y3, z3
            printf "AB: %11.6f\n", sqrt(dd12)
            printf "AC: %11.6f\n", sqrt(dd13)
            printf "BC: %11.6f\n", sqrt(dd23)
            printf "\n"
            printf "No solutions.\n"
            exit(1)
        }

        BB = BxBx + BxBx + ByBy + ByBy
        if (BB <= eps*eps) {
            printf "A:   %7.3f %7.3f %7.3f\n", x1, y1, z1
            printf "B:   %7.3f %7.3f %7.3f\n", x2, y2, z2
            printf "C:   %7.3f %7.3f %7.3f\n", x3, y3, z3
            printf "AB: %11.6f\n", sqrt(dd12)
            printf "AC: %11.6f\n", sqrt(dd13)
            printf "BC: %11.6f\n", sqrt(dd23)
            printf "\n"
            printf "Cannot solve, A and B differ only by their z coordinates.\n"
            exit(1)
        }

        BL = ACAC + ABAB - BCBC - 2*Bz*Cz

        if (Bx > 0) {
            if (By > 0)
                BS = +1;
            else
            if (By < 0)
                BS = -1;
            else
                BS = +1;
        } else
        if (Bx < 0) {
            if (By > 0)
                BS = -1;
            else
            if (By < 0)
                BS = +1;
            else
                BS = -1;
        } else {
             if (By > 0)
                BS = +1;
            else
            if (By < 0)
                BS = -1;
            else
                BS = 0;
        }

        Cx1 = x1 + (Bx * BL - BS * sqrt( ByBy * S )) / BB
        Cy1 = y1 + (By * BL + sqrt( BxBx * S )) / BB
        Cz1 = z3
        Cx2 = x1 + (Bx * BL + BS * sqrt( ByBy * S )) / BB
        Cy2 = y1 + (By * BL - sqrt( BxBx * S )) / BB
        Cz2 = z3

        err11 = distance(Cx1,Cy1,Cz1,x1,y1,z1) - sqrt(dd13)
        err12 = distance(Cx1,Cy1,Cz1,x2,y2,z2) - sqrt(dd23)
        err21 = distance(Cx2,Cy2,Cz2,x1,y1,z1) - sqrt(dd13)
        err22 = distance(Cx2,Cy2,Cz2,x2,y2,z2) - sqrt(dd23)

        if (err11 > eps || err11 < -eps ||
            err12 > eps || err12 < -eps ||
            err21 > eps || err21 < -eps ||
            err22 > eps || err22 < -eps) {
            printf "A:   %7.3f %7.3f %7.3f\n", x1, y1, z1
            printf "B:   %7.3f %7.3f %7.3f\n", x2, y2, z2
            printf "C:   %7.3f %7.3f %7.3f\n", x3, y3, z3
            printf "AB: %11.6f\n", sqrt(dd12)
            printf "AC: %11.6f\n", sqrt(dd13)
            printf "BC: %11.6f\n", sqrt(dd23)
            printf "\n"
            printf "Found: %7.3f %7.3f %7.3f errors: %+.6f %+.6f %+.6f\n", Cx1, Cy1, Cz1, distance(Cx1,Cy1,Cz1,x3,y3,z3), distance(Cx1,Cy1,Cz1,x1,y1,z1) - sqrt(dd13), distance(Cx1,Cy1,Cz1,x2,y2,z2) - sqrt(dd23)
            printf "Found: %7.3f %7.3f %7.3f errors: %+.6f %+.6f %+.6f\n", Cx2, Cy2, Cz2, distance(Cx2,Cy2,Cz2,x3,y3,z3), distance(Cx2,Cy2,Cz2,x1,y1,z1) - sqrt(dd13), distance(Cx2,Cy2,Cz2,x2,y2,z2) - sqrt(dd23)
            exit(1)
        }

    }

    printf "All %d tests successful.\n", n
}

In the output, in cases of an issue with the algorithm, the errors reported for the found points are the distance to the initially generated point, followed by the errors in distances to the two initial points.

Testing a few hundred million such points indicates the algorithm should work: that is, it finds the valid solution. If it finds two solutions, one of them is the initially generated point, but the other fulfills the requirements (distances to the other two points and the $z$ coordinate).

If $A$ and $B$ differ only in their $z$ coordinates, there is zero, one (degenerate triangle) or an infinite number of solutions. — amd, Jun 29 '17 at 23:47
The fully-expanded solutions might look messy, true, but to my mind these equations aren’t all that difficult to solve. After substituting for $z$ in the equations of the spheres, subtract one from the other, which gives you a linear equation in $x$ and $y$. From there, it’s a matter of back-substitution into either of the circle equations. One can also work up a fairly simple vector-based solution once the problem’s been reduced to intersecting a pair of coplanar circles. — amd, Jun 29 '17 at 23:56
Wow, thanks so much! I've got some real world points to try out tomorrow. I was getting stuck because points A and B can have the same X value so I was getting stuck dividing by 0. — Lance Miller, Jun 30 '17 at 04:54
@amd: Difficult, no. Straightforward, yes; but "horrible" (complicated), and tricky, for non-mathematicians. Let me show how a non-mathematician like me trips up here. The linear equation (in coordinates relative to point $A$) is $$Cx = \frac{AC^2 + AB^2 - BC^2 - 2(By Cy + Bz Cz)}{2 Bx}$$or$$Cy = \frac{AC^2 + AB^2 - BC^2 - 2(Bx Cx + Bz Cz)}{2 By}$$which look nice and simple. Yet, if $Bx = 0$ you cannot use the above to find $Cx$ when you already know $Cy$, or if $By = 0$ you cannot find $Cy$ when you already know $Cx$, due to the division. (This is what tripped Lance Miller, I believe?) [...] — Nominal Animal, Jun 30 '17 at 13:01
[...] In practice, we'd solve that by splitting the solution into two sub-paths, solving $Cy$ first if $\lvert Bx \rvert \ge \lvert By \rvert$, $Cx$ first otherwise. However, substituting the expression back to the original equations yields the "horrible" form, and because of that split, we have to implement two of those. (The two forms are relatively symmetric, and both have a divisor $2 (Bx^2 + By^2)$, so the only pitfall left in an implementation is to check that the argument of the square root is not negative: if it is, there are no real solutions.) [...] — Nominal Animal, Jun 30 '17 at 13:10
[...] What I did, instead, is to ask Maple for the "horrible" direct forms for $Cx$ and $Cy$, found their common terms, sorted and optimized (by hand) the common terms using temporary variables, and "inserted" check points at division(s) or square root(s) (and similar operations with singularities or limited argument ranges). To me, this is almost a "procedure" now, because I've implemented quite a few solutions using the same steps. The result is what you see above. Thing is, the above is what every programmer hopes they'd receive from a mathematician. [...] — Nominal Animal, Jun 30 '17 at 13:17
[...] Unfortunately, to most mathematicians, these are implementation details, and not at all interesting or useful in any way. (Also, most mathematicians are often not interested in how a solution is found, and tend to ignore that step completely, satisfied by proving that a solution exists, leaving the task of finding a solution for the reader/implementer. Which is infuriating for a programmer, because that is exactly what we need help with; we already know the solution exists.) [...] — Nominal Animal, Jun 30 '17 at 13:19
[...] What I try to show, in this answer, and in general with my answers, is that gray middle ground: how to find the solution in an efficient form that can be implemented in a computer program. I've done this several times after you or someone else has already shown the correct "mathematical" answer, or at least covered the major hurdles. However, I know from experience that that is not enough for non-mathematicians; that covering the middle ground, bridges the two approaches ("mathematical" and "programming") in a very fruitful way. — Nominal Animal, Jun 30 '17 at 13:23
@amd: If you, or anyone else, has a suggestion how to mark my answers (perhaps with a statement at the beginning, indicating that the answer only exists to help programmers implement the math, or find a particular solution?) in a way that mathematicians see the answer as being an ancillary note or guide for implementers, I would be delighted to use one. I can definitely say that there are many members (like myself) that need these, and find them useful, even if they are ancillary or pore into details not interesting to the core member group here, the mathematicians and math enthusiasts. — Nominal Animal, Jun 30 '17 at 13:27
@NominalAnimal I'm testing out your solution and I get the correct Y value of C. For x, since my points A and B have the same X value, that makes Bx and BS 0 which simplifies Cx to Cx = Ax. My test points (given by work spec) (x,y,z) = A(2445.7,-216.6,275.15), B(2445.7,216.6,275.15), C(2326.5,-109.5,256.25). Cx and Cy are passive axes on a tool that needs to pitch, roll, and yaw. — Lance Miller, Jun 30 '17 at 16:11
@LanceMiller: Thanks for the note: I'll recheck those cases (Bx=0 or By=0), and edit the answer to fix. — Nominal Animal, Jun 30 '17 at 16:36
@NominalAnimal This is really great! Thank you so much for your help! Your programmers breakdown is extremely useful too. I can understand the theory, but getting to a step by step process that a PLC can do can be a challenge at times. — Lance Miller, Jun 30 '17 at 18:41
@LanceMiller: Thanks. Many of the math answers I've seen here (by amd and others) are nothing short of brilliant, but there is still that gap between pure math and applying it in practice (especially for computer programs). My own math-fu is not that strong, but I seem to be able to help others bridge that gap somewhat. It pains me that I'm not able to express the fact that I see these answers as "implementation notes" for the proper answers; complementary, not competing, and that a "pure math" answer is still the proper one to the question. — Nominal Animal, Jun 30 '17 at 19:05
A direct vector computation might result in even simpler code. See the update to my answer. — amd, Jun 30 '17 at 20:20

3D Triangle with 2 Points, All Lengths, and Z of the 3rd Point Given

2 Answers2