From 8ed20d16bd546b9c6101711f2ee6fc95630cf372 Mon Sep 17 00:00:00 2001 From: Steve Hollasch Date: Thu, 24 Aug 2023 17:45:16 -0700 Subject: [PATCH] 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 --- CHANGELOG.md | 6 +- books/RayTracingTheNextWeek.html | 172 ++++++++++++++++++------------- src/TheNextWeek/aabb.h | 28 +++-- src/TheNextWeek/quad.h | 2 +- src/TheRestOfYourLife/aabb.h | 28 +++-- src/TheRestOfYourLife/quad.h | 2 +- 6 files changed, 142 insertions(+), 96 deletions(-) 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; }