# 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. (11Oct13) 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
@RahulNarain: Re "no?": Yes. :-) But one must actually code this. Patience, patience! –  Joseph O'Rourke Oct 12 '13 at 1:05

[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 &)] :=
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 :=


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],
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

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},
s2 = Sign[
]

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 tet integer vertices chosen at random from 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. –  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
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). –  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