diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2e68b47f..fd02b8b6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -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)
diff --git a/books/RayTracingTheNextWeek.html b/books/RayTracingTheNextWeek.html
index ef7641bb..1330f665 100644
--- a/books/RayTracingTheNextWeek.html
+++ b/books/RayTracingTheNextWeek.html
@@ -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)
@@ -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
@@ -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.
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)
-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} $$
@@ -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} $$
@@ -476,12 +478,17 @@
$$ t_1 = \frac{x_1 - A_x}{b_x} $$
-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.
+
@@ -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 ← compute_intersection_x (ray, x0, x1)
+ interval_y ← compute_intersection_y (ray, y0, y1)
+ return overlaps(interval_x, interval_y)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-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 ← compute_intersection_x (ray, x0, x1)
+ interval_y ← compute_intersection_y (ray, y0, y1)
+ interval_z ← compute_intersection_z (ray, z0, z1)
+ return overlaps(interval_x, interval_y, interval_z)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-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},
@@ -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 ← max(t_interval1.min, t_interval2.min)
+ t_max ← 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).
+
+
+To accomplish this, we'll first add a new `interval` function `expand`, which pads an interval by a
+given amount:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
class interval {
@@ -577,6 +585,10 @@
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[Listing [interval-expand]: [interval.h] interval::expand() method]
+
+
+
+Now we have everything we need to implment the new AABB class.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
#ifndef AABB_H
@@ -591,7 +603,10 @@
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
@@ -599,6 +614,8 @@
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 {
@@ -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]: [aabb.h] Axis-aligned bounding box class]
+
+
An Optimized AABB Hit Method
-----------------------------
@@ -659,7 +689,7 @@
...
};
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- [Listing [aabb-hit]: [aabb.h] Axis-aligned bounding box hit function]
+ [Listing [aabb-hit]: [aabb.h] Optional optimized AABB hit function]
Constructing Bounding Boxes for Hittables
diff --git a/src/TheNextWeek/aabb.h b/src/TheNextWeek/aabb.h
index 5ffa38c7..64034497 100644
--- a/src/TheNextWeek/aabb.h
+++ b/src/TheNextWeek/aabb.h
@@ -21,7 +21,10 @@ 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
@@ -29,6 +32,8 @@ class aabb {
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) {
@@ -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;
@@ -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);
diff --git a/src/TheNextWeek/quad.h b/src/TheNextWeek/quad.h
index 5eb73ce6..73f33640 100644
--- a/src/TheNextWeek/quad.h
+++ b/src/TheNextWeek/quad.h
@@ -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; }
diff --git a/src/TheRestOfYourLife/aabb.h b/src/TheRestOfYourLife/aabb.h
index 5ffa38c7..64034497 100644
--- a/src/TheRestOfYourLife/aabb.h
+++ b/src/TheRestOfYourLife/aabb.h
@@ -21,7 +21,10 @@ 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
@@ -29,6 +32,8 @@ class aabb {
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) {
@@ -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;
@@ -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);
diff --git a/src/TheRestOfYourLife/quad.h b/src/TheRestOfYourLife/quad.h
index 41f77b29..ae883981 100644
--- a/src/TheRestOfYourLife/quad.h
+++ b/src/TheRestOfYourLife/quad.h
@@ -30,7 +30,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; }