8.6 Fourier Basis BSDFs

While reflection models like Torrance–Sparrow and Oren–Nayar can accurately represent many materials, some materials have BRDF shapes that don’t match these models well. (Examples include layered materials like metals with smooth or rough coatings or fabrics, which are often partially retro-reflective.) One option for materials like these is to store their BSDF values in a large 3D or 4D lookup table, though this approach can require an unacceptable amount of storage—for example, if and are sampled in spherical angles with 1-degree spacing, then over one billion sample points are needed to represent the corresponding anisotropic BSDF in the form of a 4D lookup table.

Therefore, having a more compact representation that still represents the BSDF accurately is highly desirable. This section introduces the FourierBSDF, which represents BSDFs with sums of scaled cosine terms using the Fourier basis. This representation is accurate, space-efficient, and works well with Monte Carlo integration (see Section 14.1.4.) Figure 8.24 shows two instances of the dragon model rendered using this representation.

Here, we won’t discuss how BSDFs are transformed into this representation, but we will focus on its use in rendering. See the “Further Reading” section at the end of this chapter for pointers to more details on those issues and the brdfs/ directory in the pbrt distribution for a variety of BSDFs represented in this format.

The FourierBSDF represents isotropic BSDFs by parameterizing the BSDF by a pair of spherical coordinates for the incident and outgoing directions, where and denote the cosines of the incident and outgoing zenith angles, respectively, and and are the azimuth angles:

The assumption of isotropy means that the function can be rewritten with a simpler dependence on the zenith angle cosines and the azimuth difference angle :

Isotropic BSDFs are generally also even functions of the azimuth difference, i.e.:

Given these properties, the product of the BSDF with the cosine falloff is then represented as a Fourier series in the azimuth angle difference.

Note how only cosine terms and no sine terms are needed due to Equation (8.20). The function evaluations denote the Fourier coefficients for a specific pair of zenith angle cosines.

Next, the functions are discretized over their input arguments. We choose a set of zenith angle cosines and store the values of for every pair . Thus, we can think of each as a matrix, and the entire BRDF representation then consists of a set of such matrices. Each describes a different azimuthal oscillation frequency in the material’s response to incident illumination.

The maximum order needed to evaluate Equation (8.21) to satisfactory accuracy varies: it depends on the particular zenith angles, so it’s worth adapting the number of coefficients to the complexity of the BSDF for a given pair of directions. Doing this is very important for the compactness of this representation.

To see the value of being able to vary the number of coefficients, consider nearly perfect specular reflection: when , many coefficients are necessary to accurately represent the specular lobe, which is zero for almost all azimuth angle differences and then very large for a small set of directions around , where the incident and outgoing directions are nearly opposite. However, when and aren’t aligned, only a single term is needed to represent that the BSDF is zero (or has negligible value). ^† For smoother BSDFs, most or all pairs of and angles require multiple coefficients to represent their distribution accurately, but their smoothness means that not too many coefficients are generally needed for each . The FourierBSDF representation exploits this property and only stores the sparse set of coefficients that is needed to achieve a desired accuracy. Thus, for most types of realistic BSDF data, the representation of Equation (8.21) is fairly compact; a few megabytes is typical.

FourierBSDFTable is a helper structure that holds all of the data for a BSDF represented in this manner. It’s mostly a simple struct that collects data that’s directly accessible to calling code, though it does provide a few utility methods.

The Read() method initializes the structure for the BSDF stored in the given file. It returns true on success or false if an error was encountered while reading the file.

If the BSDF represents scattering at the boundary between two different media, then the FourierBSDFTable::eta member variable gives the relative index of refraction over the surface boundary (Section 8.2.3), mMax gives the maximum order for any pair of , directions; this upper limit is useful when allocating temporary buffers to store coefficients, for example.

nChannels gives the number of spectral channels available; in this implementation, it is either 1, representing a monochromatic BSDF, or 3, representing a BSDF with RGB colors. Here, the three-channel variant actually stores luminance, red, and blue, rather than red, green, and blue—representing luminance directly turns out to be useful for the Monte Carlo sampling routines defined in Section 14.1.4, since it provides aggregate information about the function over all color channels. The corresponding green color is easily computed from luminance, red, and blue, as we’ll see shortly.

The zenith angles are discretized into nMu directions, which are stored in the mu array. mu is sorted from low to high, so that binary search can be used to find the entry that’s closest to a given or angle.

To evaluate Equation (8.21), we need to know the target Fourier order and all coefficients corresponding to the directions and . For simplicity now, we’ll present the basic ideas as if only the coefficients for the closest mu directions less than or equal to and are used, though the implementation to follow interpolates between coefficients from multiple mu values around the directions.

