Ray-triangle intersection

Question

I saw that Fast Minimum Storage Ray/Triangle Intersection by Moller and Trumbore is frequently recommended.

The thing is, I don't mind pre-computing and storing any amounts of data, as long as it speeds-up the intersection.

So my question is, not caring about memory, what are the fastest methods of doing ray-triangle intersection?

Edit: I wont move the triangles, i.e. it is a static scene.

Answer 1

As others have mentioned, the most effective way to speed things up is to use an acceleration structure to reduce the number of ray-triangle intersections needed. That said, you still want your ray-triangle intersections to be fast. If you're happy to precompute stuff, you can try the following:

Convert your ray lines and your triangle edges to Plücker coordinates. This allows you to determine if your ray line passes through a triangle at 6 multiply/add's per edge. You will still need to compare your ray start and end points with the triangle plane (at 4 multiply/add's per point) to make sure it actually hits the triangle.

Worst case runtime expense is 26 multiply/add's total. Also, note that you only need to compute the ray/edge sign once per ray/edge combination, so if you're evaluating a mesh, you may be able to use each edge evaluation twice.

Also, these numbers assume everything is being done in homogeneous coordinates. You may be able to reduce the number of multiplications some by normalizing things ahead of time.

Answer 2

I have done a lot of benchmarks, and I can confidently say that the fastest (published) method is the one invented by Havel and Herout and presented in their paper Yet Faster Ray-Triangle Intersection (Using SSE4). Even without using SSE it is about twice as fast as Möller and Trumbore's algorithm.

My C implementation of Havel-Herout:

    typedef struct {
        vec3 n0; float d0;
        vec3 n1; float d1;
        vec3 n2; float d2;
    } isect_hh_data;

    void
    isect_hh_pre(vec3 v0, vec3 v1, vec3 v2, isect_hh_data *D) {
        vec3 e1 = v3_sub(v1, v0);
        vec3 e2 = v3_sub(v2, v0);
        D->n0 = v3_cross(e1, e2);
        D->d0 = v3_dot(D->n0, v0);

        float inv_denom = 1 / v3_dot(D->n0, D->n0);

        D->n1 = v3_scale(v3_cross(e2, D->n0), inv_denom);
        D->d1 = -v3_dot(D->n1, v0);

        D->n2 = v3_scale(v3_cross(D->n0, e1), inv_denom);
        D->d2 = -v3_dot(D->n2, v0);
    }

    inline bool
    isect_hh(vec3 o, vec3 d, float *t, vec2 *uv, isect_hh_data *D) {
        float det = v3_dot(D->n0, d);
        float dett = D->d0 - v3_dot(o, D->n0);
        vec3 wr = v3_add(v3_scale(o, det), v3_scale(d, dett));
        uv->x = v3_dot(wr, D->n1) + det * D->d1;
        uv->y = v3_dot(wr, D->n2) + det * D->d2;
        float tmpdet0 = det - uv->x - uv->y;
        int pdet0 = ((int_or_float)tmpdet0).i;
        int pdetu = ((int_or_float)uv->x).i;
        int pdetv = ((int_or_float)uv->y).i;
        pdet0 = pdet0 ^ pdetu;
        pdet0 = pdet0 | (pdetu ^ pdetv);
        if (pdet0 & 0x80000000)
            return false;
        float rdet = 1 / det;
        uv->x *= rdet;
        uv->y *= rdet;
        *t = dett * rdet;
        return *t >= ISECT_NEAR && *t <= ISECT_FAR;
    }