Ray intersections: spheres, triangles, boxes

Every ray tracer, path tracer and hardware RT core spends most of its time answering one question: where does this ray first meet the scene? Primary rays find what a pixel sees, shadow rays ask whether a light is blocked, bounce rays gather indirect light, and NeRF-style volume renderers march along rays through a bounding box. This lesson builds the four tests everything else is made of — plane, sphere, triangle, box — with the numerical care that separates a toy from a renderer.

1. A ray is a half-line with an interval

A ray is a point and a direction, plus the range of the parameter we care about:

p(t) = o + t·d        for t in (tmin, tmax)

• d need not be unit length. t is then measured in units of d. Camera rays are usually normalised; rays transformed into an instance's object space (a scaled mesh) are not, and renormalising them would change every t. Write every test so it works for any d ≠ 0.
• tmax shrinks as hits are found ("only accept something nearer than what I already have") and is the light's distance for a shadow ray.
• tmin is slightly above 0 for rays that start on a surface (section 7).

#include <limits>
#include <print>

struct vec3 { float x = 0, y = 0, z = 0; };
constexpr vec3 operator+(vec3 a, vec3 b) { return {a.x + b.x, a.y + b.y, a.z + b.z}; }
constexpr vec3 operator*(vec3 a, float s) { return {a.x * s, a.y * s, a.z * s}; }

struct Ray {
    vec3 o, d;
    float tmin = 0.0f;
    float tmax = std::numeric_limits<float>::infinity();
    constexpr vec3 at(float t) const { return o + d * t; }
};

int main() {
    Ray r{{0, 1, 0}, {0, 0, 2}};          // d has length 2
    vec3 p = r.at(1.5f);                   // 1.5 "d-units" = 3 units of distance
    std::println("p = ({}, {}, {}), tmax = {}", p.x, p.y, p.z, r.tmax);
}

The recipe for every primitive: write the surface as an equation F(p) = 0, substitute p = o + t·d, solve for t, keep the solutions inside the interval.

2. Ray–plane

A plane is n·p = k (normal n, offset k). Substituting:

n·(o + t·d) = k   ⇒   t = (k − n·o) / (n·d)

n·d is the cosine-like term: zero when the ray runs parallel to the plane. Then the ray is either inside the plane (every t works) or never touches it; either way, report no hit. The denominator's sign tells you which side the ray sees: n·d < 0 means the ray comes from the side the normal points to.

#include <cmath>
#include <optional>
#include <print>

struct vec3 { float x = 0, y = 0, z = 0; };
constexpr float dot(vec3 a, vec3 b) { return a.x * b.x + a.y * b.y + a.z * b.z; }

std::optional<float> intersect_plane(vec3 o, vec3 d, vec3 n, float k, float tmin, float tmax) {
    float denom = dot(n, d);
    if (std::fabs(denom) < 1e-8f) return std::nullopt;    // parallel
    float t = (k - dot(n, o)) / denom;
    if (t <= tmin || t >= tmax) return std::nullopt;
    return t;
}

int main() {
    vec3 up{0, 1, 0};                                       // the floor y = 0
    auto a = intersect_plane({0, 2, 0}, {0, -1, 0}, up, 0, 0, 100);
    auto b = intersect_plane({0, 2, 0}, {1, 0, 0}, up, 0, 0, 100);
    std::println("down: {}  sideways: {}", a ? *a : -1.0f, b ? *b : -1.0f);
}

3. Ray–sphere: a quadratic, solved carefully

A sphere of centre C and radius r is |p − C|² = r². With oc = o − C:

|oc + t·d|² = r²   ⇒   (d·d)·t² + 2(oc·d)·t + (oc·oc − r²) = 0
                          a            b = 2h          c

Writing b = 2h (the half-b form) cancels the 2s and 4s: disc = h² − a·c, roots t = (−h ± √disc) / a. Negative disc: the ray misses.

