Bounding boxes and hierarchies

This is an online bonus chapter for The Ray Tracer Challenge, by Jamis Buck. To be successful, you:

If you don't have a copy of the book, grab one today and start writing a 3D renderer of your own!

Your ray tracer probably performs great with the simpler of your test scenes. A plane or two, a handful of spheres, perhaps even a cube or cylinder thrown in for variety—these aren't going to tax your computer terribly. But as soon as your scene grows much bigger than that—and especially if it gets into the hundreds or thousands of elements—you'll find that your renderer begins to bog down, and render times start to drag on into the realm of hours, or even days.

Why? It's because your renderer has to iterate over every object in your scene for each ray it casts, and you'll probably have at least two such rays for each pixel in your image—one for the initial ray, and one for the shadow ray at the point of intersection. Thus, the more objects you have in your scene, the more objects have to be tested for every pixel you render.

Well, then. If you want your renderer to go faster, it becomes a matter of reducing the number of objects that have to be examined at each pixel. Right? The question now is: just how are you supposed to do that?

In this bonus chapter, you'll learn how to optimize large scenes by using bounding boxes and bounding volume hierarchies. By the end, you should see significant improvements when rendering scenes with large numbers of objects, like triangle meshes.

Here we go!

Bounding boxes

Imagine a scene containing a 3D grid of spheres, with 10 spheres on a side.

That's, like, a thousand spheres, right? And the way your ray tracer is currently implemented, every ray you cast is going to have to be tested against every single one of those spheres, which means you can expect your render to take roughly...let's see...multiply...carry the one... Um. Let's just call it an eternity.

BUT.

What if you were to put all of those spheres inside a big, invisible box? And then, what if you reworked your ray tracer so that it only tried to intersect the shapes if a ray happened to intersect the box that contained them?

This is what is called a bounding box. Technically, you could use a sphere instead, or a cylinder, or really any other shape, in which case you would talk about bounding volumes instead of bounding boxes. The goal is to pick a relatively simple shape that forms the tightest possible bound on its contents. Boxes (and especially the axis-aligned bounding boxes from chapter 12) work well as a general solution for this, because they are simple to construct and inexpensive to intersect.

So, that's what you're going to build here. The process will look something like this:

  1. Implement a BoundingBox structure and some associated logic to manipulate it.
  2. Add a function to each primitive shape for querying the shape's bounds (in object space).
  3. Implement a function for recursively finding the bounds of a Group object.
  4. Implement a function for recursively finding the bounds of a CSG object.
  5. Detect if a ray intersects a bounding box.
  6. Add a guard to the local_intersect() function of your aggregate shapes (Group and CSG).

Once you've nailed that down, you'll move on to the second half of this optimization: bounding volume hierarchies, but let's not get ahead of ourselves. First things first: bounding boxes. Onward!

The BoundingBox structure

A BoundingBox is a simple structure describing an axis-aligned bounding box, or AABB. It contains just two points, one identifying the minimum point on the box, and the other identifying the maximum point. The minimum point is where the x, y, and z coordinates are all smallest, and the maximum point is where those coordinates are all largest.

Initially, though, the box will be empty, which means its min and max points are a bit wonky by default. The following test demonstrates this:

bounds.feature
Scenario: Creating an empty bounding box Given box bounding_box(empty) Then box.min = point(infinity, infinity, infinity) And box.max = point(-infinity, -infinity, -infinity)

What? Yeah, like I said, wonky. An empty bounding box has it's smallest point at positive infinity, and it's largest point at negative infinity. Basically, what this is saying is that the empty box is invalid. It contains no space whatsoever.

You may be wondering, "why bother with an empty bounding box at all?" It does seem pretty pointless, but hang in there. When you start dealing with adding points to these boxes, or combining multiple boxes together, having your empty boxes start out "wonky" actually makes things easier.

To describe a bounding box that encapsulates some volume, you instantiate the box with explicit minimum and maximum points, like this:

bounds.feature
Scenario: Creating a bounding box with volume Given box bounding_box(min=point(-1, -2, -3) max=point(3, 2, 1)) Then box.min = point(-1, -2, -3) And box.max = point(3, 2, 1)

Later in this chapter you'll define several operations on these bounding boxes, but for now, you'll add just one: the ability to resize a bounding box so that it includes a given point. It looks like this:

bounds.feature
Scenario: Adding points to an empty bounding box Given box bounding_box(empty) And p1 point(-5, 2, 0) And p2 point(7, 0, -3) When p1 is added to box And p2 is added to box Then box.min = point(-5, 0, -3) And box.max = point(7, 2, 0)

Every time you add a point to a bounding box, the box adjusts its min and max points, resizing itself so that it can contain the given point. It does this by checking the x, y, and z components against the x, y, and z components of the min and max points. If any component is less than the corresponding component of min, that component of min is replaced with the new value. Similarly, if any component is greater than the corresponding component of max, that component gets replaced. Here's some pseudocode:

# adding "point" to "box" box.min.x point.x if point.x < box.min.x box.min.y point.y if point.y < box.min.y box.min.z point.z if point.z < box.min.z box.max.x point.x if point.x > box.max.x box.max.y point.y if point.y > box.max.y box.max.z point.z if point.z > box.max.z

With that implementation in mind, recall how your empty box was defined: min was positive infinity, and max was negative infinity. This means that any point you try to add to that empty box it will be smaller than min, and larger than max, allowing the pseudocode above to automatically snap your box to the given point. See? Wonky is wonderful!

Once you have that much defined, you can begin describing the bounds for each of the primitive shapes in your ray tracer.

Bounds for primitives

By the time you reach chapter 14, "Groups", of The Ray Tracer Challenge, your ray tracer ought to include the following primitives:

In this section, you'll define bounds for each of those primitive types, as well as triangles. (If you haven't yet reached chapter 15, "Triangles", you can skip that part and come back to it when you're ready.)

For each of these primitive shapes, you'll define a function called bounds_of(shape), which returns a bounding box in object space. This means that, for now, you won't need to worry about any transformations on the objects.

Spheres

The spheres in your ray tracer are always at the origin, and have a radius of 1. This means that the corresponding bounding boxes will also be centered at the origin, and extend from -1 to 1 along each axis.

spheres.feature
Scenario: A sphere has a bounding box Given shape sphere() When box bounds_of(shape) Then box.min = point(-1, -1, -1) And box.max = point(1, 1, 1)

Planes

Planes are trickier, because they stretch to infinity in x and z, and have zero thickness in y.

planes.feature
Scenario: A plane has a bounding box Given shape plane() When box bounds_of(shape) Then box.min = point(-infinity, 0, -infinity) And box.max = point(infinity, 0, infinity)

Any bounds you wrap around a plane are effectively useless, because even the tiniest rotation will cause it's bounding box to suddenly stretch to infinity along every axis. Still, having a bounding box defined for it helps avoid special cases, so it's worth doing.

Cubes

Cubes are great, because they can be perfectly bound by their bounding box. There is literally no other shape that bounds a cube more effectively than another cube of the same shape and size.

cubes.feature
Scenario: A cube has a bounding box Given shape cube() When box bounds_of(shape) Then box.min = point(-1, -1, -1) And box.max = point(1, 1, 1)

Neat, huh?

Cylinders

Like planes, an infinite cylinder (or cone) is kind of worthless when it comes to trying to bound it, because any rotation will cause it to expand to infinity. For completeness, though, it's worth defining anyway.

In object space, the cylinder extends from -1 to 1 in both x and z, and from -infinity to infinity in y:

cylinders.feature
Scenario: An unbounded cylinder has a bounding box Given shape cylinder() When box bounds_of(shape) Then box.min = point(-1, -infinity, -1) And box.max = point(1, infinity, 1)

Once you limit the cylinder (and especially if you limit it at both ends), it becomes much easier to bound:

cylinders.feature
Scenario: A bounded cylinder has a bounding box Given shape cylinder() And shape.minimum -5 And shape.maximum 3 When box bounds_of(shape) Then box.min = point(-1, -5, -1) And box.max = point(1, 3, 1)

Cones

Bounding boxes around infinite cones are perhaps the most useless of all, because cones extends to infinity along every axis.

cones.feature
Scenario: An unbounded cone has a bounding box Given shape cone() When box bounds_of(shape) Then box.min = point(-infinity, -infinity, -infinity) And box.max = point(infinity, infinity, infinity)

Limiting the cones makes them easier to bound, though the logic is more complicated than you might think at first glance.

cones.feature
Scenario: A bounded cone has a bounding box Given shape cone() And shape.minimum -5 And shape.maximum 3 When box bounds_of(shape) Then box.min = point(-5, -5, -5) And box.max = point(5, 3, 5)

