The bottleneck: intersecting rays with objects.

Our scenes have consisted of a handful of objects (spheres, triangles) so performance hasn't really been an issue yet. But more complex objects require thousands, or even millions of triangles. A for-loop over all of the objects in our scene is very slow. For $w$ pixels along the width and $h$ pixels along the height of the image, $n$ objects (maybe $n$ triangles) will take $O(whn)$ time to render with this approach. The main question we want to anser today is: how can we improve this? In other words, how can we make our ray tracer faster?

An idea: first check if ray intersects sphere enclosing the model.

There are a lot of camera rays that never intersect our model, and the pixel color we set to them is just the background sky color. Here's an idea. We know how to intersect rays with spheres, so why don't we enclose our entire model within a sphere and first check if a camera ray intersects this sphere. We'll call this a bounding sphere. If a ray doesn't intersect this sphere, then there's no way it can intersect the model. This means we can already filter out the rays that don't intersect the model and save a loop over many triangles.

This is a good start, but it still lacks some efficiency. What if we are trying to render a giraffe? The model geometry is very long in one direction and very short in another. The bounding sphere of the giraffe will be huge because of the long neck of the giraffe. And there will be a lot of empty space in the sphere between the sides of the giraffe and the surface of the sphere.

A better idea: check if ray intersects bounding box enclosing the model.

To make this idea work better, let's try a box instead. We won't use a cube (equal lengths of all sides of the box) because that will have the same issues as the sphere (and probably worse). Instead, we'll find the tightest box around the giraffe and use that as our filter to know if a ray has any chance of intersecting our model.

This type of box is called an axis-aligned bounding box which we will just refer to as an AABB in the rest of the lecture. To define the AABB of our model, we can loop through all the points (each vertex of every triangle) in the model and find the minimum and maximum coordinates, which will define two corners of the box.

Let's define our AABB class that takes in these two corner points:

class AABB {
  constructor(min, max) {
    this.min = min;
    this.max = max;
  }
}

Intersection of a ray with an Axis-Aligned Bounding Box (AABB).

We want to use our AABB of the model to filter out the rays that definitely don't intersect the model. This means we need to know how to intersect rays with boxes. And, actually, we just need to know if a ray intersects a box (not the full intersection information since the box is not something we're going to shade). We know how to do ray-plane intersections (see the Lecture 2 notes), and a box can be interpreted as the intersection of six planes. First, let's represent our ray equation as:

$$ \vec{p}(t) = \vec{o} + \vec{r}t \quad \rightarrow \quad \left[\begin{array}{c}x(t) \\ y(t) \\ z(t)\end{array}\right] = \left[ \begin{array}{c}o_x\\ o_y\\ o_z \end{array}\right] + \left[ \begin{array}{c}r_xt\\ r_yt\\ r_zt \end{array}\right]. $$

The only difference between this and our previous description is that we're more generally describing the ray origin as $\vec{o}$ (not necessarily the eye location). We're also writing it in terms of each component since it'll be friendlier to work with below.

We'll focus on the two planes of the AABB whose normal is aligned with the x-direction: $\vec{n} = [1, 0, 0]^T$. The idea will be the same for the y and z planes. Assuming the ray is not parallel to the plane, then it will intersect the first plane for which all points have the lower x-componenent of the AABB ($x_l$) at:

$$ x_l = o_x + r_x t_{l,x}\quad \rightarrow \quad t_{l,x} = \frac{x_l - o_x}{r_x}. $$

Similarly for the other constant-x plane (where all points have an x-coordinate of $x_u$) of the box,

$$ x_u = o_x + r_x t_{u,x} \quad \rightarrow \quad t_{u,x} = \frac{x_u - o_x}{r_x}. $$

We then have the interval within which the ray passes between these two planes: $[t_{l,x}, t_{u,x}]$. Note that if $r_x \lt 0$, then the interval is reversed, so we'll use min and max to get the correct interval:

$$ t_{l,x} = \min\left(\frac{x_l - o_x}{r_x}, \frac{x_u - o_x}{r_x}\right), \quad t_{u,x} = \max\left(\frac{x_l - o_x}{r_x}, \frac{x_u - o_x}{r_x}\right). $$

It's possible that $r_x$ is equal to zero since a valid ray direction might have a zero in the x-component of its direction (ray is parallel to the constant-x plane). In this case, both $x_l$ and $x_u$ will be either $+\infty$ or $-\infty$ ($\pm$ Infinity in JavaScript). This is actually fine, and our min/max trick will handle this.