Which root? Sort the roots, then return the smaller one if it lies in (tmin, tmax), else the larger one if it does. The second case is not exotic: it is every ray whose origin is inside the sphere — refraction rays in a glass ball, a camera inside a sky dome — and every ray leaving the surface inward.

Catastrophic cancellation

For a ray heading towards the sphere, h < 0, so −h > 0 and the near root (−h − √disc)/a is a difference of two positive numbers. When c is small compared with h² — the origin is close to the surface, as for every secondary ray off a big sphere — √disc ≈ |h| and the subtraction cancels the leading digits. The fix (Numerical Recipes §5.6; Haines et al., Ray Tracing Gems ch. 7) never subtracts:

q  = −(h + sign(h)·√disc)        two numbers of the same sign: no cancellation
t0 = q / a,   t1 = c / q         Vieta: t0·t1 = c/a

#include <cmath>
#include <print>

int main() {
    // Origin 0.5 outside a sphere of radius 1024 centred at 0, ray d = (0.6, 0, 0.8).
    float oz = -1024.5f, dx = 0.6f, dz = 0.8f, r = 1024.0f;
    float a = dx * dx + dz * dz;
    float h = oz * dz;                                // oc·d
    float c = oz * oz - r * r;                        // exact here: 1024.25
    float disc = h * h - a * c;

    float naive = (-h - std::sqrt(disc)) / a;         // cancels
    float q = -(h + std::copysign(std::sqrt(disc), h));
    float stable = c / q;                             // the small root, no cancellation

    double A = double(dx) * dx + double(dz) * dz, H = double(oz) * dz, Cc = double(oz) * oz - double(r) * r;
    double truth = Cc / -(H + std::copysign(std::sqrt(H * H - A * Cc), H));
    std::println("naive  {:.9f}  rel err {:.1e}", naive, std::fabs(naive - truth) / truth);
    std::println("stable {:.9f}  rel err {:.1e}", stable, std::fabs(stable - truth) / truth);
}

The naive root is off by 4·10⁻⁵ relative — 25 µm on a 1 m object, enough to put a reflection ray on the wrong side of a thin layer. Ray Tracing Gems ch. 7 also fixes the opposite case, a small sphere far away, where h² − ac itself cancels: compute disc = a·(r² − |oc − (h/a)·d|²) instead. Production code (PBRT-v4) carries error bounds on t as well.

In Python (NumPy). The same test vectorises over a million rays with no loop: h = np.einsum('ij,ij->i', oc, d), disc = h*h - a*c, q = -(h + np.copysign(np.sqrt(np.maximum(disc, 0)), h)), then np.where picks the root. That's how NeRF codebases compute ray–sphere and ray–box "near/far" bounds for whole batches.

4. Ray–triangle: Möller–Trumbore

A point of the triangle (a, b, c) is a + u·e1 + v·e2 with e1 = b − a, e2 = c − a, u, v ≥ 0, u + v ≤ 1 — barycentric coordinates: the hit point is (1−u−v)·a + u·b + v·c, and the same weights interpolate normals, UVs and colours. Setting it equal to the ray gives a 3×3 linear system:

o + t·d = a + u·e1 + v·e2   ⇔   [ −d  e1  e2 ] · (t, u, v)ᵀ = o − a = s

Möller & Trumbore (1997) solve it with Cramer's rule, writing each determinant as a scalar triple product and reusing the cross products:

p = d × e2        det = e1 · p
u = (s · p) / det
q = s × e1        v = (d · q) / det        t = (e2 · q) / det

#include <cmath>
#include <optional>
#include <print>

struct vec3 { float x = 0, y = 0, z = 0; };
constexpr vec3 operator-(vec3 a, vec3 b) { return {a.x - b.x, a.y - b.y, a.z - b.z}; }
constexpr float dot(vec3 a, vec3 b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
constexpr vec3 cross(vec3 a, vec3 b) { return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x}; }

struct TriHit { float t, u, v; };

