# Student-friendly ray-triangle intersection

I'm teaching a computer graphics course, and would like to give my students a function for calculating ray-triangle intersections (not just the point of intersection, but also the uv coordinates within the triangle). This seems like a good opportunity to review some key ideas in 3d vector geometry, so I'd like to pick an algorithm that is intuitive and helps review key vector topics.

My original approach to this problem was to write out the equations and solve the system using a matrix inverse: Three equations (the x, y, and z coordinates of the intersection point) and three unknowns (how far along the ray, and the uv coordinates in the triangle). It looks like this is similar to Möller–Trumbore. This would be good for reviewing how to solve systems of linear equations, but doesn't seem to have much vector intuition behind it.

The textbook I'm using as a reference suggests using a ray-plane intersection (a nice dot product, which is good for reviewing dot products!) then solving for u and v as a system of two equations in two unknowns.

I'm not interested in the most efficient algorithm, but something that gives me an opportunity to teach something. Any suggestions? Any resources I can link my students to?

• The Scratchapixel tutorials build up the necessary background on barycentric coordinates and provides a framework building up to Möller–Trumbore. The chapter begins with less efficient - though perhaps more immediately intuitive ray-triangle intersection, using edge test equations. Commented Mar 7, 2022 at 12:01

If you want something that looks more vectorial, start with doing a ray-plane intersection; then you can get the barycentric coordinates as follows:

For each edge of the triangle, construct the plane through that edge that's perpendicular to the plane of the triangle (use cross products, etc).

Then you can use the plane equation to measure (a) the distance from the edge to the opposite vertex of the triangle, and (b) the distance from the edge to the intersection point. Then divide b by a to obtain the barycentric coordinate associated with that opposite vertex.

(UVs in the triangle are just picking two of the three barycentric coordinates; since they add up to 1, you can reconstruct the third once you have two.)

The University of Washington, Computer Science department teaches a ray, triangle intersection algorithm that is fairly intuitive, has good documentation (they did their own write up here), and covers a nice range of concepts. It also reviews several important concepts with nice explanations. I thought about doing a write up here but the original paper is better so I will just leave the link.

I am aware of two main ways for deriving the intersection between a ray and a triangle. I believe most ray-tri intersection routines are variants of those. In fact Möller–Trumbore is not something special from a purely mathematical perspective, it's just a specific grouping of terms when solving a 3x3 system.

The first derivation involves finding barycentric and ray coordinates simultaneously which results in a 3x3 system. The second chooses to first intersect the ray with the plane, in which the triangle lies, and then finds the barycentric coordinates.

0.1) Parametric Ray Definition

A parametric definition of a line given an origin $$\vec{o}$$ and a direction $$\vec{d}$$ is:

$$\vec{r}(t) = \vec{o} + t\vec{d}.$$

A ray is the set of points on one side of the origin of the line:

$$\mathcal{R} = \{\vec{r}(t) = \vec{o} + t\vec{d}\,:\,t\in[0,\infty)\}.$$

You should convince your students that this is indeed a ray. A good approach would be to draw several point on the ray for different $$t$$ vaues (e.g. 0, 1, 0.5, but also note that -1 is on the line but not part of the ray).

0.2) Parametric Triangle Definition

A parametric definition of a plane given a point on it $$\vec{v}_0$$ and two vectors $$\vec{f}_1 = \vec{v}_1 - \vec{v}_0$$, $$\vec{f}_2 = \vec{v}_2-\vec{v}_0$$ is:

$$\vec{v}(\beta, \gamma) = \vec{v}_0 + \beta\vec{f}_1 + \gamma\vec{f}_2$$

The pair $$(\beta, \gamma)$$ are the coordinates of $$\vec{v}$$ in the 2D coordinate system with origin $$\vec{v}_0$$ and basis vectors $$\vec{f}_1, \vec{f}_2$$. For those to be basis vectors we additionally require that $$\vec{f}_1, \vec{f}_2$$ are not parallel/linearly dependent (i.e. the triangle is non-degenerate).

