Intersection between line and triangle in 3D

Question

I have a line and a triangle somewhere in 3D space. In other words, I have 3 points ([x,y,z] each) for the triangle, and two points (also [x,y,z]) for the line.

I need to figure out a way, hopefully using C++, to figure out if the line ever crosses the triangle. A line parallel to the triangle, and with more than one point in common, should be counted as "does not intersect".

I already made some code, but it doesn't work, and I always get false even when a visual representation clearly shows an intersection.

ofVec3f P1, P2;
P1 = ray.s;
P2 = ray.s + ray.t;

ofVec3f p1, p2, p3;
p1 = face.getVertex(0);
p2 = face.getVertex(1);
p3 = face.getVertex(2);

ofVec3f v1 = p1 - p2;
ofVec3f v2 = p3 - p2;

float a, b, c, d;

a = v1.y * v2.z - v1.z * v2.y;
b = -(v1.x * v2.z - v1.z * v2.x);
c = v1.x * v2.y - v1.y * v2.x;
d = -(a * p1.x + b * p1.y + c * p1.z);

ofVec3f O = P1;
ofVec3f V = P2 - P1;

float t;

t = -(a * O.x + b * O.y + c * O.z + d) / (a * V.x + b * V.y + c * V.z);

ofVec3f p = O + V * t;

float xmin = std::min(P1.x, P2.x);
float ymin = std::min(P1.y, P2.y);
float zmin = std::min(P1.z, P2.z);

float xmax = std::max(P1.x, P2.x);
float ymax = std::max(P1.y, P2.y);
float zmax = std::max(P1.z, P2.z);


if (inside(p, xmin, xmax, ymin, ymax, zmin, zmax)) {
    *result = p.length();
    return true;
}
return false;

And here is the definition of inside()

bool primitive3d::inside(ofVec3f p, float xmin, float xmax, float ymin, float ymax, float zmin, float zmax) const {
    if (p.x >= xmin && p.x <= xmax && p.y >= ymin && p.y <= ymax && p.z >= zmin && p.z <= zmax)
        return true;

    return false;
}

Answer 1

1) If you just want to know whether the line intersects the triangle (without needing the actual intersection point):

Let p1,p2,p3 denote your triangle

Pick two points q1,q2 on the line very far away in both directions.

Let SignedVolume(a,b,c,d) denote the signed volume of the tetrahedron a,b,c,d.

If SignedVolume(q1,p1,p2,p3) and SignedVolume(q2,p1,p2,p3) have different signs AND SignedVolume(q1,q2,p1,p2), SignedVolume(q1,q2,p2,p3) and SignedVolume(q1,q2,p3,p1) have the same sign, then there is an intersection.

SignedVolume(a,b,c,d) = (1.0/6.0)*dot(cross(b-a,c-a),d-a)

2) Now if you want the intersection, when the test in 1) passes

write the equation of the line in parametric form: p(t) = q1 + t*(q2-q1)

Write the equation of the plane: dot(p-p1,N) = 0 where

N = cross(p2-p1, p3-p1)

Inject p(t) into the equation of the plane: dot(q1 + t*(q2-q1) - p1, N) = 0

Expand: dot(q1-p1,N) + t dot(q2-q1,N) = 0

Deduce t = -dot(q1-p1,N)/dot(q2-q1,N)

The intersection point is q1 + t*(q2-q1)

3) A more efficient algorithm

We now study the algorithm in:

Möller and Trumbore, « Fast, Minimum Storage Ray-Triangle Intersection », Journal of Graphics Tools, vol. 2,‎ 1997, p. 21–28

(see also:)

https://en.wikipedia.org/wiki/M%C3%B6ller%E2%8%93Trumbore_intersection_algorithm

The algorithm is in the end simpler (less instructions than what we did in 1) and 2)), but sightly more complicated to understand. Let us derive it step by step.

Notation:

• O = origin of the ray,

• D = direction vector of the ray,

• A,B,C = vertices of the triangle

An arbitrary point P on the ray can be written as P = O + tD

An arbitrary point P on the triangle can be written as P = A + uE1 + vE2 where E1 = B-A and E2 = C-A, u>=0, v>=0 and (u+v)<=1

Writing both expressions of P gives:

O + tD = A + uE1 + vE2 

or:

uE1 + vE2 -tD = O-A

in matrix form:

            [u]
 [E1|E2|-D] [v] = O-A
            [t]

(where [E1|E2|-D] is the 3x3 matrix with E1,E2,-D as its columns)

Using Cramer's formula for the solution of:

   [a11 a12 a13][x1]   [b1]
   [a12 a22 a23][x2] = [b2]
   [a31 a32 a33][x3]   [b3]

gives:

       |b1 a12 a13|   |a11 a12 a13|
  x1 = |b2 a22 a23| / |a21 a22 a23|
       |b3 a32 a33|   |a31 a32 a33|

       |a11 b1 a13|   |a11 a12 a13|
  x2 = |a21 b2 a23| / |a21 a22 a23|
       |a31 b3 a33|   |a31 a32 a33|

       |a11 a12 b1|   |a11 a12 a13|
  x3 = |a21 a22 b2| / |a21 a22 a23|
       |a31 a32 b3|   |a31 a32 a33|

Now we get:

  u = (O-A,E2,-D) / (E1,E2,-D)
  v = (E1,O-A,-D) / (E1,E2,-D)
  t = (E1,E2,O-A) / (E1,E2,-D)

where (A,B,C) denotes the determinant of the 3x3 matrix with A,B,C as its column vectors.

Now we use the following identities:

  (A,B,C) = dot(A,cross(B,C))  (develop the determinant w.r.t. first column)

  (B,A,C) = -(A,B,C)           (swapping two vectors changes the sign)

  (B,C,A) =  (A,B,C)           (circular permutation does not change the sign)

Now we get:

u = -(E2,O-A,D)  / (D,E1,E2)
v =  (E1,O-A,D)  / (D,E1,E2)
t = -(O-A,E1,E2) / (D,E1,E2)  

Using:

N=cross(E1,E2);

AO = O-A; 

DAO = cross(D,AO)

We obtain finally the following code (here in GLSL, easy to translate to other languages):

bool intersect_triangle(
    in Ray R, in vec3 A, in vec3 B, in vec3 C, out float t, 
    out float u, out float v, out vec3 N
) { 
   vec3 E1 = B-A;
   vec3 E2 = C-A;
         N = cross(E1,E2);
   float det = -dot(R.Dir, N);
   float invdet = 1.0/det;
   vec3 AO  = R.Origin - A;
   vec3 DAO = cross(AO, R.Dir);
   u =  dot(E2,DAO) * invdet;
   v = -dot(E1,DAO) * invdet;
   t =  dot(AO,N)  * invdet; 
   return (det >= 1e-6 && t >= 0.0 && u >= 0.0 && v >= 0.0 && (u+v) <= 1.0);
}
 

When the function returns true, the intersection point is given by R.Origin + t * R.Dir. The barycentric coordinates of the intersection in the triangle are u, v, 1-u-v (useful for Gouraud shading or texture mapping). The nice thing is that you get them for free !

Note that the code is branchless. It is used by some of my shaders on ShaderToy

• https://www.shadertoy.com/view/tl3XRN
• https://www.shadertoy.com/view/3ltSzM