3.9 Managing Rounding Error

Thus far, we’ve been discussing ray–shape intersection algorithms purely with respect to idealized arithmetic operations based on the real numbers. This approach has gotten us far, although the fact that computers can only represent finite quantities and therefore can’t actually represent all of the real numbers is important. In place of real numbers, computers use floating-point numbers, which have fixed storage requirements. However, error may be introduced each time a floating-point operation is performed, since the result may not be representable in the designated amount of memory.

The accumulation of this error has a few implications for the accuracy of intersection tests. First, it’s possible that it will cause valid intersections to be missed completely—for example, if a computed intersection’s t value is negative even though the precise value is positive. Furthermore, computed ray–shape intersection points may be above or below the actual surface of the shape. This leads to a problem: when new rays are traced starting from computed intersection points for shadow rays and reflection rays, if the ray origin is below the actual surface, we may find an incorrect re-intersection with the surface. Conversely, if the origin is too far above the surface, shadows and reflections may appear detached. (See Figure 3.39.)

Figure 3.39: Geometric Settings for Rounding-Error Issues That Can Cause Visible Errors in Images. The incident ray on the left intersects the surface. On the left, the computed intersection point (black circle) is slightly below the surface and a too-low “epsilon” offsetting the origin of the shadow ray leads to an incorrect self-intersection, as the shadow ray origin (white circle) is still below the surface; thus the light is incorrectly determined to be occluded. On the right a too-high “epsilon” causes a valid intersection to be missed as the ray’s origin is past the occluding surface.

Typical practice to address this issue in ray tracing is to offset spawned rays by a fixed “ray epsilon” value, ignoring any intersections along the ray normal p Subscript Baseline plus t bold d closer than some t Subscript normal m normal i normal n value. Figure 3.40 shows why this approach requires fairly high t Subscript normal m normal i normal n values to work effectively: if the spawned ray is fairly oblique to the surface, incorrect ray intersections may occur quite some distance from the ray origin. Unfortunately, large t Subscript normal m normal i normal n values cause ray origins to be relatively far from the original intersection points, which in turn can cause valid nearby intersections to be missed, leading to loss of fine detail in shadows and reflections.

Figure 3.40: If the computed intersection point (filled circle) is below the surface and the spawned ray is oblique, incorrect re-intersections may occur some distance from the ray origin (open circle). If a minimum t value along the ray is used to discard nearby intersections, a relatively large t Subscript normal m normal i normal n is needed to handle oblique rays well.

In this section, we’ll introduce the ideas underlying floating-point arithmetic and describe techniques for analyzing the error in floating-point computations. We’ll then apply these methods to the ray–shape algorithms introduced earlier in this chapter and show how to compute ray intersection points with bounded error. This will allow us to conservatively position ray origins so that incorrect self-intersections are never found, while keeping ray origins extremely close to the actual intersection point so that incorrect misses are minimized. In turn, no additional “ray epsilon” values are needed.

3.9.1 Floating-Point Arithmetic

Computation must be performed on a finite representation of numbers that fits in a finite amount of memory; the infinite set of real numbers just can’t be represented on a computer. One such finite representation is fixed point, where given a 16-bit integer, for example, one might map it to positive real numbers by dividing by 256. This would allow us to represent the range left-bracket 0 comma 65535 slash 256 right-bracket equals left-bracket 0 comma 255 plus 255 slash 256 right-bracket with equal spacing of 1 slash 256 between values. Fixed-point numbers can be implemented efficiently using integer arithmetic operations (a property that made them popular on early PCs that didn’t support floating-point computation), but they suffer from a number of shortcomings; among them, the maximum number they can represent is limited, and they aren’t able to accurately represent very small numbers near zero.

An alternative representation for real numbers on computers is floating-point numbers. These are based on representing numbers with a sign, a significand, and an exponent: essentially, the same representation as scientific notation but with a fixed number of digits devoted to significand and exponent. (In the following, we will assume base-2 digits exclusively.) This representation makes it possible to represent and perform computations on numbers with a wide range of magnitudes while using a fixed amount of storage.

Programmers using floating-point arithmetic are generally aware that floating point is imprecise; this understanding sometimes leads to a belief that floating-point arithmetic is unpredictable. In this section we’ll see that floating-point arithmetic has a carefully designed foundation that in turn makes it possible to compute conservative bounds on the error introduced in a particular computation. For ray-tracing calculations, this error is often surprisingly small.

Modern CPUs and GPUs nearly ubiquitously implement a model of floating-point arithmetic based on a standard promulgated by the Institute of Electrical and Electronics Engineers (1985, 2008). (Henceforth when we refer to floats, we will specifically be referring to 32-bit floating-point numbers as specified by IEEE 754.) The IEEE 754 technical standard specifies the format of floating-point numbers in memory as well as specific rules for precision and rounding of floating-point computations; it is these rules that make it possible to reason rigorously about the error present in a given floating-point value.

Floating-Point Representation

The IEEE standard specifies that 32-bit floats are represented with a sign bit, 8 bits for the exponent, and 23 bits for the significand. With 8 bits, the exponent e ranges from 0 to 255; the actual exponent used, e Subscript normal b , is computed by biasing e :

e Subscript normal b Baseline equals e minus 127 period

The significand actually has 24 bits of precision when a normalized floating-point value is stored. When a number expressed with significand and exponent is normalized, there are no leading zeros in the significand. In binary, this means that the leading digit of the significand must be one; in turn, there’s no need to store this value explicitly. Thus, the implicit leading one digit with the 23 digits encoding the fractional part of the significand gives a total of 24 bits of precision.

Given a sign s equals plus-or-minus 1 , significand m , and exponent e , the corresponding floating-point value is

s times 1 period m times 2 Superscript e minus 127 Baseline period

For example, with a normalized significand, the floating-point number 6.5 is written as 1.101 Subscript 2 Baseline times 2 squared , where the 2 subscript denotes a base-2 value. (If binary decimals aren’t immediately intuitive, note that the first number to the right of the decimal contributes 2 Superscript negative 1 Baseline equals 1 slash 2 , and so forth.) Thus, we have

left-parenthesis 1 times 2 Superscript 0 Baseline plus 1 times 2 Superscript negative 1 Baseline plus 0 times 2 Superscript negative 2 Baseline plus 1 times 2 Superscript negative 3 Baseline right-parenthesis times 2 squared equals 1.625 times 2 squared equals 6.5 period

e Subscript normal b Baseline equals 2 , so e equals 129 equals 1000001 Subscript 2 and m equals 10100000000000000000000 Subscript 2 .

Floats are laid out in memory with the sign bit at the most significant bit of the 32-bit value (with negative signs encoded with a one bit), then the exponent, and the significand. Thus, for the value 6.5 the binary in-memory representation of the value is

0 10000001 10100000000000000000000 equals 40 normal d Baseline 00000 Subscript 16 Baseline period

Similarly the floating-point value 1.0 has m equals 0 ellipsis 0 Subscript 2 and e Subscript normal b Baseline equals 0 , so e equals 127 equals 01111111 Subscript 2 and its binary representation is:

0 01111111 00000000000000000000000 equals 3 normal f Baseline 800000 Subscript 16 Baseline period

This hexadecimal number is a value worth remembering, as it often comes up in memory dumps when debugging.

An implication of this representation is that the spacing between representable floats between two adjacent powers of two is uniform throughout the range. (It corresponds to increments of the significand bits by one.) In a range left-bracket 2 Superscript e Baseline comma 2 Superscript e plus 1 Baseline right-parenthesis , the spacing is

2 Superscript e minus 23 Baseline period
(3.6)

Thus, for floating-point numbers between 1 and 2, e equals 0 , and the spacing between floating-point values is 2 Superscript negative 23 Baseline almost-equals 1.19209 ellipsis times 10 Superscript negative 7 . This spacing is also referred to as the magnitude of a unit in last place (“ulp”); note that the magnitude of an ulp is determined by the floating-point value that it is with respect to—ulps are relatively larger at numbers with bigger magnitudes than they are at numbers with smaller magnitudes.

As we’ve described the representation so far, it’s impossible to exactly represent zero as a floating-point number. This is obviously an unacceptable state of affairs, so the minimum exponent e equals 0 , or e Subscript normal b Baseline equals negative 127 , is set aside for special treatment. With this exponent, the floating-point value is interpreted as not having the implicit leading one bit in the significand, which means that a significand of all zero bits results in

s times 0.0 midline-horizontal-ellipsis 0 Subscript 2 Baseline times 2 Superscript negative 127 Baseline equals 0 period

