Finding the coordinates of the fourth vertex of tetrahedron, given coordinates of "base" vertices and the distances to them

Question

I have a tetrahedron defined as:

• "base" vertices $P$, $Q$, $R$ are given.
• length of "remain" edges $L_P$, $L_Q$, and $L_R$ are also given.

I need to find the 4th vertex coordinates $(x, y, z)$. The image below describes my problem:

I understand that there are two symmetric solutions, one where the vertex is up and another when the vertex is below the plane defined by $(P, Q, R)$.

I tried to solve this problem by considering 3 spheres $S_P$, $S_Q$, $S_R$ with center on $P$, $Q$, $R$ and radius $L_P$, $L_Q$, and $L_R$, respectively. I'm wondering if there is an easier straighforward way to solve this.

Answer 1

Finding $(x, y, z)$ as the intersection of three spheres of radius $L_P$, $L_Q$, and $L_R$, centered at $P = (P_x, P_y, P_z)$, $Q = (Q_x, Q_y, Q_z)$, and $R = (R_x, R_y, R_z)$, respectively, is the solution.

However, if you rotate and translate the coordinate system, you can simplify the math a lot. (This is nothing special; it's just that when most of the coordinates are zeroes, the expressions simplify a lot.)

Rotate and translate the coordinate system (we'll use $(u, v, w)$ for the rotated and translated coordinates for clarity; note that distances are unchanged), $P$ is at origin $(0, 0, 0)$, $Q$ is at $(U_Q, 0, 0)$, and $R$ at $(U_R, V_R, 0)$. Then, the fourth vertex is at $$\begin{aligned} u &= \frac{L_P^2 - L_Q^2 + U_Q^2}{2 U_Q} \\ v &= \frac{L_P^2 - L_R^2 + U_R^2 + V_R^2 - 2 U_R u}{2 V_R} \\ w &= \pm\sqrt{L_P^2 - u^2 - v^2} \\ \end{aligned}$$

Rotating and translating the coordinate system is not difficult, either: $$\begin{aligned} U_Q &= \left\lVert Q - P \right\rVert \\ \hat{u} &= \frac{Q - P}{U_Q} \\ \vec{t} &= (R - P) - \hat{u}\bigr(\hat{u} \cdot (R - P)\bigr) \\ \hat{v} &= \frac{\vec{t}}{\left\lVert \vec{t} \right\rVert} \\ \hat{w} &= \hat{u} \times \hat{v} \\ U_R &= (R - P) \cdot \hat{u} \\ V_R &= (R - P) \cdot \hat{v} \\ \end{aligned}$$ Conversion back to original coordinates is similarly trivial: $$\vec{p} = P + u \hat{u} + v \hat{v} + w \hat{w}$$

---

Here is a Python 3 implementation:

# SPDX-License-Identifier: CC0-1.0
# This file is in Public Domain.

from vector import Vector, sqrt

def find_fourth_vertex(vertex1, vertex2, vertex3, distance1, distance2, distance3):
    # Use Vector type for the vertices
    p1 = Vector(vertex1[0], vertex1[1], vertex1[2])
    p2 = Vector(vertex2[0], vertex2[1], vertex2[2])
    p3 = Vector(vertex3[0], vertex3[1], vertex3[2])

    # Use float type for the distances
    r1 = float(distance1)
    r2 = float(distance2)
    r3 = float(distance3)

    u_axis = (p2 - p1).unit
    v_axis = (p3 - p1).perp(u_axis).unit
    w_axis = u_axis ^ v_axis

    u2 = (p2 - p1) | u_axis
    u3 = (p3 - p1) | u_axis
    v3 = (p3 - p1) | v_axis

    u = (r1*r1 - r2*r2 + u2*u2) / (2*u2)
    v = (r1*r1 - r3*r3 + u3*u3 + v3*v3 - 2*u*u3) / (2*v3)
    w = sqrt(r1*r1 - u*u - v*v)

    return (p1 + u*u_axis + v*v_axis + w*w_axis,
            p1 + u*u_axis + v*v_axis - w*w_axis)

if __name__ == '__main__':
    from math import sin, cos, pi
    from random import Random

    prng = Random()

    while True:
        # Generate four random vectors in (-9.9,-9.9,-9.9) - (+9.9,+9.9,+9.9)
        v = [ Vector(prng.uniform(-9.9, 9.9), prng.uniform(-9.9, 9.9), prng.uniform(-9.9, 9.9)),
              Vector(prng.uniform(-9.9, 9.9), prng.uniform(-9.9, 9.9), prng.uniform(-9.9, 9.9)),
              Vector(prng.uniform(-9.9, 9.9), prng.uniform(-9.9, 9.9), prng.uniform(-9.9, 9.9)),
              Vector(prng.uniform(-9.9, 9.9), prng.uniform(-9.9, 9.9), prng.uniform(-9.9, 9.9)) ]
        # Find their minimum pairwise distance
        rmin = (v[1] - v[0]).norm
        for i in range(0, len(v) - 1):
            for j in range(i+1, len(v)):
                rmin = min(rmin, (v[j] - v[i]).norm)
        # If they're all least 1 unit from each other, accept.
        if rmin >= 1:
            break

    v1 = v[0]
    v2 = v[1]
    v3 = v[2]

    r1 = (v[3] - v[0]).norm
    r2 = (v[3] - v[1]).norm
    r3 = (v[3] - v[2]).norm

    print("v1 = %s, distance %f" % (v1, r1))
    print("v2 = %s, distance %f" % (v2, r2))
    print("v3 = %s, distance %f" % (v3, r3))

    v4a, v4b = find_fourth_vertex(v1, v2, v3, r1, r2, r3)

    print("v4 == %s" % (v[3],))
    print("v4a = %s" % (v4a,))
    print("v4b = %s" % (v4b,))

    print("v4a distances: %f, %f, %f" % ((v4a-v1).norm, (v4a-v2).norm, (v4a-v3).norm))
    print("v4b distances: %f, %f, %f" % ((v4b-v1).norm, (v4b-v2).norm, (v4b-v3).norm))

where a.perp(b) is $\vec{a} - \vec{b}(\vec{a}\cdot\vec{b})$, a | b is $\vec{a} \cdot \vec{b}$ and a ^ b is $\vec{a} \times \vec{b}$.

When run, it generates a test tetrahedron, and displays the results when find_fourth_vertex is given three of the vertices and their distances to the fourth.

The helper Vector class is implemented by vector.py:

# SPDX-License-Identifier: CC0-1.0
# This file is in Public Domain.

from math import sqrt

class Vector(tuple):
    """Tuple subclass implementing basic 3D vectors"""

    def __new__(cls, x, y, z):
        return tuple.__new__(cls, (float(x), float(y), float(z)))

    def perp(self, other):
        """Part perpendicular to other vector"""
        dp = self[0]*other[0] + self[1]*other[1] + self[2]*other[2]
        return Vector(self[0] - dp*other[0],
                      self[1] - dp*other[1],
                      self[2] - dp*other[2])

    @property
    def unit(self):
        """Scaled to unit length"""
        n = sqrt(self[0]*self[0] + self[1]*self[1] + self[2]*self[2])
        return Vector(self[0]/n, self[1]/n, self[2]/n)

    @property
    def norm(self):
        """Euclidean length"""
        return sqrt(self[0]*self[0] + self[1]*self[1] + self[2]*self[2])

    @property
    def normsqr(self):
        """Euclidean length squared"""
        return self[0]*self[0] + self[1]*self[1] + self[2]*self[2]

    @property
    def x(self):
        """Vector x coordinate"""
        return self[0]

    @property
    def y(self):
        """Vector y coordinate"""
        return self[1]

    @property
    def z(self):
        """Vector z coordinate"""
        return self[2]

    def __bool__(self):
        """Nonzero vector"""
        return (self[0]*self[0] + self[1]*self[1] + self[2]*self[2] > 0)

    def __abs__(self):
        """abs(a): Euclidean length of vector a"""
        return sqrt(self[0]*self[0] + self[1]*self[1] + self[2]*self[2])

    def __add__(self, other):
        """a + b: Vector addition"""
        if isinstance(other, (tuple, list)) and len(other) >= 3:
            return Vector(self[0]+other[0], self[1]+other[1], self[2]+other[2])
        else:
            return NotImplemented

    def __radd__(self, other):
        """b + a: Vector addition"""
        if isinstance(other, (tuple, list)) and len(other) >= 3:
            return Vector(other[0]+self[0], other[1]+self[1], other[2]+self[2])
        else:
            return NotImplemented

    def __mul__(self, other):
        """a * b: Scalar multiplication"""
        if isinstance(other, (int, float)):
            return Vector(self[0]*other, self[1]*other, self[2]*other)
        else:
            return NotImplemented

    def __rmul__(self, other):
        """b * a: Scalar multiplication"""
        if isinstance(other, (int, float)):
            return Vector(other*self[0], other*self[1], other*self[2])
        else:
            return NotImplemented

    def __neg__(self):
        """-a: Negation"""
        return Vector(-self[0], -self[1], -self[2])

    def __or__(self, other):
        """a | b: Dot product"""
        if isinstance(other, (tuple, list)) and len(other) >= 3:
            return self[0]*other[0] + self[1]*other[1] + self[2]*other[2]
        else:
            return NotImplemented

    def __ror__(self, other):
        """b | a: Dot product"""
        if isinstance(other, (tuple, list)) and len(other) >= 3:
            return other[0]*self[0] + other[1]*self[1] + other[2]*self[2]
        else:
            return NotImplemented

    def __sub__(self, other):
        """a - b: Vector subtraction"""
        if isinstance(other, (tuple, list)) and len(other) >= 3:
            return Vector(self[0]-other[0], self[1]-other[1], self[2]-other[2])
        else:
            return NotImplemented

    def __rsub__(self, other):
        """b - a: Vector subtraction"""
        if isinstance(other, (tuple, list)) and len(other) >= 3:
            return Vector(other[0]-self[0], other[1]-self[1], other[2]-self[2])
        else:
            return NotImplemented

    def __truediv__(self, other):
        """a / b: Scalar division"""
        if isinstance(other, (int, float)):
            return Vector(self[0]/other, self[1]/other, self[2]/other)
        else:
            return NotImplemented

    def __xor__(self, other):
        """a ^ b: Vector cross product"""
        if isinstance(other, (tuple, list)) and len(other) >= 3:
            return Vector(self[1]*other[2] - self[2]*other[1],
                          self[2]*other[0] - self[0]*other[2],
                          self[0]*other[1] - self[1]*other[0])
        else:
            return NotImplemented

    def __rxor__(self, other):
        """b ^ a: Vector cross product"""
        if isinstance(other, (tuple, list)) and len(other) >= 3:
            return Vector(other[1]*self[2] - other[2]*self[1],
                          other[2]*self[0] - other[0]*self[2],
                          other[0]*self[1] - other[1]*self[0])
        else:
            return NotImplemented

that you can just put in the same directory as the example Python file. Run pydoc3 vector in that directory to see the API description for it.

Note that vector.py defines a generic 3D Euclidean vector class with basic vector algebra operations, and is in no way specific to this particular problem.

Answer 2

After some work on my own question, I think I have found an alternative way to solve this problem.

The objective is to find the vertex $E$ of a tetrahedron defined as:

• Points $P$, $Q$ and $R$
• Distances $||\vec{PE}||$, $||\vec{QE}||$ and $||\vec{RE}||$

In this solution, $E$ can be achieved by finding the angles $\sigma$ and $\theta$ in order to construct a vector $\vec{PE}$.

Note that $\theta$ is the angle between the unknow vector $\vec{PE}$ and the plane defined by the points $P$, $Q$ and $R$. $\sigma$ is the angle between the projection of $\vec{PE}$ on the same plane $PQR$ and the vector $\vec{PR}$.

As the image suggests, $\sigma$ and $\theta$ can be obtained in a straightforward way from the tetrahedron height and elementary trigonometric properties, as shown below.

Finding $\vec{PE}$ angles $\sigma$ and $\theta$

1. Find the tetrahedron $Volume$ using Calyer-Menger determinant:

$$288 Volume^2 = \left|\begin{matrix}0 & 1 & 1 & 1 & 1\cr 1 & 0 & ||\vec{RE}||^{2} & ||\vec{PE}||^{2} & ||\vec{QE}||^{2}\cr 1 & ||\vec{RE}||^{2} & 0 & \tilde||\vec{QE}||^{2} & \tilde||\vec{PE}||^{2}\cr 1 & ||\vec{PE}||^{2} & \tilde||\vec{QE}||^{2} & 0 & \tilde||\vec{RE}||^{2}\cr 1 & ||\vec{QE}||^{2} & \tilde||\vec{PE}||^{2} & \tilde||\vec{RE}||^{2} & 0\end{matrix}\right|$$

2. Find the $Area$ of triangle $P$, $Q$, $R$ using Heron's formula:

$$Area = \frac{1}{4}\sqrt{4||\vec{PE}||^2||\vec{QE}||^2-(||\vec{PE}||^2+||\vec{QE}||^2-||\vec{RE}||^2)^2}$$

3. Find the tetrahedron height $H$ using the relationship between $Volume$ and $Area$:

$$H = \frac{3\times Volume}{Area}$$

4. Find $\theta$:

$$\theta = arcsin\left (\frac{H}{||\vec{PE}||}\right )$$

Once we have $\theta$ the next step is to find the length of the projections $\vec{PE'}$ and $\vec{RE'}$ onto the plane defined by $P$, $Q$ and $R$:

$$||\vec{PE'}|| = \sqrt{||\vec{PE}||^2 - H^2}$$ $$||\vec{RE'}|| = \sqrt{||\vec{RE}||^2 - H^2}$$

5. Thus, using the cosine law, $\sigma$ is given by:

$$\sigma = arccos\left (\frac{||\vec{PE'}||^2 - ||\vec{RE'}||^2 + ||\vec{PR}||^2}{2 ||\vec{PE'}|| \times ||\vec{PR}||}\right )$$

Once we have $P$, $||\vec{PE}||$, $\sigma$ and $\theta$ we know everything we need to find $E$.

Finding $E$ given $\sigma$, $\theta$, $P$ and $||\vec{PE}||$

There are several ways to obtain $E(x, y, z)$, one of them is rotating $\vec{PR}$ by $\sigma$ and then rotating again by $\theta$, as demonstrated below.

1. Find the triangle $PQR$ normal $\vec{n}$:

$$\vec{n} = \frac{\vec{PR}\times\vec{PQ}}{||\vec{PR}|| \times ||\vec{PQ}||}$$

2. Rotate $\vec{PR}$ about $\vec{n}$ by $-\sigma$ using Rodrigues' formula:

$$\vec{PE'} = \vec{PR}cos(-\sigma) + (\vec{n} \times \vec{PR})\sin(-\sigma) + \vec{n}(\vec{n} \cdot \vec{PR}) (1 - cos(-\sigma))$$

3. Find the normal $\vec{m}$ from $\vec{PE'}$ and $\vec{n}$:

$$\vec{m} = \frac{\vec{PE'}\times\vec{n}}{||\vec{PE'}|| \times ||\vec{n}||}$$

4. Rotate $\vec{PE'}$ by $-\theta$ about $\vec{m}$:

$$\vec{PE_{dir}} = \vec{PE'}cos(-\theta) + (\vec{m} \times \vec{PE'})\sin(-\theta) + \vec{m}(\vec{m} \cdot \vec{PE'}) (1 - cos(-\theta))$$

5. Get the unit vector from $\vec{PE_{dir}}$ and multiply it by $||\vec{PE}||$ in order to obtain $\vec{PE}$:

$$\vec{PE} = \frac{\vec{PE_{dir}}}{||\vec{PE_{dir}}||} \times ||\vec{PE}||$$

Finally, $E$ is given by

$$E = \vec{PE} + P$$

It is noteworthy that the symmetric solution $E_2$ can be find by rotating $\vec{PE'}$ about $\vec{m}$ by $+\theta$ (instead of $-\theta$):

One of my future work is checking out if this approach is less computational intensive than others.

Follow some images from an experiment where $E$ is obtained by the procedure described here. This program can be visualized here: https://doleron.github.io/tetrahedron-4th-vertex/ and the source code is here: https://github.com/doleron/tetrahedron-4th-vertex

Note that the spheres are there only for comparision purposes.