3. Bounding Volume Hierarchies

Ray Tracing: The Next Week (v3.2.3): 3 Bounding Volume Hierarchies / 2.3 Bounding Volume Hierarchies

Ray-object hit testing is the largest bottleneck in a ray tracer. A naive implementation takes $O (N)$ time for $N$ objects. Because millions of rays are repeatedly fired into the same scene, an idea similar to binary search can reduce this cost to $O (\log N)$ . A bounding volume hierarchy (BVH) is the data structure that makes this possible, and it is also the most intricate component in this ray tracer.

Four new elements are required:

Aabb, which represents an axis-aligned bounding box (AABB) and tests it for intersection with a ray using Andrew Kensler's optimized method
surrounding_box, which returns the smallest AABB containing two other AABBs
An additive bounding_box method on the Hittable trait that does not break the existing r1XX crates
BvhNode, a BVH node that is itself a Hittable

Basic Idea

A bounding volume for a set of objects is a solid that completely encloses every object in the set. For example, if a bounding sphere encloses ten objects, a ray that misses it cannot hit any of those ten objects. If the ray intersects the bounding sphere, it may intersect one or more of them. The code takes the following form:

text

if the ray intersects the bounding volume
    test each object inside for intersection
otherwise
    return no hit

The important point is that the implementation partitions the set of objects rather than the image or space. Each object belongs to exactly one bounding volume, but the bounding volumes may overlap.

A Hierarchy of Bounding Volumes

To achieve sublinear^[sublinear]A complexity smaller than linear $O (N)$ , such as $O (\log N)$ . The opposite is superlinear complexity, such as $O (N^{2})$ . performance, the bounding volumes are arranged into a hierarchy. For example, divide a set of objects into red and blue groups and enclose each group in a rectangular bounding volume. These two volumes then determine a parent bounding volume that contains both.

text

if the ray intersects the parent
    hit_blue = test the objects in the blue group
    hit_red  = test the objects in the red group
    if either side is hit, return the nearer hit record
otherwise
    return no hit

The red and blue volumes may overlap inside their parent, and neither has an ordering relationship with the other. Viewed as a tree, the left and right children have no semantic distinction; both are simply contained by their parent. Each ray traverses the tree from its root and prunes an entire subtree whenever it misses the parent's AABB. At a leaf, the renderer finally performs the ordinary hit test.

Axis-Aligned Bounding Boxes

Several shapes can serve as bounding volumes, but an axis-aligned bounding box (AABB) works particularly well for many models. The standard interpretation of an $n$ -dimensional AABB is the intersection of $n$ axis-aligned intervals called slabs. In two dimensions, the intersection of an interval on $x$ and an interval on $y$ forms a rectangle.

For the ray $P (t) = A + t b$ , the parameter values $t$ at which it intersects the slab $x \in [x_{0}, x_{1}]$ are

t_{x 0} = \frac{x_{0} - A_{x}}{b_{x}}, t_{x 1} = \frac{x_{1} - A_{x}}{b_{x}}

A ray intersects a three-dimensional AABB when the intervals of $t$ obtained for $x$ , $y$ , and $z$ have a nonempty intersection. In other words,

max (t_{x 0}, t_{y 0}, t_{z 0}, t_{min}) < min (t_{x 1}, t_{y 1}, t_{z 1}, t_{max})

The test only needs to determine whether these intervals overlap. No intersection point or surface normal is required, so it can be implemented very efficiently.

There are two details to consider:

If the ray travels in the negative direction along an axis, such as when $b_{x} < 0$ , the interval is reversed and $t_{x 0} > t_{x 1}$ .
If $b_{x} = 0$ , division by zero occurs. Under IEEE 754, the result is $\pm \infty$ , and a ray outside the interval produces infinities of the same sign at both ends. A NaN may also occur.

The standard solution is to absorb these cases with fmin and fmax, together with a swap when the interval is reversed.

AABB Hit Testing

I use the optimized method proposed by Pixar's Andrew Kensler. It performs one division per axis and branches only to swap the endpoints when their order is reversed.

