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

@BrunoLevi: your algorithm does not seem to work, see the following python implementation:

def intersect_line_triangle(q1,q2,p1,p2,p3):
    def signed_tetra_volume(a,b,c,d):
        return np.sign(np.dot(np.cross(b-a,c-a),d-a)/6.0)

    s1 = signed_tetra_volume(q1,p1,p2,p3)
    s2 = signed_tetra_volume(q2,p1,p2,p3)

    if s1 != s2:
        s3 = signed_tetra_volume(q1,q2,p1,p2)
        s4 = signed_tetra_volume(q1,q2,p2,p3)
        s5 = signed_tetra_volume(q1,q2,p3,p1)
        if s3 == s4 and s4 == s5:
            n = np.cross(p2-p1,p3-p1)
            t = -np.dot(q1,n-p1) / np.dot(q1,q2-q1)
            return q1 + t * (q2-q1)
    return None

My test code is:

q0 = np.array([0.0,0.0,1.0])
q1 = np.array([0.0,0.0,-1.0])
p0 = np.array([-1.0,-1.0,0.0])
p1 = np.array([1.0,-1.0,0.0])
p2 = np.array([0.0,1.0,0.0])

print(intersect_line_triangle(q0,q1,p0,p1,p2))

gives:

[ 0.  0. -3.] 

instead of the expected

[ 0.  0. 0.]

looking at the line

t = np.dot(q1,n-p1) / np.dot(q1,q2-q1)

Subtracting p1 from the normal doesn't make sense to me, you want to project from q1 onto the plane of the triangle, so you need to project along the normal, with a distance that is proportional to the ratio of the distance from q1 to the plane and q1-q2 along the normal, right?

The following code fixes this:

n = np.cross(p2-p1,p3-p1)
t = np.dot(p1-q1,n) / np.dot(q2-q1,n)
return q1 + t * (q2-q1)

Answer 2

To find the intersection between a line and a triangle in 3D, follow this approach:

• Compute the plane supporting the triangle,

• Intersect the line with the plane supporting the triangle:

  • If there is no intersection, then there is no intersection with the triangle.

  • If there is an intersection, verify that the intersection point indeed lies in the triangle:

    • Each edge of the triangle together with the normal of the plane supporting the triangle determines a half-space bounding the inside of the triangle (the corresponding bounding plane can be derived from the normal and the edge vertices),
    • Verify that the intersection point lies on the inside of all the edge half-spaces.

Here is some sample code with detailed computations that should work:

// Compute the plane supporting the triangle (p1, p2, p3)
//     normal: n
//     offset: d
//
// A point P lies on the supporting plane iff n.dot(P) + d = 0
//
ofVec3f v21 = p2 - p1;
ofVec3f v31 = p3 - p1;

ofVec3f n = v21.getCrossed(v31);
float d = -n.dot(p1);

// A point P belongs to the line from P1 to P2 iff
//     P = P1 + t * (P2 - P1)
//
// Find the intersection point P(t) between the line and
// the plane supporting the triangle:
//     n.dot(P) + d = 0
//                  = n.dot(P1 + t (P2 - P1)) + d
//                  = n.dot(P1) + t n.dot(P2 - P1) + d
//
//     t = -(n.dot(P1) + d) / n.dot(P2 - P1)
//
ofVec3f P21 = P2 - P1;
float nDotP21 = n.dot(P21);

// Ignore line parallel to (or lying in) the plane
if (fabs(nDotP21) < Epsilon)
    return false;

float t = -(n.dot(P1) + d) / nDotP21;
ofVec3f P = P1 + t * P21;

// Plane bounding the inside half-space of edge (p1, p2): 
//     normal: n21 = n x (p2 - p1)
//     offset: d21 = -n21.dot(p1)
//
// A point P is in the inside half-space iff n21.dot(P) + d21 > 0
//

// Edge (p1, p2)
ofVec3f n21 = n.cross(v21);
float d21 = -n21.dot(p1);

if (n21.dot(P) + d21 <= 0)
    return false;

// Edge (p2, p3)
ofVec3f v32 = p3 - p2;
ofVec3f n32 = n.cross(v32);
float d32 = -n32.dot(p2);

if (n32.dot(P) + d32 <= 0)
    return false;

// Edge (p3, p1)
ofVec3f n13 = n.cross(-v31);
float d13 = -n13.dot(p3);

if (n13.dot(P) + d13 <= 0)
    return false;

return true;

Some comments on the code posted with the question:

• Predefined operations of ofVec3f (.dot() and .cross() for geometric products, etc...) should be preferred when available (more readable, avoids implementation mistakes, etc...),
• The code initially follows the approach above but then only checks that the intersection point is in the 3D axis-aligned bounding box of the line segment [P1, P2]. This combined with possible other errorscould explain why the results are incorrect.
• One can verify that the intersection point is in the 3D axis-aligned bounding box of the (whole) triangle. While this is not enough to guarantee intersection, it can however be used to cull points clearly not intersecting and avoid further complex computations.