Ray-Tracing: Rendering a Triangle

Möller-Trumbore algorithm

Möller-Trumbore algorithm

The Möller-Trumbore (or MT for the remainder of this lesson) algorithm is a fast ray-triangle intersection algorithm that was introduced in 1997 by Tomas Möller and Ben Trumbore in a paper titled "Fast, Minimum Storage Ray/Triangle Intersection". It is still considered today a fast algorithm that is often used in benchmarks to compare the performances of other methods. However, a fair comparison of ray-triangle intersection algorithms is difficult because speed can depend on many factors, such as how algorithms are implemented, the type of test scene used, whether values are precomputed, etc.

The Möller-Trumbore algorithm takes advantage of the parameterization of P, the intersection point in barycentric coordinates (which we discussed in the previous chapter). We learned in the previous chapter to calculate P, the intersection point, using the following equation:

We also learned that \(w=1-u-v\); thus, we can write:

If we develop, we get (equation 1):

Note that \( (B - A) \) and \( (C - A) \) are the edges \( AB \) and \( AC \) of the triangle \( \triangle ABC \). The intersection point \( P \) can also be written using the ray's parametric equation (Equation 2).

Where \(t\) is the distance from the ray's origin to the intersection \(P\). If we replace \(P\) in equation 1 with the ray's equation, we get (equation 3):

On the right side of the equal sign, we have:

• \( t \), \( u \), and \( v \) as three unknowns,

• respectively multiplied by \( B - A \), \( C - A \), and \( D \), which are three known variables.

We can rearrange these terms and present equation 3 using the following notation (equation 4):

The left side of the equation has been rearranged into a row-column vector multiplication (1x3*3x1). This is the simplest possible form of matrix multiplication. You take the first element of the row matrix \(-D\), \(B-A\), \(C-A\) and multiply it by the first element of the column vector (\(t\), \(u\), \(v\)). Then you add the raw matrix's second element multiplied by the column vector's second element. Then finally, add the third element of the raw matrix multiplied by the third element of the column vector (which gives you equation 3 again).

In equation 3, the term \(O-A\) on the left side of the equal sign is a vector. \(B-A\), \(C-A\), and \(D\) are vectors as well, and \(t\), \(u\), and \(v\) (unknown) are scalars (floats or doubles if you prefer to think in terms of programming). This equation is about vectors. It combines three vectors in quantities defined by \(t\), \(u\), and \(v\). It gives another vector as a result of this operation (be aware that in their paper, Möller and Trumbore use the term x-, y- and z-axis instead of \(t\), \(u\) and \(v\) when they explain the geometrical meaning of this equation).

Let's say that \( P \) is where the ray intersects the triangle \( \triangle ABC \). You will (hopefully) agree that the Cartesian coordinates of point \( P \) (x, y, z) change as we move or rotate the triangle and the ray together. However, the barycentric coordinates of \( P \) are invariant. If the triangle is rotated, scaled, stretched, or translated, the coordinates \( u \) and \( v \), which define the position of \( P \) with respect to vertices \( A \), \( B \), and \( C \), will not change (see Figure 1).

The Möller-Trumbore algorithm takes advantage of this property. Instead of solving the ray-triangle intersection equation using the \(x\), \(y\), and \(z\) coordinates for the intersection point, it expresses the intersection point in terms of the triangle's barycentric coordinates \(u\) and \(v\). A point lying on the triangle (our intersection point, if it exists) can thus be expressed either in Cartesian coordinates \((x, y, z)\) or in barycentric coordinates \((u, v)\).

Finally, as we learned in the previous chapter, let's recall that the \(u\) and \(v\) coordinates of the point lying on the surface of a triangle can't be greater than one nor lower than 0. Their sum can't be greater than one either \(u + v <= 1\). They express coordinates of points defined inside a unit triangle (this is the triangle defined in \(u\), \(v\) space by the vertices (0, 0), (1, 0), (0, 1) as shown in Figure 1).

Let's summarize. What do we have so far? We have reinterpreted the three-dimensional \(x\), \(y\), \(z\) position of point \(P\) in terms of \(u\) and \(v\) barycentric coordinates, effectively moving from one parameter space to another: from xyz-space to uv-space. We have also learned that points defined in uv-space lie within a unit triangle.

