Logo_iSurf_org

Subscribe
to eAlerts

for Geom Site
Updates

CLICK HERE

 

Intersections of Rays and Triangles (3D)

by Dan Sunday

 

 

The intersection of the most basic geometric primitives was presented in the Algorithm 5 about Intersections of Lines and Planes. 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.


Intersection of a Ray/Segment with a Plane

Assume we have a ray R (or segment S) from P0 to P1, and a plane P through V0 with normal n. The intersection of the parametric line L: P(r)=P0+r(P1-P0) and the plane P occurs at the point P(rI) with parameter value:

Eqn_line-plane2


When the denominator n_dot_(P1-P0)=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 P0 or P1 is in P ). Otherwise, when the denominator is nonzero and rI is a real number, then the ray R intersects the plane P only when rI.ge.0. A segment S intersects P only if rI.ge-0.le-1. In all algorithms, the additional test rI.le.1 is the only difference for a segment instead of a ray.

Intersection of a Ray/Segment with a Triangle

Consider a ray R (or a segment S) from P0 to P1, and a triangle T with vertices V0, V1 and V2. The triangle T lies in the plane P  through V0 with normal vector n=(V1-V0)-cross-(V2-V0). To get the intersection of R (or S) with T, one first determines the intersection of R (or S) and P . If it does not intersect, then it also does not intersect T and we are done. However, if they intersect in the point PI = P(rI), we need to determine if this point is inside 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 here an alternate method that also uses direct 3D computations to determine inclusion, avoiding the 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=V1-V0 and v=V2-V0 are edge vectors of T. Then, P=V(s,t) is in the triangle T when s.ge-0, t.ge-0, and s+t.le-1. So, given PI, one just has to find the (sI, tI) coordinate for it, and then check these inequalities to verify inclusion in T. Further, 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). And, the three vertices are given by: V0=V(0,0), V1=V(1,0), and V2=V(0,1).

Pic_line-triangle

To solve for sI and tI, we apply the method described in Algorithm 4 for Barycentric Coordinate Computation using the "generalized perp operator on P " that we defined there. Then, with w=PI-V0, which is a vector in P  (that is, n_dot_w=0), we solve the equation: w=su+tv for s and t. The final result is:

Eqn_param-Ps_inter

Eqn_param-Pt_inter

which has only 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: intersect3D_RayTriangle() ). Based on a count of the operations done up to the first rejection test, this algorithm is a bit less efficient than the MT algorithm, although we have not done any runtime performance comparisons. 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, which is needed to compute the line parameter rI. But, when the normal vectors have been precomputed and stored for all triangles in a scene (which is often the case), our algorithm would not have to compute this cross product at all. But, in this case, the MT algorithm would still compute two cross products, and be less efficient than our algorithm. So, the preferred algorithm depends on the application.

Intersection of a Triangle with a Plane

Consider a triangle T with vertices P0, P1 and P2 lying in a plane P1 with normal n1. Let P2 be a second plane through the point V0 with the normal vector n2. Unless they are parallel, the two planes P1 and P2 intersect in a line L, and when T intersects P2 it will be a segment contained in L. When T does not intersect P2 all three of its vertices must strictgly lie on the same side of the P2 plane. On the other hand, when T does intersect P2, one point of T must be on one side of (or on) P2 and the other two be on the other side of (or on) P2. We gave a test for which side of a plane a point is on (by using the signed distance from the point to the plane) in Algorithm 4 for Distance of a Point to a Plane. In fact, to determine the sideness of a point P, one only has to find the sign of n_dot_(P-V0)=0.

Suppose that P0 is on one side of P2 and that P1 and P2 are on the other side. Then the two segments P0P1 and P0P2 intersect P2 in two points I1 and I2 which are on the intersection line of P1 and P2. Then, the segment I1I2 is the intersection of triangle T and the plane P2.

Pic_triangle-plane

Intersection of a Triangle with a Triangle

Consider two triangles T1 and T2. They each lie in a plane, respectively P1 and P2, and their intersection must be on the line of intersection L for the two planes. Let the intersection of T1 and P2 be the segment S1 = I11I12, and the intersection of T2 and P1 be S2 = I21I22. If either S1 or S2 doesn’t exist (that is, one triangle does not intersect the plane of the other), then T1 and T2 do not intersect. Otherwise their intersection is equal to the intersection of the two segments S1 and S2 on the line L. This can be easily computed by projecting them onto an appropriate coordinate axis, and determining their intersection on it.

Pic_triangle-triangle


Implementations

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

// Copyright 2001 softSurfer, 2012 Dan Sunday
// 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)
 


// intersect3D_RayTriangle(): find the 3D intersection of a ray with a 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
intersect3D_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)

Tomas Moller & Eric Haines, "Intersection Test Methods" in Real-Time Rendering (3rd Edition) (2008)

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

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

J.P. Snyder and A.H. Barr, "Ray  Tracing Complex Models Containing Surface Tessellations", ACM Comp Graphics 21,  (1987)

 

© Copyright 2012 Dan Sunday, 2001 softSurfer