std::optional<TriHit> moller_trumbore(vec3 o, vec3 d, vec3 a, vec3 b, vec3 c) {
    vec3 e1 = b - a, e2 = c - a;
    vec3 p = cross(d, e2);
    float det = dot(e1, p);
    if (std::fabs(det) < 1e-8f) return std::nullopt;   // parallel to the plane
    float inv = 1.0f / det;
    vec3 s = o - a;
    float u = dot(s, p) * inv;
    if (u < 0 || u > 1) return std::nullopt;
    vec3 q = cross(s, e1);
    float v = dot(d, q) * inv;
    if (v < 0 || u + v > 1) return std::nullopt;
    float t = dot(e2, q) * inv;
    if (t <= 0) return std::nullopt;
    return TriHit{t, u, v};
}

int main() {
    vec3 a{0, 0, 0}, b{1, 0, 0}, c{0, 1, 0};
    if (auto h = moller_trumbore({0.25f, 0.5f, 2}, {0, 0, -1}, a, b, c))
        std::println("t = {}, u = {}, v = {}, w = {}", h->t, h->u, h->v, 1 - h->u - h->v);
    std::println("outside: {}", moller_trumbore({0.8f, 0.8f, 2}, {0, 0, -1}, a, b, c).has_value());
}

Three details that matter in practice:

• Parallel rays. det = 0 when d lies in the triangle's plane. Then 1/det = inf, and with the origin in the plane s·p = 0, so u = 0·inf = NaN. Every comparison with NaN is false, so NaN passes u < 0 || u > 1 and the function "hits" at t = NaN. The |det| < ε check is not an optimisation. (det scales with |d|·|e1|·|e2|: an absolute ε suits unit-scale scenes; scale-aware code compares det² with ε²·|e1|²·|p|².)
• Back faces. det > 0 exactly when the ray sees the counter-clockwise side ((b−a) × (c−a) points back at the origin). Culling back faces is if (det < ε) return miss; — free, and right for closed opaque meshes; wrong for glass, foliage and anything two-sided.
• Watertightness. Two triangles sharing an edge each evaluate "is the hit inside my edge?" with their own rounding, so a ray through the shared edge can pass both tests as outside and slip through: a speckle of background pixels along mesh edges. Woop, Benthin & Wald (2013), "Watertight Ray/Triangle Intersection", shear the triangle into a ray-aligned space where the edge functions are computed identically for both neighbours (with a double-precision fallback for exact zeros); Embree, PBRT-v4 and RT hardware do this.

5. Ray–box: the slab test

An axis-aligned box [lo, hi] is the intersection of three slabs, the regions between two parallel planes. Along the ray, the slab lo.x ≤ x ≤ hi.x is occupied for t between (lo.x − o.x)/d.x and (hi.x − o.x)/d.x. The box is the intersection of three such intervals and [tmin, tmax]: the latest entry and the earliest exit. It is a hit when entry ≤ exit.

for (int k = 0; k < 3; ++k) {
    float inv = 1.0f / r.d[k];             // precomputed once per ray in a BVH
    float tn = (box.lo[k] - r.o[k]) * inv;
    float tf = (box.hi[k] - r.o[k]) * inv;
    if (tn > tf) std::swap(tn, tf);        // d[k] < 0: the ray enters at hi
    t0 = tn > t0 ? tn : t0;                // written so a NaN tn is ignored
    t1 = tf < t1 ? tf : t1;
    if (t0 > t1) return std::nullopt;
}
return Interval{t0, t1};

Negative directions. With d.x < 0 the ray reaches hi.x first, so (lo − o)/d is the exit. Forget the swap and every ray travelling in −x, −y or −z misses.

Zero components are fine — almost. 1.0f / 0.0f is +inf in IEEE 754. A ray parallel to the x slab inside it gets (lo − o)·inf = −inf and (hi − o)·inf = +inf: "always inside", exactly right. Outside it, both are +inf (or both −inf): "never", also right. No special case needed. The one trap is an origin lying exactly on a face plane: (lo − o) = 0 and 0·inf = NaN. What happens next depends on operand order:

