A Beautiful Ray/Triangle Intersection Method

Tavian Barnes Comments

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 AA, BB, and CC. If your triangle is not degenerate, those three points define a plane. Writing the vector from AA to BB as AB=BA\overrightarrow{AB} = \overrightarrow{B} - \overrightarrow{A}, and similarly for AC\overrightarrow{AC}, the normal vector perpendicular to that plane is given by

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

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

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

We can make this space even nicer by translating AA to the origin, giving

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

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

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

P=[ABACNA0001]1.P = \begin{bmatrix} \overrightarrow{AB} & \overrightarrow{AC} & \overrightarrow{N} & \overrightarrow{A} \\ 0 & 0 & 0 & 1 \end{bmatrix}^{-1}.

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

[ounuovnvownw1]=P[oxnxoynyoznz1][nunvnw0]=P[nxnynz0]\begin{aligned} \begin{bmatrix} \mathrlap{o_u}\hphantom{n_u} \\ \mathrlap{o_v}\hphantom{n_v} \\ \mathrlap{o_w}\hphantom{n_w} \\ 1 \end{bmatrix} & = P\, \begin{bmatrix} \mathrlap{o_x}\hphantom{n_x} \\ \mathrlap{o_y}\hphantom{n_y} \\ \mathrlap{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{aligned}

Solving for the intersection point is easy:

t=ow/nwu=ou+tnuv=ov+tnv\begin{aligned} t & = -o_w/n_w \\ u & = o_u + t\,n_u \\ v & = o_v + t\,n_v \end{aligned}

And as above, this can be quickly checked with u,v0u, v \ge 0, u+v1u + v \le 1 (and t0t \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 PP, 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.



I'm a bit confused about the "cost of 1 division, 2 multiplications, 2 or 3 additions, and some comparisons". Don't you still have to multiply each ray by the matrix P? So even if P is precomputed and stored ahead of time, that's still 16 multiplications and 12 additions, plus the computations of t,u, and v.

Tavian Barnes

That's in the case that "objects are already associated with a transformation matrix," in which case you have to do a matrix-ray multiplication anyway, so it doesn't count against you. Also I think you've underestimated the cost of a matrix-ray multiplication, it ought to be 18 multiplications and 15 additions.


Hi! I wonder what is nw is 0? I watched sources from the dimension repo, that case is not handled at all:

Could you please explain what to do if nw is 0?

Tavian Barnes

Geometrically, if nwn_w is zero, then the line is exactly parallel to the plane of the triangle. For this degenerate case, it's easiest to just return "no intersection." There are a few sub-cases:

  • ow=0o_w = 0: Then we'll have t=NaNt = \mathrm{NaN}, so t0t \ge 0 fails and we return false
  • ow>0o_w > 0: We'll have t=t = -\infty, so again t0t \ge 0 fails and we return false
  • ow<0o_w < 0: We get t=t = \infty. Depending on nun_u and nvn_v,
    • nu<0n_u < 0 or nv<0n_v < 0: Then u=u = -\infty, v=v = -\infty, so at least one of u0u \ge 0 or v0v \ge 0 fails and we return false
    • nu=0n_u = 0 or nv=0n_v = 0: Results in u=NaNu = \mathrm{NaN} or v=NaNv = \mathrm{NaN}, so again one of u0u \ge 0 or v0v \ge 0 fails and we return false
    • nu>0n_u > 0 and nv>0n_v > 0: That means u=v=u = v = \infty, so u+v1u + v \le 1 fails and we return false

Without explicitly handling the nw=0n_w = 0 case, the properties of IEEE floating point arithmetic have given us the desired behaviour in all cases.