Eliminating the leading one significand bit also makes it possible to represent denormalized numbers: if the leading one was always present, then the smallest 32-bit float would be

1.0 midline-horizontal-ellipsis 0 Subscript 2 Baseline times 2 Superscript negative 127 Baseline almost-equals 5.8774718 times 10 Superscript negative 39 Baseline period

Without the leading one bit, the minimum value is

0.00 midline-horizontal-ellipsis 1 Subscript 2 Baseline times 2 Superscript negative 126 Baseline equals 2 Superscript negative 126 Baseline times 2 Superscript negative 23 Baseline almost-equals 1.4012985 times 10 Superscript negative 45 Baseline period

Providing some capability to represent these small values can make it possible to avoid needing to round very small values to zero.

Note that there is both a “positive” and “negative” zero value with this representation. This detail is mostly transparent to the programmer. For example, the standard guarantees that the comparison -0.0 == 0.0 evaluates to true, even though the in-memory representations of these two values are different.

The maximum exponent, e equals 255 , is also reserved for special treatment. Therefore, the largest regular floating-point value that can be represented has e equals 254 (or e Subscript normal b Baseline equals 127 ) and is approximately

3.402823 ellipsis times 10 Superscript 38 Baseline period

With e equals 255 , if the significand bits are all zero, the value corresponds to positive or negative infinity, according to the sign bit. Infinite values result when performing computations like 1 slash 0 in floating point, for example. Arithmetic operations with infinity result in infinity. For comparisons, positive infinity is larger than any non-infinite value and similarly for negative infinity.

The MaxFloat and Infinity constants are initialized to be the largest representable and “infinity” floating-point values, respectively. We make them available in separate constants so that code that uses these values doesn’t need to use the wordy C++ standard library calls to get their values.

<<Global Constants>>= 
static constexpr Float MaxFloat = std::numeric_limits<Float>::max(); static constexpr Float Infinity = std::numeric_limits<Float>::infinity();

With e Subscript normal b Baseline equals 255 , non-zero significand bits correspond to special NaN values, which result from operations like taking the square root of a negative number or trying to compute 0 slash 0 . NaNs propagate through computations: any arithmetic operation where one of the operands is a NaN itself always returns NaN. Thus, if a NaN emerges from a long chain of computations, we know that something went awry somewhere along the way. In debug builds, pbrt has many Assert() statements that check for NaN values, as we almost never expect them to come up in the regular course of events. Any comparison with a NaN value returns false; thus, checking for !(x == x) serves to check if a value is not a number. For clarity, we use the C++ standard library function std::isnan() to check for not-a-number values.

Utility Routines

For certain low-level operations, it can be useful to be able to interpret a floating-point value in terms of its constituent bits and to convert the bits representing a floating-point value to an actual float or double.

One natural approach to this would be to take a pointer to a value to be converted and cast it to a pointer to the other type:

float f = ...; uint32_t bits = *((uint32_t *)&f);

However, modern versions of C++ specify that it’s illegal to cast a pointer of one type, float, to a different type, uint32_t. (This restriction allows the compiler to optimize more aggressively in its analysis of whether two pointers may point to the same memory location, which can inhibit storing values in registers.)

Another common approach is to use a union with elements of both types, assigning to one type and reading from the other:

union FloatBits { float f; uint32_t ui; }; FloatBits fb; fb.f = ...; uint32_t bits = fb.ui;

This, too, is illegal: the C++ standard says that reading an element of a union different from the one last one assigned to is undefined behavior.

These conversions can be properly made using memcpy() to copy from a pointer to the source type to a pointer to the destination type:

<<Global Inline Functions>>= 
inline uint32_t FloatToBits(float f) { uint32_t ui; memcpy(&ui, &f, sizeof(float)); return ui; }

<<Global Inline Functions>>+=  
inline float BitsToFloat(uint32_t ui) { float f; memcpy(&f, &ui, sizeof(uint32_t)); return f; }

While a call to the memcpy() function may seem gratuitously expensive to avoid these issues, in practice good compilers turn this into a no-op and just reinterpret the contents of the register or memory as the other type. (Versions of these functions that convert between double and uint64_t are also available in pbrt but are similar and are therefore not included here.)

These conversions can be used to implement functions that bump a floating-point value up or down to the next greater or next smaller representable floating-point value. They are useful for some conservative rounding operations that we’ll need in code to follow. Thanks to the specifics of the in-memory representation of floats, these operations are quite efficient.

<<Global Inline Functions>>+=  
inline float NextFloatUp(float v) { <<Handle infinity and negative zero for NextFloatUp()>> 
if (std::isinf(v) && v > 0.) return v; if (v == -0.f) v = 0.f;
<<Advance v to next higher float>> 
uint32_t ui = FloatToBits(v); if (v >= 0) ++ui; else --ui; return BitsToFloat(ui);
}

There are two important special cases: if v is positive infinity, then this function just returns v unchanged. Negative zero is skipped forward to positive zero before continuing on to the code that advances the significand. This step must be handled explicitly, since the bit patterns for negative 0.0 and 0.0 aren’t adjacent.

<<Handle infinity and negative zero for NextFloatUp()>>= 
if (std::isinf(v) && v > 0.) return v; if (v == -0.f) v = 0.f;

Conceptually, given a floating-point value, we would like to increase the significand by one, where if the result overflows, the significand is reset to zero and the exponent is increased by one. Fortuitously, adding one to the in-memory integer representation of a float achieves this: because the exponent lies at the high bits above the significand, adding one to the low bit of the significand will cause a one to be carried all the way up into the exponent if the significand is all ones and otherwise will advance to the next higher significand for the current exponent. Note also that when the highest representable finite floating-point value’s bit representation is incremented, the bit pattern for positive floating-point infinity is the result.

For negative values, subtracting one from the bit representation similarly advances to the next value.

<<Advance v to next higher float>>= 
uint32_t ui = FloatToBits(v); if (v >= 0) ++ui; else --ui; return BitsToFloat(ui);

The NextFloatDown() function, not included here, follows the same logic but effectively in reverse. pbrt also provides versions of these functions for doubles.

Arithmetic Operations

IEEE 754 provides important guarantees about the properties of floating-point arithmetic: specifically, it guarantees that addition, subtraction, multiplication, division, and square root give the same results given the same inputs and that these results are the floating-point number that is closest to the result of the underlying computation if it had been performed in infinite-precision arithmetic. It is remarkable that this is possible on finite-precision digital computers at all; one of the achievements in IEEE 754 was the demonstration that this level of accuracy is possible and can be implemented fairly efficiently in hardware.

Using circled operators to denote floating-point arithmetic operations and monospace s monospace q monospace r monospace t for floating-point square root, these precision guarantees can be written as:

StartLayout 1st Row 1st Column a circled-plus b 2nd Column equals normal r normal o normal u normal n normal d left-parenthesis a plus b right-parenthesis 2nd Row 1st Column a minus b 2nd Column equals normal r normal o normal u normal n normal d left-parenthesis a minus b right-parenthesis 3rd Row 1st Column a circled-times b 2nd Column equals normal r normal o normal u normal n normal d left-parenthesis a asterisk b right-parenthesis 4th Row 1st Column a circled-division-slash b 2nd Column equals normal r normal o normal u normal n normal d left-parenthesis a slash b right-parenthesis 5th Row 1st Column monospace s monospace q monospace r monospace t monospace left-parenthesis monospace a monospace right-parenthesis 2nd Column equals normal r normal o normal u normal n normal d left-parenthesis StartRoot a EndRoot right-parenthesis EndLayout
(3.7)

where normal r normal o normal u normal n normal d left-parenthesis x right-parenthesis indicates the result of rounding a real number to the closest floating-point value.

This bound on the rounding error can also be represented with an interval of real numbers: for example, for addition, we can say that the rounded result is within an interval

StartLayout 1st Row 1st Column a circled-plus b 2nd Column equals normal r normal o normal u normal n normal d left-parenthesis a plus b right-parenthesis subset-of left-parenthesis a plus b right-parenthesis left-parenthesis 1 plus-or-minus epsilon right-parenthesis 2nd Row 1st Column Blank 2nd Column equals left-bracket left-parenthesis a plus b right-parenthesis left-parenthesis 1 minus epsilon right-parenthesis comma left-parenthesis a plus b right-parenthesis left-parenthesis 1 plus epsilon right-parenthesis right-bracket EndLayout
(3.8)

for some epsilon . The amount of error introduced from this rounding can be no more than half the floating-point spacing at a plus b —if it was more than half the floating-point spacing, then it would be possible to round to a different floating-point number with less error (Figure 3.41).