common/src/aabb.rs

rust

pub struct Aabb {
    pub minimum: Point3,
    pub maximum: Point3,
}

impl Aabb {
    pub fn new(a: Point3, b: Point3) -> Self {
        Aabb { minimum: a, maximum: b }
    }
    pub fn min(&self) -> Point3 { self.minimum }
    pub fn max(&self) -> Point3 { self.maximum }

    pub fn hit(&self, r: &Ray, mut t_min: f64, mut t_max: f64) -> bool {
        for a in 0..3 {
            let inv_d = 1.0 / r.direction()[a];
            let mut t0 = (self.minimum[a] - r.origin()[a]) * inv_d;
            let mut t1 = (self.maximum[a] - r.origin()[a]) * inv_d;
            if inv_d < 0.0 {
                std::mem::swap(&mut t0, &mut t1);
            }
            if t0 > t_min { t_min = t0; }
            if t1 < t_max { t_max = t1; }
            if t_max <= t_min {
                return false;
            }
        }
        true
    }
}

The loop over the three axes progressively updates the interval $[t_{min}, t_{max}]$ and returns false as soon as that interval becomes empty.

surrounding_box is a free function that returns the smallest AABB containing two other AABBs. It is used to calculate the AABB of a parent node during BVH construction.

common/src/aabb.rs

rust

pub fn surrounding_box(box0: Aabb, box1: Aabb) -> Aabb {
    let small = Point3::new(
        box0.min().x().min(box1.min().x()),
        box0.min().y().min(box1.min().y()),
        box0.min().z().min(box1.min().z()),
    );
    let big = Point3::new(
        box0.max().x().max(box1.max().x()),
        box0.max().y().max(box1.max().y()),
        box0.max().z().max(box1.max().z()),
    );
    Aabb::new(small, big)
}

Extending the `Hittable` Trait

Add a bounding_box method to the trait so that an AABB can be obtained from any Hittable. A moving object must return an AABB that contains the object at every time in $[t_{0}, t_{1}]$ , so the method takes time0 and time1 as arguments.

common/src/hittable.rs

rust

pub trait Hittable {
    fn hit(&self, r: &Ray, t_min: f64, t_max: f64) -> Option<HitRecord>;

    /// Returns `None` for objects without a defined AABB, such as infinite planes.
    /// The default implementation returns `None`, so objects that are not added
    /// to a BVH, including spheres created with `Sphere::new` in Book 1, require
    /// no changes.
    fn bounding_box(&self, _time0: f64, _time1: f64) -> Option<Aabb> {
        None
    }
}

Differences Between C++ and Rust

The original C++ declares bounding_box as a pure virtual function, requiring every derived class to implement it. This implementation instead defines a method with a default implementation that returns Option<Aabb>. There are two reasons:

I want this to remain an additive extension so that the Hittable implementations defined by the r1XX crates in Book 1, including the basic spheres created through Sphere::new, continue to build unchanged.
Representing the absence of an AABB with Option is more idiomatic in Rust than using a Boolean result and an output parameter.

Some objects, such as infinite planes, genuinely have no AABB. Returning Option<Aabb> is therefore a natural type for the trait as a whole.

Sphere, MovingSphere, and HittableList override the default implementation. The AABB for a Sphere extends from its center by $\pm (r, r, r)$ .

common/src/sphere.rs

rust

impl Hittable for Sphere {
    fn hit(...) -> Option<HitRecord> { /* Existing implementation. */ }

    fn bounding_box(&self, _time0: f64, _time1: f64) -> Option<Aabb> {
        let r = self.radius.abs();
        let radius_vec = Vec3::new(r, r, r);
        Some(Aabb::new(self.center - radius_vec, self.center + radius_vec))
    }
}

For MovingSphere, surrounding_box encloses the sphere's AABB at time0 together with its AABB at time1.

common/src/moving_sphere.rs

rust

