Skip to content

Commit

Permalink
Book2: Rework the AABB chapter
Browse files Browse the repository at this point in the history 
In addition, I've changed the way that AABBs are padded. Originally, the
primitive constructing an AABB needed to immediately call the
aabb::pad() function after construction. Now, the AABB constructors call
aabb::pad_to_minimum() automatically to pad any dimension smaller than
some delta in size.

Resolves #1236
  • Loading branch information
@hollasch
hollasch committed Aug 30, 2023
1 parent 2ee8f9a commit 8ed20d1
Show file tree
Hide file tree
Showing 6 changed files with 142 additions and 96 deletions.
6 changes: 5 additions & 1 deletion CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,12 +16,16 @@ then.

### In One Weekend
- Change - Update reference to "Fundamentals of Interactive Computer Graphics" to "Computer
Graphics: Principles and Practice". This is the name used by newer editions of the book.
Graphics: Principles and Practice". This is the name used by newer editions of the
book.
- Change - New BVH optimization splits the bounds according to the longest bounding box dimension,
yielding a 15-20% speedup (#1007)

### The Next Week
- Change - `perlin::turb()` no longer defaults the value for the depth parameter.
- Change - AABB automatically pads to mininmum size for any dimension; no longer requires
primitives to call aabb::pad() function.
- Change - Reworked the AABB chapter (#1236)
- New - add section on alternative 2D primitives such as triangle, ellipse and annulus (#1204,
#1205)

Expand Down
172 changes: 101 additions & 71 deletions books/RayTracingTheNextWeek.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -374,22 +374,22 @@
sticking it in this chapter so the code can run faster, and because it refactors `hittable` a
little, and when I add rectangles and boxes we won't have to go back and refactor them.

The ray-object intersection is the main time-bottleneck in a ray tracer, and the time is linear with
the number of objects. But it’s a repeated search on the same model, so we ought to be able to make
Ray-object intersection is the main time-bottleneck in a ray tracer, and the run time is linear with
the number of objects. But it’s a repeated search on the same scene, so we ought to be able to make
it a logarithmic search in the spirit of binary search. Because we are sending millions to billions
of rays on the same model, we can do an analog of sorting the model, and then each ray intersection
can be a sublinear search. The two most common families of sorting are to 1) divide the space, and
2) divide the objects. The latter is usually much easier to code up and just as fast to run for most
models.
of rays into the same scene, we can sort the objects in the scene, and then each ray intersection
can be a sublinear search. The two most common methods of sorting are to 1) subdivide the space, and
2) subdivide the objects. The latter is usually much easier to code up, and just as fast to run for
most models.


The Key Idea
-------------
The key idea of a bounding volume over a set of primitives is to find a volume that fully encloses
(bounds) all the objects. For example, suppose you computed a sphere that bounds 10 objects. Any ray
that misses the bounding sphere definitely misses all ten objects inside. If the ray hits the
bounding sphere, then it might hit one of the ten objects. So the bounding code is always of the
form:
The key idea of creating bounding volumes for a set of primitives is to find a volume that fully
encloses (bounds) all the objects. For example, suppose you computed a sphere that bounds 10
objects. Any ray that misses the bounding sphere definitely misses all ten objects inside. If the
ray hits the bounding sphere, then it might hit one of the ten objects. So the bounding code is
always of the form:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (ray hits bounding object)
Expand All @@ -398,8 +398,9 @@
return false
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

A key thing is we are dividing objects into subsets. We are not dividing the screen or the volume.
Any object is in just one bounding volume, but bounding volumes can overlap.
Note that we will use these bounding volumes to group the objects in the scene into subgroups. We
are *not* dividing the screen or the scene space. We want any given object to be in just one
bounding volume, though bounding volumes can overlap.


Hierarchies of Bounding Volumes
Expand Down Expand Up @@ -432,33 +433,34 @@
To get that all to work we need a way to make good divisions, rather than bad ones, and a way to
intersect a ray with a bounding volume. A ray bounding volume intersection needs to be fast, and
bounding volumes need to be pretty compact. In practice for most models, axis-aligned boxes work
better than the alternatives, but this design choice is always something to keep in mind if you
encounter unusual types of models.
better than the alternatives (such as the spherical bounds mentioned above), but this design choice
is always something to keep in mind if you encounter other types of bounding models.

