Book2: Rework the AABB chapter

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

Author: hollasch

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)
 

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.
 
 <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} $$
 
@@ -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} $$
 
 <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>
 
 
@@ -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},
@@ -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 {
@@ -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
@@ -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 {
@@ -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
 -----------------------------
@@ -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

src/TheNextWeek/aabb.h

@@ -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) {
@@ -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);

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; }