For BSP generation, How to intersect or locate a triangle with a plane defined by another triangle? I've hit a stumbling block in my project to draw the Utah teapot. I want to generate a binary space partition tree of a set of 3D triangles.
The decision step in the recursive tree-construction algorithm is determining whether one triangle is in front of or behind the plane defined by another triangle.
So, given two triangles, $T_0 = (P_1, P_2, P_3), T_1 = (P_4, P_5, P_6), P_1 = (x_1, y_1, z_1), P_2 = (x_2, y_2, z_2)$, ..., we take a normal vector of $T_0$, $N_0 = (P_2 - P_1) \times (P_3 - P_1)$.
And then I get confused. Help?
I think I need to describe the coordinates of $T_1$ in terms of a Plane [ a b c d ] + $N_0$. The magnitude of the resulting vector is the distance. Positive or negative determines whether it's in front of or behind the plane. But I'm fuzzy on step 2.
 A: Using determinants
You can use the sign of the determinant
$$\begin{vmatrix}
x_1 & x_2 & x_3 & x_4 \\
y_1 & y_2 & y_3 & y_4 \\
z_1 & z_2 & z_3 & z_4 \\
1 & 1 & 1 & 1
\end{vmatrix}$$
to decide on which side of the plane spanned by $P_1, P_2, P_3$ the last point $P_4$ lies. So take the three corners of one triangle, and perhaps the center of the other triangle, and compute the sign of that determinant to make your decisions.
You could even move the division from the computation of the center into the determinant like this:
$$\begin{vmatrix}
x_1 & x_2 & x_3 & x_4 + x_5 + x_6 \\
y_1 & y_2 & y_3 & y_4 + y_5 + y_6 \\
z_1 & z_2 & z_3 & z_4 + z_5 + z_6 \\
1 & 1 & 1 & 3
\end{vmatrix}$$
The last column is the center of the triangle $P_4,P_5,P_6$, scaled by 3. This will scale the value of the determinant by 3 as well, but won't affect its sign.
Normal form
And in case you have to compare multiple points against a given plane, you can do Laplace expansion along the fourth column and end up with 4 numbers, each computed as a $3\times 3$ determinant, which represent something like the normal form of the plane and can be used to perform the subsequent checks using only 4 multiplications per point. So you'd compute
\begin{align*}
a &= -\begin{vmatrix}
y_1 & y_2 & y_3 \\
z_1 & z_2 & z_3 \\
1 & 1 & 1
\end{vmatrix} &
b &= \begin{vmatrix}
x_1 & x_2 & x_3 \\
z_1 & z_2 & z_3 \\
1 & 1 & 1
\end{vmatrix} \\[2ex]
c &= -\begin{vmatrix}
x_1 & x_2 & x_3 \\
y_1 & y_2 & y_3 \\
1 & 1 & 1
\end{vmatrix} &
d &= \begin{vmatrix}
x_1 & x_2 & x_3 \\
y_1 & y_2 & y_3 \\
z_1 & z_2 & z_3 \\
\end{vmatrix}
\end{align*}
After this precomputation, you can use the sign of $ax+by+cz+d$ to decide on which side of the plane a point $(x,y,z)$ lies. If you want to use centroids without division, you'd multiply $d$ by $3$ and the others by the coordinate sums. The plane itself is described as
$$ \left\{(x,y,z)\in\mathbb R^3\;\middle\vert\;ax+by+cz+d=0\right\} $$
Notice that this is not quite the Hesse normal form of the plane, because the normal vector $(a,b,c)$ isn't normalized to unit length, and $d$ isn't neccessarily positive but instead indicates on which side of the plane the origin lies. So for computing distances you'd do some more computations, but to decide sidedness, these numbers $a$ through $d$ are all you need.
