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 omega Subscript normal i and omega Subscript normal o 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.

Figure 8.24: Dragon models rendered using the 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 mu Subscript normal i and mu Subscript normal o denote the cosines of the incident and outgoing zenith angles, respectively, and phi Subscript normal i and phi Subscript normal o are the azimuth angles:

f left-parenthesis omega Subscript normal i Baseline comma omega Subscript normal o Baseline right-parenthesis equals f left-parenthesis mu Subscript normal i Baseline comma phi Subscript normal i Baseline comma mu Subscript normal o Baseline comma phi Subscript normal o Baseline right-parenthesis period

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 phi equals phi Subscript normal i Baseline minus phi Subscript normal o :

f left-parenthesis omega Subscript normal i Baseline comma omega Subscript normal o Baseline right-parenthesis equals f left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline comma phi Subscript normal i Baseline minus phi Subscript normal o Baseline right-parenthesis equals f left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline comma phi right-parenthesis period

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

f left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline comma phi right-parenthesis equals f left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline comma negative phi right-parenthesis period

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

f left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline comma phi Subscript normal i Baseline minus phi Subscript normal o Baseline right-parenthesis StartAbsoluteValue mu Subscript normal i Baseline EndAbsoluteValue equals sigma-summation Underscript k equals 0 Overscript m minus 1 Endscripts a Subscript k Baseline left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline right-parenthesis cosine left-parenthesis k left-parenthesis phi Subscript normal i Baseline minus phi Subscript normal o Baseline right-parenthesis right-parenthesis

Note how only cosine terms and no sine terms are needed due to Equation (8.20). The function evaluations a 0 left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline right-parenthesis comma ellipsis comma a Subscript m minus 1 Baseline left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline right-parenthesis denote the Fourier coefficients for a specific pair of zenith angle cosines.

Next, the functions a Subscript k are discretized over their input arguments. We choose a set of zenith angle cosines mu overbar equals StartSet mu 0 comma ellipsis comma mu Subscript n minus 1 Baseline EndSet and store the values of a Subscript k Baseline left-parenthesis mu Subscript i Baseline comma mu Subscript j Baseline right-parenthesis for every pair 0 less-than-or-equal-to i comma j less-than n . Thus, we can think of each a Subscript k as a n times n matrix, and the entire BRDF representation then consists of a set of m such matrices. Each describes a different azimuthal oscillation frequency in the material’s response to incident illumination.

The maximum order m 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 a Subscript k 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 mu Subscript normal i Baseline almost-equals mu Subscript normal o , many coefficients are necessary to accurately represent the specular lobe, which is zero for almost all azimuth angle differences phi equals phi Subscript normal i Baseline minus phi Subscript normal o and then very large for a small set of directions around phi equals pi , where the incident and outgoing directions are nearly opposite. However, when mu Subscript normal i and mu Subscript normal o 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 mu Subscript normal i and mu Subscript normal o angles require multiple coefficients a Subscript k to represent their phi distribution accurately, but their smoothness means that not too many coefficients are generally needed for each a Subscript k . 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.

<<BSDF Declarations>>+=  
struct FourierBSDFTable { <<FourierBSDFTable Public Data>> 
Float eta; int mMax; int nChannels; int nMu; Float *mu; int *m; int *aOffset; Float *a; Float *a0; Float *cdf; Float *recip;
<<FourierBSDFTable Public Methods>> 
static bool Read(const std::string &filename, FourierBSDFTable *table); const Float *GetAk(int offsetI, int offsetO, int *mptr) const { *mptr = m[offsetO * nMu + offsetI]; return a + aOffset[offsetO * nMu + offsetI]; } bool GetWeightsAndOffset(Float cosTheta, int *offset, Float weights[4]) const;
};

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.

<<FourierBSDFTable Public Methods>>= 
static bool Read(const std::string &filename, FourierBSDFTable *table);

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 m for any pair of mu Subscript normal i , mu Subscript normal o directions; this upper limit is useful when allocating temporary buffers to store a Subscript k coefficients, for example.

<<FourierBSDFTable Public Data>>= 
Float eta; int mMax;

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.

