Ray Triangle Intersection #

Möller, Trumbore triangle intersection #

Tomas Möller and Ben Trumbore wrote a tiny paper with the title “Fast, Minimum Storage Ray/Triangle Intersection”. It describes an algorithm that can compute the intersection of a ray with a triangle. The paper was published in 1997. It contains one sentence that should draw your attention to this paper: “…, we believe it is the fastest ray/triangle intersection routine for triangles …”. I call the paper tiny because it consists only of 7 pages. Two pages of the paper consist mainly of a C implementation of the described algorithm. One page of the paper contains only some nice pictures. The paper is very brief. According to Google Scholar, the paper was cited more than 1300 times. I suggest you to read this paper now and then come back to this page.

My implementation of the algorithm looks like this:

// Möller–Trumbore intersection algorithm
// Fast Minimum Storage Ray/Triangle Intersection
template <typename ScalarType>
bool intersect_triangle_moeller_trumbore(const RayType<3, ScalarType>& ray,
                    const Point3<ScalarType>& v0,
                    const Point3<ScalarType>& v1,
                    const Point3<ScalarType>& v2,
                    ScalarType& t,
                    ScalarType& u,
                    ScalarType& v) {
    const ScalarType EPSILON = ScalarType{0.000001};

    // find vectors for two edges sharing vert0
    auto edge1 = v1 - v0;
    auto edge2 = v2 - v0;

    // begin calculating determinant - also used to calculate U parameter
    auto pvec = ray.direction.cross(edge2);

    // if determinant is near zero, ray lies in plane of triangle
    auto det = edge1.dot(pvec);

    if (det > -EPSILON && det < EPSILON)
        return false;
    auto inv_det = ScalarType{1.0} / det;

    // calculate distance from vert to ray origin
    auto tvec = ray.origin - v0;

    // calculate U parameter and test bounds
    u = tvec.dot(pvec) * inv_det;
    if(u < ScalarType{0.0} || u > ScalarType{1.0}) {
        return false;
    }

    // prepare to test V parameter
    auto qvec = tvec.cross(edge1);

    // calculate V parameter and test bounds
    v = ray.direction.dot(qvec) * inv_det;
    if(v < ScalarType{0.0} || u + v > ScalarType{1.0}) {
        return false;
    }

    // calculate t, ray intersects triangle
    t = edge2.dot(qvec) * inv_det;

    return true;
}

I just replaced the doubles with a vector and ray type and tried to stick as close as possible to the original implementation. Instead of using an int as a return code, I switched to bool. Since I do not need support for back face culling I removed this part.

The interesting part of my implementation is my unit tests.

TEST(intersect_triangle_moeller_trumbore, trivial_hit) {
    // Arrange

    // setup ray
    Point3f ray_origin{0.f,0.f,100.f};
    Vector3f ray_direction{0.f,0.f,-1.f};
    Ray3f ray{ray_origin, ray_direction, 0.f, 1000.f};

    // counter clockwise is front facing
    std::vector<Point3f> vertices = {
        {-1.f, -1.f, 0.f},
        {1.f, -1.f, 0.f},
        {0.f, 1.f, 0.f},
    };

    float t = -1.f;
    float u = -1.f;
    float v = -1.f;

    // Act
    bool hit = intersect_triangle_moeller_trumbore(ray, 
                                                    vertices[0], 
                                                    vertices[1], 
                                                    vertices[2], 
                                                    t, u, v);

    // Assert
    EXPECT_THAT(hit, true);
    EXPECT_THAT(t, 100);
    EXPECT_THAT(u, 0.25);
    EXPECT_THAT(v, 0.5);
}

Maybe you ask yourself why is \(u = 0.25\) and \(v = 0.5\) . Take a look at the paper again. It says: \(T(u,v) = (1-u-v)V_0+uV_1+uV_2\) . The total area of the triangle in the above unit test is \(2\) . The intersection point splits the triangle into three sub-triangles. Two of them have the same size the other one has the area \(1\) . This gives us the barycentric coordinates \((0.25, 0.25, 0.5)\) .

Here is another interesting unit test. It is almost the same test, but this time only the ray direction is inverted. The idea here is to miss the triangle (no hit).

TEST(intersect_triangle_moeller_trumbore, hit_even_it_is_a_miss) {
    // Arrange

    // setup ray
    Point3f ray_origin{0.f,0.f,100.f};
    Vector3f ray_direction{0.f,0.f,1.f};
    Ray3f ray{ray_origin, ray_direction, 0.f, 1000.f};

    // counter clockwise is front facing
    std::vector<Point3f> vertices = {
            {-1.f, -1.f, 0.f},
            {1.f, -1.f, 0.f},
            {0.f, 1.f, 0.f},
    };

    float t = -1.f;
    float u = -1.f;
    float v = -1.f;

    // Act
    bool hit = intersect_triangle_moeller_trumbore(ray, vertices[0], vertices[1], vertices[2], t, u, v);

    // Assert
    EXPECT_THAT(hit, true);
    EXPECT_THAT(t, -100);
    EXPECT_THAT(u, 0.25);
    EXPECT_THAT(v, 0.5);
}

To my surprise, the algorithm still reports a hit (returns true), but also gives a negative distance. When just copy & pasting the implementation without recognizing this this can cause long and complicated debug sessions.

Ray Triangle Intersection #

Möller, Trumbore triangle intersection #

Further reading #