#include <algorithm>
#include <cmath>
#include <limits>
#include <print>

int main() {
    const float inf = std::numeric_limits<float>::infinity();
    float tn = 0.0f * inf;                     // origin on the face, d[k] == 0
    float t0 = 1.5f;
    std::println("tn                  = {}", tn);
    std::println("std::max(tn, t0)    = {}", std::max(tn, t0));   // (tn < t0) ? t0 : tn -> NaN
    std::println("std::max(t0, tn)    = {}", std::max(t0, tn));   // (t0 < tn) ? tn : t0 -> 1.5
    std::println("tn > t0 ? tn : t0   = {}", tn > t0 ? tn : t0);
    std::println("std::fmax(tn, t0)   = {}", std::fmax(tn, t0));  // fmax drops NaN
}

std::max(a, b) is (a < b) ? b : a: a NaN in a comes back out. Let it into t0 and t0 > t1 is false forever after — the box test returns an interval containing NaN, and traversal goes wherever. Write the update as tn > t0 ? tn : t0 (false for NaN, so the old bound stays) — PBRT does — or handle d[k] == 0 explicitly. x86's minss/maxss instructions have the same "return the second operand" rule, which branchless SIMD slab tests exploit (Tavian Barnes, "Fast, Branchless Ray/Bounding Box Intersections, Part 2").

6. Closest hit and any hit

Over a list of primitives, keep the best hit and shrink tmax to it, so each later test only accepts something nearer. With mixed primitive types a std::variant and std::get_if (or std::visit) dispatch without virtual calls:

std::optional<Hit> best;
for (std::size_t i = 0; i < prims.size(); ++i) {
    std::optional<float> t;
    if (auto* s = std::get_if<Sphere>(&prims[i])) t = intersect_sphere(r, *s);
    else if (auto h = intersect_triangle(r, std::get<Triangle>(prims[i]))) t = h->t;
    if (t) { best = Hit{*t, i}; r.tmax = *t; }
}

A shadow ray needs only any hit in (tmin, tmax) and can return at the first one; that's why ray-tracing APIs have separate closest-hit and any-hit (or "terminate on first hit") paths. Linear loops are O(n) per ray; the next skill, gfx-bvh, makes it O(log n) — and the shrinking tmax becomes its pruning rule.

7. Shadow acne: rays that start on a surface

A hit point o + t·d computed in float is off the true surface by a rounding error — above or below it. A shadow ray from a point just below the surface hits the surface itself at t ≈ 10⁻⁶, and the point is wrongly in shadow. The result is the speckled "shadow acne" every first ray tracer produces.

Two standard fixes:

• tmin = ε: ignore hits closer than ε. Simple, and what the build task asks for. But ε is scale-dependent: 10⁻⁴ is huge in a scene measured in micrometres and invisible in one measured in kilometres, and grazing rays can still re-hit.
• Offset the origin along the geometric normal, to the side the new ray travels into (+n for reflection, −n for refraction). PBRT-v4 computes a conservative bound on the hit point's rounding error and offsets by exactly that much; Wächter & Binder, Ray Tracing Gems ch. 6, give a cheap integer-ULP version ("A Fast and Robust Method for Avoiding Self-Intersection"). Both work at any scale.

8. In an interview

• "Intersect a ray with a sphere." Derive the quadratic, say which root you return when the origin is inside, and mention the stable form unprompted.
• "How do you test a ray against a triangle?" Möller–Trumbore, barycentrics, the determinant's sign for culling, the parallel case — and that it isn't watertight, with Woop et al. as the fix.
• "Why do my shadows have speckles?" Self-intersection; tmin or normal-offset origins.
• "Your BVH occasionally misses axis-aligned rays." NaN from 0·inf in the slab test and std::min/std::max operand order.