Circumsphere of a tetrahedron

Question

I'm trying to implement Bowyer-Watson algorithm. Currently performance is stuck behind calculating the circumsphere of a tetrahedron.

I tried using this http://mathworld.wolfram.com/Circumsphere.html. It works but it is slow.

Next i tried this https://www2.mps.mpg.de/homes/daly/CSDS/t4h/tetra.htm#Q1-1-7. But this has the problem of not always returning a solution. For example vertices (-1000, 0, -1000), (1000,0,-1000), (0,0,1000), (0,1000,0). I mostly keep getting -Infinity or NaN.

Is there a way to make the first one faster (leave some equations out, make some equations simpler), or make the second one always work (so that the linear equation system can handle vectors that can have 2 coordinates 0)? Or is there some third way that is even better?

EDIT

THe way i tried using to solve the linear system.

2d1x*x + 2d1y*y + 2d1z*z = d1^2   (1)
2d2x*x + 2d2y*y + 2d2z*z = d2^2   (2)
2d3x*x + 2d3y*y + 2d3z*z = d3^2   (3)


(2) - d2x/d1x *(1)
(3) - d3x/d1x *(1)

2d1x*x + 2d1y*y + 2d1z*z = d1^2   (1)
 0    + d2x/d1x * 2d2y*y + d2x/d1x * 2d2z*z = d2x/d1x * d2^2   (2)
 0    + d3x/d1x * 2d3y*y + d3x/d1x * 2d3z*z = d3x/d1x * d3^2   (3)

(3) - (d2y*d3x)/(d2x*d3y) *(2)

2d1x*x + 2d1y*y + 2d1z*z = d1^2   (1)
0    + d2x/d1x * 2d2y*y + d2x/d1x * 2d2z*z = d2x/d1x * d2^2   (2)
0    +    0             + (d3x*d3y)/(d2x*d2y) * 2d3z*z = dx3/d1x * 2d2y * d3^2   (3)

Answer 1

Following are some alternatives.

Alternative I - compute barycentric coordinates of circumcenter $O$ directly.

Given a tetrahedron $T$ with vertices $v_0, v_1, v_2, v_3$. Let

• $V$ and $R$ be the volume and circumradius of $T$.
• $S_i$ and $R_i$ be the area and circumradius of the face opposite to vertex $v_i$.
• $\ell_{ij}$ be the distance between vertex $v_i$ and $v_j$.

Once you have computed $V$ and $S_i$, the circumradius $R$ is given by the formula${}^{\color{blue}{[1]}}$.

$$6VR = \sqrt{p(p-aa_1)(p-bb_1)(p-cc_1)} \quad\text{ where }\quad \begin{cases} a &= \ell_{01}, a_1 = \ell_{23}\\ b &= \ell_{02}, b_1 = \ell_{31}\\ c &= \ell_{03}, c_1 = \ell_{12}\\ p &= \frac12(aa_1 + bb_1 + cc_1) \end{cases} $$ The circumradii $R_i$ for the faces can be computed by corresponding $2$-d counterparts: $$4S_i R_i = \prod_{j < k;\\ j \ne i,k\ne i}\ell_{jk} \quad\iff\quad \begin{cases} 4 S_0 R_0 &= \ell_{12}\ell_{13}\ell_{23} = a_1b_1c_1\\ 4 S_1 R_1 &= \ell_{02}\ell_{03}\ell_{23} = a_1 b c \\ 4 S_2 R_2 &= \ell_{01}\ell_{03}\ell_{13} = ab_1 c\\ 4 S_3 R_3 &= \ell_{01}\ell_{02}\ell_{13} = ab c_1 \end{cases} $$ The barycentric coordinates for circumcenter $O$, i.e. those four numbers $\lambda_0,\ldots,\lambda_3$ such that

$$O = \sum_{i=0}^3 \lambda_i v_i\quad\text{ with }\quad \sum_{i=0}^3 \lambda_i = 1$$

is then given by the formula

$$3V\lambda_i = S_i\sqrt{R^2 - R_i^2}$$

