View source for Möller–Trumbore intersection algorithm

{{Short description|Method of calculating ray-triangle intersections in 3D space}}

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 (Euclidean geometry)|intersection]] of a [[ray (geometry)|ray]] and a [[triangle]] in [[Three-dimensional space|three dimensions]] without needing precomputation of the plane equation of the [[Plane (geometry)|plane]] containing the triangle.<ref>{{cite journal | doi = 10.1080/10867651.1997.10487468 | volume=2 | title=Fast, Minimum Storage Ray-Triangle Intersection | year=1997 | journal=Journal of Graphics Tools | pages=21–28 | last1 = Möller | first1 = Tomas | last2 = Trumbore | first2 = Ben}}</ref> Among other uses, it can be used in [[computer graphics]] to implement [[ray tracing (graphics)|ray tracing]] computations involving [[triangle mesh]]es.<ref>{{cite web|url=http://www.lighthouse3d.com/tutorials/maths/ray-triangle-intersection|title=Ray-Triangle Intersection|date=26 March 2011 |publisher=lighthouse3d|accessdate=2017-09-10}}</ref>

== Calculation ==
=== Definitions ===
The ray is defined by an origin point <math>O</math> and a direction vector <math>\vec{v}</math>.
Every point on the ray can be expressed by <math>\vec{r}(t) = O + t\vec{v}</math>, where the parameter <math>t</math> ranges from zero to infinity.
The triangle is defined by three vertices, named <math>v_1</math>, <math>v_2</math>, <math>v_3</math>.
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 <math>v_1</math>, and a vector that is orthogonal to every point on that plane, such as the cross product between the vector from <math>v_1</math> to <math>v_2</math> and the vector from <math>v_1</math> to <math>v_3</math>:

<math>\vec{n} \cdot (P_1 - P_2) = 0</math>, where <math>\vec{n} = (v_2 - v_1) \times (v_3 - v_1)</math>, and <math>P_1</math> and <math>P_2</math> are any points on the plane.

=== Check if the ray is parallel to the triangle ===
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.<ref>Note: If the ray's origin is itself on the plane, in addition to the ray's direction vector being parallel to the plane, than the entire ray is technically on the plane. However, since theoretical planes are infinitely thin, the ray would still be considered to not intersect the plane in that scenario.</ref> 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.

=== Check if the ray-plane intersection lies outside the triangle ===
Using [[Barycentric coordinate system|barycentric coordinates]], any point on the triangle can be expressed as a [[convex combination]] of the triangle's vertices: 
:<math>P = w v_1 + u v_2 + v v_3</math>
<!-- pictures of examples are needed -->

The coefficients must be non-negative and sum to 1, so <math>w</math> can be replaced with <math>1 - u - v</math>:
:<math>\begin{align}
P &= (1 - u - v) v_1 + u v_2 + v v_3 \\
P &= v_1 + u (v_2 - v_1) + v (v_3 - v_1)
\end{align}</math>
where <math>P</math> is any point on the plane.
Observe that <math>\vec{e_1} = v_2 - v_1</math> and <math>\vec{e_2} = v_3 - v_1</math> 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 <math>u e_1 + v e_2</math> and can be translated by <math>v_1</math> to "move" that point onto the plane that the triangle is on.

To find <math>u</math> and <math>v</math> 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.
:<math>\begin{align}
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{align}</math>

This is a system of linear equations with three equations (one each for <math>x</math>, <math>y</math>, <math>z</math>) and three unknowns (<math>t</math>, <math>u</math>, and <math>v</math>), and can be represented as a matrix-vector multiplication.
:<math>
\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
</math>
This equation will always have a solution when the matrix has three linearly independent column vectors in <math>\mathbb{R}^3</math> 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 <math>t</math>, <math>u</math>, and <math>v</math> values for an intersection, and if it lies within the triangle, the exact coordinates of the intersection can be found by plugging in <math>t</math> to the ray's equation.

==Implementations==

===C++ implementation===
The following is an implementation of the algorithm in [[C++]]:
<syntaxhighlight lang="cpp">
std::optional<vec3> ray_intersects_triangle( const vec3 &ray_origin,
    const vec3 &ray_vector,
    const triangle3& triangle)
{
    constexpr float epsilon = std::numeric_limits<float>::epsilon();

    vec3 edge1 = triangle.b - triangle.a;
    vec3 edge2 = triangle.c - triangle.a;
    vec3 ray_cross_e2 = cross(ray_vector, edge2);
    float det = dot(edge1, ray_cross_e2);

    if (det > -epsilon && det < epsilon)
        return {};    // This ray is parallel to this triangle.

    float inv_det = 1.0 / det;
    vec3 s = ray_origin - triangle.a;
    float u = inv_det * dot(s, ray_cross_e2);

    if ((u < 0 && abs(u) > epsilon) || (u > 1 && abs(u-1) > epsilon))
        return {};

    vec3 s_cross_e1 = cross(s, edge1);
    float v = 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.
    float t = inv_det * dot(edge2, s_cross_e1);

    if (t > epsilon) // ray intersection
    {
        return  vec3(ray_origin + ray_vector * t);
    }
    else // This means that there is a line intersection but not a ray intersection.
        return {};
}
</syntaxhighlight>

===Rust implementation===
The following is an implementation of the algorithm in [[Rust (programming language)|Rust]] using the glam crate:
<syntaxhighlight lang="rust">
fn moller_trumbore_intersection (origin: Vec3, direction: Vec3, triangle: Triangle) -> Option<Vec3> {
	let e1 = triangle.b - triangle.a;
	let e2 = triangle.c - triangle.a;

	let ray_cross_e2 = direction.cross(e2);
	let det = e1.dot(ray_cross_e2);

	if det > -f32::EPSILON && det < f32::EPSILON {
		return None; // This ray is parallel to this triangle.
	}

	let inv_det = 1.0 / det;
	let s = origin - triangle.a;
	let u = inv_det * s.dot(ray_cross_e2);
	if u < 0.0 || u > 1.0 {
		return None;
	}

	let s_cross_e1 = s.cross(e1);
	let v = inv_det * direction.dot(s_cross_e1);
	if v < 0.0 || u + v > 1.0 {
		return None;
	}
	// At this stage we can compute t to find out where the intersection point is on the line.
	let t = inv_det * e2.dot(s_cross_e1);

	if t > f32::EPSILON { // ray intersection
		let intersection_point = origin + direction * t;
		return Some(intersection_point);
	}
	else { // This means that there is a line intersection but not a ray intersection.
		return None;
	}
}
</syntaxhighlight>

===Java implementation===
The following is an implementation of the algorithm in [[Java (programming language)|Java]] using <code>javax.vecmath</code> from [[Java 3D]] API:
<syntaxhighlight lang="java">
public class MollerTrumbore {

    private static final double EPSILON = 0.0000001;

    public static boolean rayIntersectsTriangle(Point3d rayOrigin, 
                                                Vector3d rayVector,
                                                Triangle inTriangle,
                                                Point3d outIntersectionPoint) {
        Point3d vertex0 = inTriangle.getVertex0();
        Point3d vertex1 = inTriangle.getVertex1();
        Point3d vertex2 = inTriangle.getVertex2();
        Vector3d edge1 = new Vector3d();
        Vector3d edge2 = new Vector3d();
        Vector3d h = new Vector3d();
        Vector3d s = new Vector3d();
        Vector3d q = new Vector3d();
        double a, f, u, v;
        edge1.sub(vertex1, vertex0);
        edge2.sub(vertex2, vertex0);
        h.cross(rayVector, edge2);
        a = edge1.dot(h);

        if (a > -EPSILON && a < EPSILON) {
            return false;    // 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) {
            return false;
        }

        q.cross(s, edge1);
        v = f * rayVector.dot(q);

        if (v < 0.0 || u + v > 1.0) {
            return false;
        }

        // At this stage we can compute t to find out where the intersection point is on the line.
        double t = f * edge2.dot(q);
        if (t > EPSILON) // ray intersection
        {
            outIntersectionPoint.set(0.0, 0.0, 0.0);
            outIntersectionPoint.scaleAdd(t, rayVector, rayOrigin);
            return true;
        } else // This means that there is a line intersection but not a ray intersection.
        {
            return false;
        }
    }
}
</syntaxhighlight>

== See also ==
* [[Badouel intersection algorithm]]

== References ==
{{reflist}}

== External links ==
* [http://www.graphics.cornell.edu/pubs/1997/MT97.pdf Fast Minimum Storage Ray-Triangle Intersection]
* [https://github.com/erich666/jgt-code/tree/master/Volume_02/Number_1/Moller1997a Optimizations on the basic algorithm by Möller &amp; Trumbore], code from ''journal of graphics tools''
* [https://www.scratchapixel.com/lessons/3d-basic-rendering/ray-tracing-rendering-a-triangle/moller-trumbore-ray-triangle-intersection.html Ray-Tracing: Rendering a Triangle]
* [http://www.mathworks.com/matlabcentral/fileexchange/33073 MATLAB version] of this algorithm (highly vectorized)
* [http://jcgt.org/published/0005/03/03/ Baldwin-Weber ray-triangle intersection algorithm]
* Schlick–Subrenat algorithm<ref>''[http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.44.6955 Ray Intersection of Tessellated Surfaces: Quadrangles versus Triangles]'', Schlick C., Subrenat G.  Graphics Gems 1993</ref> for ray-quadrilateral intersection

{{DEFAULTSORT:Moller-Trumbore intersection algorithm}}
[[Category:Geometric algorithms]]
[[Category:Geometric intersection]]
[[Category:Triangle geometry]]