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

should have an account at the Ray Tracer Challenge forum. Ask all your questions there, and see how others have managed!

must have a copy of The Ray Tracer Challenge.

must have implemented your ray tracer through chapter 14 (Groups).

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!

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:

- Implement a
`BoundingBox`

structure and some associated logic to manipulate it. - Add a function to each primitive shape for querying the shape's bounds (in object space).
- Implement a function for recursively finding the bounds of a
`Group`

object. - Implement a function for recursively finding the bounds of a
`CSG`

object. - Detect if a ray intersects a bounding box.
- 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!

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.featureScenario: 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.featureScenario: 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.featureScenario: 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.

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

- spheres
- planes
- cubes
- cylinders
- cones

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.

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.featureScenario: 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 are trickier, because they stretch to infinity in `x`

and `z`

, and have zero thickness in `y`

.

planes.featureScenario: 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 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.featureScenario: 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?

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.featureScenario: 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.featureScenario: 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)

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

cones.featureScenario: 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.featureScenario: 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

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.featureScenario: 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!

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.featureScenario: 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.

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.featureScenario: 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.featureScenario 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.featureScenario 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*.

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.featureScenario: 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!

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.featureScenario: 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.featureScenario: 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.featureScenario: 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.

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.featureScenario 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.featureScenario 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.</