<<FourierBSDFTable Public Data>>+=  
int nChannels;

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 mu Subscript normal i or mu Subscript normal o angle.

<<FourierBSDFTable Public Data>>+=  
int nMu; Float *mu;

To evaluate Equation (8.21), we need to know the target Fourier order m and all coefficients a 0 comma ellipsis comma a Subscript m minus 1 Baseline corresponding to the directions omega Subscript normal i and omega Subscript normal o . For simplicity now, we’ll present the basic ideas as if only the coefficients for the closest mu directions less than or equal to mu Subscript normal i and mu Subscript normal o are used, though the implementation to follow interpolates between coefficients from multiple mu values around the directions.

The order m of the Fourier representation is always bounded by mMax but varies with respect to the incident and outgoing zenith angle cosine mu Subscript normal i and mu Subscript normal o : how many orders are needed can be determined by querying an monospace n monospace upper M monospace u times monospace n monospace upper M monospace u integer matrix m.

<<FourierBSDFTable Public Data>>+=  
int *m;

To find m 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

StartLayout 1st Row 1st Column monospace m monospace u monospace left-bracket monospace o monospace i monospace right-bracket 2nd Column less-than-or-equal-to mu Subscript normal i Baseline less-than monospace m monospace u monospace left-bracket monospace o monospace i monospace plus monospace 1 monospace right-bracket 2nd Row 1st Column monospace m monospace u monospace left-bracket monospace o monospace o monospace right-bracket 2nd Column less-than-or-equal-to mu Subscript normal o Baseline less-than monospace m monospace u monospace left-bracket monospace o monospace o monospace plus monospace 1 monospace right-bracket EndLayout

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

All of the a Subscript k 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 m coefficients starting at a[offset] give the a Subscript k values for the corresponding pair of directions. For the three color channel case, the first m coefficients after a[offset] encode coefficients for luminance, the next m correspond to the red channel, and then blue follows.

<<FourierBSDFTable Public Data>>+=  
int *aOffset; Float *a;

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

<<FourierBSDFTable Public Methods>>+= 
const Float *GetAk(int offsetI, int offsetO, int *mptr) const { *mptr = m[offsetO * nMu + offsetI]; return a + aOffset[offsetO * nMu + offsetI]; }

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

<<BxDF Declarations>>+= 
class FourierBSDF : public BxDF { public: <<FourierBSDF Public Methods>> 
Spectrum f(const Vector3f &wo, const Vector3f &wi) const; FourierBSDF(const FourierBSDFTable &bsdfTable, TransportMode mode) : BxDF(BxDFType(BSDF_REFLECTION | BSDF_TRANSMISSION | BSDF_GLOSSY)), bsdfTable(bsdfTable), mode(mode) { } Spectrum Sample_f(const Vector3f &wo, Vector3f *wi, const Point2f &u, Float *pdf, BxDFType *sampledType) const; Float Pdf(const Vector3f &wi, const Vector3f &wo) const;
private: <<FourierBSDF Private Data>> 
const FourierBSDFTable &bsdfTable; const TransportMode mode;
};

<<FourierBSDF Public Methods>>= 
FourierBSDF(const FourierBSDFTable &bsdfTable, TransportMode mode) : BxDF(BxDFType(BSDF_REFLECTION | BSDF_TRANSMISSION | BSDF_GLOSSY)), bsdfTable(bsdfTable), mode(mode) { }

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

<<FourierBSDF Private Data>>= 
const FourierBSDFTable &bsdfTable; const TransportMode mode;

Evaluating the BSDF is a matter of computing the cosines mu Subscript normal i and mu Subscript normal o , finding the corresponding coefficients a Subscript k and maximum order, and then evaluating Equation (8.21).

