From c6f803f6f9d117005940744e640535ef5a076a75 Mon Sep 17 00:00:00 2001 From: Steve Hollasch Date: Wed, 21 Aug 2024 14:15:42 -0700 Subject: [PATCH] Improve random_unit_vector() The old method had a floating-point weakness in which all three vector components, when small enough, can yield a vector length that underflows to zero, leading to a bogus [+/- infinity, +/- infinity, +/- infinity] result. This change also eliminates the `random_in_unit_sphere()` function, and does everything inside the `random_unit_vector()` function, which allows us to compute the vector length only once and then re-use it for normalization. Resolves #1606 --- CHANGELOG.md | 1 + books/RayTracingInOneWeekend.html | 70 +++++++++++++++++-------------- src/InOneWeekend/vec3.h | 11 ++--- src/TheNextWeek/vec3.h | 11 ++--- src/TheRestOfYourLife/vec3.h | 11 ++--- 5 files changed, 51 insertions(+), 53 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f5a86a99..bf6149ca 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ Change Log / Ray Tracing in One Weekend - Fix -- Big improvement to print version listing font size (#1595) and more compact line height for code listings in both print and browser. - Change -- Include hittable.h from material.h; drop `hit_record` forward declaration (#1609) + - Fix -- Fixed possible bogus values from `random_unit_vector()` due to underflow (#1606) ### In One Weekend - Fix -- Fixed usage of the term "unit cube" for a cube of diameter two (#1555, #1603) diff --git a/books/RayTracingInOneWeekend.html b/books/RayTracingInOneWeekend.html index fb78509e..75d93fec 100644 --- a/books/RayTracingInOneWeekend.html +++ b/books/RayTracingInOneWeekend.html @@ -2292,24 +2292,29 @@ surprisingly complicated to understand, and quite a bit complicated to implement. Instead, we'll use what is typically the easiest algorithm: A rejection method. A rejection method works by repeatedly generating random samples until we produce a sample that meets the desired criteria. In other words, -keep rejecting samples until you find a good one. +keep rejecting bad samples until you find a good one. There are many equally valid ways of generating a random vector on a hemisphere using the rejection method, but for our purposes we will go with the simplest, which is: -1. Generate a random vector inside of the unit sphere -2. Normalize this vector +1. Generate a random vector inside the unit sphere +2. Normalize this vector to extend it to the sphere surface 3. Invert the normalized vector if it falls onto the wrong hemisphere
First, we will use a rejection method to generate the random vector inside the unit sphere (that is, a sphere of radius 1). Pick a random point inside the cube enclosing the unit sphere (that is, where -$x$, $y$, and $z$ are all in the range $[-1,+1]$). If this point lies outside (or on) the unit -sphere, then generate a new one until we find one that lies inside the unit sphere. +$x$, $y$, and $z$ are all in the range $[-1,+1]$). If this point lies outside the unit +sphere, then generate a new one until we find one that lies inside or on the unit sphere. - ![Figure [sphere-vec]: Two vectors were rejected before finding a good one + ![Figure [sphere-vec]: Two vectors were rejected before finding a good one (pre-normalization) ](../images/fig-1.11-sphere-vec.jpg) + ![Figure [sphere-vec]: The accepted random vector is normalized to produce a unit vector + ](../images/fig-1.12-sphere-unit-vec.jpg) + +Here's our first draft of the function: + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ ... @@ -2319,49 +2324,45 @@ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight - inline vec3 random_in_unit_sphere() { + inline vec3 random_unit_vector() { while (true) { auto p = vec3::random(-1,1); - if (p.length_squared() < 1) - return p; + auto lensq = p.length_squared(); + if (lensq <= 1) + return p / sqrt(lensq); } } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - [Listing [random-in-unit-sphere]: [vec3.h] The random_in_unit_sphere() function] + [Listing [random-in-unit-sphere]: [vec3.h] The random_unit_vector() function, version one]
-
-Once we have a random vector in the unit sphere we need to normalize it to get a vector _on_ the -unit sphere. +Sadly, we have a small floating-point abstraction leak to deal with. Since floating-point numbers +have finite precision, a very small value can underflow to zero when squared. So if all three +coordinates are small enough (that is, very near the center of the sphere), the norm of the vector +will be zero, and thus normalizing will yield the bogus vector $[\pm\infty, \pm\infty, \pm\infty]$. +To fix this, we'll also reject points that lie inside this "black hole" around the center. With +double precision (64-bit floats), we can safely support values greater than $10^{-160}$. - ![Figure [sphere-vec]: The accepted random vector is normalized to produce a unit vector - ](../images/fig-1.12-sphere-unit-vec.jpg) +Here's our more robust function: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ - ... - - inline vec3 random_in_unit_sphere() { + inline vec3 random_unit_vector() { while (true) { auto p = vec3::random(-1,1); - if (p.length_squared() < 1) - return p; - } - } - - + auto lensq = p.length_squared(); ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight - inline vec3 random_unit_vector() { - return unit_vector(random_in_unit_sphere()); + if (1e-160 < lensq && lensq <= 1) + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ + return p / sqrt(lensq); + } } ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - [Listing [random-unit-vec]: [vec3.h] Random vector on the unit sphere] - -
+ [Listing [random-in-unit-sphere]: [vec3.h] The random_unit_vector() function, version one]
-And now that we have a random vector on the surface of the unit sphere, we can determine if it is on -the correct hemisphere by comparing against the surface normal: +Now that we have a random vector on the surface of the unit sphere, we can determine if it is on the +correct hemisphere by comparing against the surface normal: ![Figure [normal-hor]: The normal vector tells us which hemisphere we need ](../images/fig-1.13-surface-normal.jpg) @@ -2377,7 +2378,12 @@ ... inline vec3 random_unit_vector() { - return unit_vector(random_in_unit_sphere()); + while (true) { + auto p = vec3::random(-1,1); + auto lensq = p.length_squared(); + if (1e-160 < lensq && lensq <= 1) + return p / sqrt(lensq); + } } diff --git a/src/InOneWeekend/vec3.h b/src/InOneWeekend/vec3.h index 3afc89d8..bdf4dcf7 100644 --- a/src/InOneWeekend/vec3.h +++ b/src/InOneWeekend/vec3.h @@ -121,18 +121,15 @@ inline vec3 random_in_unit_disk() { } } -inline vec3 random_in_unit_sphere() { +inline vec3 random_unit_vector() { while (true) { auto p = vec3::random(-1,1); - if (p.length_squared() < 1) - return p; + auto lensq = p.length_squared(); + if (1e-160 < lensq && lensq <= 1.0) + return p / sqrt(lensq); } } -inline vec3 random_unit_vector() { - return unit_vector(random_in_unit_sphere()); -} - inline vec3 random_on_hemisphere(const vec3& normal) { vec3 on_unit_sphere = random_unit_vector(); if (dot(on_unit_sphere, normal) > 0.0) // In the same hemisphere as the normal diff --git a/src/TheNextWeek/vec3.h b/src/TheNextWeek/vec3.h index 3afc89d8..bdf4dcf7 100644 --- a/src/TheNextWeek/vec3.h +++ b/src/TheNextWeek/vec3.h @@ -121,18 +121,15 @@ inline vec3 random_in_unit_disk() { } } -inline vec3 random_in_unit_sphere() { +inline vec3 random_unit_vector() { while (true) { auto p = vec3::random(-1,1); - if (p.length_squared() < 1) - return p; + auto lensq = p.length_squared(); + if (1e-160 < lensq && lensq <= 1.0) + return p / sqrt(lensq); } } -inline vec3 random_unit_vector() { - return unit_vector(random_in_unit_sphere()); -} - inline vec3 random_on_hemisphere(const vec3& normal) { vec3 on_unit_sphere = random_unit_vector(); if (dot(on_unit_sphere, normal) > 0.0) // In the same hemisphere as the normal diff --git a/src/TheRestOfYourLife/vec3.h b/src/TheRestOfYourLife/vec3.h index c6c2497c..0ab08e31 100644 --- a/src/TheRestOfYourLife/vec3.h +++ b/src/TheRestOfYourLife/vec3.h @@ -121,18 +121,15 @@ inline vec3 random_in_unit_disk() { } } -inline vec3 random_in_unit_sphere() { +inline vec3 random_unit_vector() { while (true) { auto p = vec3::random(-1,1); - if (p.length_squared() < 1) - return p; + auto lensq = p.length_squared(); + if (1e-160 < lensq && lensq <= 1.0) + return p / sqrt(lensq); } } -inline vec3 random_unit_vector() { - return unit_vector(random_in_unit_sphere()); -} - inline vec3 random_on_hemisphere(const vec3& normal) { vec3 on_unit_sphere = random_unit_vector(); if (dot(on_unit_sphere, normal) > 0.0) // In the same hemisphere as the normal