Find a tetrahedron circumcenter

Tetrahdron's circumcenter

Find the center (a.k.a circumcenter) of the circumscribed sphere of a Tetrahedron.

Definition

For all tetrahedra, there exists a sphere called the circumsphere which completely encloses the tetrahedron. The tetrahedron's vertices all lie on the surface of its circumsphere. The point at the centre of the circumsphere is called the circumcentre.

Formula

The circumcenter of a tetrahedron can be computed as intersection of three bisector planes. A bisector plane is defined as the plane centered on, and orthogonal to an edge of the tetrahedron.

With this definition, the circumcenter \( \mathbf c \) of a tetrahedron with vertices \(\mathbf v_0, \mathbf v_1, \mathbf v_2, \mathbf v_3\) can be formulated as matrix product:

\[ \mathbf c = \mathbf A^{-1} \mathbf B  \text{  with  }
\mathbf A =
\left(\begin{matrix}
\left[\mathbf v_1 - \mathbf v_0\right]^T \\
\left[\mathbf v_2 - \mathbf v_0\right]^T \\
\left[\mathbf v_3 - \mathbf v_0\right]^T
\end{matrix}\right)   \text{  and  }
\mathbf B =
\frac{1}{2}
\left(\begin{matrix}
\| \mathbf v_1 - \mathbf v_0 \|^2 \\
\| \mathbf v_2 - \mathbf v_0\|^2 \\
\| \mathbf v_3 - \mathbf v_0\|^2
\end{matrix}\right) 
\]

Where:
\[
\vec e_{1} = \left[\mathbf v_1 - \mathbf v_0\right]^T  =  \{ v_{1x} - v_{0x}, v_{1y} - v_{0y}, v_{1z} - v_{0z} \} \\
\vec e_{2} = \left[\mathbf v_2 - \mathbf v_0\right]^T  =  \{ v_{2x} - v_{0x}, v_{2y} - v_{0y}, v_{2z} - v_{0z} \} \\
\vec e_{3} = \left[\mathbf v_3 - \mathbf v_0\right]^T  =  \{ v_{3x} - v_{0x}, v_{3y} - v_{0y}, v_{3z} - v_{0z} \}
\]

are the row vectors representing an edge of the tetrahedron and:

\[ \| \mathbf p \|^2 = p_x^2 + p_y^2 + p_z^2 \]

is the squared Euclidean norm.

In contrast to the centroid, the circumcenter may not always lay on the inside of a tetrahedron.
Analogously to an obtuse triangle, the circumcenter is outside of the object for an obtuse tetrahedron.


Expanded expression

Let's expand the above matrix formula in more details, here we use \( \vec e_{i} = \mathbf v_i - \mathbf v_0 \) for brievity and \( \times \) is the usual vector cross-product:
\[
c = \mathbf v_0 + \frac{ 1 } {
2 ~ det \left| \mathbf A \right |
} \| \vec e_{3} \|^2 \big (  \vec e_{1}  \times \vec e_{2} \big ) + \| \vec e_{2} \|^2 \big (  \vec e_{3}  \times  \vec e_{1}  \big ) + \| \vec e_{1} \|^2 \big (  \vec e_{2}  \times \vec e_{3} \big )

\]

The determinant of the \( 3 \times 3 \) \( \mathbf A \) matrix can be expanded as well recursively:

\[
\begin{align*}
det \left| \mathbf A \right | &= \left| \begin{matrix}
v_{1x} - v_{0x} & v_{1y} - v_{0y} & v_{1z} - v_{0z} \\
v_{2x} - v_{0x} & v_{2y} - v_{0y} & v_{2z} - v_{0z} \\
v_{3x} - v_{0x} & v_{3y} - v_{0y} & v_{3z} - v_{0z}
\end{matrix} \right | \\ ~\\
&=
  \begin{vmatrix}
    e_{1x} & e_{1y} & e_{1z} \\
    e_{2x} & e_{2y} & e_{2z} \\
    e_{3x} & e_{3y} & e_{3z}
  \end{vmatrix} \\~\\
  &=
  e_{1x}
  \begin{vmatrix}
    e_{2y} & e_{2z} \\
    e_{3y} & e_{3z}
  \end{vmatrix}
   - e_{1y}
   \begin{vmatrix}
      e_{2x} & e_{2z} \\
      e_{3x} & e_{3z}
   \end{vmatrix}
   + e_{1z}
   \begin{vmatrix}
      e_{2x} & e_{2y} \\
      e_{3x} & e_{3y}
   \end{vmatrix}\\ ~\\
&=  e_{1x}e_{2y}e_{3z} + e_{1y}e_{2z}e_{3x} + e_{1z}e_{2x}e_{3y} - e_{1z}e_{2y}e_{3x} - e_{1y}e_{2x}e_{3z} - e_{1x}e_{2z}e_{3y}
\end{align*}
\]

