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.

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.

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.

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

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.

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.

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

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

The Transform class also provides an operator() that takes a point and its own absolute error and returns the absolute error in the result, applying Equation (3.17). The definition is straightforward, so isn’t included in the text here.

The Transform class also provides methods to transform vectors and rays, returning the resulting error. The vector error bound derivations (and thence, implementations) are very similar to those for points, and so also aren’t included here.

This method is used to transform the intersection point and its error bounds in the Transform::operator() method for SurfaceInteractions.

3.9.5 Robust Spawned Ray Origins

Computed intersection points and their error bounds give us a small 3D box that bounds a region of space. We know that the precise intersection point must be somewhere inside this box and that thus the surface must pass through the box (at least enough to present the point where the intersection is). (Recall Figure 3.43.) Having these boxes makes it possible to position the origins of rays leaving the surface so that they are always on the right side of the surface so that they don’t incorrectly reintersect it. When tracing spawned rays leaving the intersection point , we offset their origins enough to ensure that they are past the boundary of the error box and thus won’t incorrectly re-intersect the surface.

In order to ensure that the spawned ray origin is definitely on the right side of the surface, we move far enough along the normal so that the plane perpendicular to the normal is outside the error bounding box. To see how to do this, consider a computed intersection point at the origin, where the plane equation for the plane going through the intersection point is just

the plane is implicitly defined by , and the normal is .

For a point not on the plane, the value of the plane equation gives the offset along the normal that gives a plane that goes through the point. We’d like to find the maximum value of for the eight corners of the error bounding box; if we offset the plane plus and minus this offset, we have two planes that don’t intersect the error box that should be (locally) on opposite sides of the surface, at least at the computed intersection point offset along the normal (Figure 3.45).

If the eight corners of the error bounding box are given by , then the maximum value of is easily computed:

Computing spawned ray origins by offsetting along the surface normal in this way has a few advantages: assuming that the surface is locally planar (a reasonable assumption, especially at the very small scale of the intersection point error bounds), moving along the normal allows us to get from one side of the surface to the other while moving the shortest distance. In general, minimizing the distance that ray origins are offset is desirable for maintaining shadow and reflection detail.

We also must handle round-off error when computing the offset point: when offset is added to p, the result will in general need to be rounded to the nearest floating-point value. In turn, it may be rounded down toward p such that the resulting point is in the interior of the error box rather than in its boundary (Figure 3.46). Therefore, the offset point is rounded away from p here to ensure that it’s not inside the box. ^†

Alternatively, the floating-point rounding mode could have been set to round toward plus or minus infinity (based on the sign of the value). Changing the rounding mode is generally fairly expensive, so we just shift the floating-point value by one ulp here. This will sometimes cause a value already outside of the error box to go slightly farther outside it, but because the floating-point spacing is so small, this isn’t a problem in practice.

Given the OffsetRayOrigin() function, we can now implement the Interaction methods that generate rays leaving intersection points.