Möller–Trumbore intersection algorithm

The Möller–Trumbore ray-triangle intersection algorithm, named after its inventors Tomas Möller and Ben Trumbore, is a fast method for calculating the intersection of a ray and a triangle in three dimensions without needing precomputation of the plane equation of the plane containing the triangle.^[1] Among other uses, it can be used in computer graphics to implement ray tracing computations involving triangle meshes.^[2]

The ray is defined by an origin point ${\displaystyle O}${\displaystyle O} and a direction vector ${\displaystyle {\vec {v}}}${\displaystyle {\vec {v}}}. Every point on the ray can be expressed by ${\displaystyle {\vec {r}}(t)=O+t{\vec {v}}}${\displaystyle {\vec {r}}(t)=O+t{\vec {v}}}, where the parameter ${\displaystyle t}${\displaystyle t} ranges from negative infinity to infinity. The triangle is defined by three vertices, named ${\displaystyle v_{1}}${\displaystyle v_{1}}, ${\displaystyle v_{2}}${\displaystyle v_{2}}, ${\displaystyle v_{3}}${\displaystyle v_{3}}. The plane that the triangle is on, which is needed to calculate the ray-triangle intersection, is defined by a point on the plane, such as ${\displaystyle v_{1}}${\displaystyle v_{1}}, and a vector that is orthogonal to every point on that plane, such as the cross product between the vector from ${\displaystyle v_{1}}${\displaystyle v_{1}} to ${\displaystyle v_{2}}${\displaystyle v_{2}} and the vector from ${\displaystyle v_{1}}${\displaystyle v_{1}} to ${\displaystyle v_{3}}${\displaystyle v_{3}}:

${\displaystyle {\vec {n}}\cdot (P_{1}-P_{2})=0}${\displaystyle {\vec {n}}\cdot (P_{1}-P_{2})=0}, where ${\displaystyle {\vec {n}}=(v_{2}-v_{1})\times (v_{3}-v_{1})}${\displaystyle {\vec {n}}=(v_{2}-v_{1})\times (v_{3}-v_{1})}, and ${\displaystyle P_{1}}${\displaystyle P_{1}} and ${\displaystyle P_{2}}${\displaystyle P_{2}} are any points on the plane.

First, find out if the line produced by the ray intersects with the plane that the triangle is on, and if it does, find the coordinates of that intersection. The only way that the line will not intersect the plane is if the ray's direction vector is parallel to the plane.^[3] When this happens, the dot product between the ray's direction vector and the plane's normal vector will be zero. Otherwise, the line does intersect the plane somewhere, but not necessarily within the triangle.

Using barycentric coordinates, any point on the triangle can be expressed as a convex combination of the triangle's vertices:

${\displaystyle P=wv_{1}+uv_{2}+vv_{3}}${\displaystyle P=wv_{1}+uv_{2}+vv_{3}}

The coefficients must be non-negative and sum to 1, so ${\displaystyle w}${\displaystyle w} can be replaced with ${\displaystyle 1-u-v}${\displaystyle 1-u-v}:

${\displaystyle {\begin{aligned}P&=(1-u-v)v_{1}+uv_{2}+vv_{3}\\P&=v_{1}+u(v_{2}-v_{1})+v(v_{3}-v_{1})\end{aligned}}}${\displaystyle {\begin{aligned}P&=(1-u-v)v_{1}+uv_{2}+vv_{3}\\P&=v_{1}+u(v_{2}-v_{1})+v(v_{3}-v_{1})\end{aligned}}}

where ${\displaystyle P}${\displaystyle P} is any point on the plane. Observe that ${\displaystyle {\vec {e_{1}}}=v_{2}-v_{1}}${\displaystyle {\vec {e_{1}}}=v_{2}-v_{1}} and ${\displaystyle {\vec {e_{2}}}=v_{3}-v_{1}}${\displaystyle {\vec {e_{2}}}=v_{3}-v_{1}} are vectors on the edge of the triangle, and together, they span a plane (which goes through the origin). Each point on that plane can be written as ${\displaystyle ue_{1}+ve_{2}}${\displaystyle ue_{1}+ve_{2}} and can be translated by ${\displaystyle v_{1}}${\displaystyle v_{1}} to "move" that point onto the plane that the triangle is on.

To find ${\displaystyle u}${\displaystyle u} and ${\displaystyle v}${\displaystyle v} for a particular intersection, set the ray expression equal to the plane expression, and put the variables on one side and the constants on the other.

${\displaystyle {\begin{aligned}O+tD&=v_{1}+u(v_{2}-v_{1})+v(v_{3}-v_{1})\\O-v_{1}&=-tD+u(v_{2}-v_{1})+v(v_{3}-v_{1})\end{aligned}}}${\displaystyle {\begin{aligned}O+tD&=v_{1}+u(v_{2}-v_{1})+v(v_{3}-v_{1})\\O-v_{1}&=-tD+u(v_{2}-v_{1})+v(v_{3}-v_{1})\end{aligned}}}

This is a system of linear equations with three equations (one each for ${\displaystyle x}${\displaystyle x}, ${\displaystyle y}${\displaystyle y}, ${\displaystyle z}${\displaystyle z}) and three unknowns (${\displaystyle t}${\displaystyle t}, ${\displaystyle u}${\displaystyle u}, and ${\displaystyle v}${\displaystyle v}), and can be represented as a matrix-vector multiplication.

${\displaystyle {\begin{bmatrix}\vert &\vert &\vert \\-D&(v2-v1)&(v3-v1)\\\vert &\vert &\vert \end{bmatrix}}{\begin{bmatrix}t\\u\\v\end{bmatrix}}=O-v_{1}}${\displaystyle {\begin{bmatrix}\vert &\vert &\vert \\-D&(v2-v1)&(v3-v1)\\\vert &\vert &\vert \end{bmatrix}}{\begin{bmatrix}t\\u\\v\end{bmatrix}}=O-v_{1}}