You should convince the students that $$\vec{v}(\beta, \gamma)$$ indeed spans the whole plane for different $$(\beta, \gamma)$$. For instance you could draw the points $$(0,0), (1,0), (0,1), (1,1)$$ which are the vertices of a parallelogram. Additionally, the points forming the 3 edges of the triangle are: $$(0,\gamma); \gamma \in[0,1]$$, $$(\beta, 0); \beta \in [0,1]$$, and $$(\beta, \gamma); \beta \geq 0, \gamma \geq 0, \beta + \gamma = 1$$. Then the set of points forming the triangle are:

$$\mathcal{T} = \{\vec{v}(\beta, \gamma)\,:\,0\leq \beta, 0\leq \gamma, \beta+\gamma \leq 1\}.$$

This should definitely be accompanied with a drawing of the geometry so that the students may grasp its meaning.

0.2.1) Barycentric Coorindates

This is something that I would personally leave for the very end, since it is not integral to any of the derivation, is less general, and more complex.

The coordinates $$\beta, \gamma$$ are in fact the barycentric coordinates of $$\vec{v}$$ associated with $$\vec{v}_1, \vec{v}_2$$. In fact we can rewrite our formulation as:

\begin{align}\vec{v}(\beta, \gamma) &= \vec{v}_0 + \beta \vec{f}_1 + \gamma \vec{f}_2 \\ \vec{v}(\beta, \gamma) &= \alpha\vec{v}_0 + \beta\vec{v}_1 + \gamma\vec{v}_2, \alpha = 1-\beta -\gamma,\end{align}

where $$\alpha, \beta, \gamma$$ are the barycentric coordinates of point $$\vec{v}$$ with respect to the triangle $$\vec{v}_0, \vec{v}_1, \vec{v}_2$$. The choice of $$(\beta, \gamma)$$ over $$(\alpha, \beta)$$ or $$(\alpha, \gamma)$$ was arbitrary. Indeed using $$\alpha+\beta+\gamma=1$$ I could have equivalently written:

$$\vec{v}(\alpha, \beta) = \vec{v}_2 + \alpha (\vec{v}_0-\vec{v}_2) + \beta (\vec{v}_1-\vec{v}_2) \\ \text{or} \\ \vec{v}(\alpha, \gamma) = \vec{v}_1 + \alpha (\vec{v}_0-\vec{v}_1) + \gamma (\vec{v}_2-\vec{v}_1).$$

However, the fact that $$\alpha, \beta, \gamma$$ are the barycentric coordinates is not really important, and can be kept as an additional fact in order to not overload the students. The 2D basis interpretation is simpler and much more general.

0.3) Problem Formulation

We want to find the closest intersection (the one with smallest $$t$$) between $$\mathcal{R}$$ and $$\mathcal{T}$$. Let $$\vec{p}$$ be an intersection point, i.e. $$\vec{p} \in \mathcal{R}$$ and $$\vec{p} \in \mathcal{T}$$, then there exist $$t\geq 0, \beta\geq 0, \gamma\geq 0, \beta+\gamma\leq 1$$ such that:

$$\vec{o}+t\vec{d} = \vec{p} = \vec{v}_0 + \beta\vec{f}_1 + \gamma\vec{f}_2.$$

Conversely, if those do not exist then we'll say that there is no intersection. Beyond this point the only difference is how one groups the terms when solving the above $$3\times 3$$ system. For example Möller–Trumbore is just a special case of a specific solution order for the unknowns $$t, \beta, \gamma$$.

1) Simultaneous Solution

Let us rewrite the equations in matrix form:

$$\vec{o}+t\vec{d} = \vec{v}_0 + \beta\vec{f}_1 + \gamma\vec{f}_2$$ $$\begin{bmatrix}\vec{d} & \vec{f}_1 & \vec{f}_2 \end{bmatrix} \begin{bmatrix} -t \\ \beta \\ \gamma \end{bmatrix} = \vec{o} - \vec{v}_0.$$

It should be made clear for the students that $$M = \begin{bmatrix}\vec{d} & \vec{f}_1 & \vec{f}_2 \end{bmatrix}$$ is a 3x3 matrix, didactically, being explicit there may help:

$$M = \begin{bmatrix}d_x & f_{1,x} & f_{2,x} \\ d_y & f_{1,y} & f_{2,y} \\ d_z & f_{1,z} & f_{2,z}\end{bmatrix}.$$