Equation 3 (or 4, as they are identical) has three unknowns: \(u\), \(v\), and \(t\). Geometrically, we have already explained the meaning of \(u\) and \(v\). But what about \(t\)? Simply put, we will consider \(t\) as the third axis of the \(u\) and \(v\) coordinate system we just introduced. We now have a coordinate system defined by three axes: \(u\), \(v\), and \(t\).

Geometrically, let's explain what \(t\) represents. We know that \(t\) expresses the distance from the ray origin to \(P\), the intersection point. If you look at Figure 1, you can see that we have created a coordinate system that allows us to fully describe the position of the intersection point \(P\) in terms of barycentric coordinates and the distance from the ray origin to that point on the triangle.

Cramer's Rule

What we need to do now is solve this equation, which has three unknowns: \(t\), \(u\), and \(v\). To solve Equation 4, Möller and Trumbore use a technique known in mathematics as Cramer's rule. Cramer's rule provides the solution to a system of linear equations in terms of determinants. It states that if the multiplication of a matrix \(M\) (a set of three numbers laid out horizontally) by a column vector \(X\) (three numbers laid out vertically) is equal to a column vector \(C\), then it is possible to find \(X_i\) (the \(i\)th element of the column vector \(X\)) by dividing the determinant of \(M_i\) by the determinant of \(M\). Here, \(M_i\) is the matrix formed by replacing the \(i\)th column of \(M\) with the column vector \(C\).

If you are unfamiliar with Cramer's rule, it is explained in the gray box below.

In linear algebra, the determinant is denoted by two vertical bars, and it reads: \( \begin{vmatrix} I & J & K \end{vmatrix} \) is the determinant of the matrix having \( I \), \( J \), and \( K \) as its columns. Following the conventions used in the gray box above, you can now see that matrix \( M \) is \( [-D, B - A, C - A] \), \( X \) is \( [t, u, v] \), and \( C \) is \( O - A \) in our system of linear equations. We are interested in finding values for \( t \), \( u \), and \( v \). To understand how we apply Cramer's rule let's recap:

Cramer's rule states that for a system of linear equations \( MX = C \), where \( M \) is a matrix of coefficients, \( X \) is a vector of unknowns, and \( C \) is the result vector, the solution for each component \( X_i \) can be found by replacing the \( i \)-th column of the matrix \( M \) with the vector \( C \), then dividing the determinant of the resulting matrix by the determinant of \( M \).

In this case, the matrix \( M \) is \( [-D, B - A, C - A] \), and the right-hand side vector \( C \) is \( O - A \). So, the solution using Cramer's rule for the unknowns \( t \), \( u \), and \( v \) is:

Where:

• \( M_t \) is the matrix \( M \) with its first column replaced by \( O - A \),

• \( M_u \) is the matrix \( M \) with its second column replaced by \( O - A \),

• \( M_v \) is the matrix \( M \) with its third column replaced by \( O - A \).

Applying this technique to our problem gives (Equation 5):

where, for brevity, we denote \( T = O - A \), \( E_1 = B - A \), and \( E_2 = C - A \). The equation for each unknown (\( t \), \( u \), and \( v \)) is expressed as the ratio of a determinant.

• The determinant of \( M \) is written as \( \left| \begin{array}{r} -D & E_1 & E_2 \end{array} \right| \), where \( E_1 = B - A \) and \( E_2 = C - A \).

• \( M_t \), \( M_u \), \( M_v \) are matrices where you replace the corresponding columns of \( M \) with \( T = O - A \). These are used to compute the determinants for \( t \), \( u \), and \( v \), respectively.

The next step is to find the values of these four determinants. In the lesson on Geometry, we explained that the determinant of a \(1 \times 3\) matrix \( \begin{vmatrix} A & B & C \end{vmatrix} \) can be calculated as follows (refer to the grey box below for more details on where these equations come from):

The determinant of a \(3 \times 3\) matrix is a scalar triple product, which combines a cross product and a dot product. Therefore, we can rewrite Equation 5 as:

where \( P = (D \times E_2) \) and \( Q = (T \times E_1) \).

Implementation