fn bounding_box(&self, time0: f64, time1: f64) -> Option<Aabb> {
    let radius_vec = Vec3::new(self.radius, self.radius, self.radius);
    let box0 = Aabb::new(
        self.center(time0) - radius_vec,
        self.center(time0) + radius_vec,
    );
    let box1 = Aabb::new(
        self.center(time1) - radius_vec,
        self.center(time1) + radius_vec,
    );
    Some(surrounding_box(box0, box1))
}

HittableList returns the result of folding its children's AABBs with surrounding_box. An empty list returns None, and None from any child is propagated.

common/src/hittable_list.rs

rust

fn bounding_box(&self, time0: f64, time1: f64) -> Option<Aabb> {
    if self.objects.is_empty() {
        return None;
    }
    let mut output: Option<Aabb> = None;
    for object in &self.objects {
        let temp = object.bounding_box(time0, time1)?;
        output = Some(match output {
            None => temp,
            Some(prev) => surrounding_box(prev, temp),
        });
    }
    output
}

BVH Nodes

A BVH is itself a kind of Hittable. Although it is only a container for other Hittable objects, it can answer whether a given ray hits it. Following the original, I represent the entire tree with a single node class rather than separating trees from nodes. The root is simply another node.

common/src/bvh.rs

rust

pub struct BvhNode {
    left: Arc<dyn Hittable>,
    right: Arc<dyn Hittable>,
    bbox: Aabb,
}

The child pointers use Arc<dyn Hittable> rather than a concrete type. A child may be another BvhNode, a Sphere, or a MovingSphere. I use Arc because BVH construction sometimes places the same leaf object in both child positions for a subtree containing a single element.

The hit implementation is straightforward: if the ray misses the parent AABB, it never descends to the children.

common/src/bvh.rs

rust

impl Hittable for BvhNode {
    fn hit(&self, r: &Ray, t_min: f64, t_max: f64) -> Option<HitRecord> {
        if !self.bbox.hit(r, t_min, t_max) {
            return None;
        }
        let hit_left = self.left.hit(r, t_min, t_max);
        // Use the t value of a left-side hit as the upper bound for the right,
        // further narrowing the right-side search interval.
        let new_t_max = match &hit_left {
            Some(rec) => rec.t,
            None => t_max,
        };
        let hit_right = self.right.hit(r, t_min, new_t_max);
        hit_right.or(hit_left)
    }

    fn bounding_box(&self, _time0: f64, _time1: f64) -> Option<Aabb> {
        Some(self.bbox)
    }
}

hit_right.or(hit_left) is an idiomatic way to return the right-side hit, which is nearer if it exists, and otherwise return the left-side result. Reducing new_t_max also enables the right-side traversal to terminate early.

Partitioning the BVH

Construction is the most complicated part of an acceleration data structure. The BVH is built inside BvhNode::new. Once the set of children has been divided, hit works correctly, so any partition is sufficient for correctness. Smaller child AABBs will naturally be faster than large ones, but that affects performance rather than correctness.

Following the original, I use the following algorithm as a compromise between simplicity and performance:

Randomly select one axis
Sort the objects along that axis
Split the range at its midpoint and recurse on the left and right halves

For one element, place the same object in both child positions and terminate the recursion. Thanks to Arc::clone, this duplicates only the pointer and is nearly free. For two elements, the comparator determines whether to swap them. For three or more elements, use sort_by and split the range.

common/src/bvh.rs

rust

fn build(
    objects: &mut Vec<Arc<dyn Hittable>>,
    start: usize, end: usize,
    time0: f64, time1: f64,
) -> Self {
    let axis = random_int(0, 2) as usize;
    let span = end - start;

    let (left, right): (Arc<dyn Hittable>, Arc<dyn Hittable>) = if span == 1 {
        (objects[start].clone(), objects[start].clone())
    } else if span == 2 {
        if box_compare(objects[start].as_ref(),
                       objects[start + 1].as_ref(), axis) == Ordering::Less {
            (objects[start].clone(), objects[start + 1].clone())
        } else {
            (objects[start + 1].clone(), objects[start].clone())
        }
    } else {
        objects[start..end]
            .sort_by(|a, b| box_compare(a.as_ref(), b.as_ref(), axis));
        let mid = start + span / 2;
        let l = Arc::new(Self::build(objects, start, mid, time0, time1));
        let r = Arc::new(Self::build(objects, mid, end, time0, time1));
        (l as Arc<dyn Hittable>, r as Arc<dyn Hittable>)
    };

    let box_left = left.bounding_box(time0, time1)
        .expect("No bounding box in BvhNode constructor");
    let box_right = right.bounding_box(time0, time1)
        .expect("No bounding box in BvhNode constructor");
    let bbox = surrounding_box(box_left, box_right);

    BvhNode { left, right, bbox }
}