Alternative II - compute circumcenter $O$ using dot and cross products.

Same notation as previous alternative. Let $u_i = v_i - v_0$ for $i = 1,2,3$ and $p = O - v_0$.
We have $$R^2 = \|p\|^2 = \| u_1 - p \|^2 = \| u_2 - p \|^2 = \| u_3 - p \|^2$$ Cancelling the common term $R^2$ from these equalities, we obtain $$ u_1 \cdot p = \frac12\ell_{01}^2,\quad u_2 \cdot p = \frac12\ell_{02}^2,\quad\text{ and }\quad u_3 \cdot p = \frac12\ell_{03}^2, $$ We can invert these equations using the dual basis for the three vectors $u_1, u_2, u_3$. The end result is following expression of $O$ which only need dot and cross products among vectors $u_k$.

$$O = v_0 + \frac{ \ell_{01}^2 ( u_2 \times u_3 ) + \ell_{02}^2 ( u_3 \times u_1 ) + \ell_{03}^2 ( u_1 \times u_2 )}{ 2 u_1\cdot(u_2 \times u_3)} $$

Notes/References

• $\color{blue}{[1]}$ - I. Todhunter, Spherical Trigonometry: For the Use of Colleges and Schools, $\S163$ (1886). An online copy is available under Project Gutenberg.

Answer 2

I'm fond of Miroslav Fiedler's method for determining the circumsphere. The starting point is the usual Cayley-Menger matrix:

$$\mathbf C=\begin{pmatrix} 0&1&1&1&1\\ 1&0&d_{12}&d_{13}&d_{14}\\ 1&d_{12}&0&d_{23}&d_{24}\\ 1&d_{13}&d_{23}&0&d_{34}\\ 1&d_{14}&d_{24}&d_{34}&0\end{pmatrix}$$

where $d_{jk}$ is the squared distance between points $\mathbf p_j$ and $\mathbf p_k$ of the tetrahedron.

From there, one computes $\mathbf M=-2\mathbf C^{-1}$, and then gets the first row of the inverse. The circumradius is given by $\frac12\sqrt{m_{1,1}}$, and the circumcenter is given by the weighted average $\dfrac{m_{1,2}\mathbf p_1+m_{1,3}\mathbf p_2+m_{1,4}\mathbf p_3+m_{1,5}\mathbf p_4}{m_{1,2}+m_{1,3}+m_{1,4}+m_{1,5}}$. (Alternatively, one can interpret the $m_{1,k}$ as unnormalized barycentric coordinates.)

---

P.S.

I gave a Mathematica demonstration of Fiedler's method here; the routine there also uses Cayley-Menger to compute the insphere. Here is the result of using the methods there on the OP's example tetrahedron:

Answer 3

Layman here. After looking at the tons of square roots, inversion of 5x5 matrix, I went back to Wikipedia and read it again. While there are also tons of square roots, somewhere in the lower-middle of the Tetrahedron page, there is a super-simple expression with a barycentric approach (at least that's what it looks like to me).
Of course, it's not Wikipedia, but Lévy, Bruno; Liu, Yang (2010). "Lp Centroidal Voronoi Tessellation and its applications". ACM: 119.

$c=M^{-1}b$ , where $M=\begin{pmatrix} [x_1-x_0]^T\\ [x_2-x_0]^T\\ [x_3-x_0]^T\end{pmatrix}$ and  $b=\frac1 2\begin{bmatrix} \|x_1\|^2-\|x_0\|^2\\ \|x_2\|^2-\|x_0\|^2\\ \|x_3\|^2-\|x_0\|^2\end{bmatrix}$,
$M$ is a 3x3 matrix, everything else is a column vector ($c$ will be the circumcenter, $x_0...x_3$ are the corners of the tetrahedron).

So far it appears to work for me. In fact I also needed this for Bowyer-Watson, and at the end I did not calculate a single square root, as storing $r^2=\|c-x_0\|^2$ is more efficient for an in/out comparison with Cartesian coordinates.

Ps: it's probably the same thing as "Alternative II" from the accepted answer, but that leaves a bit to the imagination.