Figure 3.41: The IEEE standard specifies that floating-point calculations must be implemented as if the calculation was performed with infinite-precision real numbers and then rounded to the nearest representable float. Here, an infinite precision result in the real numbers is denoted by a filled dot, with the representable floats around it denoted by ticks on a number line. We can see that the error introduced by rounding to the nearest float, delta , can be no more than half the spacing between floats.

For 32-bit floats, we can bound the floating-point spacing at a plus b from above using Equation (3.6) (i.e., an ulp at that value) by left-parenthesis a plus b right-parenthesis 2 Superscript negative 23 , so half the spacing is bounded from above by left-parenthesis a plus b right-parenthesis 2 Superscript negative 24 and so StartAbsoluteValue epsilon EndAbsoluteValue less-than-or-equal-to 2 Superscript negative 24 . This bound is the machine epsilon. For 32-bit floats, epsilon Subscript normal m Baseline equals 2 Superscript negative 24 Baseline almost-equals 5.960464 ellipsis times 10 Superscript negative 8 .

<<Global Constants>>+=  
static constexpr Float MachineEpsilon = std::numeric_limits<Float>::epsilon() * 0.5;

Thus, we have

StartLayout 1st Row 1st Column a circled-plus b 2nd Column equals normal r normal o normal u normal n normal d left-parenthesis a plus b right-parenthesis subset-of left-parenthesis a plus b right-parenthesis left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column equals left-bracket left-parenthesis a plus b right-parenthesis left-parenthesis 1 minus epsilon Subscript normal m Baseline right-parenthesis comma left-parenthesis a plus b right-parenthesis left-parenthesis 1 plus epsilon Subscript normal m Baseline right-parenthesis right-bracket period EndLayout

Analogous relations hold for the other arithmetic operators and the square root operator.

A number of useful properties follow directly from Equation (3.7). For a floating-point number x ,

  • 1 circled-times x equals x .
  • x circled-division-slash x equals 1 .
  • x circled-plus 0 equals x .
  • x minus x equals 0 .
  • 2 circled-times x and x circled-division-slash 2 are exact; no rounding is performed to compute the final result. More generally, any multiplication by or division by a power of two gives an exact result (assuming there’s no overflow or underflow).
  • x circled-division-slash 2 Superscript i Baseline equals x circled-times 2 Superscript negative i for all integer i , assuming 2 Superscript i doesn’t overflow.

All of these properties follow from the principle that the result must be the nearest floating-point value to the actual result; when the result can be represented exactly, the exact result must be computed.

Error Propagation

Using the guarantees of IEEE floating-point arithmetic, it is possible to develop methods to analyze and bound the error in a given floating-point computation. For more details on this topic, see the excellent book by Higham (2002), as well as Wilkinson’s earlier classic (1994).

Two measurements of error are useful in this effort: absolute and relative. If we perform some floating-point computation and get a rounded result a overTilde , we say that the magnitude of the difference between a overTilde and the result of doing that computation in the real numbers is the absolute error, delta Subscript normal a :

delta Subscript normal a Baseline equals StartAbsoluteValue a overTilde minus a EndAbsoluteValue period

Relative error, delta Subscript normal r , is the ratio of the absolute error to the precise result:

delta Subscript normal r Baseline equals StartAbsoluteValue StartFraction a overTilde minus a Over a EndFraction EndAbsoluteValue equals StartAbsoluteValue StartFraction delta Subscript normal a Baseline Over a EndFraction EndAbsoluteValue comma
(3.9)

as long as a not-equals 0 . Using the definition of relative error, we can thus write the computed value a overTilde as a perturbation of the exact result a :

a overTilde equals a plus-or-minus delta Subscript normal a Baseline equals a left-parenthesis 1 plus-or-minus delta Subscript normal r Baseline right-parenthesis period

As a first application of these ideas, consider computing the sum of four numbers, a , b , c , and d , represented as floats. If we compute this sum as r = (((a + b) + c) + d), Equation (3.8) gives us

StartLayout 1st Row 1st Column left-parenthesis left-parenthesis left-parenthesis a circled-plus b right-parenthesis circled-plus c right-parenthesis circled-plus d right-parenthesis 2nd Column subset-of left-parenthesis left-parenthesis left-parenthesis left-parenthesis a plus b right-parenthesis left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis right-parenthesis plus c right-parenthesis left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis plus d right-parenthesis left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column equals left-parenthesis a plus b right-parenthesis left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis cubed plus c left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis squared plus d left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis period EndLayout

Because epsilon Subscript normal m is small, higher order powers of epsilon Subscript normal m can be bounded by an additional epsilon Subscript normal m term, and so we can bound the left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis Superscript n terms with

left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis Superscript n Baseline less-than-or-equal-to left-parenthesis 1 plus-or-minus left-parenthesis n plus 1 right-parenthesis epsilon Subscript normal m Baseline right-parenthesis period

(As a practical matter, left-parenthesis 1 plus-or-minus n epsilon Subscript normal m Baseline right-parenthesis almost bounds these terms, since higher powers of epsilon Subscript normal m get very small very quickly, but the above is a fully conservative bound.)

This bound lets us simplify the result of the addition to:

StartLayout 1st Row 1st Column left-parenthesis a plus b right-parenthesis left-parenthesis 1 plus-or-minus 4 epsilon Subscript normal m Baseline right-parenthesis 2nd Column plus c left-parenthesis 1 plus-or-minus 3 epsilon Subscript normal m Baseline right-parenthesis plus d left-parenthesis 1 plus-or-minus 2 epsilon Subscript normal m Baseline right-parenthesis equals 2nd Row 1st Column Blank 2nd Column a plus b plus c plus d plus left-bracket plus-or-minus 4 epsilon Subscript normal m Baseline left-parenthesis a plus b right-parenthesis plus-or-minus 3 epsilon Subscript normal m Baseline c plus-or-minus 2 epsilon Subscript normal m Baseline d right-bracket period EndLayout

The term in square brackets gives the absolute error: its magnitude is bounded by

4 epsilon Subscript normal m Baseline StartAbsoluteValue a plus b EndAbsoluteValue plus 3 epsilon Subscript normal m Baseline StartAbsoluteValue c EndAbsoluteValue plus 2 epsilon Subscript normal m Baseline StartAbsoluteValue d EndAbsoluteValue period

Thus, if we add four floating-point numbers together with the above parenthesization, we can be certain that the difference between the final rounded result and the result we would get if we added them with infinite-precision real numbers is bounded by Equation (3.10); this error bound is easily computed given specific values of a , b , c , and  d .

This is a fairly interesting result; we see that the magnitude of a plus b makes a relatively large contribution to the error bound, especially compared to d . (This result gives a sense for why, if adding a large number of floating-point numbers together, sorting them from small to large magnitudes generally gives a result with a lower final error than an arbitrary ordering.)

Our analysis here has implicitly assumed that the compiler would generate instructions according to the expression used to define the sum. Compilers are required to follow the form of the given floating-point expressions in order to not break carefully crafted computations that may have been designed to minimize round-off error. Here again is a case where certain transformations that would be valid on expressions with integers can not be safely applied when floats are involved.

What happens if we change the expression to the algebraically equivalent float r = (a + b) + (c + d)? This corresponds to the floating-point computation

left-parenthesis left-parenthesis a circled-plus b right-parenthesis circled-plus left-parenthesis c circled-plus d right-parenthesis right-parenthesis period

If we apply the same process of applying Equation (3.8), expanding out terms, converting higher-order left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis Superscript n terms to left-parenthesis 1 plus-or-minus left-parenthesis n plus 1 right-parenthesis epsilon Subscript normal m Baseline right-parenthesis , we get absolute error bounds of

3 epsilon Subscript normal m Baseline StartAbsoluteValue a plus b EndAbsoluteValue plus 3 epsilon Subscript normal m Baseline StartAbsoluteValue c plus d EndAbsoluteValue comma

which are lower than the first formulation if StartAbsoluteValue a plus b EndAbsoluteValue is relatively large, but possibly higher if StartAbsoluteValue d EndAbsoluteValue is relatively large.

This approach to computing error is known as forward error analysis; given inputs to a computation, we can apply a fairly mechanical process that provides conservative bounds on the error in the result. The derived bounds in the result may overstate the actual error—in practice, the signs of the error terms are often mixed, so that there is cancellation when they are added. An alternative approach is backward error analysis, which treats the computed result as exact and provides bounds on perturbations on the inputs that give the same result. This approach can be more useful when analyzing the stability of a numerical algorithm but is less applicable to deriving conservative error bounds on the geometric computations we’re interested in here.