To find the bounds of a limited cone, choose the largest of the absolute values of cone.minimum and cone.maximum. Call it limit. The cone's bounding box then goes from -limit to limit in both x and z, and from cone.minimum to cone.maximum in y.

In pseudocode:

function bounds_of(cone) let a abs(cone.minimum) let b abs(cone.maximum) let limit max(a, b) return bounding_box(min=point(-limit, cone.minimum, -limit) max=point(limit, cone.maximum, limit)) end function

Triangles

Finding the bounding box for a triangle is a matter of finding the smallest and largest x, y, and z components from its three points.

triangles.feature
Scenario: A triangle has a bounding box Given p1 point(-3, 7, 2) And p2 point(6, 2, -4) And p3 point(2, -1, -1) And shape triangle(p1, p2, p3) When box bounds_of(shape) Then box.min = point(-3, -1, -4) And box.max = point(6, 7, 2)

The (conceptually) simplest way to implement this is to just declare an empty bounding box, and then add each of the triangle's points to it:

function bounds_of(triangle) let box bounding_box(empty) add triangle.p1 to box add triangle.p2 to box add triangle.p3 to box return box end function

And there you have it!

The Test Shape

Lastly, it may seem odd to put bounds on an abstract shape that you only use in tests, but it'll be useful for some of the tests later in this chapter. For the sake of test_shape(), assume its bounds extend from -1 to 1 on each axis.

shapes.feature
Scenario: Test shape has (arbitrary) bounds Given shape test_shape() When box bounds_of(shape) Then box.min = point(-1, -1, -1) And box.max = point(1, 1, 1)

There's really no reason for those particular bounds. They're arbitrary, but they're also predictable and easy to work with.

More bounding box operations

Before you can add bounding boxes to groups and CSG objects, you're going to need to add a few more operations to your bounding boxes. Specifically, you need to be able to merge two bounding boxes, check to see if a box contains a given point, and check to see if a box contains another bounding box.

Let's start with merging one bounding box into another. This works very much like adding a point did, where the first box is resized sufficiently to contain the new box. Here's a test:

bounds.feature
Scenario: Adding one bounding box to another Given box1 bounding_box(min=point(-5, -2, 0) max=point(7, 4, 4)) And box2 bounding_box(min=point(8, -7, -2) max=point(14, 2, 8)) When box2 is added to box1 Then box1.min = point(-5, -7, -2) And box1.max = point(14, 4, 8)

Make this pass by adding the min and max points from one box to the other, causing the box to resize until it contains both points. In pseudocode:

# adding "box2" to "box1" add box2.min to box1 add box2.max to box1

Next, you'll implement the ability to check whether a box contains a given point, like this:

bounds.feature
Scenario Outline: Checking to see if a box contains a given point Given box bounding_box(min=point(5, -2, 0) max=point(11, 4, 7)) And p <point> Then box_contains_point(box, p) is <result> Examples: | point | result | | point(5, -2, 0) | true | | point(11, 4, 7) | true | | point(8, 1, 3) | true | | point(3, 0, 3) | false | | point(8, -4, 3) | false | | point(8, 1, -1) | false | | point(13, 1, 3) | false | | point(8, 5, 3) | false | | point(8, 1, 8) | false |

A box contains a point if each of the point's components lie between the corresponding min and max components, inclusive, like this:

function box_contains_point(box, point) return point.x is between box.min.x and box.max.x and point.y is between box.min.y and box.max.y and point.z is between box.min.z and box.max.z end function

Taking this one step further, you'll also need to check to see whether a box contains another box:

bounds.feature
Scenario Outline: Checking to see if a box contains a given box Given box bounding_box(min=point(5, -2, 0) max=point(11, 4, 7)) And box2 bounding_box(min=<min> max=<max>) Then box_contains_box(box, box2) is <result> Examples: | min | max | result | | point(5, -2, 0) | point(11, 4, 7) | true | | point(6, -1, 1) | point(10, 3, 6) | true | | point(4, -3, -1) | point(10, 3, 6) | false | | point(6, -1, 1) | point(12, 5, 8) | false |

A box lies within another box if both its min and max points lie within that box.

Make those tests all pass. Once you've got those nailed down, there's one more operation to tackle before you can put bounding boxes around groups and CSG objects: transformations.

Transforming bounding boxes

Remember that the bounds reported by a shape will be in object space. This is just fine, but when finding the bounds of an aggregate shape like a group or a CSG object, you need to convert each child's bounding box from its original object space, to parent space—the object space of the parent group or CSG shape.

