3D ray/triangle intersections are obviously an important part of much of computer graphics. The Möller–Trumbore algorithm, for example, computes these intersections very quickly. But there is another method that I believe is more elegant, and in some cases allows you to compute the intersection for “free.”

Say your triangle has corners at \(A\), \(B\), and \(C\). If your triangle is not degenerate, those three points define a plane. Writing the vector from \(A\) to \(B\) as \(\overrightarrow{AB} = \overrightarrow{B} - \overrightarrow{A}\), and similarly for \(\overrightarrow{AC}\), the normal vector perpendicular to that plane is given by

$$

\overrightarrow{N} = \overrightarrow{AB} \times \overrightarrow{AC}.

$$

Since \(\overrightarrow{AB}\), \(\overrightarrow{AC}\), and \(\overrightarrow{N}\) are all linearly independent, they form an alternative basis for all of 3-space. That is, any point \(\vec{p}\) can be represented by a unique triple \(\langle u,v,w \rangle\) such that

$$

\vec{p} = u\,\overrightarrow{AB} + v\,\overrightarrow{AC} + w\,\overrightarrow{N}.

$$

We can make this space even nicer by translating \(A\) to the origin, giving

$$

\vec{p} = u\,\overrightarrow{AB} + v\,\overrightarrow{AC} + w\,\overrightarrow{N} + \overrightarrow{A}.

$$

In this space, the three corners of the triangle are at \(\langle 0,0,0 \rangle\), \(\langle 1,0,0 \rangle\), and \(\langle 0,1,0 \rangle\). It is easy to see that a point is within the triangle if and only if \(u,v \ge 0\), \(u + v \le 1\), and \(w = 0\).

To transform a point \((x,y,z)\) into this nicer space, we can use an affine change of basis matrix

$$

P =

\begin{bmatrix}

\overrightarrow{AB} & \overrightarrow{AC} & \overrightarrow{N} & \overrightarrow{A} \\

0 & 0 & 0 & 1

\end{bmatrix}^{-1}.

$$

A line defined by \(\overrightarrow{o} + t\,\overrightarrow{n}\) can thus be transformed into the new space by

$$

\begin{align*}

\begin{bmatrix}

\rlap{o_u}\hphantom{n_u} \\

\rlap{o_v}\hphantom{n_v} \\

\rlap{o_w}\hphantom{n_w} \\

1

\end{bmatrix}

& =

P\,

\begin{bmatrix}

\rlap{o_x}\hphantom{n_x} \\

\rlap{o_y}\hphantom{n_y} \\

\rlap{o_z}\hphantom{n_z} \\

1

\end{bmatrix} \\

\begin{bmatrix}

n_u \\

n_v \\

n_w \\

0

\end{bmatrix}

& =

P\,

\begin{bmatrix}

n_x \\

n_y \\

n_z \\

0

\end{bmatrix}

\end{align*}

$$

Solving for the intersection point is easy:

$$

\begin{align*}

t & = -o_w/n_w \\

u & = o_u + t\,n_u \\

v & = o_v + t\,n_v

\end{align*}

$$

And as above, this can be quickly checked with \(u, v \ge 0\), \(u + v \le 1\) (and \(t \ge 0\) to only count intersections “in front” of the ray).

The beauty of this method is that often, objects are already associated with a transformation matrix, so if you left-multiply that matrix by \(P\), the change of basis is performed for free at the same time as other transformations, and only the short block above is needed to actually perform the intersection test. At a cost of 1 division, 2 multiplications, 2 or 3 additions, and some comparisons, that's about as fast as possible without special-casing triangles versus other objects.

Of course, if you're working with lots of triangles (a mesh, for example), you can save memory and time by not associating a transformation matrix with each triangle. In that case, other algorithms such as the one linked above may be faster as they can avoid the two matrix multiplications.

An implementation of this method can be seen here in my ray tracer Dimension.