Ray Triangle Intersection issue

Question

I'm trying to perform ray triangle intersection on a mesh made of triangles. The below code seems to work fine but only about 50% of the time. The ray often gets into positions where no intersection is detected when there definitely should be. Is there something wrong in my function that isn't checking for intersection under all circumstances.

bool rayTriangleIntersect(Vector3 rayOrigin, 
                          Vector3 rayDirection, 
                          Vector3 vert0, 
                          Vector3 vert1, 
                          Vector3 vert2, 
                          ref Vector3 intersectionPosition)
{
    // this originally comes from http://www.graphics.cornell.edu/pubs/1997/MT97.pdf

     float u, v, t, det, inv_det;          
     Vector3 edge1, edge2, pvec, tvec, qvec;

    Vector3 direction = rayDirection;
    direction.Normalize();
    t = v = u = 0f;
    edge1 = vert1 - vert0;
    edge2 = vert2 - vert0;
    pvec = Vector3.Cross(direction, edge2);
    det = Vector3.Dot(edge1, pvec);

    if (det > -0.00001f)
        return false;

    inv_det = 1.0f / det;
    tvec = rayOrigin - vert0;
    u = Vector3.Dot(tvec, pvec) * inv_det;

    if (u < -0.001f || u > 1.001f)
        return false;

    qvec = Vector3.Cross(tvec, edge1);
    v = Vector3.Dot(direction, qvec) * inv_det;

    if (v < -0.001f || u + v > 1.001f)
        return false;

    t = Vector3.Dot(edge2, qvec) * inv_det;

    float rayMagnitude = rayDirection.magnitude;

    if (t <= 0f || t >= rayMagnitude)
        return false;

    intersectionPosition = rayOrigin + (rayDirection * (t / rayMagnitude)); // percent along ray 

    return true;
}

and im calling it with

for (z = 0; z < LENGTH - 1; z++)
{
     for (x = 0; x < WIDTH - 1; x++)
     {
           v0 = verts[(z * WIDTH) + x];
           v1 = verts[((z + 1) * WIDTH) + x];
           v2 = verts[(z * WIDTH) + x + 1];
           v3 = verts[((z + 1) * WIDTH) + x + 1];

           if(rayTriangleIntersect(rayOrigin, rayDirection, v0, v1, v2, ref intersectPosition) ||   
              rayTriangleIntersect(rayOrigin, rayDirection, v1, v2, v3, ref intersectPosition))
           {  
                 Debug.log(Intersection found at : " intersectPosition);  
           }
    }
}

I'm really stuck on this one so any guidance would really be appreciated.

Edit

The triangle are defined like this:

 for (z = 0; z < LENGTH - 1; z++)
   {
        for (x = 0; x < WIDTH - 1; x++)
        {
            triangles[index++] = (z * WIDTH) + x;
            triangles[index++] = ((z + 1) * WIDTH) + x;
            triangles[index++] = (z * WIDTH) + x + 1;

            triangles[index++] = ((z + 1) * WIDTH) + x;
            triangles[index++] = ((z + 1) * WIDTH) + x + 1;
            triangles[index++] = (z * WIDTH) + x + 1;
        }
    }

Answer 1

Try changing

 if (det > -0.00001f)
        return false;

to,

 if (det > -0.00001f && det < 0.00001f)
        return false;

Currently, your code only checks one side of the triangle

Answer 2

Your code seems to use quads and it seems that you don't use the two triangles from the quad but two overlapping ones (say the second one is wrong), try this instead (see the "v2, v3, v0" part) :

          if(rayTriangleIntersect(rayOrigin, rayDirection, v0, v1, v2, ref intersectPosition) ||   
          rayTriangleIntersect(rayOrigin, rayDirection, v2, v3, v0, ref intersectPosition))

If that doesn't work, you could try v0,v3,v1 or v1,v3,v2 or check out how your quads are actually built and figure it out.

Answer 3

I don't know what's the problem with your implementation, but I have used the method below (taken directly from this sample) on a raytracer and it worked correctly and performed well too. So you might consider swapping to this implementation:

/// <summary>
/// Checks whether a ray intersects a triangle. This uses the algorithm
/// developed by Tomas Moller and Ben Trumbore, which was published in the
/// Journal of Graphics Tools, volume 2, "Fast, Minimum Storage Ray-Triangle
/// Intersection".
/// 
/// This method is implemented using the pass-by-reference versions of the
/// XNA math functions. Using these overloads is generally not recommended,
/// because they make the code less readable than the normal pass-by-value
/// versions. This method can be called very frequently in a tight inner loop,
/// however, so in this particular case the performance benefits from passing
/// everything by reference outweigh the loss of readability.
/// </summary>
static void RayIntersectsTriangle(ref Ray ray,
                                  ref Vector3 vertex1,
                                  ref Vector3 vertex2,
                                  ref Vector3 vertex3, out float? result)
{
    // Compute vectors along two edges of the triangle.
    Vector3 edge1, edge2;

    Vector3.Subtract(ref vertex2, ref vertex1, out edge1);
    Vector3.Subtract(ref vertex3, ref vertex1, out edge2);

    // Compute the determinant.
    Vector3 directionCrossEdge2;
    Vector3.Cross(ref ray.Direction, ref edge2, out directionCrossEdge2);

    float determinant;
    Vector3.Dot(ref edge1, ref directionCrossEdge2, out determinant);

    // If the ray is parallel to the triangle plane, there is no collision.
    if (determinant > -float.Epsilon && determinant < float.Epsilon)
    {
        result = null;
        return;
    }

    float inverseDeterminant = 1.0f / determinant;

    // Calculate the U parameter of the intersection point.
    Vector3 distanceVector;
    Vector3.Subtract(ref ray.Position, ref vertex1, out distanceVector);

    float triangleU;
    Vector3.Dot(ref distanceVector, ref directionCrossEdge2, out triangleU);
    triangleU *= inverseDeterminant;

    // Make sure it is inside the triangle.
    if (triangleU < 0 || triangleU > 1)
    {
        result = null;
        return;
    }

    // Calculate the V parameter of the intersection point.
    Vector3 distanceCrossEdge1;
    Vector3.Cross(ref distanceVector, ref edge1, out distanceCrossEdge1);

    float triangleV;
    Vector3.Dot(ref ray.Direction, ref distanceCrossEdge1, out triangleV);
    triangleV *= inverseDeterminant;

    // Make sure it is inside the triangle.
    if (triangleV < 0 || triangleU + triangleV > 1)
    {
        result = null;
        return;
    }

    // Compute the distance along the ray to the triangle.
    float rayDistance;
    Vector3.Dot(ref edge2, ref distanceCrossEdge1, out rayDistance);
    rayDistance *= inverseDeterminant;

    // Is the triangle behind the ray origin?
    if (rayDistance < 0)
    {
        result = null;
        return;
    }

    result = rayDistance;
}

Answer 4

Without looking at the code in much detail, I'm going to take a guess that your triangle intersection routine only works from one side of the triangle, but your triangle "mesh" describes triangles with opposite "handedness" (if v0,v1,v2 is clockwise, then v1,v2,v3 is anticlockwise as seen by a ray from the same direction). This could be confirmed by e.g permuting the order of the latter to v1,v3,v2 and seeing if all or none of the tests succeed. If that is indeed the problem, and it's desired to have an intersection test succeed from either side of the triangle (rather than "back face culling"), then you probably need to handle the case of a negative "det" better.