This equation will always have a solution when the matrix has three linearly independent column vectors in ${\displaystyle \mathbb {R} ^{3}}${\displaystyle \mathbb {R} ^{3}} and is thus invertible. This happens if and only if the triangle vertices aren't collinear and the ray isn't parallel to the plane.

The algorithm can use Cramer's Rule to find the ${\displaystyle t}${\displaystyle t}, ${\displaystyle u}${\displaystyle u}, and ${\displaystyle v}${\displaystyle v} values for an intersection, and if it lies within the triangle, the exact coordinates of the intersection can be found by plugging in ${\displaystyle t}${\displaystyle t} to the ray's equation.

The following is an implementation of the algorithm in C++:

std::optional<vec3>ray_intersects_triangle(constvec3&ray_origin,
constvec3&ray_vector,
consttriangle3&triangle)
{
constexprfloatepsilon=std::numeric_limits<float>::epsilon();
vec3edge1=triangle.b-triangle.a;
vec3edge2=triangle.c-triangle.a;
vec3ray_cross_e2=cross(ray_vector,edge2);
floatdet=dot(edge1,ray_cross_e2);
if(det>-epsilon&&det<epsilon)
return{};// This ray is parallel to this triangle.
floatinv_det=1.0/det;
vec3s=ray_origin-triangle.a;
floatu=inv_det*dot(s,ray_cross_e2);
if((u<0&&abs(u)>epsilon)||(u>1&&abs(u-1)>epsilon))
return{};
vec3s_cross_e1=cross(s,edge1);
floatv=inv_det*dot(ray_vector,s_cross_e1);
if((v<0&&abs(v)>epsilon)||(u+v>1&&abs(u+v-1)>epsilon))
return{};
// At this stage we can compute t to find out where the intersection point is on the line.
floatt=inv_det*dot(edge2,s_cross_e1);
if(t>epsilon)// ray intersection
{
returnvec3(ray_origin+ray_vector*t);
}
else// This means that there is a line intersection but not a ray intersection.
return{};
}

The following is an implementation of the algorithm in Rust using the glam crate:

fnmoller_trumbore_intersection(origin:Vec3,direction:Vec3,triangle:Triangle)->Option<Vec3>{
lete1=triangle.b-triangle.a;
lete2=triangle.c-triangle.a;
letray_cross_e2=direction.cross(e2);
letdet=e1.dot(ray_cross_e2);
ifdet>-f32::EPSILON&&det<f32::EPSILON{
returnNone;// This ray is parallel to this triangle.
}
letinv_det=1.0/det;
lets=origin-triangle.a;
letu=inv_det*s.dot(ray_cross_e2);
ifu<0.0||u>1.0{
returnNone;
}
lets_cross_e1=s.cross(e1);
letv=inv_det*direction.dot(s_cross_e1);
ifv<0.0||u+v>1.0{
returnNone;
}
// At this stage we can compute t to find out where the intersection point is on the line.
lett=inv_det*e2.dot(s_cross_e1);
ift>f32::EPSILON{// ray intersection
letintersection_point=origin+direction*t;
returnSome(intersection_point);
}
else{// This means that there is a line intersection but not a ray intersection.
returnNone;
}
}

The following is an implementation of the algorithm in Java using javax.vecmath from Java 3D API:

publicclass MollerTrumbore{
privatestaticfinaldoubleEPSILON=0.0000001;
publicstaticbooleanrayIntersectsTriangle(Point3drayOrigin,
Vector3drayVector,
TriangleinTriangle,
Point3doutIntersectionPoint){
Point3dvertex0=inTriangle.getVertex0();
Point3dvertex1=inTriangle.getVertex1();
Point3dvertex2=inTriangle.getVertex2();
Vector3dedge1=newVector3d();
Vector3dedge2=newVector3d();
Vector3dh=newVector3d();
Vector3ds=newVector3d();
Vector3dq=newVector3d();
doublea,f,u,v;
edge1.sub(vertex1,vertex0);
edge2.sub(vertex2,vertex0);
h.cross(rayVector,edge2);
a=edge1.dot(h);
if(a>-EPSILON&&a<EPSILON){
returnfalse;// This ray is parallel to this triangle.
}
f=1.0/a;
s.sub(rayOrigin,vertex0);
u=f*(s.dot(h));
if(u<0.0||u>1.0){
returnfalse;
}
q.cross(s,edge1);
v=f*rayVector.dot(q);
if(v<0.0||u+v>1.0){
returnfalse;
}
// At this stage we can compute t to find out where the intersection point is on the line.
doublet=f*edge2.dot(q);
if(t>EPSILON)// ray intersection
{
outIntersectionPoint.set(0.0,0.0,0.0);
outIntersectionPoint.scaleAdd(t,rayVector,rayOrigin);
returntrue;
}else// This means that there is a line intersection but not a ray intersection.
{
returnfalse;
}
}
}

• Badouel intersection algorithm

• Fast Minimum Storage Ray-Triangle Intersection
• Optimizations on the basic algorithm by Möller & Trumbore, code from journal of graphics tools
• Ray-Tracing: Rendering a Triangle
• MATLAB version of this algorithm (highly vectorized)
• Baldwin-Weber ray-triangle intersection algorithm
• Schlick–Subrenat algorithm^[1] for ray-quadrilateral intersection

• Geometric algorithms
• Geometric intersection
• Triangle geometry

• Articles with short description
• Short description matches Wikidata