The public API provides two constructors: BvhNode::new and BvhNode::from_list, which directly consumes a HittableList. from_list is syntactic sugar that converts Box<dyn Hittable> values into Arc<dyn Hittable> values using the standard library's impl<T: ?Sized> From<Box<T>> for Arc<T> implementation.

common/src/bvh.rs

rust

pub fn from_list(list: HittableList, time0: f64, time1: f64) -> Self {
    let objects: Vec<Arc<dyn Hittable>> =
        list.objects.into_iter().map(Arc::<dyn Hittable>::from).collect();
    Self::new(objects, time0, time1)
}

The random_int(min, max) helper, added to common::utils, returns a uniformly distributed integer over [min, max]. It is implemented as random_double_range(min as f64, (max + 1) as f64).floor() as i32.

AABB Comparison Function

The comparator used for sorting only compares the minimum corners of two AABBs along the selected axis. As in the original, one general function that takes an axis argument is sufficient. The C++ version separates this into three unary functions because std::sort passes function pointers, but Rust closures make that unnecessary.

common/src/bvh.rs

rust

fn box_compare(a: &dyn Hittable, b: &dyn Hittable, axis: usize) -> Ordering {
    let box_a = a.bounding_box(0.0, 0.0)
        .expect("No bounding box in box_compare");
    let box_b = b.bounding_box(0.0, 0.0)
        .expect("No bounding box in box_compare");
    box_a.min()[axis]
        .partial_cmp(&box_b.min()[axis])
        .unwrap_or(Ordering::Equal)
}

Because f64 does not implement Ord, I use partial_cmp and collapse None caused by NaN to Equal. This case does not occur for a well-formed AABB.

common/src/bvh.rs

rust

use std::cmp::Ordering;
use std::sync::Arc;

use crate::aabb::{Aabb, surrounding_box};
use crate::hittable::{HitRecord, Hittable};
use crate::hittable_list::HittableList;
use crate::ray::Ray;
use crate::utils::random_int;

pub struct BvhNode {
    left: Arc<dyn Hittable>,
    right: Arc<dyn Hittable>,
    bbox: Aabb,
}

impl BvhNode {
    pub fn new(mut objects: Vec<Arc<dyn Hittable>>, time0: f64, time1: f64) -> Self {
        let len = objects.len();
        Self::build(&mut objects, 0, len, time0, time1)
    }

    pub fn from_list(list: HittableList, time0: f64, time1: f64) -> Self {
        let objects: Vec<Arc<dyn Hittable>> =
            list.objects.into_iter().map(Arc::<dyn Hittable>::from).collect();
        Self::new(objects, time0, time1)
    }

    fn build(
        objects: &mut Vec<Arc<dyn Hittable>>,
        start: usize, end: usize,
        time0: f64, time1: f64,
    ) -> Self {
        let axis = random_int(0, 2) as usize;
        let span = end - start;

        let (left, right): (Arc<dyn Hittable>, Arc<dyn Hittable>) = if span == 1 {
            (objects[start].clone(), objects[start].clone())
        } else if span == 2 {
            if box_compare(objects[start].as_ref(),
                           objects[start + 1].as_ref(), axis) == Ordering::Less {
                (objects[start].clone(), objects[start + 1].clone())
            } else {
                (objects[start + 1].clone(), objects[start].clone())
            }
        } else {
            objects[start..end]
                .sort_by(|a, b| box_compare(a.as_ref(), b.as_ref(), axis));
            let mid = start + span / 2;
            let l = Arc::new(Self::build(objects, start, mid, time0, time1))
                as Arc<dyn Hittable>;
            let r = Arc::new(Self::build(objects, mid, end, time0, time1))
                as Arc<dyn Hittable>;
            (l, r)
        };

        let box_left = left.bounding_box(time0, time1)
            .expect("No bounding box in BvhNode constructor");
        let box_right = right.bounding_box(time0, time1)
            .expect("No bounding box in BvhNode constructor");
        let bbox = surrounding_box(box_left, box_right);

        BvhNode { left, right, bbox }
    }
}