<<BxDF Method Definitions>>+=  
Spectrum FourierBSDF::f(const Vector3f &wo, const Vector3f &wi) const { <<Find the zenith angle cosines and azimuth difference angle>> 
Float muI = CosTheta(-wi), muO = CosTheta(wo); Float cosPhi = CosDPhi(-wi, wo);
<<Compute Fourier coefficients a Subscript k for left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline right-parenthesis >> 
<<Determine offsets and weights for mu Subscript normal i and mu Subscript normal o >> 
int offsetI, offsetO; Float weightsI[4], weightsO[4]; if (!bsdfTable.GetWeightsAndOffset(muI, &offsetI, weightsI) || !bsdfTable.GetWeightsAndOffset(muO, &offsetO, weightsO)) return Spectrum(0.f);
<<Allocate storage to accumulate ak coefficients>> 
Float *ak = ALLOCA(Float, bsdfTable.mMax * bsdfTable.nChannels); memset(ak, 0, bsdfTable.mMax * bsdfTable.nChannels * sizeof(Float));
<<Accumulate weighted sums of nearby a Subscript k coefficients>> 
int mMax = 0; for (int b = 0; b < 4; ++b) { for (int a = 0; a < 4; ++a) { <<Add contribution of (a, b) to a Subscript k values>> 
Float weight = weightsI[a] * weightsO[b]; if (weight != 0) { int m; const Float *ap = bsdfTable.GetAk(offsetI + a, offsetO + b, &m); mMax = std::max(mMax, m); for (int c = 0; c < bsdfTable.nChannels; ++c) for (int k = 0; k < m; ++k) ak[c * bsdfTable.mMax + k] += weight * ap[c * m + k]; }
} }
<<Evaluate Fourier expansion for angle phi >> 
Float Y = std::max((Float)0, Fourier(ak, mMax, cosPhi)); Float scale = muI != 0 ? (1 / std::abs(muI)) : (Float)0; <<Update scale to account for adjoint light transport>> 
if (mode == TransportMode::Radiance && muI * muO > 0) { float eta = muI > 0 ? 1 / bsdfTable.eta : bsdfTable.eta; scale *= eta * eta; }
if (bsdfTable.nChannels == 1) return Spectrum(Y * scale); else { <<Compute and return RGB colors for tabulated BSDF>> 
Float R = Fourier(ak + 1 * bsdfTable.mMax, mMax, cosPhi); Float B = Fourier(ak + 2 * bsdfTable.mMax, mMax, cosPhi); Float G = 1.39829f * Y - 0.100913f * B - 0.297375f * R; Float rgb[3] = { R * scale, G * scale, B * scale }; return Spectrum::FromRGB(rgb).Clamp();
}
}

There is an important difference of convention in how directions are represented within the FourierBSDF: the incident direction omega Subscript normal i 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.

<<Find the zenith angle cosines and azimuth difference angle>>= 
Float muI = CosTheta(-wi), muO = CosTheta(wo); Float cosPhi = CosDPhi(-wi, wo);

So that the reconstructed BSDF is fairly smooth, the implementation here interpolates a Subscript k coefficients over the product of the four quantized mu directions that surround mu Subscript normal i and the four that surround mu Subscript normal o . 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 a Subscript k is thus computed by

a Subscript k Baseline equals sigma-summation Underscript a equals 0 Overscript 3 Endscripts sigma-summation Underscript b equals 0 Overscript 3 Endscripts a Subscript k Baseline left-parenthesis o Subscript normal i Baseline plus a comma o Subscript normal o Baseline plus b right-parenthesis w Subscript normal i Baseline left-parenthesis a right-parenthesis w Subscript normal o Baseline left-parenthesis b right-parenthesis comma

where a Subscript k Baseline left-parenthesis i comma j right-parenthesis gives the k th Fourier coefficient for the quantized directions mu Subscript i Baseline comma mu Subscript j Baseline and w Subscript normal i and w Subscript normal o are the spline weights. This interpolation ensures sufficient smoothness even when the discretization of directions mu Subscript i is relatively coarse; the details of how these weights are computed will be explained in a few pages.

