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.
WebGL demo
Javascript source code (Editable version)
Derivation
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. Alternatively, a perpendicular edge bisector can be defined as the set of points equidistant from two vertices of the tetrahedron.
Formula
A compact expression of 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 a 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 \|^2 - \| \mathbf v_0 \|^2 \\
\| \mathbf v_2 \|^2 - \| \mathbf v_0\|^2 \\
\| \mathbf v_3 \|^2 - \| \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 |}
\left (
\| \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 )
\right )
\]
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:
- Another implementation on gamestack
- Other tetrahedron properties (volume, height etc.)
- Alternative formulation
- matlab implementation
- Alternate C++ implementation
- Geometric Tools For Computer Graphics
No comments