impl Hittable for BvhNode {
    fn hit(&self, r: &Ray, t_min: f64, t_max: f64) -> Option<HitRecord> {
        if !self.bbox.hit(r, t_min, t_max) {
            return None;
        }
        let hit_left = self.left.hit(r, t_min, t_max);
        let new_t_max = match &hit_left {
            Some(rec) => rec.t,
            None => t_max,
        };
        let hit_right = self.right.hit(r, t_min, new_t_max);
        hit_right.or(hit_left)
    }

    fn bounding_box(&self, _time0: f64, _time1: f64) -> Option<Aabb> {
        Some(self.bbox)
    }
}

fn box_compare(a: &dyn Hittable, b: &dyn Hittable, axis: usize) -> Ordering {
    let box_a = a.bounding_box(0.0, 0.0)
        .expect("No bounding box in box_compare");
    let box_b = b.bounding_box(0.0, 0.0)
        .expect("No bounding box in box_compare");
    box_a.min()[axis]
        .partial_cmp(&box_b.min()[axis])
        .unwrap_or(Ordering::Equal)
}

Demo: `r203-bvh`

The r203-bvh crate wraps the same bouncing-sphere scene from Chapter 2 in a BVH by passing it through BvhNode::from_list.

r203-bvh/src/lib.rs

rust

let world = BvhNode::from_list(random_scene(), 0.0, 1.0);

The scene contents and camera parameters are the same as in Chapter 2. Only the complexity of hit testing differs. I keep the image dimensions at 300×169 while increasing the sample count from 30 to 100. A naive linear scan would take an impractical amount of time at 100 samples per pixel, but the BVH provides an approximate $\log_{2} (N)$ speedup and completes in an acceptable time.

The WASM demo at the top of the page shows the result. The output has much the same overall appearance as in Chapter 2, but the higher sample count visibly reduces noise.

Summary

Added Aabb, Andrew Kensler's optimized hit method, and surrounding_box to common/.
Added bounding_box to the Hittable trait with a default implementation, allowing the r1XX crates from Book 1 to continue building unchanged.
Implemented bounding_box for Sphere, MovingSphere, and HittableList.
Implemented BvhNode, constructed recursively by selecting a random axis, sorting, and splitting at the midpoint.
Demonstrated the practical effect in r203-bvh by rendering the same scene with a higher sample count.

The most difficult mechanism in Book 2 is now complete. The next chapter introduces textures, adding greater artistic control on top of the BVH.

A complexity smaller than linear $O (N)$ , such as $O (\log N)$ . The opposite is superlinear complexity, such as $O (N^{2})$ . ↩︎

3. Bounding Volume Hierarchies ​

Basic Idea ​

A Hierarchy of Bounding Volumes ​

Axis-Aligned Bounding Boxes ​

AABB Hit Testing ​

Extending the Hittable Trait ​

Differences Between C++ and Rust ​

BVH Nodes ​

Partitioning the BVH ​

AABB Comparison Function ​

Demo: r203-bvh ​

Summary ​

3. Bounding Volume Hierarchies

Basic Idea

A Hierarchy of Bounding Volumes

Axis-Aligned Bounding Boxes

AABB Hit Testing

Extending the `Hittable` Trait

Differences Between C++ and Rust

BVH Nodes

Partitioning the BVH

AABB Comparison Function

Demo: `r203-bvh`

Summary