## 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 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.)

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 closer than some value. Figure 3.40 shows why this approach requires fairly high 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 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.

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 with equal spacing of 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 ranges from 0 to 255; the actual exponent used, , is computed by biasing :

The significand actually has 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 , significand , and exponent , the corresponding floating-point value is

For example, with a normalized significand, the floating-point number 6.5 is written as , where the 2 subscript denotes a base-2 value. (If non-whole binary numbers aren’t immediately intuitive, note that the first number to the right of the radix point contributes , and so forth.) Thus, we have

, so and .

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

Similarly the floating-point value has and , so and its binary representation is:

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 , the spacing is

Thus, for floating-point numbers between 1 and 2, , and the
spacing between floating-point values is . 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 , or , 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

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

Without the leading one bit, the minimum value is

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, , is also reserved for special treatment. Therefore, the largest regular floating-point value that can be represented has (or ) and is approximately

With , 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 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.

With , 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 .
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:

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:

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:

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.

`NextFloatUp()`>>

`v`to next higher float>>

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 and aren’t adjacent.

`NextFloatUp()`>>=

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.

`v`to next higher float>>=

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

#### 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 for floating-point square root, these precision guarantees can be written as:

where 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

for some . The amount of error introduced from this rounding can be no more than half the floating-point spacing at —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).

For 32-bit floats, we can bound the floating-point spacing at from
above using Equation (3.6) (i.e., an ulp at that value) by
, so half the spacing is bounded from above by
and so . This bound is the *machine
epsilon*.
For 32-bit floats, .

Thus, we have

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 ,

- .
- .
- .
- .
- and 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).
- for all integer , assuming 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 , we say that the magnitude of the difference
between and the result of doing that computation in the real
numbers is the *absolute error*, :

*Relative error*, , is the ratio of the absolute
error to the precise result:

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

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

Because is small, higher order powers of can be bounded by an additional term, and so we can bound the terms with

(As a practical matter, almost bounds these terms, since higher powers of get very small very quickly, but the above is a fully conservative bound.)

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

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

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 , , , and .

This is a fairly interesting result; we see that the magnitude of makes a relatively large contribution to the error bound, especially compared to . (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

If we apply the same process of applying Equation (3.8), expanding out terms, converting higher-order terms to , we get absolute error bounds of

which are lower than the first formulation if is relatively large, but possibly higher if 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 by is somewhat unsatisfying since it adds a whole term purely to conservatively bound the sum of various higher powers of . Higham (2002, Section 3.1) gives an approach to more tightly bound products of error terms. If we have , it can be shown that this value is bounded by , where

as long as (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 values, so it just barely increases to achieve a conservative bound.

We will denote this bound by :

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.

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

An advantage of this approach is that quotients of terms can also be bounded with the function. Given

the interval is bounded by . Thus, can be used to collect 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 terms represent intervals, canceling them would be incorrect:

the 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, and that each carry accumulated error from earlier operations, consider their product. Using the definition of , the result is in the interval:

where we’ve used the relationship , which follows directly from Equation (3.11).

The relative error in this result is bounded by:

and so the final error is thus just roughly 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

which is in the interval and so the absolute error is bounded by .

If the signs of and are the same, then the absolute error is bounded by and the relative error is around ulps around the computed value.

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

The numerator’s magnitude is proportional to the original value 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
, though it’s now in relation to a value much smaller than .

#### 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`>>

`EFloat` maintains a computed value `v` and the absolute
error bound, `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.

`EFloat`>>=

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

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

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.

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.

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.

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

`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`.

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 along the ray where
the ray enters the bounding box and the where it exits. If , the ray passes through the box; otherwise it misses it. With
floating-point arithmetic, there may be error in the computed
values—if the computed value is greater than purely due
to round-off error, the intersection test will incorrectly return a false
result.

Recall that the computation to find the value for a ray intersection with a plane perpendicular to the axis at a point is . Expressed as a floating-point computation and applying Equation (3.7), we have

and so

The difference between the computed result and the precise result is bounded by .

If we consider the intervals around the computed values that bound the fully precise value of , 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 values. In this case, increasing by twice the error bound, , before performing the comparison ensures that we conservatively return true in this case.

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

`tFar`to ensure robust ray–bounds intersection>>=

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

First, note that if , then Equation (3.12) evaluates to exactly zero, even in floating point. We therefore just need to show that if , then is never negative. If , then must be greater than or equal to . 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.

`pbrt`compute an intersection point, shown here in the 2D setting with a filled circle. The absolute error in this point is bounded by and , 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 for a value where a ray intersects a surface and then computing the hit point with . If carries some error , then we can bound the error in the computed intersection point. Considering the coordinate, for example, we have

The error term (in square brackets) is bounded by

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 (, , and ) 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 values, which in turn means that we can expect to be at least of magnitude , with 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 .

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 value along the second ray will be quite close to zero, so that the magnitude of the absolute error in 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 term won’t introduce much additional error.

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 ), 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 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 —the main focus here is to detect errors that might cause the reprojected point to fall below the surface.

Consider a ray–sphere intersection: given a computed intersection point (e.g., from the ray equation) with a sphere at the origin with radius , 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 with

and so forth. The floating-point computation is

Because , , and are all positive, the terms in the square root can share the same term, and we have

Thus, the absolute error of the reprojected coordinate is bounded by (and similarly for and ) 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.

The error bounds follow from Equation (3.14).

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

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

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

#### 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 , , and and reports an intersection if all three have the same sign. Their values can be used to find the barycentric coordinates

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

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 reciprocal of the sum of :

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

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

Using the bounds on ,

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

which is bounded by

(Note that the term could have a factor instead of , but the difference between the two is very small so we choose a slightly simpler final expression.) Equivalent bounds hold for and .

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

#### 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 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.

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

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.

#### 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 `Shape`s 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 that is exact—i.e., without any accumulated error. Given a non-projective transformation matrix with elements denoted by , the transformed coordinate is

Thus, the absolute error in the result is bounded by

Similar bounds follow for the transformed and 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.

`pt`>>

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 coordinate of the projected point is not
one; this nit currently isn’t a problem for `pbrt`’s usage of this method.

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

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

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