# Insphere for Irregular Tetrahedron

I am looking for existing Mathematica code to compute the unique sphere inscribed inside an irregular tetrahedron. I can write it myself, but I would love to find that someone already performed this calculation and can share their code.

I want to take as input four points {a, b, c, d} in 3D and return the center, radius, and four tangency points of the insphere. If no one has already done this, then I will go ahead and write it myself.

Responding to the request for an example:

a = {0, 0, 0};
b = {1, 0, 0};
c = {2, 1, 0};
d = {1, 3/2, 1}; Updated. (11 Oct. 2013) The Wikipedia pages to which ybeltukov pointed are correct. Here the inradius is all of $0.127$: I have the insphere center and radius, but not yet the four tangency points…

And here are the four points of tangency: • Care to share the actual code for the tetrahedron graphics, too? Would be handy. You can also format code etc., have a look at the format help. – Yves Klett Oct 11 '13 at 13:30
• @YvesKlett, Graphics3D@{Opacity[0.5],Polygon[{a,b,c,d}~Subsets~{3}]} – Simon Woods Oct 11 '13 at 16:00
• There is a simple general formula in the Wikipedia. Unfortunately it is unproved. You can check it and post it as the answer. It will be your contribution to the community. – ybeltukov Oct 11 '13 at 19:38
• "not yet the four tangency points" It should be easy enough to find the orthogonal projections of the incenter onto the four face planes, no? – Rahul Oct 12 '13 at 1:03
• The MO question where this was used: "Tetrahedron insphere iteration." – Joseph O'Rourke Aug 24 '16 at 19:11

## 5 Answers

[Edit: I found this method a rather pleasing application of analytic geometry, so I rewrote the explanation hopefully to do it justice. In fact, it can be applied in any dimension, and I've updated the code to be general]

Here's my way:

faces[simplex_] := Partition[simplex, Length@simplex - 1, 1, 1];

(* outward-oriented unit normals to each of faces[a,b,c,d] *)
normals[simplex_] := With[{df = Transpose[Transpose @ Rest[#] - First[#]]},
Normalize[(-#.Last@df) # &[Cross @@ (df[[;; Length @ simplex - 2]])]]] & /@
Partition[simplex, Length @ simplex, 1, 1];

insphere[simplex_?(Subtract @@ Dimensions[#] == 1 &)] :=
Module[{x, center, radius, tangentpoints},
With[{normals = normals[simplex]},
center = Array[x, Length @ simplex - 1];
tangentpoints = center + radius # & /@ normals;
insphere[center, radius, tangentpoints] /. First @ Solve[
MapThread[
Function[{f, n}, (center + radius n  -  First[f]) . n == 0],
{faces[simplex], normals}],
Flatten[{center, radius}]]
]]

(* some interfaces to the insphere data *)
insphere[center_, radius_, tangentpoints_]["Sphere"] /; Length[center] == 2 :=
Circle[center, radius];
insphere[center_, radius_, tangentpoints_]["Sphere"] /; Length[center] == 3 :=
Sphere[center, radius];
insphere[center_, radius_, tangentpoints_]["Sphere"] := {center, radius};
insphere[center_, radius_, tangentpoints_]["center"] := center;
insphere[center_, radius_, tangentpoints_]["radius"] := radius;
insphere[center_, radius_, tangentpoints_]["tangentpoints"] := tangentpoints;


The idea is that the condition that defines the insphere is that the perpendiculars dropped from the center to the faces are all equal. This leads to a system of linear equations that is easy for Solve to deal with.

For a given face $f$, let $\bf n$ be its unit normal vector pointing out of the tetrahedron. Let $(x, y, z)$ be the unknown center of the insphere and $r$ be its radius. Then we seek a condition that $(x, y, z) + r\,{\bf n}$ lies on the face $f$. Let ${\bf v}_1=(a_1, b_1, c_1)$ and ${\bf v}_2=(a_2, b_2, c_2)$ be the displacement vectors of two edges of $f$, and let ${\bf P}_0$ be one of the vertices of $f$. Then the unit outward normal $\bf n$ will be the normalization of $\pm({\bf v}_1 \times {\bf v}_2)$. The standard point-normal form for the equation of a plane containing a point ${\bf P}_0$ and perpendicular to a vector $\bf n$ is, for an arbitrary point ${\bf P}$ on the plane, $$({\bf P}-{\bf P}_0) \cdot {\bf n} = 0\,.$$ With our variables, in which ${\bf P}$ equals center + radius n and ${\bf P}_0$ equals First[f], this is equivalent to

(center + radius n - First[f]) . n == 0


To get the orientation of the normal $\bf n$ right, we take the three vertices of a face plus the fourth vertex of the tetrahedron. Subtract the first face vertex from the other three. The first two vectors represent the edges of the face coming out of the first vertex. The last vector points toward the interior side of the face, toward the fourth vertex, and its dot product with the outward normal will be negative. Scaling the cross product of the edge vectors by the negative of dot product of the cross product and the last vector, -#.Last@df, ensures the normal points outward. Normalizing the result yields our unit outward normal $\bf n$.

The function insphere[tet] returns a database in the form

insphere[center, radius, tangentpoints]


which may be accessed through functions like those at the end of the code above. Alternatively, to get a return value like the one specified by the OP, replace it in the code for insphere by

{center, radius, tangentpoints} /. First @ Solve[...]


Example

a = {0, 0, 0};
b = {1, 0, 0};
c = {2, 1, 0};
d = {1, 3/2, 1};
midpoints[tet_] := Mean /@ faces[tet];

tet0 = {a, b, c, d};
in0 = insphere[tet0];
Graphics3D[{
Opacity[0.4], Polygon[faces[tet0]],
MapThread[Arrow[{##}] &, Accumulate @ {midpoints[tet0], normals[tet0]}],
Red, in0["Sphere"],
Black, Point @ in0["tangentpoints"],
White, Opacity, Polygon[faces @ in0["tangentpoints"]]
}] Additional examples

2D

SeedRandom;
With[{in0 = insphere[#]},
Graphics[{Point @ in0["tangentpoints"], Opacity[0.4], Line @ faces @ #, in0["Sphere"]}]] &@
RandomReal[1, {3, 2}] 4D

SeedRandom;
proj0 = N @ Orthogonalize[IdentityMatrix ~Prepend~ {1, 1, 1, -2}][[2 ;; 4]];
proj[pts_?VectorQ] := proj0.pts;
proj[pts_?MatrixQ] := Transpose[proj0.Transpose @ pts];
sim = RandomReal[1, {5, 4}];
in0 = insphere[#] &@ sim;

in0["Sphere"]

Graphics3D[{Point @ proj @ in0["tangentpoints"], Opacity[0.5],
Sphere[proj @ in0["center"], in0["radius"]],
Opacity[0.25],
MapIndexed[{Hue[First @ #2 / 4], Polygon[faces @ proj @ #1]} &, faces @ sim]
}]


{{0.335297, 0.4023, 0.634872, 0.392616}, 0.059207} • Nice explanation. Better method than mine, which has quadratics. Sorry I can't give it another upvote. – Daniel Lichtblau Oct 13 '13 at 15:25
• Very impressive that this works in 4D! – Joseph O'Rourke Oct 14 '13 at 13:07
InSphere[a_, b_, c_, d_] :=
Module[{cent, radius, tangents},
{cent, radius} = InSphere0[b - a, c - a, d - a];
cent = cent + a;(*Shift center back*)
tangents =
Map[Project[cent, #] &, {{a, b, c}, {b, c, d}, {c, d, a}, {d, a, b}}];
Return[{cent, radius, tangents}]
];

(*InSphere0[] assumes one vertex @origin*)
InSphere0[a_, b_, c_] := Module[{V6, denom, numer, cent, radius},
V6 = a.(b\[Cross]c); (*6x Volume*)
denom =
Norm[b\[Cross]c] + Norm[c\[Cross]a] + Norm[a\[Cross]b] +
Norm[(b\[Cross]c) + (c\[Cross]a) + (a\[Cross]b)];
numer =
Norm[b\[Cross]c] a + Norm[c\[Cross]a] b + Norm[a\[Cross]b] c;
cent = numer/denom;
radius = V6/denom;
Return[N[{cent, radius}]]
];

(*Project p onto plane thru {a,b,c}*)
Project[p_, {a_, b_, c_}] := Module[{v},
v = Project0[p - a, {b - a, c - a}];
Return[v + a]
];

(*Project0[] assumes one point @origin*)
Project0[p_, {b_, c_}] := Module[{nvec, d},
nvec = Normalize[b\[Cross]c];
d = nvec.p;(*dist from p to plane along nvec*)
Return[p - d nvec]
];

• You'll get my +1 if you add example code for calling and displaying the results (as in your question). – bill s Oct 12 '13 at 13:39
• You got my vote for posting what is no doubt working code. For something I thought would be straightforward, this problem was far from. (As for getting other votes, what can I say? It's a tough crowd, and one's gotta deliver.) – Daniel Lichtblau Oct 12 '13 at 23:27

This may not be as neat as the other methods posted. About the only things I can say are that it is derived from basic principles, and it is fortunate that I had my hair buzzed rather short a few days ago.

planePoly[verts_List,
vars_List] := (vars - First[verts]).Cross @@ (First[verts] - #1 &) /@
Rest[verts]

insideTet[faces_, pts_, verts_, pt_] := Module[
{s1, s2, indx},
s1 = Sign[MapThread[#1 /. Thread[#2 -> pt] &, {faces, pts}]];
s2 = Sign[
MapThread[#1 /. Thread[#2 -> #3] &, {faces, pts, verts}]];
Apply[And, Thread[s1 == s2]]
]

inSphere[verts_, vars_] := Module[
{cen, r2, p, c, pts, triads = Subsets[verts, {3}], facepolys,
distpolys, perppolys, polys, sols, centers, mverts},
pts = Array[p, {4, 3}];
cen = Array[c, 3];
mverts = Map[First[Complement[verts, #]] &, triads];
facepolys = MapThread[planePoly, {triads, pts}];
distpolys = Map[(# - cen).(# - cen) - r2 &, pts];
perppolys =
MapThread[{(#1 - cen).(#2[] - #2[]), (#1 -
cen).(#2[] - #2[])} &, {pts, triads}];
polys = Expand[Flatten[{facepolys, distpolys, perppolys}]];
sols = NSolve[polys==0, Flatten[{pts, cen, r2}]];
centers = cen /. sols;
sols = Pick[sols,
Table[insideTet[facepolys, pts, mverts, centers[[j]]], {j,
Length[centers]}]];
{pts, cen, r2, With[{diff = vars - cen}, diff.diff - r2]} /.
First[sols]
]


We set up an example with four integer vertices chosen at random from the range {-5, 5}:

SeedRandom
verts = RandomInteger[{-5, 5}, {4, 3}]

Out= {{5, -5, 3}, {4, 4, 1}, {5, -3, -1}, {3, -3, 3}}


Now do the computation:

{tangentpoints, cen, rad2, eqn} = inSphere[verts, {x, y, z}]

Out= {{{4.77144649308, -2.69409447825,
2.04504506719}, {4.31499031868, -2.62640966025,
2.57785483539}, {3.769201698, -3.17219828093,
1.80599316585}, {3.63144307326, -2.76897176288,
1.73711385348}}, {4.17242821605, -2.76897176288,
2.00760642488}, 0.365831155905, -0.365831155905 + (-4.17242821605 +
x)^2 + (2.76897176288 + y)^2 + (-2.00760642488 + z)^2}


Here is how it all looks:

Graphics3D[{Sphere[cen, Sqrt[rad2]], Point[{verts}],
Map[Line, {Subsets[verts, {2}]}],
Map[{Opacity[.3], Polygon[#]} &, {Subsets[verts, {3}]}]},
Axes -> True] I remark that one can use Solve instead to get an exact solution. It will in general be long and probably not of great interest to look at.

• re hair: Pics or it didn't happen. – Yves Klett Oct 13 '13 at 13:45
• @Yves Klett I also have new glasses. You'd hardly recognize me. – Daniel Lichtblau Oct 13 '13 at 15:21

Recently, I was finally able to read the wonderful book Matrices and Graphs in Geometry by Miroslav Fiedler. One of the things I learned from that book is that one can use the Cayley-Menger matrix not only for determining the volume of a simplex, but also for determining its circumsphere and insphere. What follows is a Mathematica implementation of Fiedler's formulae:

CayleyMengerMatrix[pts_?MatrixQ] /; Subtract @@ Dimensions[pts] == 1 :=
With[{d = Length[pts] + 1},
SparseArray[{{j_ /; j > 1, 1} :> 1, {1, k_ /; k > 1} :> 1,
{j_, k_} /; j != k :>
SquaredEuclideanDistance[pts[[j - 1]], pts[[k - 1]]]},
{d, d}]]

CircumInSphere[pts_?MatrixQ] /; Subtract @@ Dimensions[pts] == 1 :=
Module[{cv, icm, iv, rc},
icm = -2 Inverse[CayleyMengerMatrix[pts]]; cv = icm[[1, 2 ;;]];
{rc, iv} = Through[{First, Rest}[Sqrt[Tr[icm, List]]]];
{{cv.pts/Total[cv], rc/2}, {iv.pts, 1}/Total[iv]}]


I might as well throw in the formula for the volume:

SimplexVolume[pts_?MatrixQ] /; Subtract @@ Dimensions[pts] == 1 :=
With[{d = Length[pts] - 1},
Sqrt[Abs[Det[CayleyMengerMatrix[pts]]/(d! (2 d)!!)]]]


Here's the tetrahedron from the OP:

pts = {{0, 0, 0}, {1, 0, 0}, {2, 1, 0}, {1, 3/2, 1}};
tet = GraphicsComplex[pts, Polygon[Partition[Range, 3, 1, 1]]];


Compute both the circumsphere and insphere:

{cs, is} = RootReduce[CircumInSphere[pts]]
{{{1/2, 3/2, -5/8}, Sqrt/8},
{{Root[9 - 114 #1 + 191 #1^2 - 102 #1^3 + 17 #1^4 &, 2],
Root[1 + 11 #1 - 42 #1^2 + 17 #1^3 + 17 #1^4 &, 3],
Root[1 - 16 #1 + 81 #1^2 - 136 #1^3 + 17 #1^4 &, 1]},
Root[1 - 16 #1 + 81 #1^2 - 136 #1^3 + 17 #1^4 &, 1]}}


Show both along with the tetrahedron:

GraphicsRow[{Graphics3D[{{Opacity[1/2], Sphere @@ cs}, tet}, Boxed -> False],
Graphics3D[{{Opacity[1/2], tet}, Sphere @@ is}, Boxed -> False]}] Of course, the whole setup also works for triangles:

tri = {{0, 0}, {3, 1}, {1, 2}};
Graphics[{FaceForm[None], EdgeForm[Black], Polygon[tri],
Disk @@@ CircumInSphere[tri]}] I might as well add that in version 10.2 you can use the built-in functions Insphere and Circumsphere to compute the insphere and circumsphere of a tetrahedron respectively. The output of both functions are Sphere objects which can be used in Graphics3D directly. So using the tetrahedron provided by the OP we can do:

pts = {{0, 0, 0}, {1, 0, 0}, {2, 1, 0}, {1, 3/2, 1}};

GraphicsRow[{Graphics3D[{{Opacity[0.5], Circumsphere[pts]},
Tetrahedron[pts]}, Boxed -> False],
Graphics3D[{{Opacity[0.5], Tetrahedron[pts]}, Insphere[pts]},
Boxed -> False]}] Here are the actual outputs:

Insphere[pts]
Circumsphere[pts]
` • Very useful to know this---Thanks! – Joseph O'Rourke Jul 29 '15 at 21:22