<<Compute Fourier coefficients a Subscript k for left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline right-parenthesis >>= 
<<Determine offsets and weights for mu Subscript normal i and mu Subscript normal o >> 
int offsetI, offsetO; Float weightsI[4], weightsO[4]; if (!bsdfTable.GetWeightsAndOffset(muI, &offsetI, weightsI) || !bsdfTable.GetWeightsAndOffset(muO, &offsetO, weightsO)) return Spectrum(0.f);
<<Allocate storage to accumulate ak coefficients>> 
Float *ak = ALLOCA(Float, bsdfTable.mMax * bsdfTable.nChannels); memset(ak, 0, bsdfTable.mMax * bsdfTable.nChannels * sizeof(Float));
<<Accumulate weighted sums of nearby a Subscript k coefficients>> 
int mMax = 0; for (int b = 0; b < 4; ++b) { for (int a = 0; a < 4; ++a) { <<Add contribution of (a, b) to a Subscript k values>> 
Float weight = weightsI[a] * weightsO[b]; if (weight != 0) { int m; const Float *ap = bsdfTable.GetAk(offsetI + a, offsetO + b, &m); mMax = std::max(mMax, m); for (int c = 0; c < bsdfTable.nChannels; ++c) for (int k = 0; k < m; ++k) ak[c * bsdfTable.mMax + k] += weight * ap[c * m + k]; }
} }

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

<<Determine offsets and weights for mu Subscript normal i and mu Subscript normal o >>= 
int offsetI, offsetO; Float weightsI[4], weightsO[4]; if (!bsdfTable.GetWeightsAndOffset(muI, &offsetI, weightsI) || !bsdfTable.GetWeightsAndOffset(muO, &offsetO, weightsO)) return Spectrum(0.f);

The various a Subscript k vectors in the 4 times 4 extent of the directions being interpolated over may have different orders m . Therefore, the implementation here allocates storage for the a Subscript k values using the maximum possible order m 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.

<<Allocate storage to accumulate ak coefficients>>= 
Float *ak = ALLOCA(Float, bsdfTable.mMax * bsdfTable.nChannels); memset(ak, 0, bsdfTable.mMax * bsdfTable.nChannels * sizeof(Float));

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

<<Accumulate weighted sums of nearby a Subscript k coefficients>>= 
int mMax = 0; for (int b = 0; b < 4; ++b) { for (int a = 0; a < 4; ++a) { <<Add contribution of (a, b) to a Subscript k values>> 
Float weight = weightsI[a] * weightsO[b]; if (weight != 0) { int m; const Float *ap = bsdfTable.GetAk(offsetI + a, offsetO + b, &m); mMax = std::max(mMax, m); for (int c = 0; c < bsdfTable.nChannels; ++c) for (int k = 0; k < m; ++k) ak[c * bsdfTable.mMax + k] += weight * ap[c * m + k]; }
} }

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.

<<Add contribution of (a, b) to a Subscript k values>>= 
Float weight = weightsI[a] * weightsO[b]; if (weight != 0) { int m; const Float *ap = bsdfTable.GetAk(offsetI + a, offsetO + b, &m); mMax = std::max(mMax, m); for (int c = 0; c < bsdfTable.nChannels; ++c) for (int k = 0; k < m; ++k) ak[c * bsdfTable.mMax + k] += weight * ap[c * m + k]; }

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 a Subscript k 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.

<<Evaluate Fourier expansion for angle phi >>= 
Float Y = std::max((Float)0, Fourier(ak, mMax, cosPhi)); Float scale = muI != 0 ? (1 / std::abs(muI)) : (Float)0; <<Update scale to account for adjoint light transport>> 
if (mode == TransportMode::Radiance && muI * muO > 0) { float eta = muI > 0 ? 1 / bsdfTable.eta : bsdfTable.eta; scale *= eta * eta; }
if (bsdfTable.nChannels == 1) return Spectrum(Y * scale); else { <<Compute and return RGB colors for tabulated BSDF>> 
Float R = Fourier(ak + 1 * bsdfTable.mMax, mMax, cosPhi); Float B = Fourier(ak + 2 * bsdfTable.mMax, mMax, cosPhi); Float G = 1.39829f * Y - 0.100913f * B - 0.297375f * R; Float rgb[3] = { R * scale, G * scale, B * scale }; return Spectrum::FromRGB(rgb).Clamp();
}

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 y Subscript lamda from red, green, and blue color components assuming the color primaries from the sRGB standard:

y Subscript lamda Baseline equals 0.212671 r plus 0.715160 g plus 0.072169 b period