C/C++ implementation

Notes on stability: these expressions are unstable only in one case:  if the denominator is close to zero.  This instability, which arises if the tetrahedron is nearly degenerate, is unavoidable.  Depending on your application, you may want to use exact arithmetic to compute the value of the determinant.

Inputs are the vertices coordinates `a[]' `b[]' `c[]' and `d[]'. The result is returned both in terms of x y z coordinates in circumcenter[] and xi-eta-zeta coordinates. Both outputs are relative to the tetrahedron's vertex `a' (that is, `a' is  the origin of both coordinate systems).  In other words, the x y z coordinates returned are NOT absolute! One must add the coordinates of `a' to find the absolute coordinates of the circumcircle. 

The xi-eta-zeta are coordinate along the tetrahedron's edges.  Vertex `a' being the origin, xi axis represent one unit along the edge `ab'. One unit of the eta axis is the edge `ac'  and logically zeta is the coordinate along the edge `ad'. These coordinate values are useful for linear interpolation.                                                          

Note: if `xi' is NULL on input, the xi-eta-zeta coordinates will not be computed.

void tetrahedron_circumcenter(
        // In:
        const double a[3],
        const double b[3],
        const double c[3],
        const double d[3],
        // Out:
        double circumcenter[3],
        double* xi,
        double* eta,
        double* zeta)
{
    double denominator;


    // Use coordinates relative to point `a' of the tetrahedron.

    // ba = b - a
    double ba_x = b[0] - a[0];
    double ba_y = b[1] - a[1];
    double ba_z = b[2] - a[2];

    // ca = c - a
    double ca_x = c[0] - a[0];
    double ca_y = c[1] - a[1];
    double ca_z = c[2] - a[2];

    // da = d - a
    double da_x = d[0] - a[0];
    double da_y = d[1] - a[1];
    double da_z = d[2] - a[2];

    // Squares of lengths of the edges incident to `a'.
    double len_ba = ba_x * ba_x + ba_y * ba_y + ba_z * ba_z;
    double len_ca = ca_x * ca_x + ca_y * ca_y + ca_z * ca_z;
    double len_da = da_x * da_x + da_y * da_y + da_z * da_z;

    // Cross products of these edges.

    // c cross d
    double cross_cd_x = ca_y * da_z - da_y * ca_z;
    double cross_cd_y = ca_z * da_x - da_z * ca_x;
    double cross_cd_z = ca_x * da_y - da_x * ca_y;

    // d cross b
    double cross_db_x = da_y * ba_z - ba_y * da_z;
    double cross_db_y = da_z * ba_x - ba_z * da_x;
    double cross_db_z = da_x * ba_y - ba_x * da_y;

    // b cross c
    double cross_bc_x = ba_y * ca_z - ca_y * ba_z;
    double cross_bc_y = ba_z * ca_x - ca_z * ba_x;
    double cross_bc_z = ba_x * ca_y - ca_x * ba_y;

    // Calculate the denominator of the formula.
    denominator = 0.5 / (ba_x * cross_cd_x + ba_y * cross_cd_y + ba_z * cross_cd_z);

    // Calculate offset (from `a') of circumcenter.
    double circ_x = (len_ba * cross_cd_x + len_ca * cross_db_x + len_da * cross_bc_x) * denominator;
    double circ_y = (len_ba * cross_cd_y + len_ca * cross_db_y + len_da * cross_bc_y) * denominator;
    double circ_z = (len_ba * cross_cd_z + len_ca * cross_db_z + len_da * cross_bc_z) * denominator;

    circumcenter[0] = circ_x;
    circumcenter[1] = circ_y;
    circumcenter[2] = circ_z;

    if (xi != (double*) nullptr) {
        // To interpolate a linear function at the circumcenter, define a
        // coordinate system with a xi-axis directed from `a' to `b',
        // an eta-axis directed from `a' to `c', and a zeta-axis directed
        // from `a' to `d'.  The values for xi, eta, and zeta are computed
        // by Cramer's Rule for solving systems of linear equations.
        denominator *= 2.0;
        *xi   = (circ_x * cross_cd_x + circ_y * cross_cd_y + circ_z * cross_cd_z) * denominator;
        *eta  = (circ_x * cross_db_x + circ_y * cross_db_y + circ_z * cross_db_z) * denominator;
        *zeta = (circ_x * cross_bc_x + circ_y * cross_bc_y + circ_z * cross_bc_z) * denominator;
    }
}

References

Find C/C++ code for the circumcenter and circumradius of tetrahedrons, 3D and 2D triangles here:

Other:

No comments

All html tags except <b> and <i> will be removed from your comment. You can make links by just typing the url or mail-address.
Anti-spam question: