# 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. Commented Oct 11, 2013 at 13:30
• @YvesKlett, Graphics3D@{Opacity[0.5],Polygon[{a,b,c,d}~Subsets~{3}]} Commented Oct 11, 2013 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. Commented Oct 11, 2013 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?
– user484
Commented Oct 12, 2013 at 1:03
• The MO question where this was used: "Tetrahedron insphere iteration." Commented Aug 24, 2016 at 19:11

[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[
Function[{f, n}, (center + radius n  -  First[f]) . n == 0],
{faces[simplex], normals}],
]]

(* some interfaces to the insphere data *)
insphere[center_, radius_, tangentpoints_]["Sphere"] /; Length[center] == 2 :=
insphere[center_, radius_, tangentpoints_]["Sphere"] /; Length[center] == 3 :=
insphere[center_, radius_, tangentpoints_]["center"] := center;
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[1], Polygon[faces @ in0["tangentpoints"]]
}]


2D

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


4D

SeedRandom[3];
proj0 = N @ Orthogonalize[IdentityMatrix[4] ~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. Commented Oct 13, 2013 at 15:25
• Very impressive that this works in 4D! Commented Oct 14, 2013 at 13:07
InSphere[a_, b_, c_, d_] :=
{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}}];
];

(*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;
];

(*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). Commented Oct 12, 2013 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.) Commented Oct 12, 2013 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];
distpolys = Map[(# - cen).(# - cen) - r2 &, pts];
perppolys =
MapThread[{(#1 - cen).(#2[[1]] - #2[[2]]), (#1 -
cen).(#2[[1]] - #2[[3]])} &, {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[111111]
verts = RandomInteger[{-5, 5}, {4, 3}]

Out[322]= {{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[323]= {{{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. Commented Oct 13, 2013 at 13:45
• @Yves Klett I also have new glasses. You'd hardly recognize me. Commented Oct 13, 2013 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[4], 3, 1, 1]]];


Compute both the circumsphere and insphere:

{cs, is} = RootReduce[CircumInSphere[pts]]
{{{1/2, 3/2, -5/8}, Sqrt[185]/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! Commented Jul 29, 2015 at 21:22