In this case, we know y Subscript lamda , r , and b . Solving for g , we can find:

g equals 1.39829 y Subscript lamda Baseline minus 0.100913 b minus 0.297375 r period

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

<<Compute and return RGB colors for tabulated BSDF>>= 
Float R = Fourier(ak + 1 * bsdfTable.mMax, mMax, cosPhi); Float B = Fourier(ak + 2 * bsdfTable.mMax, mMax, cosPhi); Float G = 1.39829f * Y - 0.100913f * B - 0.297375f * R; Float rgb[3] = { R * scale, G * scale, B * scale }; return Spectrum::FromRGB(rgb).Clamp();

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

f left-parenthesis phi right-parenthesis equals sigma-summation Underscript k equals 0 Overscript m minus 1 Endscripts a Subscript k Baseline cosine left-parenthesis k phi right-parenthesis period

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.

<<Fourier Interpolation Definitions>>= 
Float Fourier(const Float *a, int m, double cosPhi) { double value = 0.0; <<Initialize cosine iterates>> 
double cosKMinusOnePhi = cosPhi; double cosKPhi = 1;
for (int k = 0; k < m; ++k) { <<Add the current summand and update the cosine iterates>> 
value += a[k] * cosKPhi; double cosKPlusOnePhi = 2 * cosPhi * cosKPhi - cosKMinusOnePhi; cosKMinusOnePhi = cosKPhi; cosKPhi = cosKPlusOnePhi;
} return value; }

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:

cosine left-parenthesis k phi right-parenthesis equals left-parenthesis 2 cosine phi right-parenthesis cosine left-parenthesis left-parenthesis k minus 1 right-parenthesis phi right-parenthesis minus cosine left-parenthesis left-parenthesis k minus 2 right-parenthesis phi right-parenthesis

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

The implementation starts with the declaration of two variables for the current and preceding cosine variables, corresponding to values for the indices k equals negative 1 and k equals 0 . Here, it’s important to also use double precision to compute the cosine left-parenthesis k phi right-parenthesis 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.

<<Initialize cosine iterates>>= 
double cosKMinusOnePhi = cosPhi; double cosKPhi = 1;

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.

<<Add the current summand and update the cosine iterates>>= 
value += a[k] * cosKPhi; double cosKPlusOnePhi = 2 * cosPhi * cosKPhi - cosKMinusOnePhi; cosKMinusOnePhi = cosKPhi; cosKPhi = cosKPlusOnePhi;

8.6.1 Spline Interpolation

The last detail to explain is how the spline-based interpolation used to reconstruct the a Subscript k 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 f and its derivative f prime at positions x 0 comma x 1 comma ellipsis comma x Subscript k Baseline . For each interval left-bracket x Subscript i Baseline comma x Subscript i plus 1 Baseline right-bracket , we would like to approximate the function using a cubic polynomial

p Subscript i Baseline left-parenthesis x right-parenthesis equals a x cubed plus b x squared plus c x plus d comma

which is chosen so that it matches the function and its derivative at the sample locations, i.e., p Subscript i Baseline left-parenthesis x Subscript i Baseline right-parenthesis equals f left-parenthesis x Subscript i Baseline right-parenthesis , p Subscript i Baseline left-parenthesis x Subscript i plus 1 Baseline right-parenthesis equals f left-parenthesis x Subscript i plus 1 Baseline right-parenthesis , p prime Subscript i Baseline left-parenthesis x Subscript i Baseline right-parenthesis equals f prime left-parenthesis x Subscript i Baseline right-parenthesis , and p prime Subscript i Baseline left-parenthesis x Subscript i plus 1 Baseline right-parenthesis equals f prime left-parenthesis x Subscript i plus 1 Baseline right-parenthesis . For simplicity, let us just focus on the first interval and furthermore suppose that left-bracket x 0 comma x 1 right-bracket equals left-bracket 0 comma 1 right-bracket . Solving for the coefficients a comma b comma c comma and d yields

