Ray-Triangle Intersection and the Watertight Edge Test

A ray tracer spends most of its time answering one question: where, along this parametric line, does it meet a surface. Triangles are the surface of choice because they are planar, convex, and defined by exactly three points, so the intersection reduces to a linear system with three unknowns. The system is small enough to solve directly, and that direct solution is where both the speed and the numerical troubles of the whole renderer originate.

The linear system behind the test

Represent the ray as $o + t d$ with $t \ge 0$, and the triangle as the affine combination $v_0 + u(v_1 - v_0) + v(v_2 - v_0)$ with $u \ge 0$, $v \ge 0$, $u + v \le 1$. Setting them equal gives a vector equation in three unknowns:

$$o + t d = v_0 + u e_1 + v e_2$$

with $e_1 = v_1 - v_0$ and $e_2 = v_2 - v_0$. Rearranged into matrix form, the columns are $-d$, $e_1$, $e_2$ and the right-hand side is $o - v_0$. Cramer’s rule expresses each unknown as a ratio of two scalar triple products. Writing $p = d \times e_2$ and $s = o - v_0$, the three quantities come out as

$$u = \frac{(s \cdot p)}{(e_1 \cdot p)}, \qquad v = \frac{(d \cdot (s \times e_1))}{(e_1 \cdot p)}, \qquad t = \frac{(e_2 \cdot (s \times e_1))}{(e_1 \cdot p)}.$$

Every term is a triple product, which is a determinant, which means the denominator $(e_1 \cdot p)$ is proportional to the scalar triple product of $d$, $e_1$, $e_2$ — the signed volume of the parallelepiped those three vectors span. Detection of near-parallel rays and rejection of back faces come from that single number, and no plane equation needs to be precomputed or stored. The triangle costs three vertices and nothing else, which is why the formulation in Möller and Trumbore’s 1997 note became the default in production kernels: 32 bytes of triangle data on a 4-byte-float representation, no per-triangle normal, no plane offset.

The test is also cheap in the ordered sense that matters. If u fails its bounds test, the kernel can return immediately without computing v or t. Ordering the early outs is a real decision, not a cosmetic one: testing u first rejects roughly as many triangles as testing v first on evenly distributed geometry, but the u + v > 1 test can only run after both have been computed, so the first two bounds checks must be independent.

#include <math.h>

typedef struct { float x, y, z; } vec3;