From now on we will call axis-aligned bounding rectangular parallelepiped (really, that is what they
need to be called if precise) _axis-aligned bounding boxes_, or AABBs. Any method you want to use to
intersect a ray with an AABB is fine. And all we need to know is whether or not we hit it; we don’t
need hit points or normals or any of the stuff we need to display the object.
From now on we will call axis-aligned bounding rectangular parallelepipeds (really, that is what
they need to be called if we're being precise) _axis-aligned bounding boxes_, or AABBs. (In the
code, you will also come across the naming abbreviation "bbox" for "bounding box".) Any method you
want to use to intersect a ray with an AABB is fine. And all we need to know is whether or not we
hit it; we don’t need hit points or normals or any of the stuff we need to display the object.

<div class='together'>
Most people use the “slab” method. This is based on the observation that an n-dimensional AABB is
just the intersection of $n$ axis-aligned intervals, often called “slabs”. An interval is just the
points between two endpoints, _e.g._, $x$ such that $3 < x < 5$, or more succinctly $x$ in $(3,5)$.
In 2D, two intervals overlapping makes a 2D AABB (a rectangle):
just the intersection of $n$ axis-aligned intervals, often called “slabs”. Recall that an interval
is just the points within two endpoints, for example, $x$ such that $3 \leq x \leq 5$, or more
succinctly $x$ in $[3,5]$. In 2D, an AABB (a rectangle) is defined by the overlap two intervals:

![Figure [2d-aabb]: 2D axis-aligned bounding box](../images/fig-2.02-2d-aabb.jpg)

</div>

For a ray to hit one interval we first need to figure out whether the ray hits the boundaries. For
example, again in 2D, this is the ray parameters $t_0$ and $t_1$. (If the ray is parallel to the
plane, its intersection with the plane will be undefined.)
To determine if a ray hits one interval, we first need to figure out whether the ray hits the
boundaries. For example, in 1D, ray intersection with two planes will yield the ray parameters $t_0$
and $t_1$. (If the ray is parallel to the planes, its intersection with any plane will be
undefined.)

![Figure [ray-slab]: Ray-slab intersection](../images/fig-2.03-ray-slab.jpg)

In 3D, those boundaries are planes. The equations for the planes are $x = x_0$ and $x = x_1$. Where
does the ray hit that plane? Recall that the ray can be thought of as just a function that given a
$t$ returns a location $\mathbf{P}(t)$:
How do we find the intersection between a ray and a plane? Recall that the ray is just defined by a
function that--given a parameter $t$--returns a location $\mathbf{P}(t)$:

$$ \mathbf{P}(t) = \mathbf{A} + t \mathbf{b} $$

Expand All @@ -467,7 +469,7 @@

$$ x_0 = A_x + t_0 b_x $$

Thus $t$ at that hitpoint is:
So $t$ at the intersection is given by

$$ t_0 = \frac{x_0 - A_x}{b_x} $$

Expand All @@ -476,12 +478,17 @@
$$ t_1 = \frac{x_1 - A_x}{b_x} $$

<div class='together'>
The key observation to turn that 1D math into a hit test is that for a hit, the $t$-intervals need
to overlap. For example, in 2D the green and blue overlapping only happens if there is a hit:
The key observation to turn that 1D math into a 2D or 3D hit test is this: if a ray intersects the
box bounded by all pairs of planes, then all $t$-intervals will overlap. For example, in 2D the
green and blue overlapping only happens if the ray intersects the bounded box:

![Figure [ray-slab-interval]: Ray-slab $t$-interval overlap
](../images/fig-2.04-ray-slab-interval.jpg)

In this figure, the upper ray intervals to not overlap, so we know the ray does _not_ hit the 2D box
bounded by the green and blue planes. The lower ray intervals _do_ overlap, so we know the lower ray
_does_ hit the bounded box.

</div>


Expand All @@ -490,41 +497,38 @@
The following pseudocode determines whether the $t$ intervals in the slab overlap:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
compute (tx0, tx1)
compute (ty0, ty1)
return overlap?( (tx0, tx1), (ty0, ty1))
interval_x &LeftArrow; compute_intersection_x (ray, x0, x1)
interval_y &LeftArrow; compute_intersection_y (ray, y0, y1)
return overlaps(interval_x, interval_y)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

<div class='together'>
That is awesomely simple, and the fact that the 3D version also works is why people love the slab
method:
That is awesomely simple, and the fact that the 3D version trivially extends the above is why people
love the slab method:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
compute (tx0, tx1)
compute (ty0, ty1)
compute (tz0, tz1)
return overlap ? ((tx0, tx1), (ty0, ty1), (tz0, tz1))
interval_x &LeftArrow; compute_intersection_x (ray, x0, x1)
interval_y &LeftArrow; compute_intersection_y (ray, y0, y1)
interval_z &LeftArrow; compute_intersection_z (ray, z0, z1)
return overlaps(interval_x, interval_y, interval_z)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

