How to know IF a line segment intersects a triangle in 3d space?

Question

I have a triangle defined by 3 points in 3d space. I also have a line segment defined by 2 points in 3d space. I want to know if they intersect. I don't really need to know the point of intersection.

I don't know any calculus but I know some trig. I know some about matrices but I understand vectors well (3d vectors specifically). Please keep it simple.

Can you walk me through the example problem:

triangle:

a: -4, 3, 0

b: 4, 3, 0

c: -3, -5, 4

line segment:

d: 1, -2, 0

e: -2, 6, 2

EDIT:

I am going to use this in a c++ physics engine.

One answer involved tetrahedron volume calculation from 4 vertices. Please provide formula or show it in code.

UPDATE:

As meowgoesthedog pointed out, I could try to use the Moller-Trumbore intersection algorithm. See my answer below for an alternate solution.

Answer 1

Here is one way to solve your problem. Compute the volume of the tetrahedron Td = (a,b,c,d) and Te = (a,b,c,e). If either volume of Td or Te is zero, then one endpoint of the segment de lies on the plane containing triangle (a,b,c). If the volumes of Td and Te have the same sign, then de lies strictly to one side, and there is no intersection. If Td and Te have opposite signs, then de crosses the plane containing (a,b,c).

From there there are several strategies. One is to compute the point p where de crosses that plane. Then project down to 2D, and solve the point-in-triangle problem in 2D.

Another route is to compute the volumes of the tetrahedra (a,b,d,e), (b,c,d,e), and (c,a,d,e). Then only if all three have the same sign, does de intersect the triangle (a,b,c).

How to compute the volume of a tetrahedron from its corner coordinates is all over the web, and also in Computational Geometry in C.

Answer 2

I implemented the great answer that Joseph gave in python and thought I would share. The function takes a set of line segments and triangles and computes for each line segment if it intersects any of the given triangles.

The first input to the function is a 2xSx3 array of line segments where the first index specifies the start or end point of the segment, the second index refers to the s^th line segment, and the third index points to the x, y,z coordinate of the line segment point.

