Bounding volume hierarchies

A 1080p frame at one sample per pixel with a shadow ray and a bounce is about 6 million rays. Against a million-triangle scene, testing every triangle is 6·10¹² ray–triangle tests per frame — hopeless. A bounding volume hierarchy (BVH) brings it down to a few dozen box tests and a handful of triangle tests per ray: if a ray misses a box, it misses everything inside. Every production ray tracer (Embree, OptiX, PBRT, Cycles) and every RT core in a GPU traverses a BVH; collision detection and spatial queries in physics engines use the same structure.

1. Objects in boxes, boxes in boxes

A BVH partitions objects, not space: each triangle is in exactly one leaf; each node stores the axis-aligned box (AABB) enclosing everything below it. Children's boxes may overlap — unlike a k-d tree or an octree, which split space and may have to reference a triangle from several cells. The pay-off: with leaves of at least one triangle there are at most n leaves and 2n − 1 nodes, so memory is known before building, and a node's box can be updated when geometry moves without restructuring anything.

The cost per ray grows roughly with log n (here a bit faster, because the scene grows and rays travel farther through it). That is the whole reason the structure exists.

2. Memory layout: a flat array of 32-byte nodes

Pointer-based trees scatter nodes over the heap. Renderers store the tree in one array, children allocated as a pair so one index finds both, and the node packed into 32 bytes so two fit in a cache line (Jacco Bikker, "How to build a BVH"; PBRT's LinearBVHNode is the same size):

#include <cstddef>
#include <cstdint>
#include <print>

struct vec3 { float x, y, z; };
struct AABB { vec3 lo, hi; };

struct BVHNode {
    AABB bounds;               // vec3 lo, hi: 24 bytes
    std::uint32_t left_first;  // interior: index of the left child (right = +1)
                               // leaf:     first entry in the triangle-index array
    std::uint32_t count;       // 0 = interior node, > 0 = leaf with this many triangles
};
static_assert(sizeof(BVHNode) == 32);

struct FatNode {               // the "obvious" pointer version
    AABB bounds;
    FatNode* left;
    FatNode* right;
    std::size_t first, count;
};

int main() {
    std::println("BVHNode {} bytes, {} per 64-byte line", sizeof(BVHNode), 64 / sizeof(BVHNode));
    std::println("FatNode {} bytes (offset of left: {})", sizeof(FatNode), offsetof(FatNode, left));
}

Triangles never move. The builder sorts an index array instead, so each node owns a contiguous range of it, and a leaf is just (first, count):

3. Building top-down

Start with one root owning every triangle. While a node has more than max_leaf triangles: choose a split, partition its index range in place (like quicksort's partition), append two children for the halves, recurse.

The simplest split: the longest axis of the node's centroid bounds, at its midpoint. Centroid bounds, not node bounds — a big triangle can stretch the node box far beyond where the triangles' centres are, and splitting there can leave a side empty. Two degenerate inputs break a naive builder:

• a split with an empty side makes a child identical to its parent, and the recursion never ends (stack overflow). Never accept a split that leaves a side empty.
• identical centroids (duplicated geometry, instanced meshes flattened into one) have zero extent on every axis: there is no spatial split at all. Split the range in half by index, or make an oversized leaf.

#include <algorithm>
#include <cstdint>
#include <print>
#include <vector>

// A midpoint-split builder over 1-D "triangles" (intervals) to show the mechanics:
// index array, in-place partition, children appended as a pair.
struct Node { float lo, hi; std::uint32_t left_first, count; };
struct Seg { float a, b; float centroid() const { return (a + b) / 2; } };

struct Builder {
    std::vector<Seg> segs;
    std::vector<std::uint32_t> idx;
    std::vector<Node> nodes;

    void bounds(Node& n) {
        n.lo = 1e30f, n.hi = -1e30f;
        for (auto i = n.left_first; i < n.left_first + n.count; ++i)
            n.lo = std::min(n.lo, segs[idx[i]].a), n.hi = std::max(n.hi, segs[idx[i]].b);
    }
    void subdivide(std::uint32_t ni) {
        std::uint32_t first = nodes[ni].left_first, count = nodes[ni].count;
        if (count <= 2) return;
        float cmin = 1e30f, cmax = -1e30f;
        for (auto i = first; i < first + count; ++i)
            cmin = std::min(cmin, segs[idx[i]].centroid()), cmax = std::max(cmax, segs[idx[i]].centroid());
        float split = (cmin + cmax) / 2;
        std::uint32_t i = first, j = first + count;            // half-open: [i, j) unclassified
        while (i < j) {
            if (segs[idx[i]].centroid() < split) ++i;
            else std::swap(idx[i], idx[--j]);
        }
        if (i == first || i == first + count) i = first + count / 2;   // never an empty side
        auto left = static_cast<std::uint32_t>(nodes.size());
        nodes.push_back({0, 0, first, i - first});
        nodes.push_back({0, 0, i, first + count - i});
        nodes[ni].left_first = left;
        nodes[ni].count = 0;
        bounds(nodes[left]), bounds(nodes[left + 1]);
        subdivide(left), subdivide(left + 1);
    }
};

int main() {
    Builder b;
    b.segs = {{0, 1}, {8, 9}, {1, 2}, {9, 10}, {4, 5}, {2, 3}, {5, 5.5f}, {7, 8}};
    for (std::uint32_t i = 0; i < b.segs.size(); ++i) b.idx.push_back(i);
    b.nodes.push_back({0, 0, 0, static_cast<std::uint32_t>(b.segs.size())});
    b.bounds(b.nodes[0]);
    b.subdivide(0);
    for (std::size_t n = 0; n < b.nodes.size(); ++n) {
        const Node& x = b.nodes[n];
        std::println("node {}: [{}, {}] {} {}", n, x.lo, x.hi, x.count ? "leaf, first" : "children at", x.left_first);
    }
}

4. The surface area heuristic

Midpoint splits ignore where the rays go. The surface area heuristic (SAH; Goldsmith & Salmon 1987, MacDonald & Booth 1990) estimates the cost of a split. For uniformly distributed rays, the probability that a ray hitting a convex box P also hits a convex box C inside it is the ratio of their surface areas A(C)/A(P) (a classic result of integral geometry). So the expected cost of tracing a ray through a node split into L and R is

cost = C_trav + A(L)/A(P) · N(L) · C_isect + A(R)/A(P) · N(R) · C_isect

with N the number of triangles. To compare splits of the same node, constants and A(P) don't matter: minimise A(L)·N(L) + A(R)·N(R). Big, crowded children are bad; a tight box around a dense cluster next to a big sparse box is good — exactly the "teapot in a stadium" case where midpoint splits waste work.

Binning (Wald 2007) makes it cheap. Divide the centroid bounds along an axis into ~16 equal bins; one pass over the node's triangles records each bin's count and the bounds of its triangles. A split can only fall on one of the 15 bin boundaries; sweeping from the left and from the right accumulates A and N for both sides of every boundary in O(bins). Do it for all three axes, take the cheapest.

// For one axis: bins[b] = {count, box}. left_*[b] / right_*[b] describe bins [0, b) / [b, B).
AABB acc; std::uint32_t n = 0;
for (int b = 1; b < B; ++b) { acc.grow(box[b - 1]); n += cnt[b - 1]; left_area[b] = acc.area(); left_n[b] = n; }
acc = {}; n = 0;
for (int b = B - 1; b >= 1; --b) { acc.grow(box[b]); n += cnt[b]; right_area[b] = acc.area(); right_n[b] = n; }
for (int b = 1; b < B; ++b) {
    if (left_n[b] == 0 || right_n[b] == 0) continue;           // not a split
    float cost = left_area[b] * left_n[b] + right_area[b] * right_n[b];
    if (cost < best_cost) best_cost = cost, best = {axis, b};
}

A worked comparison — a dense cluster of 100 small triangles and 4 big ground triangles far away, split two ways:

#include <print>

struct Box { float lo[3], hi[3]; };
float area(const Box& b) {
    float e[3];
    for (int k = 0; k < 3; ++k) e[k] = b.hi[k] - b.lo[k];
    return 2 * (e[0] * e[1] + e[1] * e[2] + e[2] * e[0]);
}
float sah(const Box& l, int nl, const Box& r, int nr) { return area(l) * nl + area(r) * nr; }

int main() {
    // Option A (SAH's pick): the cluster alone vs the ground.
    Box cluster{{60, 0, 60}, {62, 2, 62}}, ground{{-100, 0, -100}, {100, 0.1f, 100}};
    // Option B (midpoint of x): each side gets half the ground, one side also the cluster.
    Box left{{-100, 0, -100}, {0, 0.1f, 100}}, right_with_cluster{{0, 0, -100}, {100, 2, 100}};
    std::println("cluster | ground:      {:>9.0f}", sah(cluster, 100, ground, 4));
    std::println("midpoint of x:         {:>9.0f}", sah(left, 2, right_with_cluster, 102));
}

The midpoint split puts the cluster's 100 triangles in a box as big as half the stadium: every ray crossing that half has to descend into it. The SAH split costs about 13 times less.

Then partition with the same bin function used for counting (bin_of(centroid) < best.b), so the partition agrees exactly with the counts the cost was computed from — in particular, both sides are non-empty.

When to stop. This skill's task splits until every leaf has at most max_leaf triangles. Production builders also compare the best split with the cost of not splitting, N·C_isect, and stop early when splitting doesn't pay; they cap the depth too (fixed-size traversal stacks assume it). Leaves of 1–4 triangles suit CPUs; GPU traversal often prefers slightly bigger leaves.

5. Traversal: nearer child first, prune by the closest hit

Recursion is replaced by an explicit stack of (node, entry t) — cheaper on a CPU and the only option in a GPU shader. Two rules make it fast:

1. Prune by the current closest hit. Every hit sets r.tmax = t; box tests use [tmin, tmax], so a box that starts behind the closest hit is a miss, and a stacked node whose entry is now ≥ tmax is dropped when popped.
2. Visit the nearer child first. Push the farther child, then descend into the nearer one. The first hit found is then usually the closest, and rule 1 discards everything behind it. Without the ordering, pruning still works but kicks in late.

while (sp > 0) {
    const Entry e = stack[--sp];
    if (e.t >= r.tmax) continue;                         // a nearer hit was found meanwhile
    const BVHNode& node = nodes_[e.node];
    if (node.is_leaf()) {
        for (auto i = node.left_first; i < node.left_first + node.count; ++i)
            if (auto h = intersect_triangle(r, tris_[idx_[i]])) {
                best = Hit{h->t, h->u, h->v, idx_[i]};
                r.tmax = h->t;                           // rule 1
            }
        continue;
    }
    std::uint32_t near = node.left_first, far = near + 1;
    float tn = box_entry(nodes_[near].bounds, r.o, inv_d, r.tmin, r.tmax);
    float tf = box_entry(nodes_[far].bounds, r.o, inv_d, r.tmin, r.tmax);
    if (tf < tn) { std::swap(near, far); std::swap(tn, tf); }   // rule 2
    if (tf != kInf) stack[sp++] = {far, tf};
    if (tn != kInf) stack[sp++] = {near, tn};
}

inv_d = 1/d is computed once per ray (section 5 of gfx-ray-intersections: infinities for zero components are fine, NaN needs care). A shadow ray uses any-hit traversal: no ordering, return at the first hit — often 2–3× cheaper.

The stack holds at most one entry per level, so 32–64 entries suffice for any sane tree; GPU traversal keeps it in registers or shared memory, which is one more reason builders cap the depth.

In Python. Production NeRF and Gaussian-splatting pipelines rarely write a BVH in Python; they call one. trimesh's ray queries use Embree, and PyTorch code calls custom CUDA/OptiX kernels — but the build/traverse logic is exactly this, and a NumPy prototype of the SAH sweep is a good interview warm-up.

6. Animation: refit or rebuild

When vertices move but the mesh keeps its topology (a skinned character, cloth), there are two options:

• Refit: keep the tree, recompute every box bottom-up — leaves from their triangles, interior nodes from their two children. Because children always come after their parent in the array, one backwards sweep is a valid order: O(n), no allocation, trivially parallel per level.
• Rebuild: build from scratch. O(n log n) with binned SAH, but the tree fits the new geometry.

Refitting keeps the old grouping: triangles that were neighbours stay in the same leaf even after they move apart, so boxes grow and overlap and traversal slows. Engines refit every frame and rebuild periodically, or when the root's SAH cost has grown by some factor. The ray-tracing APIs expose this choice directly: DXR/Vulkan acceleration structures built with ALLOW_UPDATE can be updated (refit) in place.

7. Parallel builds: LBVH

Top-down SAH building is sequential at the top. Karras (2012) builds a BVH on the GPU with every internal node in parallel: compute a 30-bit Morton code for each centroid (interleave the bits of its quantised x, y, z — so sorting by code orders triangles along a Z-order curve), radix-sort the codes, and then each internal node of a binary radix tree over the sorted keys can find its range and split point independently. Boxes are filled in by a bottom-up pass. The trees are worse than SAH (splits fall on Morton-curve boundaries, not cost minima), so renderers refine them (treelet restructuring, Karras & Aila 2013) or use LBVH for dynamic geometry only. Morton codes and the octrees built from them are the subject of gfx-morton-octree.

8. BVHs in RT hardware: BLAS and TLAS

GPU ray tracing (DXR, Vulkan RT, OptiX, Metal) uses a two-level hierarchy:

• a bottom-level acceleration structure (BLAS) per mesh, in object space, built once (or refitted for deforming meshes);
• a top-level acceleration structure (TLAS) over instances: each a transform plus a reference to a BLAS. Moving an object, or drawing 10,000 copies of a tree, touches only the TLAS, which is small and rebuilt every frame.

Traversal transforms the ray into each instance's object space when it moves from TLAS to BLAS — one reason ray directions are not assumed to be unit length. The driver chooses the actual node format: typically wide BVHs (4–8 children per node, boxes quantised to 8 bits relative to the parent) that the RT cores test in parallel. The builds trade quality for speed with flags such as PREFER_FAST_TRACE and PREFER_FAST_BUILD. You choose the hierarchy (what is a BLAS, what is instanced, what is refitted); the hardware does the traversal from section 5.

9. In an interview

• "How would you speed up ray casting against a big mesh?" A BVH: say why objects-in-boxes (bounded memory, refit-able), the SAH formula and its meaning, binning, and closest-first traversal with pruning.
• "Your BVH build hangs on some meshes." Empty-side splits or identical centroids: guard both.
• "Some rays miss geometry they should hit." A node box that doesn't contain its subtree (centroid bounds stored as node bounds, stale boxes after a refit in the wrong order), or NaN in the slab test.
• "Characters animate every frame." Refit vs rebuild, and BLAS/TLAS.