Implementing the Möller-Trumbore algorithm is straightforward. The scalar triple product \( (A \times B) \cdot C \) is the same as \( A \cdot (B \times C) \). Therefore, the denominator \( (D \times E2) \cdot E1 \) in Equation 5 is the same as \( D \cdot (E1 \times E2) \). The cross product of \( E1 \) and \( E2 \) gives the normal of the triangle. We know that if the dot product of the ray direction \( D \) and the triangle normal is 0, then the triangle and the ray are parallel (and therefore, there is no intersection).

Remember that the user might want to cull (discard) back-facing triangles (the back side of the triangle, which typically will never face the camera). If the triangle is front-facing, the determinant is positive; otherwise, it is negative. In the code, we will first compute this determinant. There is no intersection if the result is negative and the triangle is single-sided, or if the determinant is close to 0. If the triangle is double-sided, we will still need to check if the determinant is close to 0, which is more complex than in the previous case, as the value is signed.

bool rayTriangleIntersect(
    const Vec3f &orig, const Vec3f &dir,
    const Vec3f &v0, const Vec3f &v1, const Vec3f &v2,
    float &t, float &u, float &v)
{
    Vec3f v0v1 = v1 - v0;
    Vec3f v0v2 = v2 - v0;
    Vec3f pvec = dir.crossProduct(v0v2);
    float det = v0v1.dotProduct(pvec);
#ifdef CULLING
    // if the determinant is negative, the triangle is 'back facing.'
    // if the determinant is close to 0, the ray misses the triangle
    if (det < kEpsilon) return false;
#else
    // ray and triangle are parallel if det is close to 0
    if (fabs(det) < kEpsilon) return false;
#endif
    ...
}

