Csdrhrt

I'm doing some 3D graphics stuff and in order to add some simple lighting, the program needs to determine whether or not a point is below a triangle. (For now just the plane defined by the triangle. I'll tackle the point-in-triangle problem later.) Y is up. Thanks!

HINT

We can proceed as follow

• find the plane containing the triangle


• find the intersection point $Q$ between the plane and the vertical line through the point $P$
  


• calculate $z_P-z_Q$
  

Suppose the triangle is $ABC$ and the point is $P$.

Calculate the vectors $B-A$, $C-A$ and $P-A$. Then the sign of the vector triple product
$$
(P-A) cdot ((B-A) times (C-A))
$$
determines which side of the plane of the triangle contains $P$.

That triple product is just a determinant, which should make the coding easy.

Write the equation of the plane of the triangle as $$z=ax+by+d$$

If your point is $p= (x_0, y_0,z_0)$ compare the $z=ax_0+by_0+d$ with $z_0$.

If $z_0 > ax_0+by_0+d$ then $p$ is above otherwise it is either below or on the plane ( in case of equality).

Although others have already answered the stated question, I'd like to give some advice to those doing similar to what OP is doing.

1. 
   

   Precalculate the unit normal vector to the triangle, and store it.

   
   
   

   If the triangle vertices are $vec{v}_0$, $vec{v}_1$, and $vec{v}_2$, then the triangle normal is usually defined as
   $$vec{n} = bigl(vec{v}_1 - vec{v}_0bigr)timesbigl(vec{v}_2 - vec{v}_0bigr)$$
   where $times$ denotes vector cross product. This is scaled to unit length,
   $$hat{n} = frac{vec{n}}{leftlVertvec{n}rightrVert} = frac{vec{n}}{sqrt{vec{n} cdot vec{n}}}$$

   
   
   

   This is best precalculated and stored, because the square root operation tends to be slow, and it is just one vector per triangle.
    

   
   


2. 
   

   A point is "above" the triangle if it is in the same halfspace as the triange normal vector, with respect to the plane formed by the triangle.

   
   
   

   (If you extend the triangle plane to infinity, it splits space into two halves. If you look at the triangle sideways, the point is "above" the triangle if it is in the same side as the normal is; and "below" if it is on the other side. If the point is on the triangle plane exactly, I recommend you use a clear and unambigous term for it; I use "exactly on the triangle" myself. I'm sure mathematicians have more precise terms, but for those of us who only use math as a tool for solving problems, making sure the terms are unambiguous suffices.)

   
   
   

   This means that point $vec{p}$ has signed distance
   $$d = bigl(vec{p} - vec{v}_0bigr) cdot hat{n}$$
   from the plane of the triangle. (Because all three triangle vertices are on the plane, you can use $vec{v}_1$ or $vec{v}_2$ in above; you'll still get the exact same result, within rounding error.)

   
   
   

   If $hat{n}$ is not an unit vector, i.e. $leftlVerthat{n}rightrVertne 1$, the distance is in units of the normal vector length.
    

   
   



3. When you know the signed distance $d$ of the point $vec{p}$ from the plane, you can trivially project the point to the plane:
   $$vec{p}_text{plane} = vec{p} - d hat{n}$$
   If $hat{n}$ is not an unit vector, i.e. $leftlVerthat{n}rightrVertne1$, the right side is actually $vec{p} - hat{n}(d / lVerthat{n}rVert)$, and you need a "slow" square root operation here. Because you probably test many points per each triangle, it is well worth precalculating and storing the $hat{n}$ beforehand.
    



4. 
   

   To determine if $vec{p}_text{plane}$ is within the triangle, find out the barycentric coordinates $(u, v)$ corresponding to $vec{p}_text{plane}$ with respect to that triangle.

   
   
   

   To do this, check which axes have the smallest values in magnitude in $hat{n}$. Their order does not matter. (This essentially picks the axis along which to look at the triangle, with the triangle being as large as possible. We do this to minimize numerical errors.) There are three possibilities: $x$ and $y$; $x$ and $z$; and $y$ and $z$. The solution is obtained using the coordinates on those two axes.

   
   
   

   The solution for $vec{p}_text{plane} = (x, y, z)$ with $hat{n}$ having $x$ and $y$ components closest to zero, is
   $$begin{cases}
   displaystyle u = frac{ (x - x_0)(y_2 - y_0) - (y - y_0)(x_2 - x_0)}{ (x - x_2)(y_1 - y_0) - (y - y_2)(x_1 - x_0) } \
   displaystyle v = frac{ (y - y_0)(x_1 - x_0) - (x - x_0)(y_1 - y_0)}{ x_0 (y_1 - y_2) + x_1 (y_2 - y_0) + x_2 (y_0 - y_1) } \
   end{cases}$$
   The solution for the other axis pairs you'll find just by permuting the coordinates ($y$, $y_0$, $y_1$, $y_2$ to $z$, $z_0$, $z_1$, $z_2$ for $x$ and $z$ axes; and also $x$, $x_0$, $x_1$, $x_2$ to $y$, $y_0$, $y_1$, $y_2$ for $y$ and $z$).

   
   
   

   Also note that you'll want to calculate the divisor for $u$ first. It will be zero at $vec{v}_2$, where $u = 0$, so when the divisor is very close to zero, just use $u = 0$.

   
   
   

   In barycentric coordinates, $vec{v}_0$ is at $(0, 0)$; $vec{v}_1$ is at $(1,0)$, and $vec{v}_2$ is at $(0,1)$. The entire triangle is $0 le u, v, u + v le 1$. The edges are at $u = 0$, $v = 0$, and $u + v = 1$.

   
   
   

   In practice, you define the triangle inclusion test as
   $$Bigl(u ge -epsilonBigr) text{ and }
   Bigl(v ge -epsilonBigr) text{ and }
   Bigl(u le 1+epsilonBigr) text{ and }
   Bigl(v le 1+epsilonBigr) text{ and }
   Bigl(u + v le 1+epsilonBigr)$$
   where $epsilon$ represents the precision of your numerical coordinates. It is the largest positive number you consider equal to zero. (Similarly, if $D$ is the divisor for $u$ earlier, if $-epsilon le D le +epsilon$, use $u = 0$.)

   
   
   

   To linearly interpolate any value $c$ within the triangle, with $c = c_0$ at $vec{v}_0$, $c_1$ at $vec{v}_1$, and $c_2$ at $vec{v}_2$, you can use
   $$c(u,v) = (1 - v)bigl((1-u) c_0 + u c_1bigr) + v c_2$$

   
   

All together, this means that you need only 15 (scalar) multiplications and less than thirty additions or subtractions, to find out if a point is within the triangle, or within a triangular prism extending from the triangle along its normal, including the barycentric coordinates needed for texture or color interpolation. (Color interpolation itself is just four multiplications per component, so typically 12 additional multiplications.) Nifty, eh?