Möller-Trumbore algorithm

Möller-Trumbore algorithm

The Möller-Trumbore (or MT for the reminding of this lesson) algorithm is a fast ray-triangle intersection algorithm which was introduced in 1997 by Tomas Möller and Ben Trumbore in a paper entitled "Fast, Minimum Storage Ray/Triangle Intersection". It is still considered today a fast algorithm which is often used in benchmarks to compare performances of other methods although, a fair comparaison of ray-triangle intersection algorithms is actually a difficult thing to do, because speed can depends on many factors such as the way algorithms are implemented, the type of test scene which is used, whether values are pre-computed, etc. 

The Möller-Trumbore algorithm takes advantage of the parametrisation of P, the intersection point in terms of barycentric coordinates (which we talked about in the previous chapter). We have learned in the previous section to compute P, the intersection point, using the following equations:

$$\scriptsize P = wA + uB + vC$$

Now, we know we can also express P in terms, of u, v and w where, as we learned, we can replace w with 1-u-v which leads to the following equation:

$$\scriptsize P = (1 - u - v)A + uB + vC$$

Therefore we can write (equation 1):

$$\scriptsize P=A - uA - vA + uB + vC=A + u(B - A) + v(C - A)$$

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

$$\scriptsize P=O+tD$$

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):

$$\scriptsize\begin{array}{l} O+tD & = & A + u(B - A) + v(C - A)\\ O-A &=&-tD+u(B-A)+v(C-A)\end{array}$$

On the left side of the equal sign, we have three unknowns (t, u, v) multiplied to three known terms (B-A, C-A, D). We can re-arrange these terms and present equation 3 using the following notation (equation 4):

$$\scriptsize\begin{bmatrix}-D&(B-A)&(C-A)\end{bmatrix}\begin{bmatrix} t\\u\\v\end{bmatrix}=O-A$$

The left side of the equation has been rearranged into a row-colum vector multiplication. This is the simplest possible form of matrix multiplication. You just take the first element of the raw matrix (-D, B-A, C-A) and multiply it by the first element of the column vector (t, u, v). Then you add the second element of the raw matrix multiplied by the second element of the column vector. 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 on the right side of the equal sign is a vector (O-A). B-A, C-A and D are vectors as well, and t, u, and v (unknown) are scalars. This equation is about vectors. It combines three vectors in quantities defined by t, u and v and 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). 

Figure 1: we can express the position of P in t, u, v space. t indicates the distance from P to the ray origin and is parallel to the t-axis. if P lies in the unit triangle it means that 0 <= u <= 1, 0 <= v <= 1 and u+v <=1.

Let's imagine that you have a point P inside the triangle ABC. If we transform the triangle in any ways (for example if we translate the triangle, rotate it or scale it), the coordinates of point P expressed in the the three-dimensional coordinate system x, y, z will change. On the other hand if you express the position of P using barycentric coordinates, transformations applied to the triangle will not affect the intersection point barycentric coordinates. if the triangle is rotated, scaled, stretched or translated, the coordinates u, v defining the position of P in regards to vertices A, B and C will not change (see figure xx). That's a property upon which the MT algorithm is built. Indeed what they do, is to define a new coordinate system where the coordinates of P are not defined in terms of x, y and z but in terms of u and v. We can perfectly express the same point in different coordinate systems (x, y, z space or u, v space) which is exactly what you do if you use meters instead of inches to express the distance between two points. The distance is the same but you just use a different metrics to measure it. Now, u, v coordinates as we have learned can't be greater than 1 nor lower than 0. Their sum can't be greater than 1 either (u+v <= 1). In fact they express coordinates of points define inside the unit triangle (this is the triangle defined in u, v space by the vertices (0, 0), (1, 0), (0, 1), see figure xx). Lets summarize. What do we have so far ? We have re-interpreted the three-dimensional x, y, z position of point P in terms of u, v barycentric coordinates going from one parameter space to another: the x-, y-, z- coordinate space to the u- v- space. We also have learned that points defined in u, v space are inside the unit triangle.

