Eric and Tomas are the authors of Real-Time Rendering (A.K. Peters, 1999). They both can be contacted at http:// www.realtimerendering.com/.
An important and commonly performed operation is testing whether a point is inside a triangle. In this article, we will discuss 2D and 3D cases of this operation and present methods that also give a measure of the point's location relative to the triangles edges and vertices.
2D Point-Triangle Test
A classic method of testing whether a point is inside a triangle, defined by three points, p0, p1, and p2, is to send an imaginary ray from the point outwards to infinity. If this test ray crosses an odd number of edges, the point is inside the triangle. This test also works for polygons or any closed-curve figure. However, since this algorithm has been examined many times before (see "Point in Polygon Strategies," by Eric Haines in Graphics Gems IV, edited by Paul S. Heckbert, AP Professional, 1994, http://www.acm.org/tog/GraphicsGems/; and Real-Time Rendering, by Tomas Moller and Eric Haines, A.K. Peters, 1999), we will focus on a comparable method that has some additional useful features.
The algorithm we'll discuss here is based on using barycentric coordinates. A point, (u,v), on a triangle is given by the explicit formula i(u,v)=(1-u-v)p0+ up1+vp2, where (u,v) are the barycentric coordinates, which must fulfill u0, v0, and u+v1. That is, u or v is the amount by which to weigh the corresponding vertex's contribution to a particular location, with (1-u-v) being the weight for the third vertex. Figure 1 shows the coordinate system in terms of u and v. Note that (u,v) can be used for texture mapping, normal interpolation, color interpolation, and so on.
To test whether a point is inside a triangle is straightforward. If the point fulfills the conditions just given, then it must be inside the triangle. But how do you determine u and v? The answer is to transform the test point into a coordinate system based on the axes defined by the u and v values. Example 1 is pseudocode for testing some point i for inclusion in a triangle.
For simplicity of presentation, intermediate results are all computed at the start of Example 1, in steps 1-3. In the actual code, these are computed only as needed. Step 4 tests whether p0 and p1 have the same x-coordinate (a possibly common case). Step 5 checks if all three x-coordinates are the same; if so, the triangle has no area. Steps 7 and 10 check if u and v are in bounds. The u=1 test in step 7 is not strictly necessary, but often quickly helps classify points as being outside the triangle with little additional work. Our goal is to exit the routine as soon as possible to avoid additional computations. Steps 12 through 17 are for computing the common case, and involve a little more computation. Step 18 checks the third edge of the triangle against the test point. Steps 5, 8, and 13 test if the triangle has zero area; testing for closeness to zero with some very small number may be preferable for robustness.
For a full derivation of the math involved in deriving this algorithm, see Didier Badouel's presentation, "An Efficient Ray-Polygon Intersection" and Graphics Gems, edited by Andrew S. Glassner, Academic Press, 1990, http:// www.acm.org/ tog/GraphicsGems/. The code presented here is a more efficient version of Badouel's.
3D Ray-Triangle Test
A common method of testing whether a ray intersects a triangle is to intersect the ray with the plane the triangle defines, then determine whether the point of intersection within the plane is inside the triangle (see "Essential Ray Tracing Algorithms," by Eric Haines in An Introduction to Ray Tracing, edited by Andrew Glassner, Academic Press, 1989, http:// www.education.siggraph.org/materials/ HyperGraph/raytrace/rtinter0.htm). This second step is done by dropping one of the three coordinates of the vertices and point involved, thereby turning this into a 2D problem, which we solved in the previous section. (This algorithm was first presented by Moller and Trumbore in "Fast, Minimum Storage Ray-Triangle Intersection," Journal of Graphics Tools, 2(1), 1997.) There is a more elegant solution, however, which has the added advantage of not needing to precompute and store the plane's normal. In addition, barycentric coordinates are generated for the intersection point.
To compute the intersection between a ray and a triangle in three dimensions using this newer method, we need to review a little linear algebra. Assume that we have three vectors, m0, m1, and m2, that make up a basis (a frame of reference). To transform the vectors of this basis so that they become mutually perpendicular and align with the x-, y- and z-axes, we first form a matrix M=(m0 m1 m2), where the basis vectors thus form columns in the matrix. Then the inverse of this matrix can be used to transform as we wanted. Essentially, this is just how the multiplication of a matrix by its inverse (where the result of the multiplication is equal to the identity matrix) is interpreted geometrically. A fine thing about this property is that it works as long as the vectors do not all lie in the same plane. The basis transform is illustrated in Figure 2.
Now to the problem that we want to solve. We have a triangle defined by three points: p0, p1, and p2 (in counter-clockwise order if seen from the front), and we have a ray r(t)=o+td, where o is the origin of the ray, d is the normalized direction vector of the ray, and t is the parameter of the ray equation (that is, the distance along the ray). This situation is illustrated in the upper left of Figure 3. The task is to compute whether the ray intersects the triangle, and if so, also compute the distance to the intersection point and the barycentric coordinates of it.
To solve this problem we first translate the triangle vertices and the ray by -p0 so that the first triangle vertex (p0) is at the origin. This is depicted in the upper right of Figure 3. Now we form a basis from -d, p1-p0, and p2-p0, and transform the translated origin, o-p0, with the inverse of that matrix. This is shown in the lower left of Figure 3 and in Example 2.
So now we only have to interpret what the resulting vector w is. Looking at the lower right of Figure 3, this is quite simple. The x-component of w is the positive distance from the ray origin to the front facing triangle, which means that it is actually the value of the t parameter of the ray equation. The y- and the z-components of w are the barycentric coordinates for the intersection point; that is, it must hold that wy0, wz0, and wy+wz1 if the intersection point is to be inside the triangle. wy and wz are denoted here instead as u and v. This means that w=(t,u,v)T, where t is the distance to the intersection point, and (u,v) are the barycentric coordinates of the intersection.
At this point you might argue that this is all very well in theory, but solving Example 2 involves the inversion of a 3X3 matrix, which is not considered a fast operation. However, it turns out that there are many factors that can be reused if the matrix is solved with Cramer's rule for matrix inversion. If we denote e1=p1-p0, e2=p2-p0, and s=o-p0, then the solution to Example 2 is in Example 3, where det() is the determinant of three vectors. Knowing that det(a,b,c)=(aX b)c=-(aX c)b-(cX b)a, Example 3 is restructured as in Example 4, where r=dX e2 and q=sX e1. This last expression is what we need in order to make an efficient implementation of the ray-triangle intersection test. Recall that it is important to jump out of the code as soon as we can find out that we do not have an intersection. So the main idea of the code is to compute a value, check if it is within its valid bounds, continue if it is, and otherwise exit.
Before showing the pseudocode, we will present a few more optimizations. First, if we care only about intersecting front-facing triangles, we can test if the determinant is less than a very small number, . If it is, then exit, since a negative determinant means that the basis is negatively oriented, which in turn means that the ray hits the back face of the triangle. On some modern CPU architectures, division is done in parallel with other operations, so dividing early and using the result much later can be advantageous.
When we care about both front- and back-facing triangles, we can do the division quite early in the code, but use its result as late as possible -- hoping that the division is computed by then; see Example 5.
Lines 1 and 2 compute two of the edges of the triangle, while lines 3 and 4 compute the determinant a. Line 5 computes the difference between the ray origin and the first triangle vertex. The inverse of the determinant f is computed at line 6. This division usually takes many cycles, but because we only need its result at lines 23-25 (if we reach that point in the code), we may hope that it is finished or at least close to finished at that point. This turned out to be about 5-10 percent faster than the original implementation, but this depends on the CPU architecture and on the percentage of intersections. Line 7 computes q, which later is used to compute v. What follows at line 8 is a test on the determinant, and if that branch is taken the triangle is front facing. The else branch is taken for a negative determinant (back-facing triangle), and thus we have changed from less than to greater than and vice versa in the inner if statements. At line 22 we know that the determinant was close to zero, and thus the triangle was close to or parallel to the ray direction. If we have not exited from the code, lines 23-25 compute the final result (t,u,v). All that is needed to get the values of u and v at this point is to multiply with the inverse of the determinant.
The source code is available electronically; see "Resource Center," page 5. Some lines can be swapped without altering the final result, but your CPU architecture may favor a certain order.
Further Reading and Resources
In Real-Time Rendering (A.K. Peters, 1999), we describe related tests in more detail. The book covers a wide variety of other intersection tests and methods; see http://realtimerendering.com/. For code for other intersection tests, we recommend visiting the Graphics Gems web site at http://www.acm.org/tog/GraphicsGems/ and Magic software at http://www.magic- software.com/.