Here's a test to start the discussion. It defines a bounding box, and then rotates it around the y and x axes.

bounds.feature
Scenario: Transforming a bounding box Given box bounding_box(min=point(-1, -1, -1) max=point(1, 1, 1)) And matrix rotation_x(π / 4) * rotation_y(π / 4) When box2 transform(box, matrix) Then box2.min = point(-1.4142, -1.7071, -1.7071) And box2.max = point(1.4142, 1.7071, 1.7071)

You might think transforming a bounding box is "just" a matter of multiplying the box's min and max points by some transformation matrix...and you'd almost be right, so long as the matrix only translates or scales.

Rotation, though, throws a wrench into the works. To see why, consider the following illustration of a bounding box being rotated 30º around the y axis.

If you construct a new bounding box using only the transformed min and max points, check out what happens:

The behavior we expect is that the transformed box ought to be able to contain the volume of the rotated cube, but as you can see, that's not what happens here. It's gone all skinny, and definitely no longer contains all the volume of the original cube. This is because rotating an AABB means it is no longer axis-aligned! A new axis-aligned bounding box needs to be computed for that rotated box.

The right way to do this is to actually transform the points at all eight corners of the cube, and then find a new bounding box that contains all eight transformed points. In pseudocode, it might look something like this:

function transform(bbox, matrix) let p1 bbox.min let p2 point(bbox.min.x, bbox.min.y, bbox.max.z) let p3 point(bbox.min.x, bbox.max.y, bbox.min.z) let p4 point(bbox.min.x, bbox.max.y, bbox.max.z) let p5 point(bbox.max.x, bbox.min.y, bbox.min.z) let p6 point(bbox.max.x, bbox.min.y, bbox.max.z) let p7 point(bbox.max.x, bbox.max.y, bbox.min.z) let p8 bbox.max let new_bbox bounding_box(empty) for p in (p1, p2, p3, p4, p5, p6, p7, p8) add (matrix * p) to new_bbox end for return new_bbox end function

Doing this ensures that the new bounding box is sufficiently large to hold the original cube after rotation.

Got that all passing? Awesome! You're finally ready to wrap some bounding boxes around aggregate shapes like groups and CSG objects. Read on!

Bounds for groups and CSG

So, now that you can transform a bounding box, you can implement a function for reporting a shape's bounds in the space of the shape's parent.

shapes.feature
Scenario: Querying a shape's bounding box in its parent's space Given shape sphere() And set_transform(shape, translation(1, -3, 5) * scaling(0.5, 2, 4)) When box parent_space_bounds_of(shape) Then box.min = point(0.5, -5, 1) And box.max = point(1.5, -1, 9)

This is accomplished by transforming the shape's bounding box by the shape's transformation matrix.

function parent_space_bounds_of(shape) return transform(bounds_of(shape), shape.transform) end function

With that function in hand, you're ready to find the bounds of a Group object. The following test sets up a group with a transformed sphere and cylinder, and shows that the bounds of the group are sufficient to include the child shapes.

groups.feature
Scenario: A group has a bounding box that contains its children Given s sphere() And set_transform(s, translation(2, 5, -3) * scaling(2, 2, 2)) And c cylinder() And c.minimum -2 And c.maximum 2 And set_transform(c, translation(-4, -1, 4) * scaling(0.5, 1, 0.5)) And shape group() And add_child(shape, s) And add_child(shape, c) When box bounds_of(shape) Then box.min = point(-4.5, -3, -5) And box.max = point(4, 7, 4.5)

CSG objects, being aggregate shapes like groups, should behave the same way, by reporting a bounding box sufficiently large to inculde their children.

csg.feature
Scenario: A CSG shape has a bounding box that contains its children Given left sphere() And right sphere() with: | transform | translation(2, 3, 4) | And shape csg("difference", left, right) When box bounds_of(shape) Then box.min = point(-1, -1, -1) And box.max = point(3, 4, 5)

Make these two tests pass by implementing the bounds_of(shape) function for each of Group and CSG. Both functions must find the parent-space bounds of each child object, and then merge them all together into a single bounding box. Here's the pseudocode for the Group version, to get you started:

function bounds_of(group) let box bounding_box(empty) for each child of group let cbox parent_space_bounds_of(child) add cbox to box end for return box end function

Got that? Great! Next up, you'll do a handy bit of refactoring to be able to intersect rays with these bounding boxes.