Equation 3 (or 4, they are identical) has three unknowns: u, v and t. Geometrically, we have just explained the meaning of u and v. But what about t? Simple: we will consider that t is the third axis of the u, v coordinate system we just introduced. Now we have a coordinate system defined by three axes, u, v and t. But again, geometrically let's explain what it means. In fact we know that t expresses the distance from the ray origin to the intersection P. Now if you look at figure xx, we can see that we have developed a coordinate system that allows us to fully express the position of the intersection P in terms of barycentric coordinates and distance from the ray origin to that point on the triangle. Note that the third axis t is a coordinate. To find the position of the ray origin in t, u, v space, you trace a line parallel to the t-axis of length t and starting at P (figure 1. You can find the same figure in the original paper).

Möller and Trumbore explain that the first part of equation 4 (the term O-A) can be looked at as a transformation moving the triangle from its original world space position to the origin (the first vertex of the triangle coincides with the origin). The other side of the equation has for effect to transform the intersection point from x, y, z space to t, u, v space as explained above.

Cramer's Rule

What we have to do now is solve this equation which has three unknowns: t, u and v. To solve equation 4, Möller and Trumbore chose a technique which is known in mathematics, as the Cramer's rule. Cramer's rule gives the solution to a system of linear equations in terms of determinant.  It says that if the multiplication of a matrix M (three numbers disposed in a horizontal manner) by a colum vector X (three number disposed in vertical manner) is equal to a column vector C, then it is possible to find Xi (the ith element of the column vector X) by dividing the determinant of Mi by the determinant of M. Mi is the matrix formed by replacing the ith column of M by the column vector C.

Imaging we have to solve for x y and z using this set of linear equations:

$$ \scriptsize \begin{array}{l} 2x + 1y + 1z = 3 \\ 1x - 1y - 1z = 0 \\ 1x + 2y + 1z = 0 \end{array}$$

We can rewrite these equations using a matrix representation, where the coefficients of the equations on the left will become the coefficients of the matrix M and the numbers on the right of the equal sign will become the matrix C's coefficient:

$$ \scriptsize \begin{array}{l} M=\left[ \begin{array}{rrr}2 & 1 & 1\\1 & -1 & -1\\1 & 2 & 1\end{array}\right] \\ C=\begin{bmatrix}\color{\red}{3 \\ 0 \\ 0}\end{bmatrix}\\ M*\begin{bmatrix}x \\ y \\ z \end{bmatrix} = C \end{array} $$

By replacing the first column in M with C we can create a matrix Mx:

$$ \scriptsize M_x=\left[ \begin{array}{rrr}\color{\red}{3} & 1 & 1 \\ \color{\red}{0} & -1 & -1 \\ \color{\red}{0} & 2 & 1\end{array} \right] $$

Similarly, we can create My and Mz:

$$ \scriptsize \begin{array}{l} M_y=\left[ \begin{array}{rrr}2 & \color{\red}{3} & 1 \\ 1 & \color{\red}{0} & -1 \\ 1 & \color{\red}{0} & 1\end{array} \right]\\ M_z=\left[ \begin{array}{rrr}2 & 1 & \color{\red}{3} \\ 1 & -1 & \color{\red}{0} \\ 1 & 2 & \color{\red}{0}\end{array} \right] \end{array} $$

Evaluating each determinant, we get:

$$ \scriptsize \begin{array}{l} M=\left| \begin{array}{rrr} 2 & 1 & 1\\1 & -1 & -1\\1 & 2 & 1 \end{array} \right| = (-2)+(-1)+(2)-(-1)-(-4)-(1)=3\\ M_x=\left| \begin{array}{rrr} \color{\red}{3} & 1 & 1 \\ \color{\red}{0} & -1 & -1 \\ \color{\red}{0} & 2 & 1 \end{array} \right| = (-3) + (0) + (0) - (0)-(-6) - (0) = 3\\ M_y=\left| \begin{array}{rrr} 2 & \color{\red}{3} & 1 \\ 1 & \color{\red}{0} & -1 \\ 1 & \color{\red}{0} & 1 \end{array} \right| = (0) + (-3) + (0) - (0) - (0) -(3) = -6\\ M_z=\left| \begin{array}{rrr} 2 & 1 & \color{\red}{3} \\ 1 & -1 & \color{\red}{0} \\ 1 & 2 & \color{\red}{0} \end{array} \right| = (0) + (0) + (6) -(-3) - (0) - (0) = 9\\ \end{array} $$