Solving the system is as simple as:

$$\begin{bmatrix} -t \\ \beta \\ \gamma \end{bmatrix} = M^{-1}(\vec{o}-\vec{v}_0)$$

Either $$3\times 3$$ matrix inversion followed by matrix-vector multiplication can be used, or the equivalent Cramer's rule formulation. There's also a trick involving the cross product that works for $$3\times 3$$ matrix inversion that some people like to use (which is elegant, but didactically not very useful):

$$M^{-1} = \frac{1}{det[M]}adj(M) = \frac{1}{\vec{d}\cdot(\vec{f}_1\times\vec{f}_2)}\begin{bmatrix} (\vec{f}_1\times\vec{f}_2)^T \\ (\vec{f}_2 \times \vec{d})^T \\ (\vec{d}\times \vec{f}_1)^T\end{bmatrix}$$

You can identify the terms in Möller–Trumbore's algorithm with the terms from above, albeit I would altogether just recommend skipping Möller–Trumbore, and giving your students a 3x3 matrix inverse procedure and a matrix-vector multiplication procedure which should be enough to find the vector $$\begin{bmatrix} -t & \beta & \gamma \end{bmatrix}^T$$.

Once $$t, \beta, \gamma$$ are computed one needs to check that this is indeed a valid intersection (i.e. it belongs to both $$\mathcal{R}$$ and $$\mathcal{T}$$):

$$(0\leq t < \infty) \land (0\leq \beta) \land (0\leq \gamma) \land (\beta+\gamma\leq 1).$$

In C++ the above check would fail also when the matrix is not invertible, since NaNs and infinities will fail this check, so no additional code is necessary for the special cases. It is instructive however to consider what the special cases represent. One option is that $$\vec{f}_1$$ and $$\vec{f}_2$$ are parallel/linearly dependent - that means that we have a degenerate triangle (it's rather a line segment in this case). Another option is that $$\vec{d}, \vec{f}_1, \vec{f}_2$$ are coplanar/linearly dependent. This may occur if $$\vec{d}$$ is parallel to the plane spanned by $$\vec{f}_1,\vec{f}_2$$. Then if $$\vec{o}$$ doesn't lie in the plane there is no intersection. If $$\vec{o}$$ lies in the plane, then there are infinitely many intersections, and the closest one is at $$t=0$$. For convenience one can treat this as a non-intersection however, since the ray is inside the plane, which is infinitely thin.

2) Split Solution

The other solution method that I like and which has a nice geometric interpretation is the one where we first intersect the plane in which the triangle lies, and only after do we find the $$\beta, \gamma$$ coordinates. This is loosely equivalent to solving for the $$t$$ variable in the system, and then substituting back in order to get a $$2\times 2$$ system. We know that the plane has a normal $$\vec{n} = \vec{f}_1\times\vec{f}_2$$ and that the set of points on the plane are:

$$\mathcal{P} = \{\vec{p}\,:\, (\vec{p}-\vec{v}_0)\cdot \vec{n} = 0\}.$$

Then plugging in $$\vec{r}(t)$$ for $$\vec{p}$$ yields:

$$(\vec{o}+t\vec{d}-\vec{v}_0)\cdot \vec{n} = 0 \implies t = \frac{(\vec{v}_0-\vec{o})\cdot\vec{n}}{\vec{d}\cdot \vec{n}}.$$

Note also that this is exactly what you'll get using Cramer's rule for $$t$$ using the $$3\times 3$$ system. Provided that $$t \in [0,\infty)$$ we can use the fact that we know the intersection point on the plane: $$\vec{p} = \vec{o}+t\vec{d}$$. We are now left with 3 equations but two unknowns:

$$\vec{p} = \vec{v}_0 + \beta\vec{f}_1 + \gamma\vec{f}_2.$$

Ideally we would like to reduce this to a system of $$2$$ equations, this can be done by dotting the equations with $$\vec{f}_1$$ and $$\vec{f}_2$$:

$$\vec{f}_1\cdot(\vec{p}-\vec{v}_0) = \beta\vec{f}_1\cdot\vec{f}_1 + \gamma\vec{f}_1\cdot\vec{f}_2 \\ \vec{f}_2\cdot(\vec{p}-\vec{v}_0) = \beta\vec{f}_2\cdot\vec{f}_1 + \gamma\vec{f}_2\cdot\vec{f}_2.$$

