<p>The University of Washington, Computer Science department teaches a ray, triangle intersection algorithm that is fairly intuitive, has good documentation <a href="https://courses.cs.washington.edu/courses/csep557/10au/lectures/triangle_intersection.pdf" rel="nofollow noreferrer">(they did their own write up here)</a>, 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.</p>

CC BY-SA 4.0
2022-03-05T13:42:01.183Z
2022-03-05T13:42:01.183Z
pmw1234
<p>If you want something that looks more vectorial, start with doing a ray-plane intersection; then you can get the barycentric coordinates as follows:</p>
<p>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).</p>
<p>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.</p>
<p>(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.)</p>

CC BY-SA 4.0
2022-03-04T18:19:08.967Z
2022-03-04T18:19:08.967Z
Nathan Reed
<p>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.</p>
<p>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.</p>
<p><strong>0.1) Parametric Ray Definition</strong></p>
<p>A parametric definition of a line given an origin <span class="math-container">$\vec{o}$</span> and a direction <span class="math-container">$\vec{d}$</span> is:</p>
<p><span class="math-container">$$\vec{r}(t) = \vec{o} + t\vec{d}.$$</span></p>
<p>A ray is the set of points on one side of the origin of the line:</p>
<p><span class="math-container">$$\mathcal{R} = \{\vec{r}(t) = \vec{o} + t\vec{d}\,:\,t\in[0,\infty)\}.$$</span></p>
<p>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 <span class="math-container">$t$</span> vaues (e.g. 0, 1, 0.5, but also note that -1 is on the line but not part of the ray).</p>
<p><strong>0.2) Parametric Triangle Definition</strong></p>
<p>A parametric definition of a plane given a point on it <span class="math-container">$\vec{v}_0$</span> and two vectors <span class="math-container">$\vec{f}_1 = \vec{v}_1 - \vec{v}_0$</span>, <span class="math-container">$\vec{f}_2 = \vec{v}_2-\vec{v}_0$</span> is:</p>
<p><span class="math-container">$$\vec{v}(\beta, \gamma) = \vec{v}_0 + \beta\vec{f}_1 + \gamma\vec{f}_2$$</span></p>
<p>The pair <span class="math-container">$(\beta, \gamma)$</span> are the coordinates of <span class="math-container">$\vec{v}$</span> in the 2D coordinate system with origin <span class="math-container">$\vec{v}_0$</span> and basis vectors <span class="math-container">$\vec{f}_1, \vec{f}_2$</span>. For those to be basis vectors we additionally require that <span class="math-container">$\vec{f}_1, \vec{f}_2$</span> are not parallel/linearly dependent (i.e. the triangle is non-degenerate).</p>
<p>You should convince the students that <span class="math-container">$\vec{v}(\beta, \gamma)$</span> indeed spans the whole plane for different <span class="math-container">$(\beta, \gamma)$</span>. For instance you could draw the points <span class="math-container">$(0,0), (1,0), (0,1), (1,1)$</span> which are the vertices of a parallelogram. Additionally, the points forming the 3 edges of the triangle are: <span class="math-container">$(0,\gamma); \gamma \in[0,1]$</span>, <span class="math-container">$(\beta, 0); \beta \in [0,1]$</span>, and <span class="math-container">$(\beta, \gamma); \beta \geq 0, \gamma \geq 0, \beta + \gamma = 1$</span>. Then the set of points forming the triangle are:</p>
<p><span class="math-container">$$\mathcal{T} = \{\vec{v}(\beta, \gamma)\,:\,0\leq \beta, 0\leq \gamma, \beta+\gamma \leq 1\}.$$</span></p>
<p>This should definitely be accompanied with a drawing of the geometry so that the students may grasp its meaning.</p>
<p><strong>0.2.1) Barycentric Coorindates</strong></p>
<p>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.</p>
<p>The coordinates <span class="math-container">$\beta, \gamma$</span> are in fact the barycentric coordinates of <span class="math-container">$\vec{v}$</span> associated with <span class="math-container">$\vec{v}_1, \vec{v}_2$</span>. In fact we can rewrite our formulation as:</p>
<p><span class="math-container">$$\begin{align}\vec{v}(\beta, \gamma) &amp;= \vec{v}_0 + \beta \vec{f}_1 + \gamma \vec{f}_2 \\ \vec{v}(\beta, \gamma) &amp;= \alpha\vec{v}_0 + \beta\vec{v}_1 + \gamma\vec{v}_2, \alpha = 1-\beta -\gamma,\end{align}$$</span></p>
<p>where <span class="math-container">$\alpha, \beta, \gamma$</span> are the barycentric coordinates of point <span class="math-container">$\vec{v}$</span> with respect to the triangle <span class="math-container">$\vec{v}_0, \vec{v}_1, \vec{v}_2$</span>. The choice of <span class="math-container">$(\beta, \gamma)$</span> over <span class="math-container">$(\alpha, \beta)$</span> or <span class="math-container">$(\alpha, \gamma)$</span> was arbitrary. Indeed using <span class="math-container">$\alpha+\beta+\gamma=1$</span> I could have equivalently written:</p>
<p><span class="math-container">$$\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).
$$</span></p>
<p>However, the fact that <span class="math-container">$\alpha, \beta, \gamma$</span> 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.</p>
<p><strong>0.3) Problem Formulation</strong></p>
<p>We want to find the closest intersection (the one with smallest <span class="math-container">$t$</span>) between <span class="math-container">$\mathcal{R}$</span> and <span class="math-container">$\mathcal{T}$</span>. Let <span class="math-container">$\vec{p}$</span> be an intersection point, i.e. <span class="math-container">$\vec{p} \in \mathcal{R}$</span> and <span class="math-container">$\vec{p} \in \mathcal{T}$</span>, then there exist <span class="math-container">$t\geq 0, \beta\geq 0, \gamma\geq 0, \beta+\gamma\leq 1$</span> such that:</p>
<p><span class="math-container">$$\vec{o}+t\vec{d} = \vec{p} = \vec{v}_0 + \beta\vec{f}_1 + \gamma\vec{f}_2.$$</span></p>
<p>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 <span class="math-container">$3\times 3$</span> system. For example Möller–Trumbore is just a special case of a specific solution order for the unknowns <span class="math-container">$t, \beta, \gamma$</span>.</p>
<p><strong>1) Simultaneous Solution</strong></p>
<p>Let us rewrite the equations in matrix form:</p>
<p><span class="math-container">$$\vec{o}+t\vec{d} = \vec{v}_0 + \beta\vec{f}_1 + \gamma\vec{f}_2$$</span>
<span class="math-container">$$\begin{bmatrix}\vec{d} &amp; \vec{f}_1 &amp; \vec{f}_2 \end{bmatrix} \begin{bmatrix} -t \\ \beta \\ \gamma \end{bmatrix} = \vec{o} - \vec{v}_0.$$</span></p>
<p>It should be made clear for the students that <span class="math-container">$M = \begin{bmatrix}\vec{d} &amp; \vec{f}_1 &amp; \vec{f}_2 \end{bmatrix}$</span> is a 3x3 matrix, didactically, being explicit there may help:</p>
<p><span class="math-container">$$M = \begin{bmatrix}d_x &amp; f_{1,x} &amp; f_{2,x} \\  d_y &amp; f_{1,y} &amp; f_{2,y} \\ d_z &amp; f_{1,z} &amp; f_{2,z}\end{bmatrix}.$$</span></p>
<p>Solving the system is as simple as:</p>
<p><span class="math-container">$$\begin{bmatrix} -t \\ \beta \\ \gamma \end{bmatrix} = M^{-1}(\vec{o}-\vec{v}_0)$$</span></p>
<p>Either <span class="math-container">$3\times 3$</span> 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 <span class="math-container">$3\times 3$</span> matrix inversion that some people like to use (which is elegant, but didactically not very useful):</p>
<p><span class="math-container">$$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}$$</span></p>
<p>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 <span class="math-container">$\begin{bmatrix} -t &amp; \beta &amp; \gamma \end{bmatrix}^T$</span>.</p>
<p>Once <span class="math-container">$t, \beta, \gamma$</span> are computed one needs to check that this is indeed a valid intersection (i.e. it belongs to both <span class="math-container">$\mathcal{R}$</span> and <span class="math-container">$\mathcal{T}$</span>):</p>
<p><span class="math-container">$$(0\leq t &lt; \infty) \land (0\leq \beta) \land (0\leq \gamma) \land (\beta+\gamma\leq 1).$$</span></p>
<p>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 <span class="math-container">$\vec{f}_1$</span> and <span class="math-container">$\vec{f}_2$</span> 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 <span class="math-container">$\vec{d}, \vec{f}_1, \vec{f}_2$</span> are coplanar/linearly dependent. This may occur if <span class="math-container">$\vec{d}$</span> is parallel to the plane spanned by <span class="math-container">$\vec{f}_1,\vec{f}_2$</span>. Then if <span class="math-container">$\vec{o}$</span> doesn't lie in the plane there is no intersection. If <span class="math-container">$\vec{o}$</span> lies in the plane, then there are infinitely many intersections, and the closest one is at <span class="math-container">$t=0$</span>. For convenience one can treat this as a non-intersection however, since the ray is inside the plane, which is infinitely thin.</p>
<p><strong>2) Split Solution</strong></p>
<p>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 <span class="math-container">$\beta, \gamma$</span> coordinates. This is loosely equivalent to solving for the <span class="math-container">$t$</span> variable in the system, and then substituting back in order to get a <span class="math-container">$2\times 2$</span> system. We know that the plane has a normal <span class="math-container">$\vec{n} = \vec{f}_1\times\vec{f}_2$</span> and that the set of points on the plane are:</p>
<p><span class="math-container">$$\mathcal{P} = \{\vec{p}\,:\, (\vec{p}-\vec{v}_0)\cdot \vec{n} = 0\}.$$</span></p>
<p>Then plugging in <span class="math-container">$\vec{r}(t)$</span> for <span class="math-container">$\vec{p}$</span> yields:</p>
<p><span class="math-container">$$(\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}}.$$</span></p>
<p>Note also that this is exactly what you'll get using Cramer's rule for <span class="math-container">$t$</span> using the <span class="math-container">$3\times 3$</span> system. Provided that <span class="math-container">$t \in [0,\infty)$</span> we can use the fact that we know the intersection point on the plane: <span class="math-container">$\vec{p} = \vec{o}+t\vec{d}$</span>. We are now left with 3 equations but two unknowns:</p>
<p><span class="math-container">$$\vec{p} = \vec{v}_0 + \beta\vec{f}_1 + \gamma\vec{f}_2.$$</span></p>
<p>Ideally we would like to reduce this to a system of <span class="math-container">$2$</span> equations, this can be done by dotting the equations with <span class="math-container">$\vec{f}_1$</span> and <span class="math-container">$\vec{f}_2$</span>:</p>
<p><span class="math-container">$$\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.$$</span></p>
<p>It's also instructive to look at what this looks like in terms of matrix notation:</p>
<p><span class="math-container">$$F = \begin{bmatrix} \vec{f}_1 &amp; \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}.$$</span></p>
<p>Then the solution reads:</p>
<p><span class="math-container">$$\begin{bmatrix} \beta \\ \gamma \end{bmatrix} = (F^TF)^{-1}F^T(\vec{p}-\vec{v}_0),$$</span></p>
<p>and one must check that <span class="math-container">$(0\leq \beta) \land (0\leq \gamma) \land (\beta+\gamma\leq 1)$</span> as before.</p>
<p><strong>2.1) Algebraic and Geometric Meaning</strong></p>
<p>These are additional mathematically interesting details, that are probably best left for the very end.</p>
<p>The matrix <span class="math-container">$F^TF$</span> is a <span class="math-container">$2\times 2$</span> matrix that is trivial to invert and has the following form:</p>
<p><span class="math-container">$$F^TF = \begin{bmatrix} \vec{f}_1\cdot\vec{f}_1 &amp; \vec{f}_1\cdot\vec{f}_2 \\ \vec{f}_1\cdot\vec{f}_2 &amp; \vec{f}_2\cdot\vec{f}_2\end{bmatrix}.$$</span></p>
<p>The matrix <span class="math-container">$G=F^TF$</span> 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 <span class="math-container">$\vec{u}$</span> and <span class="math-container">$\vec{w}$</span> in the plane described through <span class="math-container">$\vec{u}_P = \begin{bmatrix} \beta_u &amp; \gamma_u\end{bmatrix}^T$</span> and <span class="math-container">$\vec{w}_P = \begin{bmatrix} \beta_w &amp; \gamma_w\end{bmatrix}^T$</span> the metric tensor can be used to compute the equivalent of the dot product between those:</p>
<p><span class="math-container">$$\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}).$$</span></p>
<p>In the above <span class="math-container">$\vec{u}, \vec{w}$</span> are the representations of the vectors with respect to the canonical basis <span class="math-container">$(1,0,0),\, (0,1,0),\, (0,0,1)$</span>. It is also clear that for an orthonormal basis <span class="math-container">$\vec{f}_1, \vec{f}_2$</span>, 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: <span class="math-container">$\vec{u}^TG\vec{w}$</span> for <span class="math-container">$G$</span> 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: <span class="math-container">$\|\vec{u}_P\|^2 = \vec{u}_P^TG\vec{u}_P = \vec{u}^T\vec{u}$</span> it is also clear that the metric tensor helps with measuring length in non-orthonormal coordinate systems.</p>
<p>Another interpretation of the system is that:</p>
<p><span class="math-container">$$F^T(\vec{p}-\vec{v}_0) = F^TF\begin{bmatrix} \beta \\ \gamma \end{bmatrix}$$</span></p>
<p>are the normal equations, which are used when we have an overdetermined system (more linearly independent equations than unknowns, note that <span class="math-container">$F$</span> is rectangular with more rows than columns). In fact even if <span class="math-container">$\vec{p}$</span> was not in the plane of the triangle, the above would find the <span class="math-container">$\beta, \gamma$</span> for its projection on this plane, since the solution of the normal equations is a minimizer of:</p>
<p><span class="math-container">$$\left\|\vec{v}_0 + F\begin{bmatrix} \beta \\ \gamma \end{bmatrix} - \vec{p}\right\|^2.$$</span></p>
<p>Moreover in:</p>
<p><span class="math-container">$$\begin{bmatrix} \beta \\ \gamma \end{bmatrix} = (F^TF)^{-1}F^T(\vec{p}-\vec{v}_0),$$</span></p>
<p>the matrix <span class="math-container">$(F^TF)^{-1}F^T$</span> 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.</p>