Cramer's rule says that x = Mx / M, y = My / M and z = Mz / M. That is:

$$ \scriptsize\begin{array}{l} x = {3 \over 3} = 1\\ y = {-6 \over 3}= -2\\ z = {9 \over 3} = 3 \end{array} $$

In linear algebra the determinant is denoted by two vertical bars. And it reads: |A B C| is the determinant of the matrix having A B & C as its columns. You can easily see by now that the 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 and by applying Cramer's rule we can therefore write (equation 5):

$$ \scriptsize \left[ \begin{array}{r} t \\ u \\ v\end{array}\right] = {1 \over \left[ \left| \begin{array}{r} -D & E_1 & E_2 \end{array}\right| \right]} \left[ \begin{array}{c} \left| \begin{array}{c} T & E_1 & E_2\end{array}\right| \\ \left| \begin{array}{c} -D & T & E_2\end{array}\right| \\ \left| \begin{array}{c} -D & E_1 & T\end{array}\right| \\ \end{array}\right], $$

where for brevity we denote \(\scriptsize T = O - A\) \(\scriptsize E_1 = B - A\) and \(\scriptsize E_2 = C - A\). The next step is to find a value for these four deteminants. In lesson 4, we showed that the determinant of a 1x3 matrix |A B C| has the value:

$$\scriptsize |A B C| = -(A \times C) \cdot B = -(C \times B) \cdot A.$$

The determinant (of a 3x3 matrix) is nothing else than a scalar triple product which is a combination of a cross and a dot product (this is explained in details in lesson 4). Therefore we can rewrite equation 5 as:

$$ \scriptsize \left[ \begin{array}{r} t \\ u \\ v\end{array}\right] = {1 \over {(D \times E_2) \cdot E_1}} \left[ \begin{array}{c} (T \times E_1) \cdot E_2 \\ (D \times E_2) \cdot T \\ (T \times E_1) \cdot D \end{array}\right] = = {1 \over {P \cdot E_1}} \left[ \begin{array}{c} Q \cdot E_2 \\ P \cdot T \\ Q \cdot D \end{array}\right],$$

where \(\scriptsize P=(D \times E_2)\) and \(\scriptsize Q=(T \times E_1)\). As you can see it is easy for us now to find the values t, u and v. We can get them with cross and dot products between known variables: the vertices of the triangle, the origin and the ray direction.

Implementation

Implementing the MT algorithm is straightforward. The scalar triple product (AxB).C is the same as A.(BxC). Therefore the denominator (DxE2).E1 in equation 5 is the same as D.(E1xE2). The cross product of E1 and E2 gives the normal of the triangle. And we know that if the dot product of the ray direction D and the triangle normal is 0, the triangle and the ray are parallel (and therefore there is not intersection). Remember as well that the user might want to cull (discard) backfacing triangles. If the triangle is front-facing the determinant is positive otherwise it is negative. In the code we will first compute this determinant. If the result is negative and the triangle single sided or close to 0 there is not intersection. If the triangle is double sided we will need to check if the determinant is close to 0 which is a bit more complex than in the previous case as the value is signed.

#define EPSILON 1e-6 int intersectTriangle( const point va, const point vb, const point vc, int singleSided, const point rayOrig, const vector rayDir, float &t, float &u, float &v) { vector e1 = vb - va; vector e2 = vc - va; vector pvec = cross(rayDir, e2); float det = dot(e1, pvec); if (singleSided) { // if the determinant is negative the triangle is backfacing // if the determinant is close to 0, the ray misses the triangle if (det < EPSILON) return 0; } else { // ray and triangle are parallel if det is close to 0 if (det > -EPSILON && det < EPSILON) return 0; } ... }