StartLayout 1st Row 1st Column a 2nd Column equals f prime left-parenthesis x 0 right-parenthesis plus f prime left-parenthesis x 1 right-parenthesis plus 2 f left-parenthesis x 0 right-parenthesis minus 2 f left-parenthesis x 1 right-parenthesis comma 2nd Row 1st Column b 2nd Column equals 3 f left-parenthesis x 1 right-parenthesis minus 3 f left-parenthesis x 0 right-parenthesis minus 2 f prime left-parenthesis x 0 right-parenthesis minus f prime left-parenthesis x 1 right-parenthesis comma 3rd Row 1st Column c 2nd Column equals f prime left-parenthesis x 0 right-parenthesis comma 4th Row 1st Column d 2nd Column equals f left-parenthesis x 0 right-parenthesis period EndLayout

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

StartLayout 1st Row 1st Column p left-parenthesis x right-parenthesis 2nd Column equals left-parenthesis 2 x cubed minus 3 x squared plus 1 right-parenthesis f left-parenthesis x 0 right-parenthesis 2nd Row 1st Column Blank 2nd Column plus left-parenthesis minus 2 x cubed plus 3 x squared right-parenthesis f left-parenthesis x 1 right-parenthesis 3rd Row 1st Column Blank 2nd Column plus left-parenthesis x cubed minus 2 x squared plus x right-parenthesis f prime left-parenthesis x 0 right-parenthesis 4th Row 1st Column Blank 2nd Column plus left-parenthesis x cubed minus x squared right-parenthesis f prime left-parenthesis x 1 right-parenthesis period EndLayout

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 f left-parenthesis x Subscript i Baseline right-parenthesis using central differences based on the two adjacent function values f left-parenthesis x Subscript i minus 1 Baseline right-parenthesis and f left-parenthesis x Subscript i plus 1 Baseline right-parenthesis . The estimated derivative is then

f prime left-parenthesis x 0 right-parenthesis almost-equals StartFraction f left-parenthesis x 1 right-parenthesis minus f left-parenthesis x Subscript negative 1 Baseline right-parenthesis Over x 1 minus x Subscript negative 1 Baseline EndFraction equals StartFraction f left-parenthesis x 1 right-parenthesis minus f left-parenthesis x Subscript negative 1 Baseline right-parenthesis Over 1 minus x Subscript negative 1 Baseline EndFraction period

Similarly, the derivative at f left-parenthesis x 1 right-parenthesis can be estimated using the two adjacent function values:

f prime left-parenthesis x 1 right-parenthesis almost-equals StartFraction f left-parenthesis x 2 right-parenthesis minus f left-parenthesis x 0 right-parenthesis Over x 2 minus x 0 EndFraction equals StartFraction f left-parenthesis x 2 right-parenthesis minus f left-parenthesis x 0 right-parenthesis Over x 2 EndFraction period

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

StartLayout 1st Row 1st Column p left-parenthesis x right-parenthesis 2nd Column equals StartFraction x cubed minus 2 x squared plus x Over x Subscript negative 1 Baseline minus 1 EndFraction f left-parenthesis x Subscript negative 1 Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column plus left-parenthesis 2 x cubed minus 3 x squared plus 1 minus StartFraction x cubed minus x squared Over x 2 EndFraction right-parenthesis f left-parenthesis x 0 right-parenthesis 3rd Row 1st Column Blank 2nd Column plus left-parenthesis minus 2 x cubed plus 3 x squared plus StartFraction x cubed minus 2 x squared plus x Over 1 minus x Subscript negative 1 Baseline EndFraction right-parenthesis f left-parenthesis x 1 right-parenthesis 4th Row 1st Column Blank 2nd Column plus StartFraction x cubed minus x squared Over x 2 EndFraction f left-parenthesis x 2 right-parenthesis comma EndLayout

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

p left-parenthesis x right-parenthesis equals w 0 f left-parenthesis x Subscript negative 1 Baseline right-parenthesis plus w 1 f left-parenthesis x 0 right-parenthesis plus w 2 f left-parenthesis x 1 right-parenthesis plus w 3 f left-parenthesis x 2 right-parenthesis comma

with