</div>

There are some caveats that make this less pretty than it first appears. First, suppose the ray is
travelling in the negative $\mathbf{x}$ direction. The interval $(t_{x0}, t_{x1})$ as computed above
might be reversed, _e.g._ something like $(7, 3)$. Second, the divide in there could give us
infinities. And if the ray origin is on one of the slab boundaries, we can get a `NaN`. There are
many ways these issues are dealt with in various ray tracers’ AABB. (There are also vectorization
issues like SIMD which we will not discuss here. Ingo Wald’s papers are a great place to start if
you want to go the extra mile in vectorization for speed.) For our purposes, this is unlikely to be
a major bottleneck as long as we make it reasonably fast, so let’s go for simplest, which is often
fastest anyway! First let’s look at computing the intervals:
There are some caveats that make this less pretty than it first appears. Consider again the 1D
equations for $t_0$ and $t_1$:

$$ t_{x0} = \frac{x_0 - A_x}{b_x} $$
$$ t_{x1} = \frac{x_1 - A_x}{b_x} $$
$$ t_0 = \frac{x_0 - A_x}{b_x} $$
$$ t_1 = \frac{x_1 - A_x}{b_x} $$

One troublesome thing is that perfectly valid rays will have $b_x = 0$, causing division by zero.
Some of those rays are inside the slab, and some are not. Also, the zero will have a ± sign when
using IEEE floating point. The good news for $b_x = 0$ is that $t_{x0}$ and $t_{x1}$ will both be +∞
or both be -∞ if not between $x_0$ and $x_1$. So, using min and max should get us the right answers:
First, suppose the ray is traveling in the negative $\mathbf{x}$ direction. The interval $(t_{x0},
t_{x1})$ as computed above might be reversed, like $(7, 3)$ for example. Second, the denominator
$b_x$ could be zero, yielding infinite values. And if the ray origin lies on one of the slab
boundaries, we can get a `NaN`, since both the numerator and the denominator can be zero. Also, the
zero will have a ± sign when using IEEE floating point.

The good news for $b_x = 0$ is that $t_{x0}$ and $t_{x1}$ will be equal: both +∞ or -∞, if not
between $x_0$ and $x_1$. So, using min and max should get us the right answers:

$$ t_{x0} = \min(
\frac{x_0 - A_x}{b_x},
Expand All @@ -536,25 +540,29 @@
\frac{x_1 - A_x}{b_x})
$$

The remaining troublesome case if we do that is if $b_x = 0$ and either $x_0 - A_x = 0$ or $x_1 -
A_x = 0$ so we get a `NaN`. In that case we can probably accept either hit or no hit answer, but
we’ll revisit that later.
The remaining troublesome case if we do that is if $b_x = 0$ and either $x_0 - A_x = 0$ or
$x_1 - A_x = 0$ so we get a `NaN`. In that case we can arbitrarily interpret that as either hit or
no hit, but we’ll revisit that later.

Now, let’s look at that overlap function. Suppose we can assume the intervals are not reversed (so
the first value is less than the second value in the interval) and we want to return true in that
case. The boolean overlap that also computes the overlap interval $(f, F)$ of intervals $(d, D)$ and
$(e, E)$ would be:
Now, let’s look at the pseudo-function `overlaps`. Suppose we can assume the intervals are not
reversed, and we want to return true when the intervals overlap. The boolean `overlaps()` function
computes the overlap of the $t$ intervals `t_interval1` and `t_interval2`, and uses that to
determine if that overlap is non-empty:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
bool overlap(d, D, e, E, f, F)
f = max(d, e)
F = min(D, E)
return (f < F)
bool overlaps(t_interval1, t_interval2)
t_min &LeftArrow; max(t_interval1.min, t_interval2.min)
t_max &LeftArrow; min(t_interval1.max, t_interval2.max)
return t_min < t_max
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

If there are any `NaN`s running around there, the compare will return false so we need to be sure
our bounding boxes have a little padding if we care about grazing cases (and we probably should
because in a ray tracer all cases come up eventually). Here's the implementation:
If there are any `NaN`s running around there, the compare will return false, so we need to be sure
our bounding boxes have a little padding if we care about grazing cases (and we probably _should_
because in a ray tracer all cases come up eventually).

