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

`FourierBSDF`model. The surface of the dragon on the left has a BSDF that models the appearence of rough gold; the one on the right is coated copper.

*(Model courtesy of Christian Schüller.)*

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:

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

`ak`coefficients>>

`(a, b)`to values>> } }

`scale`to account for adjoint light transport>>

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.

`ak`coefficients>>

`(a, b)`to values>> } }

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.

`ak`coefficients>>=

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

`(a, b)`to values>> } }

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

`(a, b)`to values>>=

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.

`scale`to account for adjoint light transport>>

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 `float`s 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:

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

`false`if

`x`is out of bounds>>

`idx`containing

`x`>>

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

`false`if

`x`is out of bounds>>=

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.

`idx`containing

`x`>>=

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.