The order of the Fourier representation is always bounded by mMax but varies with respect to the incident and outgoing zenith angle cosine and : how many orders are needed can be determined by querying an integer matrix m.

To find for a particular set of angles, we first perform two binary searches in the discretized directions mu to give the offsets oi and oo such that

Using these indices, the requisite order can be fetched from m[oo * nMu + oi].

All of the coefficients for all of the pairs of discretized directions mu are packed into the a array. Because the maximum order (and thus, number of coefficients) varies and can even be zero depending on the characteristics of the BSDF for a given pair of directions, finding the offset into the a array is a two-step process:

1. First, the offsets oi and oo are used to index into the aOffset array to get an offset into a: offset = aOffset[oo * nMu + oi]. (The aOffset array thus has a total of nMu * nMu entries.)
2. Next, the coefficients starting at a[offset] give the values for the corresponding pair of directions. For the three color channel case, the first coefficients after a[offset] encode coefficients for luminance, the next correspond to the red channel, and then blue follows.

GetAk() is a small convenience method that, given offsets into the mu array for the incident and outgoing direction cosines, returns the order of coefficients for the directions and a pointer to their coefficients.

The FourierBSDF class provides a bridge between the FourierBSDFTable representation and the BxDF interface. Instances of this class are created by the FourierMaterial class.

The FourierBSDF class stores only a const reference to the table; the table is large enough that we definitely don’t want to make a separate copy of it for each FourierBSDF instance. Only read-access is needed here, so this approach doesn’t cause any problems. (FourierMaterial is responsible for the FourierBSDFTable storage.)

Evaluating the BSDF is a matter of computing the cosines and , finding the corresponding coefficients and maximum order, and then evaluating Equation (8.21).

There is an important difference of convention in how directions are represented within the FourierBSDF: the incident direction is negated compared to the usual approach in pbrt. This difference is helpful when performing other computations such as computing BSDFs for layered materials using this representation. ^†

So that the reconstructed BSDF is fairly smooth, the implementation here interpolates coefficients over the product of the four quantized mu directions that surround and the four that surround . The interpolation is performed with a tensor-product spline, where weights for the sampled function values are computed separately for each parameter and then multiplied together. Each final Fourier coefficient is thus computed by

where gives the th Fourier coefficient for the quantized directions and and are the spline weights. This interpolation ensures sufficient smoothness even when the discretization of directions is relatively coarse; the details of how these weights are computed will be explained in a few pages.

For each direction and , GetWeightsAndOffset() gives the offset of the first of the four mu values to be interpolated over and an array of four floating-point weights.

The various vectors in the extent of the directions being interpolated over may have different orders . Therefore, the implementation here allocates storage for the values using the maximum possible order times the number of channels for the size. For the multiple-channel case, the first bsdfTable.mMax entries of the ak array allocated here will be used for the first channel, the next mMax are for the second channel, and so forth. (Thus there is generally some unused space in the ak array for the usual case that the maximum order over the sixteen directions is less than mMax.) All of this storage is initialized to zero, so that subsequent code can add terms of Equation (8.22) to the corresponding entry in ak directly.

Given weights, offsets, and storage for the result, the interpolation can now be performed.

Given the weights and the starting offsets, adding each term of the sum in Equation (8.22) is a matter of getting the corresponding coefficients from the table for the offset and adding them to the running sum in ak.

Given the final weighted coefficients in ak, a call to Fourier() computes the BSDF value for the first color channel. Error in Fourier reconstruction can manifest itself as negative values, so the returned value must be clamped to zero.

Recall from Equation (8.21) that the coefficients represent the cosine-weighted BSDF. This cosine factor must be removed from the value returned from the f() method; the scale term encodes this factor.

As with specular transmission, radiance is scaled as it passes from a medium with one index of refraction to another, but this scaling isn’t applied to rays starting from the camera. A definition and discussion of the fragment <<Update scale to account for adjoint light transport>>, which handles this adjustment, is provided in Section 16.1.

As mentioned earlier, when there are three color channels, the first channel encodes luminance and the next two are red and blue, respectively. To see how to compute a green channel value, consider the implementation of the function RGBToXYZ(), which uses the following equation to compute from red, green, and blue color components assuming the color primaries from the sRGB standard:

In this case, we know , , and . Solving for , we can find:

As before, any color coefficients with negative values due to error in Fourier reconstruction must be clamped to zero.

We’ll now define the Fourier() function, which takes an array of coefficients , the maximum order , and the cosine of the angle . It evaluates the weighted sum of cosines in Equation (8.21), which can be written in the simpler form with now known:

The implementation of this function uses double precision for the sum of terms in order to minimize the impact of floating-point round-off error in computing the sum.

As the number of coefficients increases, a naïve evaluation of Equation (8.23) involves a correspondingly large number of trigonometric function calls. This can lead to severe performance issues: current CPU architectures require a few hundred processor cycles for a single invocation of std::cos(). Therefore, it pays to use the multiple angle formula for cosines:

This expression expresses cosine of summand in Equation (8.23) in terms of those used for the summands and .

The implementation starts with the declaration of two variables for the current and preceding cosine variables, corresponding to values for the indices and . Here, it’s important to also use double precision to compute the values; once m has values in the thousands, accumulated floating-point rounding error with 32-bit floats can become noticeable when using the multiple angle formula.

The body of the loop then adds the product of the current coefficient and cosine value to a running sum and computes the cosine values for the next iteration.

8.6.1 Spline Interpolation

The last detail to explain is how the spline-based interpolation used to reconstruct the coefficients works. The implementation here uses the Catmull–Rom spline, which can be expressed in 1D as a weighted sum over four control points, where the weight and the particular control points used depend on the parametric location along the curve’s path where its value is being computed. The spline passes through the given control points and follows a fairly smooth curve along the way.

To understand how these weights are computed, first suppose we are given a set of values of a function and its derivative at positions . For each interval , we would like to approximate the function using a cubic polynomial

which is chosen so that it matches the function and its derivative at the sample locations, i.e., , , , and . For simplicity, let us just focus on the first interval and furthermore suppose that . Solving for the coefficients and yields

Note how all of the coefficients are linear in the function and derivative values, which lets us rewrite Equation (8.25) as

This kind of interpolant is convenient but unfortunately still too restrictive, since we cannot generally expect derivative information to be available: analytic derivatives of reflectance models are often cumbersome, and measured data does not provide them at all. We therefore estimate the derivatives at each using central differences based on the two adjacent function values and . The estimated derivative is then

Similarly, the derivative at can be estimated using the two adjacent function values:

If we substitute these two expressions into Equation (8.26) and again collect terms, we have:

Note that the weights are independent of the function values: we can therefore also express this interpolation as

with

The CatmullRomWeights() function takes the variable x and the number of interpolation nodes and their positions as arguments. It does not use the function values in any way but instead computes the index offset and an array with four weights corresponding to the expressions in Equation (8.28).

The code that computes these weights is useful beyond the task of BSDF reconstruction and is thus defined in the files core/interpolation.h and core/interpolation.cpp.

The first statement returns a failure when is outside the domain of the function. Note the somewhat peculiar way of writing the conditional logic in negated form: this way, we can also catch NaN arguments, which by convention cause comparisons to evaluate to false.

The FindInterval() helper function efficiently locates the index of the interval containing x via binary search. With its result, we can now set the *offset return value to the index of the node and set variables x0 and x1 that delimit the extent of the domain of the corresponding spline segment.

Note that it’s possible that this offset would cause an out-of-bounds array access when Equation (8.27) is evaluated. (Specifically, in the case where the offset is one element before the start of the nodes array, when idx == 0, or if idx equals the size of the array minus one.) In these cases, the corresponding interpolation weights will always be set to zero for any out-of-bounds entries. Code that uses these weights in pbrt is therefore carefully written to never access the function values array for any indices where the weight is zero.

Because our derivation of the spline assumed the unit interval, we’ll define a scaled variable t in in the code here. It’s also useful to precompute some integer powers of t.

The implementation starts by initializing the second and third weights and using the results from Equation (8.28). For starters, only the terms in parenthesis are included.

There are two important details involved in computing the weights and from Equation (8.28). First, we need to introduce a scale factor of x1-x0, which corrects for the fact that the t values used in the code here incorporate a rescaling to the unit interval, while we actually want derivatives with respect to the original parameterization of the function.

Second, it’s necessary to handle an edge condition: the usual case is that idx > 0 and a previous neighbor exists; in this case, weights[0] can be initialized directly and the term can be added to weights[2], completing its initialization. If there is no previous neighbor, then the derivative is instead approximated with the forward difference . In this case, a similar algebraic process can be followed as was used to find the weights in Equation (8.28); the result is used here.

The computation for follows similarly and so the fragment that implements this part of the function, <<Compute last node weight >>, isn’t included here.

Given this machinery, we can now define the implementation of the FourierBDFTable:: GetWeightsAndOffsets() method, which just calls into CatmullRomWeights(), passing it the sampled mu array.