The second input is a 3XTX3 array of triangle vertices, where the first index specifies one of the three vertices (which don't have to be in any particular order), the second index refers to the t^th triangle, and the third index points to the x,y,z coordinates the the triangle vertex.

The output of this function is a binary array of size S which tells you whether the s^th line segment intersects any of the triangles given. If you want to know which triangles the segments intersect, then just remove the summation of the last line of the function.

def signedVolume(a, b, c, d):
    """Computes the signed volume of a series of tetrahedrons defined by the vertices in 
    a, b c and d. The ouput is an SxT array which gives the signed volume of the tetrahedron defined
    by the line segment 's' and two vertices of the triangle 't'."""

    return np.sum((a-d)*np.cross(b-d, c-d), axis=2)

def segmentsIntersectTriangles(s, t):
    """For each line segment in 's', this function computes whether it intersects any of the triangles
    given in 't'."""
    # compute the normals to each triangle
    normals = np.cross(t[2]-t[0], t[2]-t[1])
    normals /= np.linalg.norm(normals, axis=1)[:, np.newaxis]

    # get sign of each segment endpoint, if the sign changes then we know this segment crosses the
    # plane which contains a triangle. If the value is zero the endpoint of the segment lies on the 
    # plane.
    # s[i][:, np.newaxis] - t[j] -> S x T x 3 array
    sign1 = np.sign(np.sum(normals*(s[0][:, np.newaxis] - t[2]), axis=2)) # S x T
    sign2 = np.sign(np.sum(normals*(s[1][:, np.newaxis] - t[2]), axis=2)) # S x T

    # determine segments which cross the plane of a triangle. 1 if the sign of the end points of s is 
    # different AND one of end points of s is not a vertex of t
    cross = (sign1 != sign2)*(sign1 != 0)*(sign2 != 0) # S x T 

    # get signed volumes
    v1 = np.sign(signedVolume(t[0], t[1], s[0][:, np.newaxis], s[1][:, np.newaxis])) # S x T
    v2 = np.sign(signedVolume(t[1], t[2], s[0][:, np.newaxis], s[1][:, np.newaxis])) # S x T
    v3 = np.sign(signedVolume(t[2], t[0], s[0][:, np.newaxis], s[1][:, np.newaxis])) # S x T

    same_volume = np.logical_and((v1 == v2), (v2 == v3)) # 1 if s and t have same sign in v1, v2 and v3

    return (np.sum(cross*same_volume, axis=1) > 0)

Answer 3

Thanks for the help! This is an alternate solution. The question was for c++ and as meowgoesthedog pointed out, I could try to use the Moller-Trumbore intersection algorithm. This is what I came up with:

#include <math.h>

class vec3 {
public:
    float x, y, z;

    float dot(const vec3 & b) {
        return vec3::x * b.x + vec3::y * b.y + vec3::z * b.z;
    }

    vec3 cross(const vec3 & b) {
        return vec3::vec3(
            vec3::y * b.z - vec3::z * b.y,
            vec3::z * b.x - vec3::x * b.z,
            vec3::x * b.y - vec3::y * b.x
        );
    }

    vec3 normalize() {
        const float s = 1.0f / sqrtf(vec3::x * vec3::x + vec3::y * vec3::y + vec3::z * vec3::z);
        return vec3::vec3(vec3::x * s, vec3::y * s, vec3::z * s);
    }

    vec3 operator+(const vec3 & b) {
        return vec3::vec3(
            vec3::x + b.x,
            vec3::y + b.y,
            vec3::z + b.z
        );
    }
    vec3 operator+=(const vec3 & b) {
        *this = vec3::operator+(b);
        return *this;
    }

    vec3 operator-(const vec3 & b) {
        return vec3::vec3(
            vec3::x - b.x,
            vec3::y - b.y,
            vec3::z - b.z
        );
    }
    vec3 operator-=(const vec3 & b) {
        *this = vec3::operator-(b);
        return *this;
    }

    vec3 operator*(const vec3 & b) {
        return vec3::vec3(
            vec3::x * b.x,
            vec3::y * b.y,
            vec3::z * b.z
        );
    }
    vec3 operator*=(const vec3 & b) {
        *this = vec3::operator*(b);
        return *this;
    }
    vec3 operator*(float b) {
        return vec3::vec3(
            vec3::x * b,
            vec3::y * b,
            vec3::z * b
        );
    }
    vec3 operator*=(float b) {
        *this = vec3::operator*(b);
        return *this;
    }

    vec3 operator/(const vec3 & b) {
        return vec3::vec3(
            vec3::x / b.x,
            vec3::y / b.y,
            vec3::z / b.z
        );
    }
    vec3 operator/=(const vec3 & b) {
        *this = vec3::operator/(b);
        return *this;
    }
    vec3 operator/(float b) {
        return vec3::vec3(
            vec3::x * b,
            vec3::y * b,
            vec3::z * b
        );
    }
    vec3 operator/=(float b) {
        *this = vec3::operator/(b);
        return *this;
    }

    vec3(float x, float y, float z) {
        vec3::x = x;
        vec3::y = y;
        vec3::z = z;
    }
    vec3(float x) {
        vec3::x = x;
        vec3::y = x;
        vec3::z = x;
    }
    vec3() {
        //
    }
    ~vec3() {
        //
    }
};

#define EPSILON 0.000001f
bool lineSegIntersectTri(
    vec3 line[2],
    vec3 tri[3],
    vec3 * point
) {
    vec3 e0 = tri[1] - tri[0];
    vec3 e1 = tri[2] - tri[0];

    vec3 dir = line[1] - line[0];
    vec3 dir_norm = dir.normalize();

    vec3 h = dir_norm.cross(e1);
    const float a = e0.dot(h);

    if (a > -EPSILON && a < EPSILON) {
        return false;
    }

    vec3 s = line[0] - tri[0];
    const float f = 1.0f / a;
    const float u = f * s.dot(h);

    if (u < 0.0f || u > 1.0f) {
        return false;
    }

    vec3 q = s.cross(e0);
    const float v = f * dir_norm.dot(q);

    if (v < 0.0f || u + v > 1.0f) {
        return false;
    }

    const float t = f * e1.dot(q);

    if (t > EPSILON && t < sqrtf(dir.dot(dir))) { // segment intersection
        if (point) {
            *point = line[0] + dir_norm * t;
        }

        return true;
    }

    return false;
}

For running a few tests:

#include <stdio.h>

const char * boolStr(bool b) {
    if (b) {
        return "true";
    }

    return "false";
}

int main() {
    vec3 tri[3] = {
        { -1.0f, -1.0f, 0.0f },
        { 1.0f, -1.0f, 0.0f },
        { 1.0f, 1.0f, 0.0f },
    };

    vec3 line0[2] = { // should intersect
        { 0.5f, -0.5f, -1.0f },
        { 0.5f, -0.5f, 1.0f },
    };

    vec3 line1[2] = { // should not intersect
        { -0.5f, 0.5f, -1.0f },
        { -0.5f, 0.5f, 1.0f },
    };

    printf(
        "line0 intersects? : %s\r\n"
        "line1 intersects? : %s\r\n",
        boolStr(lineSegIntersectTri(line0, tri, NULL)),
        boolStr(lineSegIntersectTri(line1, tri, NULL))
    );

    return 0;
}

Answer 4

C# version:

public class AlgoritmoMollerTrumbore
  {
    private const double EPSILON = 0.0000001;

    public static bool lineIntersectTriangle(Point3D[] line,
                                             Point3D[] triangle,
                                             out Point3D outIntersectionPoint)
    {
      outIntersectionPoint = new Point3D(0, 0, 0);

      Point3D rayOrigin = line[0];
      Vector3D rayVector = Point3D.Subtract(line[1], line[0]);
      rayVector.Normalize();

      Point3D vertex0 = triangle[0];
      Point3D vertex1 = triangle[1]; 
      Point3D vertex2 = triangle[2];

      Vector3D edge1 = Point3D.Subtract(vertex1, vertex0);
      Vector3D edge2 = Point3D.Subtract(vertex2, vertex0);
      Vector3D h = Vector3D.CrossProduct(rayVector, edge2);

      double a = Vector3D.DotProduct(edge1, h);
      if (a > -EPSILON && a < EPSILON)
      {
        return false;    // This ray is parallel to this triangle.
      }

      double f = 1.0 / a;      
      Vector3D s = Point3D.Subtract(rayOrigin, vertex0);

      double u = f * (Vector3D.DotProduct(s, h));
      if (u < 0.0 || u > 1.0)
      {
        return false;
      }
      Vector3D q = Vector3D.CrossProduct(s, edge1);
      double v = f * Vector3D.DotProduct(rayVector, q);
      if (v < 0.0 || u + v > 1.0)
      {
        return false;
      }
      // At this stage we can compute t to find out where the intersection point is on the line.
      double t = f * Vector3D.DotProduct(edge2, q);
      if (t > EPSILON && t < Math.Sqrt(Vector3D.DotProduct(rayVector, rayVector))) // ray intersection
      {
        outIntersectionPoint = rayOrigin + rayVector * t;
        return true;
      }
      else // This means that there is a line intersection but not a ray intersection.
      {
        return false;
      }
    }
  }

Answer 5

Unity C# version:

public static bool RaycastTriangle(
    Vector3 origin, 
    Vector3 pointB, 
    Vector3 a, 
    Vector3 b, 
    Vector3 c, 
    out Vector3 intersection
) {
    intersection = Vector3.zero;
    var direction = pointB - origin;
    var normalizedDirection = direction.normalized;
    
    Vector3 edge1 = b - a;
    Vector3 edge2 = c - a;
    Vector3 normal = Vector3.Cross(normalizedDirection, edge2);

    float dotProduct = Vector3.Dot(edge1, normal);

    if (dotProduct > -float.Epsilon && dotProduct < float.Epsilon)
        return false; // Parallel

    float inverseDotProduct = 1f / dotProduct;
    Vector3 toStart = origin - a;
    float triangleParam = inverseDotProduct * Vector3.Dot(toStart, normal);

    if (triangleParam is < 0f or > 1f) return false;
    Vector3 edgeCross = Vector3.Cross(toStart, edge1);
    float raycastParam = inverseDotProduct * Vector3.Dot(normalizedDirection, edgeCross);

    if (raycastParam < 0f || triangleParam + raycastParam > 1f) return false;
    
    float distance = inverseDotProduct * Vector3.Dot(edge2, edgeCross);
    if (distance > float.Epsilon && distance < direction.magnitude)
    {
        intersection = origin + normalizedDirection * distance;
        return true;
    }
    return false;
}