# Confusion about order of operations with point-in-tetrahedron formula

I am not a math student, but I am attempting to roll my own GJK-based hit detection function. It would seem that most of the Internet that I've searchedd through chooses to either ignore or obfuscate the test case after the triangle - namely, the tetrahedral case.

I have found the following two checks on http://en.wikipedia.org/wiki/Tetrahedron#Interior_point which appear to be straightforward - however I'm unclear on the terminology as well as the intended order of operations.

Here's the first inequality:

Let P be any interior point of a tetrahedron of volume V for which the vertices are A, B, C, and D, and for which the areas of the opposite faces are Fa , Fb , Fc , and Fd. Then[17]:p.62,#1609

PA {dot} F_a + PB {dot} F_b + PC {dot} F_c + PD {dot} F_d >= 9V.

Again - I am not a math major, but I read this to mean "P times (A dot the area of opposite side BCD) + P times (B dot the area of opposite side CDA)" and so on. Which I think is the result of the distributive property and could be rewritten as "P times ((A dot area BCD) + (B dot area CDA) + (C dot area DAB) + (D dot area ABC))". But I am not sure that's correct.

In either event, what I do not fully understand is: how can the area of a simplex be involved in a dot product calculation? Is this my misunderstanding of the notation? Is there a different order of operations intended here? So far I have been operating on the idea that dot product requires vectors, whereas I understood area to be an integer or a float. Perhaps there is a way to express area as a vector, and I just have no idea? (that feels counter-intuitive, but then again, I'm a bit out of my depth here)

Likewise, for the following interior point check:

For vertices A, B, C, and D, interior point P, and feet J, K, L, and M of the perpendiculars from P to the faces:

PA+PB+PC+PD >= 3(PJ+PK+PL+PM).

First of all - I am not certain at all what is meant by "feet J K L M". What is meant by feet? Is it just the cross product of the triangle?

Second, if I were to plain-English this inequality, I would read it as "(P times A) + (P times B) + (P times C) + (P times D) >= 3((P times J) +...)". Or again, if I recall my distributive property correctly (not saying I do, haha), "P(A+B+C+D) >= 3(P(J+K+L+M))". Which looks really pretty! If only I were sure of what J K L and M actually were.

To summarize my questions:

• in the first inequality, what is the order of operations? And: If I am supposed to get the dot product of a point and the area of the opposite face, how do I do this?
• in the second inequality, what is meant by the term feet?
• in both inequalities: am I correctly applying the distributive property, or do I need to go back to third grade?

Finally: if this is just an "XY Problem" and I'm chasing my tail here, I would love to know about that. But, I would also appreciate gaining a better understanding of the inequalities I've posted in addition to being pointed at the more correct solution for my issue.

Thank you very much.

• For the inequality $$PA \cdot F_a + PB \cdot F_b + PC \cdot F_c + PD \cdot F_c \ge 9 V$$ to make sense, $PA$ should be interpreted as length of line segment between $P$ and $A$. Both $PA$ and $F_a$ are real numbers and $\cdot$ is scalar multiplication. In fact, if we let $\vec{PA}$ be the vector from $P$ to $A$ and $\vec{F}_a$ be the outward pointing normal on the face opposite to $A$ with magnitude $F_a$. One can show $$\vec{PA} \cdot \vec{F}_a + \vec{PB} \cdot \vec{F}_b + \vec{PC} \cdot \vec{F}_c + \vec{PD} \cdot \vec{F}_d = -9V$$ where $\cdot$ now stands for vector dot product. – achille hui Jan 20 '15 at 7:57
• @achillehui Do you have a proof that the second inequality is true in general? The reference on Wikipedia points to a source that only states this inequality for tetrahedra that have equal-area faces (and doesn't give a proof of even that). – mollyerin Jan 20 '15 at 8:11
• Hmm.... I think I misinterpret the definition of $J,K,L,M$. They are the projections of $P$ to the four faces. e.g $J$ is the point lying on the plane containing the face $BCD$ such that $PJ$ is perpendicular to that plane. Once again, $PA$, $PJ$ are lengths of corresponding line segments. – achille hui Jan 20 '15 at 8:16
• @mollyerin I don't have a proof for the second inequality but I suspect it is true. – achille hui Jan 20 '15 at 8:18
• Here's the easy way to test if a point is in a tetrahedron: steve.hollasch.net/cgindex/geometry/ptintet.html – Rahul Jan 23 '15 at 5:02

## Inside / Outside a Linear Tetrahedron

This answer is a straightforward generalization of the Inside / Outside problem in de 2-D case, for a linear triangle. The clue is provided by a Finite Element approach.

Let's consider the simplest non-trivial finite element shape in 3-D, which is a linear tetrahedron.

Function behaviour is approximated inside such a tetrahedron by a linear interpolation between the function values at the vertices, also called nodal points. Let $\,T,$ be such a function, and $\,x,y,z\,$ coordinates, then: $$T = A.x + B.y + C.z + D$$ Where the constants A, B, C, D are yet to be determined. Substitute $x=x_k$ , $y=y_k$ , $z=z_k$ with $k=0,1,2,3$ . Start with: $$T_0 = A.x_0 + B.y_0 + C.z_0 + D$$ Clearly, the first of these equations can already be used to eliminate the constant $D$, once and forever: $$T - T_0 = A.(x - x_0) + B.(y - y_0) + C.(z - z_0)$$ Then the constants $A$ , $B$ , $C$ are determined by: $$\begin{array}{ll} T_1 - T_0 = A.(x_1 - x_0) + B.(y_1 - y_0) + C.(z_1 - z_0) \\ T_2 - T_0 = A.(x_2 - x_0) + B.(y_2 - y_0) + C.(z_2 - z_0) \\ T_3 - T_0 = A.(x_3 - x_0) + B.(y_3 - y_0) + C.(z_3 - z_0) \end{array}$$ Three equations with three unknowns. A solution can be found: $$\left[ \begin{array}{c} A \\ B \\ C \end{array} \right] = \left[ \begin{array}{ccc} x_1-x_0 & y_1-y_0 & z_1-z_0 \\ x_2-x_0 & y_2-y_0 & z_2-z_0 \\ x_3-x_0 & y_3-y_0 & z_3-z_0 \end{array} \right]^{-1} \left[ \begin{array}{c} T_1-T_0 \\ T_2-T_0 \\ T_3-T_0 \end{array} \right]$$ It is concluded that $A,B,C$ and hence $(T-T_0)$ must be a linear expression in the $(T_k-T_0)$: $$T - T_0 = \xi.(T_1 - T_0) + \eta.(T_2 - T_0) + \zeta.(T_3 - T_0)$$ $$= \left[ \begin{array}{ccc} \xi & \eta & \zeta \end{array} \right] \left[ \begin{array}{c} T_1-T_0 \\ T_2-T_0 \\ T_3-T_0 \end{array} \right]$$ $$= \left[ \begin{array}{ccc} \xi & \eta & \zeta \end{array} \right] \left[ \begin{array}{ccc} x_1-x_0 & y_1-y_0 & z_1-z_0 \\ x_2-x_0 & y_2-y_0 & z_2-z_0 \\ x_3-x_0 & y_3-y_0 & z_3-z_0 \end{array} \right] \left[ \begin{array}{c} A \\ B \\ C \end{array} \right]$$ $$= T - T_0 = \left[ \begin{array}{ccc} x-x_0 & y-y_0 & z-z_0 \end{array} \right] \left[ \begin{array}{c} A \\ B \\ C \end{array} \right]$$ Hence: $$\begin{array}{ll} x - x_0 = \xi .(x_1 - x_0) + \eta.(x_2 - x_0) + \zeta.(x_3 - x_0) \\ y - y_0 = \xi .(y_1 - y_0) + \eta.(y_2 - y_0) + \zeta.(y_3 - y_0) \\ z - z_0 = \xi .(z_1 - z_0) + \eta.(z_2 - z_0) + \zeta.(z_3 - z_0) \end{array}$$ But also: $$T - T_0 = \xi.(T_1 - T_0) + \eta.(T_2 - T_0) + \zeta.(T_3 - T_0)$$ Therefore the same expression holds for the function $\,T\,$ as well as for the coordinates $\,x,y,z$ . This is called an isoparametric transformation. It is remarked without proof that the local coordinates $\,\xi,\eta,\zeta\,$ within a tetrahedron can be interpreted as sub-volumes, spanned by the vectors $\,\vec{r}_k-\vec{r}_0\,$ and $\,\vec{r}-\vec{r}_0\,$ where $\,\vec{r}=(x,y,z)\,$ and $\,k=1,2,3$ .
Reconsider the expression: $$T - T_0 = \xi.(T_1 - T_0) + \eta.(T_2 - T_0) + \zeta.(T_3 - T_0)$$ Partial differentiation to $\xi$ , $\eta$ , $\zeta$ gives: $$\partial T / \partial \xi = T_1 - T_0 \quad ; \quad \partial T / \partial \eta = T_2 - T_0 \quad ; \quad \partial T / \partial \zeta = T_3 - T_0$$ Therefore: $$T = T(0) + \xi \frac{\partial T}{\partial \xi} + \eta \frac{\partial T}{\partial \eta} + \zeta \frac{\partial T}{\partial \zeta}$$ This is part of a Taylor series expansion around the node $(0)$. Such Taylor series expansions are very common in Finite Difference analysis. Now rewrite as follows: $$T = (1 - \xi - \eta - \zeta).T_0 + \xi.T_1 + \eta.T_2 + \zeta.T_3$$ Here the functions $\,(1-\xi-\eta-\zeta)\,,\,\xi\,,\,\eta\,,\,\zeta\,$ are called the shape functions of a Finite Element. Shape functions $\,N_k\,$ have the property that they are unity in one of the nodes (k), and zero in all other nodes. In our case: $$N_0 = 1-\xi-\eta-\zeta \quad ; \quad N_1 = \xi \quad ; \quad N_2 = \eta \quad ; \quad N_3 = \zeta$$ So we have two representations, which are allmost trivially equivalent: $$\begin{array}{ll} T = T_0 + \xi.(T_1 - T_0) + \eta.(T_2 - T_0) + \zeta.(T_3 - T_0) \quad & \mbox{: Finite Difference like} \\ T = (1 - \xi - \eta - \zeta).T_0 + \xi.T_1 + \eta.T_2 + \zeta.T_3 \quad & \mbox{: Finite Element like} \end{array}$$

Let's solve $\,\xi,\,\eta,\,\zeta$ from the equations, with Cramer's rule : $$\begin{array}{ll} x - x_0 = \xi .(x_1 - x_0) + \eta.(x_2 - x_0) + \zeta.(x_3 - x_0) \\ y - y_0 = \xi .(y_1 - y_0) + \eta.(y_2 - y_0) + \zeta.(y_3 - y_0) \\ z - z_0 = \xi .(z_1 - z_0) + \eta.(z_2 - z_0) + \zeta.(z_3 - z_0) \end{array}$$ $$\xi = \left| \begin{array}{ccc} x-x_0 & x_2-x_0 & x_3-x_0 \\ y-y_0 & y_2-y_0 & y_3-y_0 \\ z-z_0 & z_2-z_0 & z_3-z_0 \end{array} \right| / D \\ \eta = \left| \begin{array}{ccc} x_1-x_0 & x-x_0 & x_3-x_0 \\ y_1-y_0 & y-y_0 & y_3-y_0 \\ z_1-z_0 & z-z_0 & z_3-z_0 \end{array} \right| / D \\ \zeta = \left| \begin{array}{ccc} x_1-x_0 & x_2-x_0 & x-x_0 \\ y_1-y_0 & y_2-y_0 & y-y_0 \\ z_1-z_0 & z_2-z_0 & z-z_0 \end{array} \right| / D$$ Where: $$D = \left| \begin{array}{ccc} x_1-x_0 & x_2-x_0 & x_3-x_0 \\ y_1-y_0 & y_2-y_0 & y_3-y_0 \\ z_1-z_0 & z_2-z_0 & z_3-z_0 \end{array} \right|$$ The Inside / Outside function IO of the tetrahedron is defined as follows. $$\mbox{IO} =\mbox{min}(N_0,N_1,N_2,N_3) = \mbox{min}(\xi,\,\eta,\,\zeta,\,1-\xi-\eta-\zeta)$$ The maximum of the IO function is reached for $\,\xi = \eta = \zeta = 1-\xi-\eta-\zeta = 1/4$ , hence at the midpoint (barycenter) of the tetrahedron. If we draw straight lines from the the midpoint towards the vertices, and further, then the whole space is subdvided into four regions, one where IO $= \xi$ , one where IO $= \eta$ , one where IO $= \zeta$ , and one where IO $= 1-\xi-\eta-\zeta$ ; hope you get the 4-D picture. The IO function is positive for points $(x,y,z)$ inside the tetrahedron and negative for points $(x,y,z)$ outside the tetrahedron, and zero for points at the boundaries. Moreover, the value of the function is sort of a measure for the distance between $(x,y,z)$ and the barycenter.

I've taken notice of the fact that the OP has some programming experience. So here comes a (Delphi Pascal) program snippet that is supposed to do the mathematics:

program Stick;
type
matrix = array of array of double;
punt = record
x,y,z : double;
end;
punten = array of punt;
function det(mat : matrix) : double;
{
Determinant of 3 x 3 matrix
}
begin
det := mat[0,0]*mat[1,1]*mat[2,2]-mat[0,0]*mat[1,2]*mat[2,1]
-mat[1,0]*mat[0,1]*mat[2,2]+mat[1,0]*mat[0,2]*mat[2,1]
+mat[2,0]*mat[0,1]*mat[1,2]-mat[2,0]*mat[0,2]*mat[1,1] ;
end ;
function Identical(M : matrix) : matrix;
{
Make copy of matrix
}
var
i,j : integer;
S : matrix;
begin
SetLength(S,3,3);
for i := 0 to 2 do
for j := 0 to 2 do
S[i,j] := M[i,j];
Identical := S;
end;
procedure test;
var
T : array[0..3] of punt;
r : punt;
M,S : matrix;
i : integer;
D,xi,eta,zeta,IO : double;
begin
{ Vertices of Tetrahedron }
for i := 0 to 3 do
begin
T[i].x := Random;
T[i].y := Random;
T[i].z := Random;
end;
SetLength(M,3,3);
for i := 0 to 2 do
begin
M[0,i] := T[i+1].x-T[0].x;
M[1,i] := T[i+1].y-T[0].y;
M[2,i] := T[i+1].z-T[0].z;
end;
D := det(M);
SetLength(S,3,3);
{ Barycenter test }
r.x := (T[0].x+T1.x+T[2].x+T[3].x)/4;
r.y := (T[0].y+T1.y+T[2].y+T[3].y)/4;
r.z := (T[0].z+T1.z+T[2].z+T[3].z)/4;
{ Point in space }
r.x := Random;
r.y := Random;
r.z := Random;
S := Identical(M);
S[0,0] := r.x-T[0].x;
S[1,0] := r.y-T[0].y;
S[2,0] := r.z-T[0].z;
xi := det(S)/D;
S := Identical(M);
S[0,1] := r.x-T[0].x;
S[1,1] := r.y-T[0].y;
S[2,1] := r.z-T[0].z;
eta := det(S)/D;
S := Identical(M);
S[0,2] := r.x-T[0].x;
S[1,2] := r.y-T[0].y;
S[2,2] := r.z-T[0].z;
zeta := det(S)/D;
IO := 1-xi-eta-zeta;
if IO > xi then IO := xi;
if IO > eta then IO := eta;
if IO > zeta then IO := zeta;
Writeln(IO);
if IO > 0 then Readln; { Inside }
end;
begin
while true do
test;
end.

Output, paused at Inside hit:

-3.57217102952215E+0000
-3.29649515391294E-0001
-6.94579094637227E-0001
-6.02022821324288E-0001
-3.00028147738175E-0001
-4.38988084865881E-0001
-4.44247003351889E-0001
-1.10578952640585E+0001
-4.84149144259875E-0001
-8.38910075699883E-0001
-4.40603226584686E-0001
6.93752218750558E-0003

Note. This approach is essentially not much different from establishing whether a point is inside / outside a sphere with equation $\,R^2-(x^2+y^2+z^2)=0$ .