The intersection of the most basic geometric primitives was presented in the Algorithm 5 about Intersections of Lines, Segments and Planes (2D and 3D).  We will now extend those algorithms to include 3D triangles which are common elements of 3D surface and polyhedron models.  We only consider transversal intersections where the two intersecting objects do not lie in the same plane.  Ray and triangle intersection computation is perhaps the most frequent nontrivial operation in computer graphics rendering using ray tracing.  Because of its importance, there are several published algorithms for this problem (see: [Badouel, 1990], [Moller & Trumbore, 1997], [O'Rourke, 1998], [Moller & Haines, 1999]).  We present an improvement of these algorithms for ray (or segment) and triangle intersection.  We also give algorithms for triangle-plane and triangle-triangle intersection.

Objects are represented by vertices: a segment by it's endpoints, a ray as an extended segment, and a triangle by its three vertices.  Lines will be given the parametric equation: P(r) = P₀ + r(P₁-P₀), which is also used for the ray starting at P₀ going in the direction (P₁-P₀), and for the segment from P₀ to P₁.   A plane will be given by a point V₀ on it and a normal vector n.

Ray/Segment-Plane Intersection

This topic was treated in Algorithm 5 (see: Line-Plane Intersection).  We recall the result here for a a ray R (or segment S) from P₀ to P₁, and a plane p through V₀ with normal n.  The intersection of the parametric line L: P(r) = P₀ + r(P₁-P₀) and the plane p occurs at the point P(r_I) with parameter value:

When the denominator n·(P₁-P₀) = 0, the line L is parallel to the plane p, and thus either does not intersect it or else lies completely in the plane (whenever either P₀ or P₁ is in p).  Otherwise, when the denominator is nonzero and r_I is a real number, then the ray R intersects the plane p only when r_I>=0.  A segment S intersects p only if 0<=r_I<=1.  In all algorithms, the additional test r_I<=1 is the only difference for a segment instead of a ray.

Ray/Segment-Triangle Intersection

Consider a ray R (or a segment S) from P₀ to P₁, and a triangle T with vertices V₀, V₁ and V₂.  The triangle T lies in the plane p through V₀ with normal vector n = (V₁-V₀)×(V₂-V₀).  To get the intersection of R (or S) and T, one first determines the intersection of R (or S) and p.  If they do not intersect, then the ray (or segment) also does not intersect T and we are done.  However, if they intersect in the point P_I=P(r_I), we need to determine if this point is in the triangle T for there to be a valid intersection.

There are a number of ways to test for the inclusion of a point inside a 3D planar triangle.  The algorithms given by [Badouel, 1990] and [O'Rourke, 1998] project the point and triangle onto a 2D coordinate plane where inclusion is tested.  To implement these algorithms, one must choose a projection coordinate plane which avoids a degenerate projection.  This is done by excluding the coordinate which has the largest component in the plane normal vector n [Synder & Barr, 1987].  The intent is to reduce the 3D problem to a simpler 2D problem which has an efficient solution.  However, there is a small overhead involved in selecting and applying the projection function.  The algorithm of [Moller-Trumbore, 1997] (MT) does not project into 2D, and finds a solution using direct 3D computations.  Testing with some complex models shows that the MT algorithm is faster than the one by Badouel.

We present an alternate method which also uses direct 3D computations to determine inclusion, avoiding projection onto a 2D coordinate plane.  As a result, the code is cleaner and more compact.  Like [Moller-Trumbore, 1997], we use the parametric equation of p relative to T, but derive a different method of solution which computes the parametric coordinates of the intersection point in the plane.  The parametric plane equation (see: Planes in Algorithm 4) is given by:

where s and t are real numbers, and u = (V₁-V₀) and v = (V₂-V₀) are edge vectors of T.  A point P = V(s,t) is in the triangle T when s>=0, t>=0, and s+t<=1.  So, given P_I, one just has to find the (s_I,t_I)-coordinate for it, and then check these inequalities to verify inclusion in T.  Further, a point P = V(s,t) is on an edge of T if one of the conditions s=0, t=0, or s+t=1 is true (each condition corresponds to one edge).  Also, the three vertices are given by: V₀=V(0,0), V₁=V(1,0), and V₂=V(0,1).

To solve for s_I and t_I, we define a 3D generalization of Hill's "perp-dot" product [Hill, 1994].  For an embedded 2D plane p with a normal vector n, and any vector a in the plane (that is, a satisfies n·a = 0), define the "generalized perp operator on p" by: a^{^} = n×a.  Then, a^{^} is another vector in the plane p (since n·a^{^} = 0), and it is perpendicular to a (since a·a^{^} = 0).  Further, this perp operator is linear on vectors in p; that is, (Aa + Bb)^{^} = Aa^{^} + Bb^{^} where a and b are vectors in p, and A and B are scalars.  Note that if p is the 2D xy-plane (z=0) with n = (0,0,1), then this is exactly the same 2D perp operator given by [Hill, 1994].

We can now use this generalized perp operator to solve the plane's parametric equation for the intersection point P_I.  Put w = (P_I-V₀) which is another vector in p.  Then, we want to solve the equation: w = su + tv for s and t.  Take the dot product of both sides with v^{^} to get: w · v^{^} = su · v^{^} + tv · v^{^} = su · v^{^}, and solve for s_I.  Similarly, taking the dot product with u^{^}, we get: w · u^{^} = su · u^{^} + tv · u^{^} = tv · u^{^}, and solve for t_I.  We get:

The denominators are nonzero whenever the triangle T is nondegenerate (that is, has a nonzero area).  When T is degenerate, it is either a segment or a point, and in either case does not uniquely define a plane and the computed normal vector is (0,0,0).  So, this case must be treated separately as a 3D segment-segment intersection which was treated in Algorithm 5 (see: Lines and Segments).

Altogether we have used 3 cross products which is a lot of computation.  However, for any three vectors a b c, the following identity holds for left association of the cross product:  (a × b) × c = (a · c) b - (b · c) a.  [Note: the cross product is not associative, and so there is a different (but similar) identity for right association].  Applying the left association identity results in the simplifications:

And we can now compute s_I and t_I using only dot products as:

with 5 distinct dot products.  We have arranged terms so that the two denominators are the same and only need to be calculated once.

This solution yields a straightforward ray/segment-triangle intersection algorithm (see our implementation: intersect_RayTriangle() ).  Based on a count of the operations done up to the first reject test, this algorithm is not as efficient as the [Moller-Trumbore, 1997] (MT) algorithm, although we have not yet done any runtime performance tests.  However, the MT algorithm uses two cross products whereas our algorithm uses only one, and the one we use computes the normal vector of the triangle's plane.  Thus, when the normal vectors have been precomputed and stored for all triangles in a scene (which is often the case), our algorithm would not compute a cross product at all, making it even more efficient.  On the other hand, in this case, the MT algorithm would still compute two cross products, and be less efficient than our algorithm.

Triangle-Plane Intersection

Consider a triangle T with vertices P₀, P₁ and P₂ lying in a plane p₁ with normal n₁.  Let p₂ be a second plane through the point V₀ with the normal vector n₂.  Unless they are parallel, the two planes p₁ and p₂ intersect in a line L, and when T intersects p₂ it will be a segment contained in L.  When T does not intersect p₂ all three of its vertices lie on the same side of the p₂ plane.  Otherwise, when T does intersect p₂, one point of T must be on one side of p₂ and the other two be on the other side.  A test for which side of a plane a point is on (by using the signed distance from the point to the plane) was given in Algorithm 4 (see: Distance of a Point to a Plane). 

Suppose that P₀ is on one side of p₂ and that P₁ and P₂ are on the other side.  Then the two segments P₀P₁ and P₀P₂ intersect p₂ in two points I₁ and I₂ which are on the intersection line of p₁ and p₂.  The segment I₁I₂ is the intersection of triangle T and the plane p₂.

Triangle-Triangle Intersection

Consider two triangles T₁ and T₂.  They each lie a plane, p₁ and p₂, and their intersection must be on the line of intersection L for the two planes.  Let the intersection of T₁ and p₂ be S₁ = I₁₁I₁₂, and the intersection of T₂ and p₁ be S₂ = I₂₁I₂₂.  If either S₁ or S₂ is empty (that is, one triangle does not intersect the plane of the other), then T₁ and T₂ do not intersect.  Otherwise their intersection is equal to the intersection of the two segments S₁ and S₂ on the line L.  This can be easily computed by projecting them onto an appropriate coordinate axis, and determining their intersection on it.

Implementations

Here are some sample "C++" implementations of these algorithms. 

// Copyright 2001, softSurfer (www.softsurfer.com)
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.

// Assume that classes are already given for the objects:
//    Point and Vector with
//        coordinates {float x, y, z;}
//        operators for:
//            == to test equality
//            != to test inequality
//            (Vector)0 = (0,0,0)         (null vector)
//            Point  = Point ± Vector
//            Vector = Point - Point
//            Vector = Scalar * Vector    (scalar product)
//            Vector = Vector * Vector    (cross product)
//    Line and Ray and Segment with defining points {Point P0, P1;}
//        (a Line is infinite, Rays and Segments start at P0)
//        (a Ray extends beyond P1, but a Segment ends at P1)
//    Plane with a point and a normal {Point V0; Vector n;}
//    Triangle with defining vertices {Point V0, V1, V2;}
//    Polyline and Polygon with n vertices {int n; Point *V;}
//        (a Polygon has V[n]=V[0])
//===================================================================

#define SMALL_NUM  0.00000001 // anything that avoids division overflow
// dot product (3D) which allows vector operations in arguments
#define dot(u,v)   ((u).x * (v).x + (u).y * (v).y + (u).z * (v).z)

// intersect_RayTriangle(): intersect a ray with a 3D triangle
//    Input:  a ray R, and a triangle T
//    Output: *I = intersection point (when it exists)
//    Return: -1 = triangle is degenerate (a segment or point)
//             0 = disjoint (no intersect)
//             1 = intersect in unique point I1
//             2 = are in the same plane
int
intersect_RayTriangle( Ray R, Triangle T, Point* I )
{
    Vector    u, v, n;             // triangle vectors
    Vector    dir, w0, w;          // ray vectors
    float     r, a, b;             // params to calc ray-plane intersect

    // get triangle edge vectors and plane normal
    u = T.V1 - T.V0;
    v = T.V2 - T.V0;
    n = u * v;             // cross product
    if (n == (Vector)0)            // triangle is degenerate
        return -1;                 // do not deal with this case

    dir = R.P1 - R.P0;             // ray direction vector
    w0 = R.P0 - T.V0;
    a = -dot(n,w0);
    b = dot(n,dir);
    if (fabs(b) < SMALL_NUM) {     // ray is parallel to triangle plane
        if (a == 0)                // ray lies in triangle plane
            return 2;
        else return 0;             // ray disjoint from plane
    }

    // get intersect point of ray with triangle plane
    r = a / b;
    if (r < 0.0)                   // ray goes away from triangle
        return 0;                  // => no intersect
    // for a segment, also test if (r > 1.0) => no intersect

    *I = R.P0 + r * dir;           // intersect point of ray and plane

    // is I inside T?
    float    uu, uv, vv, wu, wv, D;
    uu = dot(u,u);
    uv = dot(u,v);
    vv = dot(v,v);
    w = *I - T.V0;
    wu = dot(w,u);
    wv = dot(w,v);
    D = uv * uv - uu * vv;

    // get and test parametric coords
    float s, t;
    s = (uv * wv - vv * wu) / D;
    if (s < 0.0 || s > 1.0)        // I is outside T
        return 0;
    t = (uv * wu - uu * wv) / D;
    if (t < 0.0 || (s + t) > 1.0)  // I is outside T
        return 0;

    return 1;                      // I is in T
}

References

Didier Badouel, "An Efficient Ray-Polygon Intersection" in Graphics Gems (1990)

Francis Hill, "The Pleasures of 'Perp Dot' Products" in Graphics Gems IV (1994)
[Note: the first critical definition has a typo, and should be: a^{^} = (-a_y, a_x).]

Tomas Moller & Eric Haines, Real-Time Rendering, Chapter 10 "Intersection Test Methods" (1999)

Tomas Moller & Ben Trumbore, "Fast, Minimum Storage Ray-Triangle Intersection",  J. Graphics Tools 2(1), 21-28 (1997)

Joseph O'Rourke, Computational Geometry in C (2nd Edition), Section 7.3 "Segment-Triangle Intersection" (1998)