The intervals for all x, y and z axes need to overlap.

So we can determine the interval in which a ray intersects the two lower/upper planes in a particular axis direction, but we derived this for infinite planes - it doesn't tell us much about the bounding box. The key observation when applying this idea to the full box is: all three of these $t$ intervals (for each x, y and z direction) need to overlap if the ray intersects the box.

From "Ray Tracing: The Next Week" by Peter Shirley.

In the image above, there are two rays. Rays intersect the bounding box if the green and blue lines overlap. Therefore, in the image above, only the lower ray intersects the bounding box.

For our algorithm, we'll start with an assumption about the valid range of the ray intersection, like $(0, \infty)$. Then, we'll keep updating this range as we loop through each dimension. If at any point, the current maximum of the interval is less than the minimum, then there is no intersection, because there's no overlap.

So our intersect function of the AABB class looks like this:

AABB.prototype.intersect = function(ray, tmin, tmax) {
  for (let d = 0; d < 3; d++) {
    // calculate the interval for this dimension
    let tld = Math.min((this.min[d] - ray.origin[d]) / ray.direction[d], (this.max[d] - ray.origin[d]) / ray.direction[d]);
    let tud = Math.max((this.min[d] - ray.origin[d]) / ray.direction[d], (this.max[d] - ray.origin[d]) / ray.direction[d]);
    // update the interval
    tmin = Math.max(tld, tmin);
    tmax = Math.min(tud, tmax);
    // early return if no overlap condition is already satisfied
    if (tmax <= tmin) return undefined; // no intersection with AABB
  }
  return true; // the ray hits this AABB
}

Bounding Volume Hierarchy: A binary tree of bounding boxes.

Great! So we have slightly more efficient way of filtering out whether a ray intersects a model by first checking the bounding box. But it's still not as efficient as it could be. Let's expand upon this idea of bounding boxes.

What if we broke up our model geometry into different pieces, and put a bounding box over each of them. Then, maybe we split up those pieces into smaller pieces, each with a bounding box. For example, maybe we have a bounding box over the giraffe, and break up the giraffe into two pieces: upper body and lower body. Then maybe the upper body is split into one piece for the neck and another for the head, and then maybe the head piece is broken up into a nose, ears, etc. Each of these pieces would have an AABB associated with it (which we can calculate). If a ray does not intersect the upper-body AABB, there's no way it intersects the head AABB, and certainly no way it will intersect the ears. This filters out a lot of triangles to check!

We can represent this using a Bounding Volume Hierarchy, which is a binary tree of these AABBs. In the scene below, the top-level (root) BVH node is labelled as A which has an AABB enclosing all the shapes. A's left child is B which encloses 4 shapes and its right child is C which encloses two shapes.

From Wikipedia.

Building a BVH.

We'll use a top-down approach to building the BVH. For an incoming set of triangles to a particular BVH node, we'll build the two children (left and right) by splitting up the incoming array of triangles into two subarrays. We'll first sort the list according to the coordinates of each triangle along some chosen axis and then put the lower half in the left child and the upper half in the right node. We'll chose this axis randomly, but you can also alternate with each level in the tree (i.e. level % 3). You can also try to find a better split axis which will impact the efficiency of the traversal (for intersections).

Luckily, JavaScript has a built in sorting function for arrays (in-place) and we can pass a custom callback to ask it to sort an array with our own comparator. We'll write a custom comparator that compares two triangles by looking at their lower AABB corner along our random axis:

const axis = Math.floor(Math.random() * 3);
const compare = (triangleA, triangleB) => {
  return triangleA.box.min[axis] < triangleB.box.min[axis];
}

We'll then slice the incoming array of objects in half when constructing the left and right nodes:

const mid = Math.floor(objects.length / 2);
this.left = new BVHNode(objects.slice(0, mid));
this.right = new BVHNode(objects.slice(mid, objects.length));

We'll also need to handle the cases when there are either 1 or two objects in the incoming objects (triangles) array. These will be the base cases which terminate the BVH construction. When there is only one object, we'll assign it to both the left and right children.

Here is the documentation for the Array slice method. Note the end element is not included in the sliced subarray.

After the left and right children are constructed, we need to give this node an AABB as well since this will be used when we check for intersections. We could loop through all the objects to find the AABB, but this is unnecessary and inefficient. Instead, we'll use the bounding box of the left and right children to determine the bounding box of this node:

this.box = enclosingBox(this.left.box, this.right.box);

where the enclosingBox function creates a new AABB using the minimum and maximum of the two bounding box coordinates:

const enclosingBox = (boxL, boxR) => {
  const lo = vec3.fromValues(
    Math.min(boxL.min[0], boxR.min[0]),
    Math.min(boxL.min[1], boxR.min[1]),
    Math.min(boxL.min[2], boxR.min[2])
  );

  const hi = vec3.fromValues(
    Math.max(boxL.max[0], boxR.max[0]),
    Math.max(boxL.max[1], boxR.max[1]),
    Math.max(boxL.max[2], boxR.max[2])
  );
  return new AABB(lo, hi);
}

Checking for intersections.

Now we want to use our BVH to check for intersections. As with our earlier ideas, we'll start by checking for an intersection of the root's bounding box with the ray. If no intersection, return undefined as we have been doing. There is a slight difference with how we've been doing object/ray intersections so far. This time, we'll pass in the admissible interval of the ray parameter ($t$) values: tmin and tmax. This will help us filter out whether we need to check a particular branch, as well as define the admissible interval for the AABB intersect method.

If the ray hits the AABB of the current node, first check the left node:

const ixnL = this.left.intersect(ray, tmin, tmax);

If there was a hit, use the resulting t value to bound the check for the right node:

const ixnR = this.right.intersect(ray, tmin, ixnL ? ixnL.t : tmax);

This means we'll automatically ignore any intersection values that are larger than ixnL.t since that wouldn't be the closest anyway. Note that this assumes the object (triangle) eventually returns a ray t value for the intersection, which we have been doing.

If ixnR is not undefined, then we found a closer intersection point than ixnL (if there was one), so we'll prioritize ixnR over ixnL:

return ixnR || ixnL;

Note that a single ray-BVH intersection takes $O(\log n)$ time, which is the depth of the tree. This is much better than looping through all the triangles ($O(n)$)!

Resources for finding and editing 3d models.

Now that you can render your 3d models much faster, you can try finding your own models and creating your own scenes. Here are some resources for finding 3d models:

• TurboSquid: https://www.turbosquid.com/
• Sketchfab: https://sketchfab.com/
• CGTrader: https://www.cgtrader.com/

Look for files in .obj format since the scripts we have been developing use the webgl-obj-loader.

If you model still has too many triangles and takes too long to render, I would suggest reducing the number of triangles through a process called simplication/decimation. This is a technique for collapsing triangles that have a minimal effect on the overall shape. You can try the MeshLab software (download here: https://www.meshlab.net/#download):

1. Import your model (File -> Import Mesh).
2. Click on Filters -> Remeshing, Simplification and Reconstruction -> Simplification: Quadric Edge Collapse Decimation.
3. Select the number of triangles (faces) you want and click Apply.
4. Export your model (File -> Export Mesh As, select Alias Wavefront Object (.obj) ). You can uncheck all the options or click None in the bottom left corner of the options window.

Actually, Step 2 is one of the algorithms we implement in Geometric Modeling (CSCI 422).

Other options for speeding up your ray tracer.

Kd-Trees. We built a binary tree representation of our model by dividing the objects (triangles) at each node in the tree. However, we could have also divided the objects by dividing the space along some plane at each node. This is known as a Binary Space Partitioning (BSP) and a special type of BSP is one in which each plane is aligned with a coordinate axis, which is called a Kd-Tree. Kd-Trees are another popular acceleration structure for computer graphics and geometric modeling applications. For more information, please see the pbrt description here.

Parallelization. The camera ray that is sent through a particular pixel does not depend on any information from another pixel. Therefore, our ray tracing algorithm is very well-suited for parallel computation. We could use a multi-threaded approach on the CPU or even write our ray tracer in a language like CUDA or OpenCL for the GPU. Note that NVIDIA has developed the OptiX Ray Tracing Engine for developing ray tracing applications on their GPUs.

Complete Bounding Volume Hierarchy source code

Here is the complete source code for building a BVH and finding the closest intersection object. It is assumed that the incoming list of objects provides a intersect method similar to what we have been developing so far for triangles and spheres. One potential improvement is to pick a better splitting axis.

© Philip Claude Caplan, 2023 (Last updated: 2023-10-10)