It's also instructive to look at what this looks like in terms of matrix notation:

$$F = \begin{bmatrix} \vec{f}_1 & \vec{f}_2\end{bmatrix} \\ \vec{p}-\vec{v}_0 = F \begin{bmatrix} \beta \\ \gamma \end{bmatrix} \\ F^T(\vec{p}-\vec{v}_0) = F^TF\begin{bmatrix} \beta \\ \gamma \end{bmatrix}.$$

$$\begin{bmatrix} \beta \\ \gamma \end{bmatrix} = (F^TF)^{-1}F^T(\vec{p}-\vec{v}_0),$$

and one must check that $$(0\leq \beta) \land (0\leq \gamma) \land (\beta+\gamma\leq 1)$$ as before.

2.1) Algebraic and Geometric Meaning

These are additional mathematically interesting details, that are probably best left for the very end.

The matrix $$F^TF$$ is a $$2\times 2$$ matrix that is trivial to invert and has the following form:

$$F^TF = \begin{bmatrix} \vec{f}_1\cdot\vec{f}_1 & \vec{f}_1\cdot\vec{f}_2 \\ \vec{f}_1\cdot\vec{f}_2 & \vec{f}_2\cdot\vec{f}_2\end{bmatrix}.$$

The matrix $$G=F^TF$$ is in fact the metric tensor on the plane. To motivate the metric tensor you can describe it as a tool to measure angles and lengths in a non-orthonormal coordinate system. Given two vectors $$\vec{u}$$ and $$\vec{w}$$ in the plane described through $$\vec{u}_P = \begin{bmatrix} \beta_u & \gamma_u\end{bmatrix}^T$$ and $$\vec{w}_P = \begin{bmatrix} \beta_w & \gamma_w\end{bmatrix}^T$$ the metric tensor can be used to compute the equivalent of the dot product between those:

$$\vec{u}_P^TG\vec{w}_P = \vec{w}_P^TG\vec{u}_P = \vec{u}\cdot\vec{w} = \|\vec{u}\|\|\vec{w}\|\cos\angle(\vec{u},\vec{w}).$$

In the above $$\vec{u}, \vec{w}$$ are the representations of the vectors with respect to the canonical basis $$(1,0,0),\, (0,1,0),\, (0,0,1)$$. It is also clear that for an orthonormal basis $$\vec{f}_1, \vec{f}_2$$, the metric tensor is the identity matrix, recovering the usual dot product. In fact any inner product in a finite-dimensional real space has the form: $$\vec{u}^TG\vec{w}$$ for $$G$$ symmetric and positive definite. So this can be used to introduce both the concept of metric tensor and inner product with a suitable motivation. Using the definition for length: $$\|\vec{u}_P\|^2 = \vec{u}_P^TG\vec{u}_P = \vec{u}^T\vec{u}$$ it is also clear that the metric tensor helps with measuring length in non-orthonormal coordinate systems.

Another interpretation of the system is that:

$$F^T(\vec{p}-\vec{v}_0) = F^TF\begin{bmatrix} \beta \\ \gamma \end{bmatrix}$$

are the normal equations, which are used when we have an overdetermined system (more linearly independent equations than unknowns, note that $$F$$ is rectangular with more rows than columns). In fact even if $$\vec{p}$$ was not in the plane of the triangle, the above would find the $$\beta, \gamma$$ for its projection on this plane, since the solution of the normal equations is a minimizer of:

$$\left\|\vec{v}_0 + F\begin{bmatrix} \beta \\ \gamma \end{bmatrix} - \vec{p}\right\|^2.$$

Moreover in:

$$\begin{bmatrix} \beta \\ \gamma \end{bmatrix} = (F^TF)^{-1}F^T(\vec{p}-\vec{v}_0),$$

the matrix $$(F^TF)^{-1}F^T$$ is known as the Moore–Penrose (left pseudo-)inverse. So there are in fact many algebraic and geometric interpretations that can be made in order to connect this problem to topics that the students have studied or are about to study. It's a good motivational tool.