The next 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 1 and if u+v is greater than 1). At this point we know the ray intersects the triangle and we can compute t.

... float invDet = 1 / det; // prepare to compute u vector tvec = rayOrig - va; u = dot(tvec, pvec) * invDet; if (u < 0 || u > 1) return 0; // prepare to compute v vector qvec = cross(tvec, e1); v = dot(rayDir, qvec) * invDet; if (v < 0 || v > 1 || u + v > 1) return 0; // calculate t, ray intersects triangle t = dot(e2, qvec) * invDet; return 1; }

Errata: "Because we multiply u and v by invDet we can delay this multiplication until the end and test the bounds against the determinant which results in a slight speed up (if there is no intersection, we avoid a division and two multiplies)". This sentence was actually part of a previous version of the lesson but this statement is in fact partially inexact (at least confusing) and we decided not only to remove it but also make an additional comment regarding the algorithm. In fact, 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 showed 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 straighforward 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 given by:

$$\scriptsize \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:

$$\scriptsize pvec.z = raydir.x * edge.y - raydir.y * edge.x = 0*0 -(-1 * 1) = 1,$$

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

$$\scriptsize 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 simply ignore it. When culling is active, rays intersecting triangles from behind will be discarded. This can easily be checked for, using the sign of the determinant. If culling is active and the determinant is lower than 0, then we the ray doesn't intersect the triangle. However when when culling is off, we want to intersect triangles regardless of the normal orientation with respect to the ray direction.

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

However, the point here is that when culling is off and the determinant negative, it becomes important 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 also negative, then the sign of u when "normalized" (multiplied by invDet) is reversed (it becomes positive). We could delay the multiplication of u by the inverse of the determinant further down the code but it would then mean making more tests on the value of u such as:

if (((det < 0 && u >= det && u =< 0) || (det > 0 && u <= det && u >= 0)) && det > -EPSILON && det < EPSILON) ...

or

if (u <= fabs(det) && fabs(u) >= 0) && fabs(det) < EPSILON) ...

As you can guess this is not an optimisation because to potential save a couple of multiplications and a division, we have to make the code more complex in other places. In other words, it is simply best to keep the code straightforward, multiply u by the inverse of the determinant right away and test the bounds of u and v against the range [0,1]. This example shows that being obssessed with micro-optimisation is actually not a good way of using your time when you should actually be more concerned about writing elegant code (this is a note for people who think too much about micro-optimisation at the expense of code simplicity, readibility and maintainability).

Source Code

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
bool intersect(const Ray<float> &r, IsectData &isectData) const { #ifdef MOLLER_TRUMBORE Vec3f edge1 = v1 - v0; Vec3f edge2 = v2 - v0; Vec3f pvec = cross(r.dir, edge2); float det = dot(edge1, pvec); if (det == 0) return false; float invDet = 1 / det; Vec3f tvec = r.orig - v0; isectData.u = dot(tvec, pvec) * invDet; if (isectData.u < 0 || isectData.u > 1) return false; Vec3f qvec = cross(tvec, edge1); isectData.v = dot(r.dir, qvec) * invDet; if (isectData.v < 0 || isectData.u + isectData.v > 1) return false; isectData.t = dot(edge2, qvec) * invDet; #else ... #endif return true; }

Notes

These are important notes we don't want to address in this lesson which is supposed to be an introduction to the topic. But we will develop them as soon as we can. In the meatime we will briefly mention them.

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 it requires to compute them on the fly. Therefore this is not necessarily proven to be a realiable optimisation.

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 is generally enough but the ray-triangle intersection test may in some conditions suffer from precision errors. For example the intersection point might not exactly lie on the triangle plane and be either slightly above or below the surface of the triangle. We will learn more about this in the lesson on shading.

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

Chapter 6 of 6