StartLayout 1st Row 1st Column w 0 2nd Column equals StartFraction x cubed minus 2 x squared plus x Over x Subscript negative 1 Baseline minus 1 EndFraction 2nd Row 1st Column w 1 2nd Column equals 2 x cubed minus 3 x squared plus 1 minus StartFraction x cubed minus x squared Over x 2 EndFraction equals left-parenthesis 2 x cubed minus 3 x squared plus 1 right-parenthesis minus w 3 3rd Row 1st Column w 2 2nd Column equals minus 2 x cubed plus 3 x squared plus StartFraction x cubed minus 2 x squared plus x Over 1 minus x Subscript negative 1 Baseline EndFraction equals left-parenthesis minus 2 x cubed plus 3 x squared right-parenthesis minus w 0 4th Row 1st Column w 3 2nd Column equals StartFraction x cubed minus x squared Over x 2 EndFraction period EndLayout

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.

<<Spline Interpolation Definitions>>= 
bool CatmullRomWeights(int size, const Float *nodes, Float x, int *offset, Float *weights) { <<Return false if x is out of bounds>> 
if (!(x >= nodes[0] && x <= nodes[size-1])) return false;
<<Search for the interval idx containing x>> 
int idx = FindInterval(size, [&](int i) { return nodes[i] <= x; }); *offset = idx - 1; Float x0 = nodes[idx], x1 = nodes[idx+1];
<<Compute the t parameter and powers>> 
Float t = (x - x0) / (x1 - x0), t2 = t * t, t3 = t2 * t;
<<Compute initial node weights w 1 and w 2 >> 
weights[1] = 2 * t3 - 3 * t2 + 1; weights[2] = -2 * t3 + 3 * t2;
<<Compute first node weight w 0 >> 
if (idx > 0) { Float w0 = (t3 - 2 * t2 + t) * (x1 - x0) / (x1 - nodes[idx - 1]); weights[0] = -w0; weights[2] += w0; } else { Float w0 = t3 - 2 * t2 + t; weights[0] = 0; weights[1] -= w0; weights[2] += w0; }
<<Compute last node weight w 3 >> 
if (idx+2 < size) { Float w3 = (t3 - t2) * (x1 - x0) / (nodes[idx+2] - x0); weights[1] -= w3; weights[3] = w3; } else { Float w3 = t3 - t2; weights[1] -= w3; weights[2] += w3; weights[3] = 0; }
return true; }

The first statement returns a failure when monospace x 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.

<<Return false if x is out of bounds>>= 
if (!(x >= nodes[0] && x <= nodes[size-1])) return 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 x Subscript i minus 1 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.

<<Search for the interval idx containing x>>= 
int idx = FindInterval(size, [&](int i) { return nodes[i] <= x; }); *offset = idx - 1; Float x0 = nodes[idx], x1 = nodes[idx+1];

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

<<Compute the t parameter and powers>>= 
Float t = (x - x0) / (x1 - x0), t2 = t * t, t3 = t2 * t;

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

<<Compute initial node weights w 1 and w 2 >>= 
weights[1] = 2 * t3 - 3 * t2 + 1; weights[2] = -2 * t3 + 3 * t2;

There are two important details involved in computing the weights w 0 and w 3 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 w 0 term can be added to weights[2], completing its initialization. If there is no previous neighbor, then the derivative f prime left-parenthesis x 0 right-parenthesis is instead approximated with the forward difference f left-parenthesis x 1 right-parenthesis minus f left-parenthesis x 0 right-parenthesis . 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.

<<Compute first node weight w 0 >>= 
if (idx > 0) { Float w0 = (t3 - 2 * t2 + t) * (x1 - x0) / (x1 - nodes[idx - 1]); weights[0] = -w0; weights[2] += w0; } else { Float w0 = t3 - 2 * t2 + t; weights[0] = 0; weights[1] -= w0; weights[2] += w0; }

The computation for w 3 follows similarly and so the fragment that implements this part of the function, <<Compute last node weight w 3 >>, 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.

<<BxDF Method Definitions>>+=  
bool FourierBSDFTable::GetWeightsAndOffset(Float cosTheta, int *offset, Float weights[4]) const { return CatmullRomWeights(nMu, mu, cosTheta, offset, weights); }