The conservative bounding of left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis Superscript n by left-parenthesis 1 plus-or-minus left-parenthesis n plus 1 right-parenthesis epsilon Subscript normal m Baseline right-parenthesis is somewhat unsatisfying since it adds a whole epsilon Subscript normal m term purely to conservatively bound the sum of various higher powers of epsilon Subscript normal m . Higham (2002, Section 3.1) gives an approach to more tightly bound products of left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis error terms. If we have left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis Superscript n , it can be shown that this value is bounded by 1 plus theta Subscript n , where

StartAbsoluteValue theta Subscript n Baseline EndAbsoluteValue less-than-or-equal-to StartFraction n epsilon Subscript normal m Baseline Over 1 minus n epsilon Subscript normal m Baseline EndFraction comma

as long as n epsilon Subscript normal m Baseline less-than 1 (which will certainly be the case for the calculations we’re considering). Note that the denominator of this expression will be just less than one for reasonable n values, so it just barely increases n epsilon Subscript normal m to achieve a conservative bound.

We will denote this bound by gamma Subscript n :

gamma Subscript n Baseline equals StartFraction n epsilon Subscript normal m Baseline Over 1 minus n epsilon Subscript normal m Baseline EndFraction period

The function that computes its value is declared as constexpr so that any invocations with compile-time constants will generally be replaced with the corresponding floating-point return value.

<<Global Inline Functions>>+=  
inline constexpr Float gamma(int n) { return (n * MachineEpsilon) / (1 - n * MachineEpsilon); }

Using the gamma notation, our bound on the error of the sum of the four values is

StartAbsoluteValue a plus b EndAbsoluteValue gamma 3 plus StartAbsoluteValue c EndAbsoluteValue gamma 2 plus StartAbsoluteValue d EndAbsoluteValue gamma 1 period

An advantage of this approach is that quotients of left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis Superscript n terms can also be bounded with the gamma function. Given

StartFraction left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis Superscript m Baseline Over left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis Superscript n Baseline EndFraction comma

the interval is bounded by left-parenthesis 1 plus-or-minus gamma Subscript m plus n Baseline right-parenthesis . Thus, gamma can be used to collect epsilon Subscript normal m terms from both sides of an equality over to one side by dividing them through; this will be useful in some of the following derivations. (Note that because left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis terms represent intervals, canceling them would be incorrect:

StartFraction left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis Superscript m Baseline Over left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis Superscript n Baseline EndFraction not-equals left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis Superscript m minus n Baseline semicolon

the gamma Subscript m plus n bounds must be used instead.)

Given inputs to some computation that themselves carry some amount of error, it’s instructive to see how this error is carried through various elementary arithmetic operations. Given two values, a left-parenthesis 1 plus-or-minus gamma Subscript i Baseline right-parenthesis and b left-parenthesis 1 plus-or-minus gamma Subscript j Baseline right-parenthesis that each carry accumulated error from earlier operations, consider their product. Using the definition of circled-times , the result is in the interval:

a left-parenthesis 1 plus-or-minus gamma Subscript i Baseline right-parenthesis circled-times b left-parenthesis 1 plus-or-minus gamma Subscript j Baseline right-parenthesis subset-of a b left-parenthesis 1 plus-or-minus gamma Subscript i plus j plus 1 Baseline right-parenthesis comma

where we’ve used the relationship left-parenthesis 1 plus-or-minus gamma Subscript i Baseline right-parenthesis left-parenthesis 1 plus-or-minus gamma Subscript j Baseline right-parenthesis subset-of left-parenthesis 1 plus-or-minus gamma Subscript i plus j Baseline right-parenthesis , which follows directly from Equation (3.11).

The relative error in this result is bounded by:

StartAbsoluteValue StartFraction a b gamma Subscript i plus j plus 1 Baseline Over a b EndFraction EndAbsoluteValue equals gamma Subscript i plus j plus 1 Baseline comma

and so the final error is thus just roughly left-parenthesis i plus j plus 1 right-parenthesis slash 2 ulps at the value of the product—about as good as we might hope for given the error going into the multiplication. (The situation for division is similarly good.)

Unfortunately, with addition and subtraction, it’s possible for the relative error to increase substantially. Using the same definitions of the values being operated on, consider

a left-parenthesis 1 plus-or-minus gamma Subscript i Baseline right-parenthesis circled-plus b left-parenthesis 1 plus-or-minus gamma Subscript j Baseline right-parenthesis comma

which is in the interval a left-parenthesis 1 plus-or-minus gamma Subscript i plus 1 Baseline right-parenthesis plus b left-parenthesis 1 plus-or-minus gamma Subscript j plus 1 Baseline right-parenthesis comma and so the absolute error is bounded by StartAbsoluteValue a EndAbsoluteValue gamma Subscript i plus 1 plus StartAbsoluteValue b EndAbsoluteValue gamma Subscript j plus 1 .

If the signs of a and b are the same, then the absolute error is bounded by StartAbsoluteValue a plus b EndAbsoluteValue gamma Subscript i plus j plus 1 and the relative error is around left-parenthesis i plus j plus 1 right-parenthesis slash 2 ulps around the computed value.

However, if the signs of a and b differ (or, equivalently, they are the same but subtraction is performed), then the relative error can be quite high. Consider the case where a almost-equals negative b : the relative error is

StartFraction StartAbsoluteValue a EndAbsoluteValue gamma Subscript i plus 1 Baseline plus StartAbsoluteValue b EndAbsoluteValue gamma Subscript j plus 1 Baseline Over a plus b EndFraction almost-equals StartFraction 2 StartAbsoluteValue a EndAbsoluteValue gamma Subscript i plus j plus 1 Baseline Over a plus b EndFraction period

The numerator’s magnitude is proportional to the original value StartAbsoluteValue a EndAbsoluteValue yet is divided by a very small number, and thus the relative error is quite high. This substantial increase in relative error is called catastrophic cancellation. Equivalently, we can have a sense of the issue from the fact that the absolute error is in terms of the magnitude of StartAbsoluteValue a EndAbsoluteValue , though it’s now in relation to a value much smaller than a .

Running Error Analysis

In addition to working out error bounds algebraically, we can also have the computer do this work for us as some computation is being performed. This approach is known as running error analysis. The idea behind it is simple: each time a floating-point operation is performed, we also compute terms that compute intervals based on Equation (3.7) to compute a running bound on the error that has been accumulated so far. While this approach can have higher run-time overhead than deriving expressions that give an error bound directly, it can be convenient when derivations become unwieldy.

pbrt provides a simple EFloat class, which mostly acts like a regular float but uses operator overloading to provide all of the regular arithmetic operations on floats while computing these error bounds.

Similar to the Interval class from Chapter 2, EFloat keeps track of an interval that describes the uncertainty of a value of interest. In contrast to Interval, EFloat’s intervals arise due to errors in intermediate floating-point arithmetic rather than uncertainty of the input parameters.

<<EFloat Public Methods>>= 
EFloat() { } EFloat(float v, float err = 0.f) : v(v), err(err) { <<Store high-precision reference value in EFloat>> 
#ifndef NDEBUG ld = v; #endif // NDEBUG
}

EFloat maintains a computed value v and the absolute error bound, err.

<<EFloat Private Data>>= 
float v; float err;

In debug builds, EFloat also maintains a highly precise version of v that can be used as a reference value to compute an accurate approximation of the relative error. In optimized builds, we’d generally rather not pay the overhead for computing this additional value.

<<Store high-precision reference value in EFloat>>= 
#ifndef NDEBUG ld = v; #endif // NDEBUG

<<EFloat Private Data>>+= 
#ifndef NDEBUG long double ld; #endif // NDEBUG

The implementation of the addition operation for this class is essentially an implementation of the relevant definitions. We have:

StartLayout 1st Row 1st Column left-parenthesis a plus-or-minus delta Subscript a Baseline right-parenthesis circled-plus left-parenthesis b plus-or-minus delta Subscript b Baseline right-parenthesis 2nd Column equals left-parenthesis left-parenthesis a plus-or-minus delta Subscript a Baseline right-parenthesis plus left-parenthesis b plus-or-minus delta Subscript b Baseline right-parenthesis right-parenthesis left-parenthesis 1 plus-or-minus gamma 1 right-parenthesis 2nd Row 1st Column Blank 2nd Column equals a plus b plus left-bracket plus-or-minus delta Subscript a Baseline plus-or-minus delta Subscript b Baseline plus-or-minus left-parenthesis a plus b right-parenthesis gamma 1 plus-or-minus gamma 1 delta Subscript a Baseline plus-or-minus gamma 1 delta Subscript b Baseline right-bracket period EndLayout