CC BY-SA 4.0
2022-03-04T13:36:46.300Z
2022-03-04T20:06:20.200Z
2022-03-04T20:06:20.200Z
lightxbulb
<p>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.</p>
<p>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.</p>
<p>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.</p>
<p>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?</p>

CC BY-SA 4.0
2022-03-07T12:01:26.683Z
The *Scratchapixel* tutorials build up the necessary background on barycentric coordinates and provides a framework building up to [Möller–Trumbore](https://www.scratchapixel.com/lessons/3d-basic-rendering/ray-tracing-rendering-a-triangle/moller-trumbore-ray-triangle-intersection). The chapter begins with less efficient - though perhaps more immediately intuitive ray-triangle intersection, using edge test equations.
Brett Hale
<rendering><geometry><ray>
Student-friendly ray-triangle intersection
TomKern
# Student-friendly ray-triangle intersection

(Asked by: username_0 on 2022-03-03)

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?

**Tags:** rendering, geometry, ray

> **Comments:**
> **username_1**: The *Scratchapixel* tutorials build up the necessary background on barycentric coordinates and provides a framework building up to [Möller–Trumbore](https://www.scratchapixel.com/lessons/3d-basic-rendering/ray-tracing-rendering-a-triangle/moller-trumbore-ray-triangle-intersection). The chapter begins with less efficient - though perhaps more immediately intuitive ray-triangle intersection, using edge test equations.  

---

## 3 answers

### Answer by username_2 on 2022-03-05. Score: 1

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)](https://courses.cs.washington.edu/courses/csep557/10au/lectures/triangle_intersection.pdf), 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.


---

### Answer by username_3 on 2022-03-04. Score: 1

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


---

### Answer by username_4 on 2022-03-04. Score: 0

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}.$$

Then the solution reads:

$$\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.


---