Intersecting a bounding box

Having a bounding box for every shape in your ray tracer is great and all, but it won't do you any good unless you also teach your ray tracer how to see if a ray intersects those boxes. Fortunately, you've already implemented this feature!

Well, almost.

Back in chapter 12, "Cubes", of The Ray Tracer Challenge, you wrote some code that tests a ray against an axis-aligned bounding box (specifically, a cube) centered at the origin. For this section, you'll take that code and refactor it slightly so that it works with AABB's of any dimension and position.

Start with the following test, which tests the intersection between a ray and a cube-shaped AABB at the origin.

bounds.feature
Scenario Outline: Intersecting a ray with a bounding box at the origin Given box bounding_box(min=point(-1, -1, -1) max=point(1, 1, 1)) And direction normalize(<direction>) And r ray(<origin>, direction) Then intersects(box, r) is <result> Examples: | origin | direction | result | | point(5, 0.5, 0) | vector(-1, 0, 0) | true | | point(-5, 0.5, 0) | vector(1, 0, 0) | true | | point(0.5, 5, 0) | vector(0, -1, 0) | true | | point(0.5, -5, 0) | vector(0, 1, 0) | true | | point(0.5, 0, 5) | vector(0, 0, -1) | true | | point(0.5, 0, -5) | vector(0, 0, 1) | true | | point(0, 0.5, 0) | vector(0, 0, 1) | true | | point(-2, 0, 0) | vector(2, 4, 6) | false | | point(0, -2, 0) | vector(6, 2, 4) | false | | point(0, 0, -2) | vector(4, 6, 2) | false | | point(2, 0, 2) | vector(0, 0, -1) | false | | point(0, 2, 2) | vector(0, -1, 0) | false | | point(2, 2, 0) | vector(-1, 0, 0) | false |

The intersects(box, ray) function returns a boolean value here. This is because (for the purposes of this feature) you don't need to know exactly where the bounding box is intersected, just whether or not it was. However, if you care about reusing code (and you should!), and you'd like your cube's local_intersect() function to be able to share this logic, you'll need to allow the function to report where the intersections occur. For this bonus chapter, though, the tests will only worry about true and false return values.

You can make that test pass by (re)using the logic from your cube's local_intersect() function. The next test, though, will require you to "massage" that code a bit, so that it can work with a bounding box that is not centered at the origin, and is not a perfect cube.

bounds.feature
Scenario Outline: Intersecting a ray with a non-cubic bounding box Given box bounding_box(min=point(5, -2, 0) max=point(11, 4, 7)) And direction normalize(<direction>) And r ray(<origin>, direction) Then intersects(box, r) is <result> Examples: | origin | direction | result | | point(15, 1, 2) | vector(-1, 0, 0) | true | | point(-5, -1, 4) | vector(1, 0, 0) | true | | point(7, 6, 5) | vector(0, -1, 0) | true | | point(9, -5, 6) | vector(0, 1, 0) | true | | point(8, 2, 12) | vector(0, 0, -1) | true | | point(6, 0, -5) | vector(0, 0, 1) | true | | point(8, 1, 3.5) | vector(0, 0, 1) | true | | point(9, -1, -8) | vector(2, 4, 6) | false | | point(8, 3, -4) | vector(6, 2, 4) | false | | point(9, -1, -2) | vector(4, 6, 2) | false | | point(4, 0, 9) | vector(0, 0, -1) | false | | point(8, 6, -1) | vector(0, -1, 0) | false | | point(12, 5, 4) | vector(-1, 0, 0) | false |

Your current check_axis() function (from the "Cubes" chapter), assumes that the cube extends from -1 to 1 along each axis, so it hard codes -1 and 1. To make it work for this more general algorithm, you'll need to pass the minimum and maximum extents of the bounding box for the axis in question. The following pseudocode adds the min and max parameters to the function. Careful, though: these are not the same as the min and max points of your bounding box. For this function, they are floating point values, just the x or y or z components from those points.

function check_axis(origin, direction, min, max) tmin_numerator = (min - origin) tmax_numerator = (max - origin) ...

The intersect function, then, sends the minimum and maximum values with each invocation of check_axis(), like this:

xtmin, xtmax check_axis(ray.origin.x, ray.direction.x, min.x, max.x) # and so forth for y and z axes

With those changes, your tests should all be passing. You're ready now to plug this in and make it all work together.</