And so the absolute error (in brackets) is bounded by

delta Subscript a Baseline plus delta Subscript b Baseline plus gamma 1 left-parenthesis StartAbsoluteValue a plus b EndAbsoluteValue plus delta Subscript a Baseline plus delta Subscript b Baseline right-parenthesis period

<<EFloat Public Methods>>+=  
EFloat operator+(EFloat f) const { EFloat r; r.v = v + f.v; #ifndef NDEBUG r.ld = ld + f.ld; #endif // DEBUG r.err = err + f.err + gamma(1) * (std::abs(v + f.v) + err + f.err); return r; }

The implementations for the other arithmetic operations for EFloat are analogous.

Note that this implementation neglects the issue that the computation of errors will itself be affected by rounding error. If this was a concern, we could switch the floating-point rounding mode so that it always rounded the error bounds up to positive infinity, but this tends to be a fairly expensive operation since it causes a full pipeline flush on current processors. Here, we use the default rounding mode; in the following, the error bounds are expanded by one ulp when they are used to account for this issue.

The float value in an EFloat is available via a type conversion operator; it has an explicit qualifier to require the caller to have an explicit (float) cast to extract the floating-point value. The requirement to use an explicit cast reduces the risk of an unintended round trip from EFloat to Float and back, thus losing the accumulated error bounds.

<<EFloat Public Methods>>+=  
explicit operator float() const { return v; }

If a series of computations is performed using EFloat rather than float-typed variables, then at any point in the computation, the GetAbsoluteError() method can be called to find a bound on the absolute error of the computed value.

<<EFloat Public Methods>>+=  
float GetAbsoluteError() const { return err; }

The bounds of the error interval are available via the UpperBound() and LowerBound() methods. Their implementations use NextFloatUp() and NextFloatDown() to expand the returned values by one ulp, respectively, ensuring that the interval is conservative.

<<EFloat Public Methods>>+=  
float UpperBound() const { return NextFloatUp(v + err); } float LowerBound() const { return NextFloatDown(v - err); }

In debug builds, methods are available to get both the relative error as well as the precise value maintained in ld.

<<EFloat Public Methods>>+= 
#ifndef NDEBUG float GetRelativeError() const { return std::abs((ld - v)/ld); } long double PreciseValue() const { return ld; } #endif

pbrt also provides a variant of the Quadratic() function that operates on coefficients that may have error and returns error bounds with the t0 and t1 values. The implementation is the same as the regular Quadratic() function, just using EFloat.

<<EFloat Inline Functions>>= 
inline bool Quadratic(EFloat A, EFloat B, EFloat C, EFloat *t0, EFloat *t1);

With the floating-point error fundamentals in place, we’ll now focus on using these tools to provide robust intersection operations.

3.9.2 Conservative Ray–Bounds Intersections

Floating-point round-off error can cause the ray–bounding box intersection test to miss cases where a ray actually does intersect the box. While it’s acceptable to have occasional false positives from ray–box intersection tests, we’d like to never miss an actual intersection—getting this right is important for the correctness of the BVHAccel acceleration data structure in Section 4.3 so that valid ray–shape intersections aren’t missed. The ray–bounding box test introduced in Section 3.1.2 is based on computing a series of ray–slab intersections to find the parametric t Subscript normal m normal i normal n along the ray where the ray enters the bounding box and the t Subscript normal m normal a normal x where it exits. If t Subscript normal m normal i normal n Baseline less-than t Subscript normal m normal a normal x , the ray passes through the box; otherwise it misses it. With floating-point arithmetic, there may be error in the computed t values—if the computed t Subscript normal m normal i normal n value is greater than t Subscript normal m normal a normal x purely due to round-off error, the intersection test will incorrectly return a false result.

Recall that the computation to find the t value for a ray intersection with a plane perpendicular to the x axis at a point x is t equals left-parenthesis x minus normal o Subscript x Baseline right-parenthesis slash bold d Subscript x . Expressed as a floating-point computation and applying Equation (3.7), we have

t equals left-parenthesis x minus normal o Subscript x Baseline right-parenthesis circled-times left-parenthesis 1 circled-division-slash bold d Subscript x Baseline right-parenthesis subset-of StartFraction x minus normal o Subscript x Baseline Over bold d Subscript x Baseline EndFraction left-parenthesis 1 plus-or-minus epsilon right-parenthesis cubed comma

and so

t left-parenthesis 1 plus-or-minus gamma 3 right-parenthesis equals StartFraction x minus normal o Subscript x Baseline Over bold d Subscript x Baseline EndFraction period

The difference between the computed result t and the precise result is bounded by gamma 3 StartAbsoluteValue t EndAbsoluteValue .

If we consider the intervals around the computed t values that bound the fully precise value of t , then the case we’re concerned with is when the intervals overlap; if they don’t, then the comparison of computed values will give the correct result (Figure 3.42). If the intervals do overlap, it’s impossible to know the actual ordering of the t values. In this case, increasing t Subscript normal m normal a normal x by twice the error bound, 2 gamma 3 t Subscript normal m normal a normal x , before performing the comparison ensures that we conservatively return true in this case.

Figure 3.42: If the error bounds of the computed t Subscript normal m normal i normal n and t Subscript normal m normal a normal x values overlap, the comparison t Subscript normal m normal i normal n Baseline less-than t Subscript normal m normal a normal x may not actually indicate if a ray hit a bounding box. It’s better to conservatively return true in this case than to miss an actual intersection. Extending t Subscript normal m normal a normal x by twice its error bound ensures that the comparison is conservative.

We can now define the fragment for the ray–bounding box test in Section 3.1.2 that makes this adjustment.

<<Update tFar to ensure robust ray–bounds intersection>>= 
tFar *= 1 + 2 * gamma(3);

The fragments for the Bounds3::IntersectP() method, <<Update tMax and tyMax to ensure robust bounds intersection>> and <<Update tzMax to ensure robust bounds intersection>>, are similar and therefore not included here.

3.9.3 Robust Triangle Intersections

The details of the ray–triangle intersection algorithm in Section 3.6.2 were carefully designed to avoid cases where rays could incorrectly pass through an edge or vertex shared by two adjacent triangles without generating an intersection. Fittingly, an intersection algorithm with this guarantee is referred to as being watertight.

Recall that the algorithm is based on transforming triangle vertices into a coordinate system with the ray’s origin at its origin and the ray’s direction aligned along the plus z axis. Although round-off error may be introduced by transforming the vertex positions to this coordinate system, this error doesn’t affect the watertightness of the intersection test, since the same transformation is applied to all triangles. (Further, this error is quite small, so it doesn’t significantly impact the accuracy of the computed intersection points.)

Given vertices in this coordinate system, the three edge functions defined in Equation (3.1) are evaluated at the point left-parenthesis 0 comma 0 right-parenthesis ; the corresponding expressions, Equation (3.2), are quite straightforward. The key to the robustness of the algorithm is that with floating-point arithmetic, the edge function evaluations are guaranteed to have the correct sign. In general, we have

left-parenthesis a circled-times b right-parenthesis minus left-parenthesis c circled-times d right-parenthesis period

First, note that if a b equals c d , then Equation (3.12) evaluates to exactly zero, even in floating point. We therefore just need to show that if a b greater-than c d , then left-parenthesis a circled-times b right-parenthesis minus left-parenthesis c circled-times d right-parenthesis is never negative. If a b greater-than c d , then left-parenthesis a circled-times b right-parenthesis must be greater than or equal to left-parenthesis c circled-times d right-parenthesis . In turn, their difference must be greater than or equal to zero. (These properties both follow from the fact that floating-point arithmetic operations are all rounded to the nearest representable floating-point value.)

If the value of the edge function is zero, then it’s impossible to tell whether it is exactly zero or whether a small positive or negative value has rounded to zero. In this case, the fragment <<Fall back to double-precision test at triangle edges>> reevaluates the edge function with double precision; it can be shown that doubling the precision suffices to accurately distinguish these cases, given 32-bit floats as input.

The overhead caused by this additional precaution is minimal: in a benchmark with 88 million ray intersection tests, the double-precision fallback had to be used in less than 0.0000023% of the cases.

3.9.4 Bounding Intersection Point Error

We’ll now apply this machinery for analyzing rounding error to derive conservative bounds on the absolute error in computed ray-shape intersection points, which allows us to construct bounding boxes that are guaranteed to include an intersection point on the actual surface (Figure 3.43). These bounding boxes provide the basis of the algorithm for generating spawned ray origins that will be introduced in Section 3.9.5.

Figure 3.43: Shape intersection algorithms in pbrt compute an intersection point, shown here in the 2D setting with a filled circle. The absolute error in this point is bounded by delta Subscript x and delta Subscript y , giving a small box around the point. Because these bounds are conservative, we know that the actual intersection point on the surface (open circle) must lie somewhere within the box.

It’s useful to start by looking at the sources of error in conventional approaches to computing intersection points. It is common practice in ray tracing to compute 3D intersection points by first solving the parametric ray equation normal o plus t bold d for a value t Subscript normal h normal i normal t where a ray intersects a surface and then computing the hit point normal p Subscript with normal p Subscript Baseline equals normal o plus t Subscript normal h normal i normal t Baseline bold d . If t Subscript normal h normal i normal t carries some error delta Subscript t , then we can bound the error in the computed intersection point. Considering the x coordinate, for example, we have

StartLayout 1st Row 1st Column x 2nd Column equals normal o Subscript x Baseline circled-plus left-parenthesis t Subscript normal h normal i normal t Baseline plus-or-minus delta Subscript t Baseline right-parenthesis circled-times bold d Subscript x Baseline 2nd Row 1st Column Blank 2nd Column subset-of normal o Subscript x circled-plus left-parenthesis t Subscript normal h normal i normal t Baseline plus-or-minus delta Subscript t Baseline right-parenthesis bold d Subscript x Baseline left-parenthesis 1 plus-or-minus gamma 1 right-parenthesis 3rd Row 1st Column Blank 2nd Column subset-of normal o Subscript x Baseline left-parenthesis 1 plus-or-minus gamma 1 right-parenthesis plus left-parenthesis t Subscript normal h normal i normal t Baseline plus-or-minus delta Subscript t Baseline right-parenthesis bold d Subscript x Baseline left-parenthesis 1 plus-or-minus gamma 2 right-parenthesis 4th Row 1st Column Blank 2nd Column equals normal o Subscript x Baseline plus t Subscript normal h normal i normal t Baseline bold d Subscript x Baseline plus left-bracket plus-or-minus normal o Subscript x Baseline gamma 1 plus-or-minus delta Subscript t Baseline bold d Subscript x Baseline plus-or-minus t Subscript normal h normal i normal t Baseline bold d Subscript x Baseline gamma 2 plus-or-minus delta Subscript t Baseline bold d Subscript x Baseline gamma 2 right-bracket period EndLayout

The error term (in square brackets) is bounded by

gamma 1 StartAbsoluteValue normal o Subscript x Baseline EndAbsoluteValue plus delta Subscript t Baseline left-parenthesis 1 plus gamma 2 right-parenthesis StartAbsoluteValue bold d Subscript x Baseline EndAbsoluteValue plus gamma 2 StartAbsoluteValue t Subscript normal h normal i normal t Baseline bold d Subscript x Baseline EndAbsoluteValue period

There are two things to see from Equation (3.13): first, the magnitudes of the terms that contribute to the error in the computed intersection point ( normal o Subscript x , bold d Subscript x , and t Subscript normal h normal i normal t Baseline bold d Subscript x ) may be quite different from the magnitude of the intersection point. Thus, there is a danger of catastrophic cancellation in the intersection point’s computed value. Second, ray intersection algorithms generally perform tens of floating-point operations to compute t values, which in turn means that we can expect delta Subscript t to be at least of magnitude gamma Subscript n Baseline t , with n in the tens (and possibly much more, due to catastrophic cancellation). Each of these terms may be significant with respect to the magnitude of the computed point x .

Together, these factors can lead to relatively large error in the computed intersection point. We’ll develop better approaches shortly.

Reprojection: Quadrics

We’d like to reliably compute intersection points on surfaces with just a few ulps of error rather than the hundreds of ulps of error that intersection points computed with the parametric ray equation may have. Previously, Woo et al. (1996) suggested using the first intersection point computed as a starting point for a second ray–plane intersection, for ray–polygon intersections. From the bounds in Equation (3.13), we can see why the second intersection point will be much closer to the surface than the first: the t Subscript normal h normal i normal t value along the second ray will be quite close to zero, so that the magnitude of the absolute error in t Subscript normal h normal i normal t will be quite small, and thus using this value in the parametric ray equation will give a point quite close to the surface (Figure 3.44). Further, the ray origin will have similar magnitude to the intersection point, so the gamma 1 StartAbsoluteValue normal o Subscript x Baseline EndAbsoluteValue term won’t introduce much additional error.

Figure 3.44: Re-intersection to Improve the Accuracy of the Computed Intersection Point. Given a ray and a surface, an initial intersection point has been computed with the ray equation (filled circle). This point may be fairly inaccurate due to rounding error but can be used as the origin for a second ray–shape intersection. The intersection point computed from this second intersection (open circle) is much closer to the surface, though it may be shifted from the true intersection point due to error in the first computed intersection.

Although the second intersection point computed with this approach is much closer to the plane of the surface, it still suffers from error by being offset due to error in the first computed intersection. The farther away the ray origin from the intersection point (and thus, the larger the absolute error in t Subscript normal h normal i normal t ), the larger this error will be. In spite of this error, the approach has merit: we’re generally better off with a computed intersection point that is quite close to the actual surface, even if offset from the most accurate possible intersection point, than we are with a point that is some distance above or below the surface (and likely also far from the most accurate intersection point).

Rather than doing a full re-intersection computation, which may not only be computationally costly but also will still have error in the computed t value, an effective approach is to refine computed intersection points by reprojecting them to the surface. The error bounds for these reprojected points are often remarkably small.

It should be noted that these reprojection error bounds don’t capture tangential errors that were present in the original intersection normal p Subscript —the main focus here is to detect errors that might cause the reprojected point normal p prime to fall below the surface.

Consider a ray–sphere intersection: given a computed intersection point (e.g., from the ray equation) normal p Subscript with a sphere at the origin with radius r , we can reproject the point onto the surface of the sphere by scaling it with the ratio of the sphere’s radius to the computed point’s distance to the origin, computing a new point normal p prime equals left-parenthesis x prime comma y prime comma z prime right-parenthesis with

x prime equals x StartFraction r Over StartRoot x squared plus y squared plus z squared EndRoot EndFraction comma

and so forth. The floating-point computation is

StartLayout 1st Row 1st Column x prime 2nd Column equals x circled-times r circled-division-slash monospace s monospace q monospace r monospace t monospace left-parenthesis monospace left-parenthesis monospace x circled-times monospace x monospace right-parenthesis circled-plus monospace left-parenthesis monospace y circled-times monospace y monospace right-parenthesis circled-plus monospace left-parenthesis monospace z circled-times monospace z monospace right-parenthesis monospace right-parenthesis 2nd Row 1st Column Blank 2nd Column subset-of StartFraction x r left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis squared Over StartRoot x squared left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis cubed plus y squared left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis cubed plus z squared left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis squared EndRoot left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis EndFraction 3rd Row 1st Column Blank 2nd Column subset-of StartFraction x r left-parenthesis 1 plus-or-minus gamma 2 right-parenthesis Over StartRoot x squared left-parenthesis 1 plus-or-minus gamma 3 right-parenthesis plus y squared left-parenthesis 1 plus-or-minus gamma 3 right-parenthesis plus z squared left-parenthesis 1 plus-or-minus gamma 2 right-parenthesis EndRoot left-parenthesis 1 plus-or-minus gamma 1 right-parenthesis EndFraction EndLayout

Because x squared , y squared , and z squared are all positive, the terms in the square root can share the same gamma term, and we have

StartLayout 1st Row 1st Column x prime 2nd Column subset-of StartFraction x r left-parenthesis 1 plus-or-minus gamma 2 right-parenthesis Over StartRoot left-parenthesis x squared plus y squared plus z squared right-parenthesis left-parenthesis 1 plus-or-minus gamma 4 right-parenthesis EndRoot left-parenthesis 1 plus-or-minus gamma 1 right-parenthesis EndFraction 2nd Row 1st Column Blank 2nd Column equals StartFraction x r left-parenthesis 1 plus-or-minus gamma 2 right-parenthesis Over StartRoot left-parenthesis x squared plus y squared plus z squared right-parenthesis EndRoot StartRoot left-parenthesis 1 plus-or-minus gamma 4 right-parenthesis EndRoot left-parenthesis 1 plus-or-minus gamma 1 right-parenthesis EndFraction 3rd Row 1st Column Blank 2nd Column subset-of StartFraction x r Over StartRoot left-parenthesis x squared plus y squared plus z squared right-parenthesis EndRoot EndFraction left-parenthesis 1 plus-or-minus gamma 5 right-parenthesis 4th Row 1st Column Blank 2nd Column equals x prime left-parenthesis 1 plus-or-minus gamma 5 right-parenthesis period EndLayout

Thus, the absolute error of the reprojected x coordinate is bounded by gamma 5 StartAbsoluteValue x prime EndAbsoluteValue (and similarly for y prime and z prime ) and is thus no more than 2.5 ulps in each dimension from a point on the surface of the sphere.

Here is the fragment that reprojects the intersection point for the Sphere shape.

<<Refine sphere intersection point>>= 
pHit *= radius / Distance(pHit, Point3f(0, 0, 0));

The error bounds follow from Equation (3.14).

<<Compute error bounds for sphere intersection>>= 
Vector3f pError = gamma(5) * Abs((Vector3f)pHit);

Reprojection algorithms and error bounds for other quadrics can be defined similarly: for example, for a cylinder along the z axis, only the x and y coordinates need to be reprojected, and the error bounds in x and y turn out to be only gamma 3 times their magnitudes.

<<Refine cylinder intersection point>>= 
Float hitRad = std::sqrt(pHit.x * pHit.x + pHit.y * pHit.y); pHit.x *= radius / hitRad; pHit.y *= radius / hitRad;

<<Compute error bounds for cylinder intersection>>= 
Vector3f pError = gamma(3) * Abs(Vector3f(pHit.x, pHit.y, 0));

The disk shape is particularly easy; we just need to set the z coordinate of the point to lie on the plane of the disk.

<<Refine disk intersection point>>= 
pHit.z = height;

In turn, we have a point with zero error; it lies exactly on the surface on the disk.

<<Compute error bounds for disk intersection>>= 
Vector3f pError(0, 0, 0);

Parametric Evaluation: Triangles

Another effective approach to computing precise intersection points is to use the parametric representation of a shape to compute accurate intersection points. For example, the triangle intersection algorithm in Section 3.6.2 computes three edge function values e 0 , e 1 , and e 2 and reports an intersection if all three have the same sign. Their values can be used to find the barycentric coordinates

b Subscript i Baseline equals StartFraction e Subscript i Baseline Over e 0 plus e 1 plus e 2 EndFraction period

Attributes v Subscript i at the triangle vertices (including the vertex positions) can be interpolated across the face of the triangle by

v prime equals b 0 v 0 plus b 1 v 1 plus b 2 v 2 period

We can show that interpolating the positions of the vertices in this manner gives a point very close to the surface of the triangle. First consider precomputing the inverse sum of e Subscript i :

StartLayout 1st Row 1st Column d 2nd Column equals 1 circled-division-slash left-parenthesis e 0 circled-plus e 1 circled-plus e 2 right-parenthesis 2nd Row 1st Column Blank 2nd Column subset-of StartFraction 1 Over left-parenthesis e 0 plus e 1 right-parenthesis left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis squared plus e 2 left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis EndFraction left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis period EndLayout

Because all e Subscript i have the same sign if there is an intersection, we can collect the e Subscript i terms and conservatively bound d :

StartLayout 1st Row 1st Column d 2nd Column subset-of StartFraction 1 Over left-parenthesis e 0 plus e 1 plus e 2 right-parenthesis left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis squared EndFraction left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column subset-of StartFraction 1 Over e 0 plus e 1 plus e 2 EndFraction left-parenthesis 1 plus-or-minus gamma 3 right-parenthesis period EndLayout

If we now consider interpolation of the x coordinate of the position in the triangle corresponding to the edge function values, we have

StartLayout 1st Row 1st Column x prime 2nd Column equals left-parenthesis left-parenthesis e 0 circled-times x 0 right-parenthesis circled-plus left-parenthesis e 1 circled-times x 1 right-parenthesis circled-plus left-parenthesis e 2 circled-times x 2 right-parenthesis right-parenthesis circled-times d 2nd Row 1st Column Blank 2nd Column subset-of left-parenthesis e 0 x 0 left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis cubed plus e 1 x 1 left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis cubed plus e 2 x 2 left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis squared right-parenthesis d left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis 3rd Row 1st Column Blank 2nd Column subset-of left-parenthesis e 0 x 0 left-parenthesis 1 plus-or-minus gamma 4 right-parenthesis plus e 1 x 1 left-parenthesis 1 plus-or-minus gamma 4 right-parenthesis plus e 2 x 2 left-parenthesis 1 plus-or-minus gamma 3 right-parenthesis right-parenthesis d period EndLayout

Using the bounds on d ,

StartLayout 1st Row 1st Column x 2nd Column subset-of StartFraction e 0 x 0 left-parenthesis 1 plus-or-minus gamma 7 right-parenthesis plus e 1 x 1 left-parenthesis 1 plus-or-minus gamma 7 right-parenthesis plus e 2 x 2 left-parenthesis 1 plus-or-minus gamma 6 right-parenthesis Over e 0 plus e 1 plus e 2 EndFraction 2nd Row 1st Column Blank 2nd Column equals b 0 x 0 left-parenthesis 1 plus-or-minus gamma 7 right-parenthesis plus b 1 x 1 left-parenthesis 1 plus-or-minus gamma 7 right-parenthesis plus b 2 x 2 left-parenthesis 1 plus-or-minus gamma 6 right-parenthesis period EndLayout

Thus, we can finally see that the absolute error in the computed x prime value is in the interval

plus-or-minus b 0 x 0 gamma 7 plus-or-minus b 1 x 1 gamma 7 plus-or-minus b 2 x 2 gamma 7 comma

which is bounded by

gamma 7 left-parenthesis StartAbsoluteValue b 0 x 0 EndAbsoluteValue plus StartAbsoluteValue b 1 x 1 EndAbsoluteValue plus StartAbsoluteValue b 2 x 2 EndAbsoluteValue right-parenthesis period

(Note that the b 2 x 2 term could have a gamma 6 factor instead of gamma 7 , but the difference between the two is very small so we choose a slightly simpler final expression.) Equivalent bounds hold for y prime and z prime .

Equation (3.15) lets us bound the error in the interpolated point computed in Triangle::Intersect().

<<Compute error bounds for triangle intersection>>= 
Float xAbsSum = (std::abs(b0 * p0.x) + std::abs(b1 * p1.x) + std::abs(b2 * p2.x)); Float yAbsSum = (std::abs(b0 * p0.y) + std::abs(b1 * p1.y) + std::abs(b2 * p2.y)); Float zAbsSum = (std::abs(b0 * p0.z) + std::abs(b1 * p1.z) + std::abs(b2 * p2.z)); Vector3f pError = gamma(7) * Vector3f(xAbsSum, yAbsSum, zAbsSum);

Other Shapes

For shapes where we may not want to derive reprojection methods and tight error bounds, running error analysis can be quite useful: we implement all of the intersection calculations using EFloat instead of Float, compute a t Subscript normal h normal i normal t value, and use the parametric ray equation to compute a hit point. We can then find conservative bounds on the error in the computed intersection point via the EFloat GetAbsoluteError() method.

<<Compute error bounds for intersection computed with ray equation>>= 
EFloat px = ox + tShapeHit * dx; EFloat py = oy + tShapeHit * dy; EFloat pz = oz + tShapeHit * dz; Vector3f pError = Vector3f(px.GetAbsoluteError(), py.GetAbsoluteError(), pz.GetAbsoluteError());

This approach is used for cones, paraboloids, and hyperboloids in pbrt.

<<Compute error bounds for cone intersection>>= 
<<Compute error bounds for intersection computed with ray equation>> 
EFloat px = ox + tShapeHit * dx; EFloat py = oy + tShapeHit * dy; EFloat pz = oz + tShapeHit * dz; Vector3f pError = Vector3f(px.GetAbsoluteError(), py.GetAbsoluteError(), pz.GetAbsoluteError());

Because the Curve shape orients itself to face incident rays, rays leaving it must be offset by twice the curve’s width in order to not incorrectly re-intersect it when it’s reoriented to face them.

<<Compute error bounds for curve intersection>>= 
Vector3f pError(2 * hitWidth, 2 * hitWidth, 2 * hitWidth);

Effect of Transformations

The last detail to attend to in order to bound the error in computed intersection points is the effect of transformations, which introduce additional rounding error when they are applied to computed intersection points.

The quadric Shapes in pbrt transform world space rays into object space before performing ray–shape intersections and then transform computed intersection points back to world space. Both of these transformation steps introduce rounding error that needs to be accounted for in order to maintain robust world space bounds around intersection points.

If possible, it’s best to try to avoid coordinate-system transformations of rays and intersection points. For example, it’s better to transform triangle vertices to world space and intersect world space rays with them than to transform rays to object space and then transform intersection points to world space. Transformations are still useful—for example, for the quadrics and for object instancing, so we’ll show how to bound the error that they introduce.

We’ll start by considering the error introduced by transforming a point left-parenthesis x comma y comma z right-parenthesis that is exact—i.e., without any accumulated error. Given a 4 times 4 non-projective transformation matrix with elements denoted by m Subscript i comma j , the transformed coordinate x prime is

StartLayout 1st Row 1st Column x prime 2nd Column equals left-parenthesis left-parenthesis m Subscript 0 comma 0 Baseline circled-times x right-parenthesis circled-plus left-parenthesis m Subscript 0 comma 1 Baseline circled-times y right-parenthesis right-parenthesis circled-plus left-parenthesis left-parenthesis m Subscript 0 comma 2 Baseline circled-times z right-parenthesis circled-plus m Subscript 0 comma 3 Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column subset-of m Subscript 0 comma 0 Baseline x left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis cubed plus m Subscript 0 comma 1 Baseline y left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis cubed plus m Subscript 0 comma 2 Baseline z left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis cubed plus m Subscript 0 comma 3 Baseline left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis squared 3rd Row 1st Column Blank 2nd Column subset-of left-parenthesis m Subscript 0 comma 0 Baseline x plus m Subscript 0 comma 1 Baseline y plus m Subscript 0 comma 2 Baseline z plus m Subscript 0 comma 3 Baseline right-parenthesis plus gamma 3 left-parenthesis plus-or-minus m Subscript 0 comma 0 Baseline x plus-or-minus m Subscript 0 comma 1 Baseline y plus-or-minus m Subscript 0 comma 2 Baseline z plus-or-minus m Subscript 0 comma 3 Baseline right-parenthesis 4th Row 1st Column Blank 2nd Column subset-of left-parenthesis m Subscript 0 comma 0 Baseline x plus m Subscript 0 comma 1 Baseline y plus m Subscript 0 comma 2 Baseline z plus m Subscript 0 comma 3 Baseline right-parenthesis plus-or-minus gamma 3 left-parenthesis StartAbsoluteValue m Subscript 0 comma 0 Baseline x EndAbsoluteValue plus StartAbsoluteValue m Subscript 0 comma 1 Baseline y EndAbsoluteValue plus StartAbsoluteValue m Subscript 0 comma 2 Baseline z EndAbsoluteValue plus StartAbsoluteValue m Subscript 0 comma 3 Baseline EndAbsoluteValue right-parenthesis period EndLayout

Thus, the absolute error in the result is bounded by

gamma 3 left-parenthesis StartAbsoluteValue m Subscript 0 comma 0 Baseline x EndAbsoluteValue plus StartAbsoluteValue m Subscript 0 comma 1 Baseline y EndAbsoluteValue plus StartAbsoluteValue m Subscript 0 comma 2 Baseline z EndAbsoluteValue plus StartAbsoluteValue m Subscript 0 comma 3 Baseline EndAbsoluteValue right-parenthesis period

Similar bounds follow for the transformed y prime and z prime coordinates.

We’ll use this result to add a method to the Transform class that also returns the absolute error in the transformed point due to applying the transformation.

<<Transform Inline Functions>>+= 
template <typename T> inline Point3<T> Transform::operator()(const Point3<T> &p, Vector3<T> *pError) const { T x = p.x, y = p.y, z = p.z; <<Compute transformed coordinates from point pt>> 
T xp = m.m[0][0] * x + m.m[0][1] * y + m.m[0][2] * z + m.m[0][3]; T yp = m.m[1][0] * x + m.m[1][1] * y + m.m[1][2] * z + m.m[1][3]; T zp = m.m[2][0] * x + m.m[2][1] * y + m.m[2][2] * z + m.m[2][3]; T wp = m.m[3][0] * x + m.m[3][1] * y + m.m[3][2] * z + m.m[3][3];
<<Compute absolute error for transformed point>> 
T xAbsSum = (std::abs(m.m[0][0] * x) + std::abs(m.m[0][1] * y) + std::abs(m.m[0][2] * z) + std::abs(m.m[0][3])); T yAbsSum = (std::abs(m.m[1][0] * x) + std::abs(m.m[1][1] * y) + std::abs(m.m[1][2] * z) + std::abs(m.m[1][3])); T zAbsSum = (std::abs(m.m[2][0] * x) + std::abs(m.m[2][1] * y) + std::abs(m.m[2][2] * z) + std::abs(m.m[2][3])); *pError = gamma(3) * Vector3<T>(xAbsSum, yAbsSum, zAbsSum);
if (wp == 1) return Point3<T>(xp, yp, zp); else return Point3<T>(xp, yp, zp) / wp; }

The fragment <<Compute transformed coordinates from point pt>> isn’t included here; it implements the same matrix/point multiplication as in Section 2.8.

Note that the code that computes error bounds is buggy if the matrix is projective and the homogeneous w coordinate of the projected point is not one; this nit currently isn’t a problem for pbrt’s usage of this method.

<<Compute absolute error for transformed point>>= 
T xAbsSum = (std::abs(m.m[0][0] * x) + std::abs(m.m[0][1] * y) + std::abs(m.m[0][2] * z) + std::abs(m.m[0][3])); T yAbsSum = (std::abs(m.m[1][0] * x) + std::abs(m.m[1][1] * y) + std::abs(m.m[1][2] * z) + std::abs(m.m[1][3])); T zAbsSum = (std::abs(m.m[2][0] * x) + std::abs(m.m[2][1] * y) + std::abs(m.m[2][2] * z) + std::abs(m.m[2][3])); *pError = gamma(3) * Vector3<T>(xAbsSum, yAbsSum, zAbsSum);

The result in Equation (3.16) assumes that the point being transformed is exact. If the point itself has error bounded by delta Subscript x , delta Subscript y , and delta Subscript z , then the transformed x coordinate is given by:

x prime equals left-parenthesis m Subscript 0 comma 0 Baseline circled-times left-parenthesis x plus-or-minus delta Subscript x Baseline right-parenthesis circled-plus m Subscript 0 comma 1 Baseline circled-times left-parenthesis y plus-or-minus delta Subscript y Baseline right-parenthesis right-parenthesis circled-plus left-parenthesis m Subscript 0 comma 2 Baseline circled-times left-parenthesis z plus-or-minus delta Subscript z Baseline right-parenthesis circled-plus m Subscript 0 comma 3 Baseline right-parenthesis period

Applying the definition of floating-point addition and multiplication’s error bounds, we have:

StartLayout 1st Row 1st Column x prime 2nd Column equals m Subscript 0 comma 0 Baseline left-parenthesis x plus-or-minus delta Subscript x Baseline right-parenthesis left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis cubed plus m Subscript 0 comma 1 Baseline left-parenthesis y plus-or-minus delta Subscript y Baseline right-parenthesis left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis cubed plus 2nd Row 1st Column Blank 2nd Column m Subscript 0 comma 2 Baseline left-parenthesis z plus-or-minus delta Subscript z Baseline right-parenthesis left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis cubed plus m Subscript 0 comma 3 Baseline left-parenthesis 1 plus-or-minus epsilon Subscript normal m Baseline right-parenthesis squared period EndLayout

Transforming to use gamma , we can find the absolute error term to be bounded by

StartLayout 1st Row 1st Column left-parenthesis gamma 3 plus 1 right-parenthesis left-parenthesis StartAbsoluteValue m Subscript 0 comma 0 Baseline EndAbsoluteValue delta Subscript x Baseline 2nd Column plus StartAbsoluteValue m Subscript 0 comma 1 Baseline EndAbsoluteValue delta Subscript y Baseline plus StartAbsoluteValue m Subscript 0 comma 2 Baseline EndAbsoluteValue delta Subscript z Baseline right-parenthesis plus 2nd Row 1st Column Blank 2nd Column gamma 3 left-parenthesis StartAbsoluteValue m Subscript 0 comma 0 Baseline x EndAbsoluteValue plus StartAbsoluteValue m Subscript 0 comma 1 Baseline y EndAbsoluteValue plus StartAbsoluteValue m Subscript 0 comma 2 Baseline z EndAbsoluteValue plus StartAbsoluteValue m Subscript 0 comma 3 Baseline EndAbsoluteValue right-parenthesis period EndLayout