static inline vec3 sub(vec3 a, vec3 b) { return (vec3){a.x - b.x, a.y - b.y, a.z - b.z}; }
static inline float dot(vec3 a, vec3 b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
static inline vec3 cross(vec3 a, vec3 b) {
    return (vec3){a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
}

/* Ray o + t*d against triangle (v0, v1, v2). Returns 1 on a front- or back-face hit. */
int ray_triangle(vec3 o, vec3 d, vec3 v0, vec3 v1, vec3 v2,
                 float *t_out, float *u_out, float *v_out)
{
    vec3 e1 = sub(v1, v0), e2 = sub(v2, v0);
    vec3 p  = cross(d, e2);
    float det = dot(e1, p);
    if (fabsf(det) < 1e-12f) return 0;   /* ray lies in the triangle plane */
    float inv = 1.0f / det;              /* sign of det carries the facing */
    vec3 s = sub(o, v0);
    float u = dot(s, p) * inv;
    if (u < 0.0f || u > 1.0f) return 0;
    vec3 q = cross(s, e1);
    float v = dot(d, q) * inv;
    if (v < 0.0f || u + v > 1.0f) return 0;
    float t = dot(e2, q) * inv;
    if (t <= 0.0f) return 0;             /* behind the origin */
    *t_out = t; *u_out = u; *v_out = v;
    return 1;
}

The constant 1e-12f in that snippet is the part that should not be copied. A determinant threshold has units of length cubed when the ray direction is not normalised: it scales with the triangle’s area and with the length of the direction vector. A build that normalises directions can use a threshold in area units; a build that stores unnormalised directions, as most traversal kernels do because dividing by the length costs more than the error it removes, needs a threshold derived from the geometry’s own scale. The robust choice is a relative test — compare the determinant against an epsilon multiplied by the magnitudes of the three edge terms — and to keep the absolute epsilon only as a floor for degenerate triangles whose area really is zero.

Why the naive test leaks light at shared edges

Tessellated surfaces share vertices, and a shared vertex is evaluated twice: once as corner $v_0$ of one triangle, once as corner $v_2$ of its neighbour. Floating-point arithmetic is not associative across different orderings, so the two evaluations compute slightly different values of the same geometric quantity. If the ray passes within a few units in the last place of the shared edge, one triangle can compute u + v as a hair above 1 and reject, while the other computes it a hair above 1 as well but in a different parameterisation and rejects too. Nothing is hit, and the pixel or shadow ray sees through a solid surface.

The symptom is specific and recognisable: single-pixel pinholes scattered along polygon boundaries, more frequent on large triangles far from the origin, and entirely absent in the interior of any triangle. Shadow rays make it visible because a missed shadow ray means a lit point where the surface should be occluded, producing bright specks along the seams of a tessellated mesh. Adding an epsilon to the barycentric bounds is the usual first fix: accept u > -eps, u < 1 + eps, v > -eps, u + v < 1 + eps, with eps around $10^{-6}$ relative to the triangle’s scale. That closes most holes, and it is not watertight, because the epsilon band is a fixed distance in barycentric space while the disagreement between two evaluations scales with the coordinates of the vertices, not with the barycentric magnitude.

Woop, Benthin, and Wald’s 2013 formulation removes the disagreement instead of widening the tolerance. The construction translates the ray so the origin is at zero (subtracting the vertices, a single translation that is exact in the sense that it is applied identically to every triangle sharing a vertex), then permutes the axes so that the largest-magnitude component of the direction is the $z$ axis and is positive. Permutation is exact — it moves bits, it does not round them — so any two triangles that share an edge see the same vertex coordinates in the same order after the transform. The intersection then becomes a 2D problem: project onto the $xy$ plane and test the ray’s 2D line against the three 2D edges with a sign-consistent edge function.

The rotated direction is $(d_x, d_y, d_z)$ sheared so that $d_z = 1$, giving a 2D ray with direction $(d_x/d_z, d_y/d_z)$ from the sheared origin. The edge functions are

$$E_i = (x_i - o_x)(y_{i+1} - o_y) - (y_i - o_y)(x_{i+1} - o_x),$$

computed with the two endpoints in a fixed order determined only by their labels, never by which triangle is being tested. When the three $E_i$ share a sign, the point is inside. Two triangles that share an edge compute that edge’s function from the same two vertices in the same order with the same origin, so they cannot disagree — bit-identical if the compiler emits the same instruction sequence, and identical in the mathematical sense regardless. The cost is that the test is conservative rather than exact: a ray grazing exactly along an edge can be reported as hitting both triangles, so a hit can fall marginally outside the triangle’s true area. For rendering that is the correct direction of error, since a double hit at a shared edge shades the same surface point twice with a consistent normal and produces no visible artifact, whereas a miss produces a hole.

Watertightness is not exactness

It is worth being precise about what the property guarantees. A watertight test guarantees that for any ray and any closed, watertight mesh, the sequence of intersections along the ray alternates between entering and leaving the surface: the parity is correct. It does not guarantee that the reported $t$ equals the mathematically exact intersection, and it does not guarantee monotone ordering of $t$ values for nearly coincident surfaces. Those are different properties, and only the first is needed for shadow rays, for inside/outside tests in solid geometry, and for a path that must not terminate inside a wall.

The parametric $t$ reported by the watertight form is derived from the sheared coordinates, so it is computed in the permuted frame and must be un-sheared. The standard construction recovers it as a ratio of the edge function evaluated at the untransformed $z$ coordinate and the direction’s $z$ component, which keeps the division by d_z — already known to be the largest component, so the division is well conditioned. Via such care the far intersection point stays within a small relative error of the exact one, though no floating-point scheme reduces the error to zero; that residual is discussed in solving quadratic intersections without losing the far root for curved primitives, where the conditioning problem is worse.

The trade-offs against the direct Möller–Trumbore form are concrete. The watertight test replaces three triple products with three 2D edge functions, which is fewer multiplies but adds an axis permutation and a shear division that depend on the ray, so the work per ray is amortised across all triangles tested rather than per triangle. It needs the full set of vertices in a canonical order, which constrains mesh layouts: a compressed triangle format that packs vertices into 16-bit quantised coordinates must reproduce the same permutation result on both sides of every shared edge, and a format that deduplicates vertices by welding them with a tolerance can break the guarantee if the welding itself is not shared. And it changes the accepted hit set: because the test is conservative, code that depends on u and v for attribute interpolation must clamp them to the triangle rather than trusting them, since a point just outside an edge can now be reported. Clamping is one multiply and two compares, which is why the conservative behaviour is affordable.

Where the test sits in the wider kernel

Every intersection kernel is a function of a ray and a candidate, and both of those come from structures that decide which candidates are worth testing. The slab test for boxes prunes whole subtrees before any triangle is touched, and its correctness requirement — never reject a box that the enclosed triangle would hit — is the mirror image of the watertight requirement in the triangle test: one must not lose hits, the other must not lose triangles. Both constraints are about intervals, and the ray interval itself, including its lower bound and the origin offset applied after a hit, is the subject of ray parameter intervals and the epsilon offset problem.

Two performance notes follow from the derivation rather than from measurement, so they hold generally. First, the test is bound by load latency on the vertex data rather than by the arithmetic: three vertices at 12 bytes each is one and a bit cache lines if the layout is tight, and the kernel’s throughput on modern cores is limited by how many independent triangle tests can be in flight. Second, back-face culling is free in this formulation — det > 0 culls one orientation, det < 0 the other, and dropping the fabsf costs nothing — but culling changes which triangles a ray can hit, which matters for open surfaces such as single-sided cloth and for any mesh with inconsistent winding. A renderer that culls must know its meshes are consistently wound, and the watertight form’s conservative behaviour is slightly different depending on whether culling is enabled, because a degenerate triangle with zero projected area is accepted by both signs.

The remaining limitation is geometric rather than numerical: a triangle has no thickness, so a ray that starts inside a closed mesh and continues outward will hit the far side, and a point exactly on a surface needs an offset before spawning secondary rays. Both of those concerns belong to the interval arithmetic of the next layer up, not to this test.