<div class='together'>
To accomplish this, we'll first add a new `interval` function `expand`, which pads an interval by a
given amount:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class interval {
Expand All @@ -577,6 +585,10 @@
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [interval-expand]: <kbd>[interval.h]</kbd> interval::expand() method]

</div>

<div class='together'>
Now we have everything we need to implment the new AABB class.

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#ifndef AABB_H
Expand All @@ -591,14 +603,19 @@
aabb() {} // The default AABB is empty, since intervals are empty by default.

aabb(const interval& ix, const interval& iy, const interval& iz)
: x(ix), y(iy), z(iz) { }
: x(ix), y(iy), z(iz)
{
pad_to_minimums();
}

aabb(const point3& a, const point3& b) {
// Treat the two points a and b as extrema for the bounding box, so we don't require a
// particular minimum/maximum coordinate order.
x = interval(fmin(a[0],b[0]), fmax(a[0],b[0]));
y = interval(fmin(a[1],b[1]), fmax(a[1],b[1]));
z = interval(fmin(a[2],b[2]), fmax(a[2],b[2]));

pad_to_minimums();
}

const interval& axis(int n) const {
Expand All @@ -620,12 +637,25 @@
}
return true;
}

private:

void pad_to_minimums() {
// Adjust the AABB so that no side is narrower than some delta, padding if necessary.

double delta = 0.0001;
if (x.size() < delta) x = x.expand(delta);
if (y.size() < delta) y = y.expand(delta);
if (z.size() < delta) z = z.expand(delta);
}
};

#endif
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [aabb]: <kbd>[aabb.h]</kbd> Axis-aligned bounding box class]

</div>


An Optimized AABB Hit Method
-----------------------------
Expand Down Expand Up @@ -659,7 +689,7 @@
...
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [aabb-hit]: <kbd>[aabb.h]</kbd> Axis-aligned bounding box hit function]
[Listing [aabb-hit]: <kbd>[aabb.h]</kbd> Optional optimized AABB hit function]


Constructing Bounding Boxes for Hittables
Expand Down
28 changes: 17 additions & 11 deletions src/TheNextWeek/aabb.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,14 +21,19 @@ class aabb {
aabb() {} // The default AABB is empty, since intervals are empty by default.

aabb(const interval& ix, const interval& iy, const interval& iz)
: x(ix), y(iy), z(iz) { }
: x(ix), y(iy), z(iz)
{
pad_to_minimums();
}

aabb(const point3& a, const point3& b) {
// Treat the two points a and b as extrema for the bounding box, so we don't require a
// particular minimum/maximum coordinate order.
x = interval(fmin(a[0],b[0]), fmax(a[0],b[0]));
y = interval(fmin(a[1],b[1]), fmax(a[1],b[1]));
z = interval(fmin(a[2],b[2]), fmax(a[2],b[2]));

pad_to_minimums();
}

aabb(const aabb& box0, const aabb& box1) {
Expand All @@ -37,16 +42,6 @@ class aabb {
z = interval(box0.z, box1.z);
}

aabb pad() {
// Return an AABB that has no side narrower than some delta, padding if necessary.
double delta = 0.0001;
interval new_x = (x.size() >= delta) ? x : x.expand(delta);
interval new_y = (y.size() >= delta) ? y : y.expand(delta);
interval new_z = (z.size() >= delta) ? z : z.expand(delta);

return aabb(new_x, new_y, new_z);
}

const interval& axis(int n) const {
if (n == 1) return y;
if (n == 2) return z;
Expand Down Expand Up @@ -83,6 +78,17 @@ class aabb {
}

static const aabb empty, universe;

private:

void pad_to_minimums() {
// Adjust the AABB so that no side is narrower than some delta, padding if necessary.

double delta = 0.0001;
if (x.size() < delta) x = x.expand(delta);
if (y.size() < delta) y = y.expand(delta);
if (z.size() < delta) z = z.expand(delta);
}
};

const aabb aabb::empty = aabb(interval::empty, interval::empty, interval::empty);
Expand Down
2 changes: 1 addition & 1 deletion src/TheNextWeek/quad.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ class quad : public hittable {
}

virtual void set_bounding_box() {
bbox = aabb(Q, Q + u + v).pad();
bbox = aabb(Q, Q + u + v);
}

aabb bounding_box() const override { return bbox; }
Expand Down
Loading

0 comments on commit 8ed20d1

Please sign in to comment.