The following steps are simple. We first compute u and reject the triangle if u is either lower than 0 or greater than 1\. If we successfully pass the computation of \(u\), then we compute \(v\) and apply the same tests (there's no intersection if \(v\) is lower than 0 or greater than one and if u+v is greater than 1). We know the ray intersects the triangle at this point, and we can compute \(t\).

bool rayTriangleIntersect(
    const Vec3f &orig, const Vec3f &dir,
    const Vec3f &v0, const Vec3f &v1, const Vec3f &v2,
    float &t, float &u, float &v)
{
    ...
    float invDet = 1 / det;

    Vec3f tvec = orig - v0;
    u = tvec.dotProduct(pvec) * invDet;
    if (u < 0 || u > 1) return false;

    Vec3f qvec = tvec.crossProduct(v0v1);
    v = dir.dotProduct(qvec) * invDet;
    if (v < 0 || u + v > 1) return false;
    
    t = v0v2.dotProduct(qvec) * invDet;
    
    return true;
}

Note that when the ray and the normal of the triangle are facing each other (they go in opposite directions), the determinant is positive (greater than 0), as shown in the adjacent figure (top). On the other hand, when the ray hits the triangle from "behind" (the ray and the normal points in the same direction), the determinant is negative (adjacent figure, bottom). The figure on the side illustrates these two possible cases. The numbers are simple enough that cross and dot products between vectors are straightforward to compute. The coordinates of edge2 are (1,0,0), and the cross product between the ray direction (0,-1,0) and edge2 gives pvec. Remember that the cross product of two vectors is provided by:

$$ \begin{array}{l} x = v1.y v2.z - v1.z v2.y\\ y = v1.z v2.x - v1.x v2.z\\ z = v1.x v2.y - v1.y v2.x \end{array} $$

Thus:

$$ \begin{array}{l} pvec.z &=& dir.x * edge.y - dir.y * edge.x \\ &=& 0*0 -(-1 * 1) = 1, \end{array} $$

and pvec = (0,0,1). The dot product between this vector and edge1 (the determinant) is positive. When the ray has a direction (0,1,0), we get:

$$pvec.z = 0*0 - (1 * 1) = -1,$$

and pvec = (0,0,-1). The determinant is, in this case, negative.

In the original implementation of the algorithm (which you can find in the original paper, see reference below), the code can be either compiled to handle culling or ignored. When culling is active, rays intersecting triangles from behind will be discarded. This can easily be verified using the sign of the determinant. If culling is on and the determinant is lower than 0, then the ray doesn't intersect the triangle. However, when culling is off, we want to intersect triangles regardless of the normal orientation with resepct to the ray direction.

int intersect_triangle(...)
{
    ...
    det = DOT(edge1, pvec);
#ifdef CULLING        /* define CULLING if culling is desired */
    if (det < kEpsilon)
        return 0;
    ...
#else        /* the non-culling branch */
    if (det > -kEpsilon && det < kEpsilon)
        return 0;
    ...
#endif
}

The point we want to make here is that when culling is off and the determinant is negative, it becomes essential to "normalize \(u\)" by multiplying it by the inverse of the determinant. Because we later test if u is greater than 0 and lower than 1, when \(u\) is negative, but the determinant is also negative, then the sign of \(u\) when "normalized" (multiplied by invDet) is inverted (it becomes positive).

Source Code

The C++ implementation of the MT algorithm is straightforward (if you have an excellent vector library) and is also very compact, as you can see in this code snippet:

bool rayTriangleIntersect(
    const Vec3f &orig, const Vec3f &dir,
    const Vec3f &v0, const Vec3f &v1, const Vec3f &v2,
    float &t, float &u, float &v)
{
#ifdef MOLLER_TRUMBORE
    Vec3f v0v1 = v1 - v0;
    Vec3f v0v2 = v2 - v0;
    Vec3f pvec = dir.crossProduct(v0v2);
    float det = v0v1.dotProduct(pvec);
#ifdef CULLING
    // if the determinant is negative, the triangle is 'back facing'
    // if the determinant is close to 0, the ray misses the triangle
    if (det < kEpsilon) return false;
#else
    // ray and triangle are parallel if det is close to 0
    if (fabs(det) < kEpsilon) return false;
#endif
    float invDet = 1 / det;

    Vec3f tvec = orig - v0;
    u = tvec.dotProduct(pvec) * invDet;
    if (u < 0 || u > 1) return false;

    Vec3f qvec = tvec.crossProduct(v0v1);
    v = dir.dotProduct(qvec) * invDet;
    if (v < 0 || u + v > 1) return false;
    
    t = v0v2.dotProduct(qvec) * invDet;
    
    return true;
#else
    ...
#endif
}

Compile with:

c++ -o raytri raytri.cpp -O3 -std=c++11 -DMOLLER_TRUMBORE

Important note regarding the MT algorithm implementation: If you compile the code provided with this lesson, you might notice that the image produced when the MT algorithm is used is not the same as the image of the triangle displayed above. This can be confusing. That's because the image of the triangle shown above was rendered using the geometry solution. Not the MT algorithm. The image below shows what the two algorithms output. On the left is the Möller–Trumbore solution. On the right is the geometric solution.

Why are they different? Mathematically they are not. It's just that the two algorithms are not calculating u, v, and w in the same order. The coordinate u in the geometric solution is equal to w in the MT solution, etc., and u in the first solution is numerically similar to w in the second. So don't worry. It's not because they produce a different image that one or the other is incorrect. They are both correct. You'll need to potentially swap the coordinates when using them for texturing and shading.

Notes

One way of accelerating the ray triangle is to store the constants A, B, C, and D of the plane equation with the triangle's vertices (we precompute them once the mesh is loaded and store them in memory). However, storing variables in memory and accessing them can be more expansive than the time required to compute them on the fly. Therefore this is not necessarily proven to be a reliable optimization.

One way of accelerating the algorithm is to reduce the ray-triangle intersection to two dimensions. To do so, we determine the dominant axis of the triangle's normal and use this information to select the two coordinates we use for the projection. This technique is described in "Ray-Polygon Intersection. An Efficient Ray-Polygon Intersection" (see references section below).

Algorithms to ray-trace quads already exist.

Our code uses single precision floating numbers, which are generally enough, but the ray-triangle intersection test may suffer from precision errors in some conditions. For example, the intersection point might not lie on the triangle plane and be either slightly above or below the triangle's surface. We will learn more about this in the lesson on shading.

Many of these topics are or will be addressed in the advanced section.

Exercises

• Add the camera transformation back to the source code of the main program.

• Compute more than one triangle with proper visibility using the trace() function we used in the previous lesson.

• Render a triangle from the back and test the culling branch of the code.

References

Fast, Minimum Storage Ray/Triangle Intersection. Möller & Trumbore. Journal of Graphics Tools, 1997.

Ray-Polygon Intersection. An Efficient Ray-Polygon Intersection. Didier Badouel from Graphics Gems I. 1990