14.1 Sampling Reflection Functions

The BxDF::Sample_f() method chooses a direction according to a distribution that is similar to its corresponding scattering function. In Section 8.2, this method was used for finding reflected and transmitted rays from perfectly specular surfaces; later in this section, we will show how that sampling process is a special case of the sampling techniques we’ll now implement for other types of BSDFs.

BxDF::Sample_f() takes two sample values in the range left-bracket 0 comma 1 right-parenthesis squared that are intended to be used by an inversion method-based sampling algorithm (recall Section 13.3.1). The routine calling it can use stratified or low-discrepancy sampling techniques to generate these sample values, thus making it possible for the sampled directions themselves to be well distributed. Other sampling methods like rejection sampling could in theory be supported by passing a Sampler instance, though this is not done in pbrt as stratified sampling of a distribution that is similar to the BSDF generally produces superior results.

This method returns the value of the BSDF for the chosen direction as well as the sampled direction in *wi and the value of p left-parenthesis omega Subscript normal i Baseline right-parenthesis in *pdf. The PDF value returned should be measured with respect to solid angle, and both the outgoing direction omega Subscript normal o and the sampled incident direction omega Subscript normal i should be in the standard reflection coordinate system (see Section “Geometric Setting” at the start of Chapter 8).

The default implementation of this method samples the unit hemisphere with a cosine-weighted distribution. Samples from this distribution will give correct results for any BRDF that isn’t described by a delta distribution, since there is some probability of sampling all directions where the BRDF’s value is non-0: p left-parenthesis omega Subscript Baseline right-parenthesis greater-than 0 for all omega Subscript in the hemisphere. (BTDFs will thus always need to override this method but can sample the opposite hemisphere uniformly if they don’t have a better sampling method.)

<<BxDF Method Definitions>>+=  
Spectrum BxDF::Sample_f(const Vector3f &wo, Vector3f *wi, const Point2f &u, Float *pdf, BxDFType *sampledType) const { <<Cosine-sample the hemisphere, flipping the direction if necessary>> 
*wi = CosineSampleHemisphere(u); if (wo.z < 0) wi->z *= -1;
*pdf = Pdf(wo, *wi); return f(wo, *wi); }

There is a subtlety related to the orientation of the normal that must be accounted for here: the direction returned by CosineSampleHemisphere() is in the hemisphere around left-parenthesis 0 comma 0 comma 1 right-parenthesis in the reflection coordinate system. If omega Subscript normal o is in the opposite hemisphere, then omega Subscript normal i must be flipped to lie in the same hemisphere as omega Subscript normal o . This issue is a direct consequence of the fact that pbrt does not flip the normal to be on the same side of the surface as the omega Subscript normal o direction.

<<Cosine-sample the hemisphere, flipping the direction if necessary>>= 
*wi = CosineSampleHemisphere(u); if (wo.z < 0) wi->z *= -1;

While the value of the PDF returned by BxDF::Sample_f() is for the direction it chose, the BxDF::Pdf() method returns the value of the PDF for a given pair of directions. This method is useful for multiple importance sampling, where it is necessary to be able to find one sampling distribution’s PDF for directions sampled from other distributions. It is crucial that any BxDF implementation that overrides the BxDF::Sample_f() method also override the BxDF::Pdf() method so that the two return consistent results.

To actually evaluate the PDF for the cosine-weighted sampling method (which we showed earlier was p left-parenthesis omega Subscript Baseline right-parenthesis equals cosine theta slash pi ), it is first necessary to check that omega Subscript normal o and omega Subscript normal i lie on the same side of the surface; if not, the sampling probability is 0. Otherwise, the method computes StartAbsoluteValue bold n Subscript Baseline dot omega Subscript normal i Baseline EndAbsoluteValue . One potential pitfall with this method is that the order of the omega Subscript normal o and omega Subscript normal i arguments is significant. For the cosine-weighted distribution, p left-parenthesis omega Subscript normal o Baseline right-parenthesis not-equals p left-parenthesis omega Subscript normal i Baseline right-parenthesis in general. Code that calls this method must be careful to use the correct argument ordering.

<<BxDF Method Definitions>>+=  
Float BxDF::Pdf(const Vector3f &wo, const Vector3f &wi) const { return SameHemisphere(wo, wi) ? AbsCosTheta(wi) * InvPi : 0; }

<<BSDF Inline Functions>>+= 
inline bool SameHemisphere(const Vector3f &w, const Vector3f &wp) { return w.z * wp.z > 0; }

This sampling method works well for Lambertian BRDFs, and it works well for the Oren–Nayar model as well, so we will not override it for those classes.

14.1.1 Microfacet BxDFs

The microfacet-based reflection models defined in Section 8.4 were based on a distribution of microfacets upper D left-parenthesis omega Subscript normal h Baseline right-parenthesis where each microfacet exhibited perfect specular reflection and/or transmission. Because the upper D left-parenthesis omega Subscript normal h Baseline right-parenthesis function is primarily responsible for the overall shape of the Torrance–Sparrow BSDF (Section 8.4.4), approaches based on sampling from its distribution are fairly effective. With this approach, first a particular microfacet orientation is sampled from the microfacet distribution, and then the incident direction is found using the specular reflection or transmission formula.

Therefore, MicrofacetDistribution implementations must implement a method for sampling from their distribution of normal vectors.

<<MicrofacetDistribution Public Methods>>+= 
virtual Vector3f Sample_wh(const Vector3f &wo, const Point2f &u) const = 0;

The classic approach to sampling a microfacet orientation is to sample from upper D left-parenthesis omega Subscript normal h Baseline right-parenthesis directly. We’ll start by showing the derivation of this approach for the isotropic Beckmann distribution but will then describe a more effective sampling method that samples from the distribution of visible microfacets from a given outgoing direction, which can be quite different from the overall distribution.

The MicrofacetDistribution base class stores a Boolean value that determines which sampling technique will be used. In practice, the one based on sampling visible microfacet area is much more effective than the one based on sampling the overall distribution, so it isn’t possible to select between the two strategies in pbrt scene description files; the option to sample the overall distribution is only available for tests and comparisons.

<<MicrofacetDistribution Protected Methods>>= 
MicrofacetDistribution(bool sampleVisibleArea) : sampleVisibleArea(sampleVisibleArea) { }

<<MicrofacetDistribution Protected Data>>= 
const bool sampleVisibleArea;

The BeckmannDistribution’s Sample_wh() method’s implementation uses this value to determine which sampling technique to use.

<<MicrofacetDistribution Method Definitions>>+=  
Vector3f BeckmannDistribution::Sample_wh(const Vector3f &wo, const Point2f &u) const { if (!sampleVisibleArea) { <<Sample full distribution of normals for Beckmann distribution>> 
<<Compute tangent squared theta and phi for Beckmann distribution sample>> 
Float tan2Theta, phi; if (alphax == alphay) { Float logSample = std::log(1 - u[0]); if (std::isinf(logSample)) logSample = 0; tan2Theta = -alphax * alphax * logSample; phi = u[1] * 2 * Pi; } else { <<Compute tan2Theta and phi for anisotropic Beckmann distribution>> 
Float logSample = std::log(u[0]); phi = std::atan(alphay / alphax * std::tan(2 * Pi * u[1] + 0.5f * Pi)); if (u[1] > 0.5f) phi += Pi; Float sinPhi = std::sin(phi), cosPhi = std::cos(phi); Float alphax2 = alphax * alphax, alphay2 = alphay * alphay; tan2Theta = -logSample / (cosPhi * cosPhi / alphax2 + sinPhi * sinPhi / alphay2);
}
<<Map sampled Beckmann angles to normal direction wh>> 
Float cosTheta = 1 / std::sqrt(1 + tan2Theta); Float sinTheta = std::sqrt(std::max((Float)0, 1 - cosTheta * cosTheta)); Vector3f wh = SphericalDirection(sinTheta, cosTheta, phi); if (!SameHemisphere(wo, wh)) wh = -wh;
return wh;
} else { <<Sample visible area of normals for Beckmann distribution>> 
Vector3f wh; bool flip = wo.z < 0; wh = BeckmannSample(flip ? -wo : wo, alphax, alphay, u[0], u[1]); if (flip) wh = -wh; return wh;
} }

The sampling method for the Beckmann–Spizzichino distribution’s full distribution of normals returns tangent squared theta and the angle phi in spherical coordinates, which in turn are converted to a normalized direction vector wh.

<<Sample full distribution of normals for Beckmann distribution>>= 
<<Compute tangent squared theta and phi for Beckmann distribution sample>> 
Float tan2Theta, phi; if (alphax == alphay) { Float logSample = std::log(1 - u[0]); if (std::isinf(logSample)) logSample = 0; tan2Theta = -alphax * alphax * logSample; phi = u[1] * 2 * Pi; } else { <<Compute tan2Theta and phi for anisotropic Beckmann distribution>> 
Float logSample = std::log(u[0]); phi = std::atan(alphay / alphax * std::tan(2 * Pi * u[1] + 0.5f * Pi)); if (u[1] > 0.5f) phi += Pi; Float sinPhi = std::sin(phi), cosPhi = std::cos(phi); Float alphax2 = alphax * alphax, alphay2 = alphay * alphay; tan2Theta = -logSample / (cosPhi * cosPhi / alphax2 + sinPhi * sinPhi / alphay2);
}
<<Map sampled Beckmann angles to normal direction wh>> 
Float cosTheta = 1 / std::sqrt(1 + tan2Theta); Float sinTheta = std::sqrt(std::max((Float)0, 1 - cosTheta * cosTheta)); Vector3f wh = SphericalDirection(sinTheta, cosTheta, phi); if (!SameHemisphere(wo, wh)) wh = -wh;
return wh;

The isotropic Beckmann–Spizzichino distribution was defined in Equation (8.9). To derive a sampling method, we’ll consider it expressed in spherical coordinates. As an isotropic distribution, it is independent of phi , and so the PDF for this distribution, p Subscript normal h Baseline left-parenthesis theta comma phi right-parenthesis , is separable into p Subscript normal h Baseline left-parenthesis theta right-parenthesis and p Subscript normal h Baseline left-parenthesis phi right-parenthesis .

p Subscript normal h Baseline left-parenthesis phi right-parenthesis is constant with a value of 1 slash left-parenthesis 2 pi right-parenthesis , and thus phi values can be sampled by

phi equals 2 pi xi Subscript Baseline period

For p Subscript normal h Baseline left-parenthesis theta right-parenthesis , we have

p Subscript normal h Baseline left-parenthesis theta right-parenthesis equals StartFraction 2 normal e Superscript minus tangent squared theta slash alpha squared Baseline sine theta Over alpha squared cosine cubed theta EndFraction comma

where alpha is the roughness coefficient. We can apply the inversion method to find how to sample a direction theta prime from this distribution given a uniform random number xi Subscript . First, taking the PDF from Equation (14.1), and integrating to find the CDF, we have

StartLayout 1st Row 1st Column upper P Subscript normal h Baseline left-parenthesis theta prime right-parenthesis 2nd Column equals integral Subscript 0 Superscript theta Superscript prime Baseline Baseline StartFraction 2 normal e Superscript minus tangent squared theta slash alpha squared Baseline sine theta Over alpha squared cosine cubed theta EndFraction normal d theta Subscript Baseline 2nd Row 1st Column Blank 2nd Column equals 1 minus normal e Superscript minus tangent squared theta Super Superscript prime Superscript slash alpha squared Baseline period EndLayout

To find the sampling technique, we need to solve

xi Subscript Baseline equals 1 minus normal e Superscript minus tangent squared theta prime slash alpha squared

for theta prime in terms of xi Subscript . In this case, tangent squared theta prime suffices to find the microfacet orientation and is more efficient to compute, so we will compute

tangent squared theta Superscript prime Baseline equals minus alpha squared log left-parenthesis 1 minus xi Subscript Baseline right-parenthesis period

The sampling code follows directly.

<<Compute tangent squared theta and phi for Beckmann distribution sample>>= 
Float tan2Theta, phi; if (alphax == alphay) { Float logSample = std::log(1 - u[0]); if (std::isinf(logSample)) logSample = 0; tan2Theta = -alphax * alphax * logSample; phi = u[1] * 2 * Pi; } else { <<Compute tan2Theta and phi for anisotropic Beckmann distribution>> 
Float logSample = std::log(u[0]); phi = std::atan(alphay / alphax * std::tan(2 * Pi * u[1] + 0.5f * Pi)); if (u[1] > 0.5f) phi += Pi; Float sinPhi = std::sin(phi), cosPhi = std::cos(phi); Float alphax2 = alphax * alphax, alphay2 = alphay * alphay; tan2Theta = -logSample / (cosPhi * cosPhi / alphax2 + sinPhi * sinPhi / alphay2);
}

The algorithm to sample tangent squared theta and phi for anisotropic Beckmann–Spizzichino distributions can be derived following a similar process, though we won’t include the derivation or implementation in the text here.

Given tangent squared theta , we can compute cosine theta using the identities tangent squared theta equals sine squared theta slash cosine squared theta and sine squared theta plus cosine squared theta equals 1 . sine theta follows, and we have enough information to compute the microfacet orientation using the spherical coordinates formula. Because pbrt transforms the normal to left-parenthesis 0 comma 0 comma 1 right-parenthesis in the reflection coordinate system, we can almost use the computed direction from spherical coordinates directly. The last detail to handle is that if omega Subscript normal o is in the opposite hemisphere than the normal, then the half-angle vector needs to be flipped to be in that hemisphere as well.

<<Map sampled Beckmann angles to normal direction wh>>= 
Float cosTheta = 1 / std::sqrt(1 + tan2Theta); Float sinTheta = std::sqrt(std::max((Float)0, 1 - cosTheta * cosTheta)); Vector3f wh = SphericalDirection(sinTheta, cosTheta, phi); if (!SameHemisphere(wo, wh)) wh = -wh;

While sampling a microfacet orientation from the full microfacet distribution gives correct results, the efficiency of this approach is limited by the fact that only one term, upper D left-parenthesis omega Subscript normal h Baseline right-parenthesis , of the full microfacet BSDF (defined in Equation (8.18)) is accounted for. A better approach can be found using the observation that the distribution of microfacets that are visible from a given direction isn’t the same as the full distribution of microfacets; Figure 14.1 shows the geometric setting and describes why the distributions differ.

Figure 14.1: The distribution of visible microfacets from a direction oblique to the underlying geometric normal is quite different from the distribution upper D left-parenthesis omega Subscript normal h Baseline right-parenthesis . First, some microfacet orientations are backfacing and will never be seen. Others are shadowed by other microfacets. Finally, the projected area of a visible microfacet increases as its orientation approaches the viewing direction. These factors are accounted for in upper D Subscript omega Sub Subscript Baseline left-parenthesis omega Subscript normal h Baseline right-parenthesis , the distribution of visible microfacets in the direction omega Subscript .

Equation (8.12) in Section 8.4.2 defined a relationship between the visible area of microfacets from a given direction and microfacet orientation. This equation can be rearranged to get the distribution of visible normals in a direction omega Subscript :

upper D Subscript omega Sub Subscript Subscript Baseline left-parenthesis omega Subscript normal h Baseline right-parenthesis equals StartFraction upper D left-parenthesis omega Subscript normal h Baseline right-parenthesis upper G 1 left-parenthesis omega Subscript Baseline comma omega Subscript normal h Baseline right-parenthesis max left-parenthesis 0 comma omega Subscript Baseline dot omega Subscript normal h Baseline right-parenthesis Over cosine theta EndFraction period

Here, the upper G 1 term accounts for microfacet self-shadowing, and the max left-parenthesis 0 comma left-parenthesis omega Subscript Baseline dot omega Subscript normal h Baseline right-parenthesis right-parenthesis slash cosine theta term accounts for both backfacing microfacets and the interaction of microfacet orientation and projected area as a function of viewing direction that was shown in Figure 14.1.

Figure 14.2 compares the overall distribution of microfacets with the Beckmann–Spizzichino model with the visible distribution from a fairly oblique viewing direction. Note that many orientations are no longer visible at all (as they are backfacing) and that the microfacet orientations that are in the vicinity of the outgoing direction omega Subscript normal o have a higher probability of being visible than they do in the overall distribution upper D left-parenthesis omega Subscript normal h Baseline right-parenthesis .

Figure 14.2: (a) Full Beckmann–Spizzichino microfacet distribution upper D left-parenthesis omega Subscript normal h Baseline right-parenthesis for alpha equals 0.3 and (b) visible microfacet distribution upper D Subscript omega Sub Subscript normal o Baseline left-parenthesis omega Subscript normal h Baseline right-parenthesis for cosine theta Subscript normal o Baseline equals 0.1 . With this relatively oblique viewing direction, the two distributions are quite different and sampling from upper D Subscript omega Sub Subscript normal o Baseline left-parenthesis omega Subscript normal h Baseline right-parenthesis is a much more effective approach.

It turns out that samples can be taken directly from the distribution defined by Equation (14.2); because this distribution better matches the full Torrance–Sparrow model (Equation (8.18)) than upper D left-parenthesis omega Subscript normal h Baseline right-parenthesis alone, there is much less variance in resulting images (see Figure 14.3). We won’t go into the extensive details of how to directly sample this distribution or include the code in the book; see the “Further Reading” section and the source code, respectively, for more information. (The TrowbridgeReitzDistribution::Sample_wh() method similarly samples from either the full distribution of microfacet normals or from the visible distribution; see the source code for details.)

Figure 14.3: Comparison of Microfacet Sampling Techniques. With the same number of samples taken per pixel there is higher variance when (1) sampling the full microfacet distribution upper D left-parenthesis omega Subscript normal h Baseline right-parenthesis than when (2) sampling the visible microfacet distribution upper D Subscript omega Sub Subscript normal o Baseline left-parenthesis omega Subscript normal h Baseline right-parenthesis . Note in particular the variance improvement at surfaces that are oblique to the camera, where the visible distribution is quite different than the full distribution.

The implementation of the MicrofacetDistribution::Pdf() method now follows directly; it’s just a matter of returning the density from the selected sampling distribution.

<<MicrofacetDistribution Method Definitions>>+= 
Float MicrofacetDistribution::Pdf(const Vector3f &wo, const Vector3f &wh) const { if (sampleVisibleArea) return D(wh) * G1(wo) * AbsDot(wo, wh) / AbsCosTheta(wo); else return D(wh) * AbsCosTheta(wh); }

Given the ability to sample from distributions of microfacet orientations, the MicrofacetReflection::Sample_f() method can now be implemented.

<<BxDF Method Definitions>>+=  
Spectrum MicrofacetReflection::Sample_f(const Vector3f &wo, Vector3f *wi, const Point2f &u, Float *pdf, BxDFType *sampledType) const { <<Sample microfacet orientation omega Subscript normal h and reflected direction omega Subscript normal i >> 
Vector3f wh = distribution->Sample_wh(wo, u); *wi = Reflect(wo, wh); if (!SameHemisphere(wo, *wi)) return Spectrum(0.f);
<<Compute PDF of wi for microfacet reflection>> 
*pdf = distribution->Pdf(wo, wh) / (4 * Dot(wo, wh));
return f(wo, *wi); }

The implementation first uses Sample_wh() to find a microfacet orientation and reflects the outgoing direction about the microfacet’s normal. If the reflected direction is in the opposite hemisphere from omega Subscript normal o , then its direction is under the surface and no light is reflected.

<<Sample microfacet orientation omega Subscript normal h and reflected direction omega Subscript normal i >>= 
Vector3f wh = distribution->Sample_wh(wo, u); *wi = Reflect(wo, wh); if (!SameHemisphere(wo, *wi)) return Spectrum(0.f);

There’s an important detail to take care of to compute the value of the PDF for the sampled direction. The microfacet distribution gives the distribution of normals around the half-angle vector, but the reflection integral is with respect to the incoming vector. These distributions are not the same, and we must convert the half-angle PDF to the incoming angle PDF. In other words, we must change from a density in terms of omega Subscript normal h to one in terms of omega Subscript normal i using the techniques introduced in Section 13.5. Doing so requires applying the adjustment for a change of variables normal d omega Subscript normal h slash normal d omega Subscript normal i .

A simple geometric construction gives the relationship between the two distributions. Consider the spherical coordinate system oriented about omega Subscript normal o (Figure 14.4). Measured with respect to omega Subscript normal h , the differential solid angles normal d omega Subscript normal i and normal d omega Subscript normal h are sine theta Subscript normal i Baseline normal d theta Subscript normal i Baseline normal d phi Subscript normal i and sine theta Subscript normal h Baseline normal d theta Subscript normal h Baseline normal d phi Subscript normal h , respectively; thus,

StartFraction normal d omega Subscript normal h Baseline Over normal d omega Subscript normal i Baseline EndFraction equals StartFraction sine theta Subscript normal h Baseline normal d theta Subscript normal h Baseline normal d phi Subscript normal h Baseline Over sine theta Subscript normal i Baseline normal d theta Subscript normal i Baseline normal d phi Subscript normal i Baseline EndFraction period

Figure 14.4: The adjustment for change of variable from sampling from the half-angle distribution to sampling from the incident direction distribution can be derived with an observation about the relative angles involved: theta Subscript normal i Baseline equals 2 theta Subscript normal h . (Note the difference in notation: for example, here we use theta Subscript normal h as the angle between omega Subscript normal o and omega Subscript normal h ; elsewhere in the book, theta Subscript normal h denotes the angle between omega Subscript normal h and the surface normal.)

Because omega Subscript normal i is computed by reflecting omega Subscript normal o about omega Subscript normal h , theta Subscript normal i Baseline equals 2 theta Subscript normal h . Furthermore, because phi Subscript normal i Baseline equals phi Subscript normal h , we can find the desired conversion factor:

StartLayout 1st Row 1st Column StartFraction normal d omega Subscript normal h Baseline Over normal d omega Subscript normal i Baseline EndFraction 2nd Column equals StartFraction sine theta Subscript normal h Baseline normal d theta Subscript normal h Baseline normal d phi Subscript normal h Baseline Over sine 2 theta Subscript normal h Baseline 2 normal d theta Subscript normal h Baseline normal d phi Subscript normal h Baseline EndFraction 2nd Row 1st Column Blank 2nd Column equals StartFraction sine theta Subscript normal h Baseline Over 4 cosine theta Subscript normal h Baseline sine theta Subscript normal h Baseline EndFraction 3rd Row 1st Column Blank 2nd Column equals StartFraction 1 Over 4 cosine theta Subscript normal h Baseline EndFraction 4th Row 1st Column Blank 2nd Column equals StartFraction 1 Over 4 left-parenthesis omega Subscript normal i Baseline dot omega Subscript normal h Baseline right-parenthesis EndFraction equals StartFraction 1 Over 4 left-parenthesis omega Subscript normal o Baseline dot omega Subscript normal h Baseline right-parenthesis EndFraction period EndLayout

Therefore, the PDF after transformation is

p left-parenthesis theta right-parenthesis equals StartFraction p Subscript normal h Baseline left-parenthesis theta Subscript h Baseline right-parenthesis Over 4 left-parenthesis omega Subscript normal o Baseline dot omega Subscript normal h Baseline right-parenthesis EndFraction period

<<Compute PDF of wi for microfacet reflection>>= 
*pdf = distribution->Pdf(wo, wh) / (4 * Dot(wo, wh));

The same computation is implemented in the MicrofacetReflection::Pdf() method.

<<BxDF Method Definitions>>+=  
Float MicrofacetReflection::Pdf(const Vector3f &wo, const Vector3f &wi) const { if (!SameHemisphere(wo, wi)) return 0; Vector3f wh = Normalize(wo + wi); return distribution->Pdf(wo, wh) / (4 * Dot(wo, wh)); }

The approach for transmission is analogous: given a sampled omega Subscript normal h microfacet orientation, the outgoing direction is refracted about that normal direction to get the sampled incident direction. In the case of total internal reflection, Refract() returns false, and a black SPD is returned.

<<BxDF Method Definitions>>+=  
Spectrum MicrofacetTransmission::Sample_f(const Vector3f &wo, Vector3f *wi, const Point2f &u, Float *pdf, BxDFType *sampledType) const { Vector3f wh = distribution->Sample_wh(wo, u); Float eta = CosTheta(wo) > 0 ? (etaA / etaB) : (etaB / etaA); if (!Refract(wo, (Normal3f)wh, eta, wi)) return 0; *pdf = Pdf(wo, *wi); return f(wo, *wi); }

The PDF for the transmitted direction also requires an adjustment for the change of variables. This value is stored in dwh_dwi. We won’t derive this term here; the “Further Reading” section at the end of the chapter has more information.

<<BxDF Method Definitions>>+=  
Float MicrofacetTransmission::Pdf(const Vector3f &wo, const Vector3f &wi) const { if (SameHemisphere(wo, wi)) return 0; <<Compute omega Subscript normal h from omega Subscript normal o and omega Subscript normal i for microfacet transmission>> 
Float eta = CosTheta(wo) > 0 ? (etaB / etaA) : (etaA / etaB); Vector3f wh = Normalize(wo + wi * eta);
<<Compute change of variables dwh_dwi for microfacet transmission>> 
Float sqrtDenom = Dot(wo, wh) + eta * Dot(wi, wh); Float dwh_dwi = std::abs((eta * eta * Dot(wi, wh)) / (sqrtDenom * sqrtDenom));
return distribution->Pdf(wo, wh) * dwh_dwi; }

In the transmissive case, the meaning of the half-angle vector omega Subscript normal h is generalized to denote the normal of the microfacet that is responsible for refracting light from omega Subscript normal i to omega Subscript normal o . This vector can be derived following the setting in Figure 8.11.

<<Compute omega Subscript normal h from omega Subscript normal o and omega Subscript normal i for microfacet transmission>>= 
Float eta = CosTheta(wo) > 0 ? (etaB / etaA) : (etaA / etaB); Vector3f wh = Normalize(wo + wi * eta);

14.1.2 FresnelBlend

The FresnelBlend class is a mixture of a diffuse and glossy term. A straightforward approach to sampling this BRDF is to sample from both a cosine-weighted distribution as well as the microfacet distribution. The implementation here chooses between the two with equal probability based on whether xi 1 is less than or greater than 0.5 . In both cases, it remaps xi 1 to cover the range left-bracket 0 comma 1 right-parenthesis after using it to make this decision. (Otherwise, values of xi 1 used for the cosine-weighted sampling would always be less than 0.5 , for example.) Using the sample xi 1 for two purposes in this manner slightly reduces the quality of the stratification of the left-parenthesis xi 1 comma xi 2 right-parenthesis values that are actually used for sampling directions.

<<BxDF Method Definitions>>+=  
Spectrum FresnelBlend::Sample_f(const Vector3f &wo, Vector3f *wi, const Point2f &uOrig, Float *pdf, BxDFType *sampledType) const { Point2f u = uOrig; if (u[0] < .5) { u[0] = 2 * u[0]; <<Cosine-sample the hemisphere, flipping the direction if necessary>> 
*wi = CosineSampleHemisphere(u); if (wo.z < 0) wi->z *= -1;
} else { u[0] = 2 * (u[0] - .5f); <<Sample microfacet orientation omega Subscript normal h and reflected direction omega Subscript normal i >> 
Vector3f wh = distribution->Sample_wh(wo, u); *wi = Reflect(wo, wh); if (!SameHemisphere(wo, *wi)) return Spectrum(0.f);
} *pdf = Pdf(wo, *wi); return f(wo, *wi); }

The PDF for this sampling strategy is simple; it is just an average of the two PDFs used.

<<BxDF Method Definitions>>+=  
Float FresnelBlend::Pdf(const Vector3f &wo, const Vector3f &wi) const { if (!SameHemisphere(wo, wi)) return 0; Vector3f wh = Normalize(wo + wi); Float pdf_wh = distribution->Pdf(wo, wh); return .5f * (AbsCosTheta(wi) * InvPi + pdf_wh / (4 * Dot(wo, wh))); }

14.1.3 Specular Reflection and Transmission

The Dirac delta distributions that were previously used to define the BRDF for specular reflection and the BTDF for specular transmission fit into this sampling framework well, as long as a few conventions are kept in mind when using their sampling and PDF functions.

Recall that the Dirac delta distribution is defined such that

delta left-parenthesis x right-parenthesis equals 0 for all x not-equals 0

and

integral Subscript negative normal infinity Superscript normal infinity Baseline delta left-parenthesis x right-parenthesis normal d x equals 1 period

Thus, it is a probability density function, where the PDF has a value of 0 for all x not-equals 0 . Generating a sample from such a distribution is trivial; there is only one possible value for it to take on. When thought of in this way, the implementations of Sample_f() for the SpecularReflection and SpecularTransmission BxDFs can be seen to fit naturally into the Monte Carlo sampling framework.

It is not as simple to determine which value should be returned for the value of the PDF. Strictly speaking, the delta distribution is not a true function but must be defined as the limit of another function—for example, one describing a box of unit area whose width approaches 0; see Chapter 5 of Bracewell (2000) for further discussion and references. Thought of in this way, the value of delta left-parenthesis 0 right-parenthesis tends toward infinity. Certainly, returning an infinite or very large value for the PDF is not going to lead to correct results from the renderer.

However, recall that BSDFs defined with delta components also have these delta components in their f Subscript normal r functions, a detail that was glossed over when we returned values from their Sample_f() methods in Chapter 8. Thus, the Monte Carlo estimator for the scattering equation with such a BSDF is written

StartLayout 1st Row 1st Column Blank 2nd Column StartFraction 1 Over upper N EndFraction sigma-summation Underscript i Overscript upper N Endscripts StartFraction f Subscript normal r Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline comma omega Subscript i Baseline right-parenthesis upper L Subscript normal i Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript i Baseline right-parenthesis StartAbsoluteValue cosine theta Subscript i Baseline EndAbsoluteValue Over p left-parenthesis omega Subscript i Baseline right-parenthesis EndFraction equals StartFraction 1 Over upper N EndFraction sigma-summation Underscript i Overscript upper N Endscripts StartStartFraction rho Subscript normal h normal d Baseline left-parenthesis omega Subscript normal o Baseline right-parenthesis StartFraction delta left-parenthesis omega Subscript Baseline minus omega Subscript i Baseline right-parenthesis Over StartAbsoluteValue cosine theta Subscript i Baseline EndAbsoluteValue EndFraction upper L Subscript normal i Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript i Baseline right-parenthesis StartAbsoluteValue cosine theta Subscript i Baseline EndAbsoluteValue OverOver p left-parenthesis omega Subscript i Baseline right-parenthesis EndEndFraction comma EndLayout

where rho Subscript normal h normal d Baseline left-parenthesis omega Subscript normal o Baseline right-parenthesis is the hemispherical–directional reflectance and omega Subscript is the direction for perfect specular reflection or transmission.

Because the PDF p left-parenthesis omega Subscript i Baseline right-parenthesis has a delta term as well, p left-parenthesis omega Subscript i Baseline right-parenthesis equals delta left-parenthesis omega Subscript Baseline minus omega Subscript i Baseline right-parenthesis , the two delta distributions cancel out, and the estimator is

rho Subscript normal h normal d Baseline left-parenthesis omega Subscript normal o Baseline right-parenthesis upper L Subscript normal i Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis comma

exactly the quantity computed by the Whitted integrator, for example.

Therefore, the implementations here return a constant value of 1 for the PDF for specular reflection and transmission when sampled using Sample_f(), with the convention that for specular BxDFs there is an implied delta distribution in the PDF value that is expected to cancel out with the implied delta distribution in the value of the BSDF when the estimator is evaluated. The respective Pdf() methods therefore return 0 for all directions, since there is zero probability that another sampling method will randomly find the direction from a delta distribution.

<<SpecularReflection Public Methods>>+= 
Float Pdf(const Vector3f &wo, const Vector3f &wi) const { return 0; }

<<SpecularTransmission Public Methods>>+= 
Float Pdf(const Vector3f &wo, const Vector3f &wi) const { return 0; }

There is a potential pitfall with this convention: when multiple importance sampling is used to compute weights, PDF values that include these implicit delta distributions can’t be freely mixed with regular PDF values. This isn’t a problem in practice, since there’s no reason to apply MIS when there’s a delta distribution in the integrand. The light transport routines in this and the next two chapters have appropriate logic to avoid this error.

The FresnelSpecular class encapsulates both specular reflection and transmission, with the relative contributions modulated by a dielectric Fresnel term. By combining these two together, it’s able to use the value of the Fresnel term for the outgoing direction omega Subscript normal o to determine which component to sample—for example, for glancing angles where reflection is high, it’s much more likely to return a reflected direction than a transmitted direction. This approach improves Monte Carlo efficiency when rendering scenes with these sorts of surfaces, since the rays that are sampled will tend to have larger contributions to the final result.

<<BxDF Method Definitions>>+=  
Spectrum FresnelSpecular::Sample_f(const Vector3f &wo, Vector3f *wi, const Point2f &u, Float *pdf, BxDFType *sampledType) const { Float F = FrDielectric(CosTheta(wo), etaA, etaB); if (u[0] < F) { <<Compute specular reflection for FresnelSpecular>> 
<<Compute perfect specular reflection direction>> 
*wi = Vector3f(-wo.x, -wo.y, wo.z);
if (sampledType) *sampledType = BxDFType(BSDF_SPECULAR | BSDF_REFLECTION); *pdf = F; return F * R / AbsCosTheta(*wi);
} else { <<Compute specular transmission for FresnelSpecular>> 
<<Figure out which eta is incident and which is transmitted>> 
bool entering = CosTheta(wo) > 0; Float etaI = entering ? etaA : etaB; Float etaT = entering ? etaB : etaA;
<<Compute ray direction for specular transmission>> 
if (!Refract(wo, Faceforward(Normal3f(0, 0, 1), wo), etaI / etaT, wi)) return 0;
Spectrum ft = T * (1 - F); <<Account for non-symmetry with transmission to different medium>> 
if (mode == TransportMode::Radiance) ft *= (etaI * etaI) / (etaT * etaT);
if (sampledType) *sampledType = BxDFType(BSDF_SPECULAR | BSDF_TRANSMISSION); *pdf = 1 - F; return ft / AbsCosTheta(*wi);
} }

Specular reflection is chosen with probability equal to F; given that choice, the computations performed are the same as those in SpecularReflection.

<<Compute specular reflection for FresnelSpecular>>= 
<<Compute perfect specular reflection direction>> 
*wi = Vector3f(-wo.x, -wo.y, wo.z);
if (sampledType) *sampledType = BxDFType(BSDF_SPECULAR | BSDF_REFLECTION); *pdf = F; return F * R / AbsCosTheta(*wi);

Otherwise, with probability 1-F, specular transmission is selected.

<<Compute specular transmission for FresnelSpecular>>= 
<<Figure out which eta is incident and which is transmitted>> 
bool entering = CosTheta(wo) > 0; Float etaI = entering ? etaA : etaB; Float etaT = entering ? etaB : etaA;
<<Compute ray direction for specular transmission>> 
if (!Refract(wo, Faceforward(Normal3f(0, 0, 1), wo), etaI / etaT, wi)) return 0;
Spectrum ft = T * (1 - F); <<Account for non-symmetry with transmission to different medium>> 
if (mode == TransportMode::Radiance) ft *= (etaI * etaI) / (etaT * etaT);
if (sampledType) *sampledType = BxDFType(BSDF_SPECULAR | BSDF_TRANSMISSION); *pdf = 1 - F; return ft / AbsCosTheta(*wi);

14.1.4 Fourier BSDF

In addition to being able to compactly and accurately represent a variety of measured and synthetic BSDFs, the representation used by the FourierBSDF (Section 8.6) also admits a fairly efficient exact importance sampling scheme. Figure 14.5 compares the result of using this sampling approach to using uniform hemispherical sampling for the FourierBSDF.

Figure 14.5: The Effect of Accurately Sampling the FourierBSDF. Reflection from both objects is modeled using the FourierBSDF, rendered using 32 samples per pixel. (a) Uniform hemisphere sampling to compute reflection. (b) The exact sampling scheme implemented in FourierBSDF::Sample_f(). Variance is much lower and overall rendering time increased by only 20%.

Recall the BSDF representation from Equation (8.21), which was

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 comma

where the function a Subscript k was discretized over the incident and outgoing zenith angle cosines left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline right-parenthesis element-of StartSet mu 0 comma ellipsis comma mu Subscript n minus 1 Baseline EndSet squared with endpoints mu 0 equals negative 1 and mu Subscript n minus 1 Baseline equals 1 . An even Fourier expansion with real coefficients was used to model the dependence on the azimuth angle difference parameter phi equals phi Subscript normal i Baseline minus phi Subscript normal o .

The task now is to first sample mu Subscript normal i given mu Subscript normal o and then sample the angle phi Subscript normal i relative to phi Subscript normal o . A helpful property of the order 0 Fourier coefficients greatly simplifies both of these steps. The even Fourier basis functions form an orthogonal basis of the vector space of square-integrable even functions—this means that the basis coefficients of any function g satisfying these criteria can be found using an inner product between g and the individual basis functions analogous to orthogonal projections of vectors on Euclidean vector spaces. Here, we are dealing with continuous functions on left-bracket 0 comma pi right-bracket , where a suitable inner product can be defined as

mathematical left-angle g comma h mathematical right-angle equals StartFraction 1 Over pi EndFraction integral Subscript 0 Superscript pi Baseline g left-parenthesis phi right-parenthesis h left-parenthesis phi right-parenthesis normal d phi Subscript Baseline period

The Fourier basis function associated with order 0 is simply the unit constant; hence the coefficients in a 0 relate to the cosine-weighted BSDF as

a 0 left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline right-parenthesis equals StartFraction 1 Over pi EndFraction integral Subscript 0 Superscript pi Baseline f left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline comma phi right-parenthesis StartAbsoluteValue mu Subscript normal i Baseline EndAbsoluteValue normal d phi Subscript Baseline period

This quantity turns out to be very helpful in constructing a method for importance sampling the BSDF: disregarding normalization factors, this average over phi can be interpreted as the marginal distribution of the cosine-weighted BSDF with respect to pairs of mu angle cosines (Section 13.6 discussed marginal density functions).

It will be useful to be able to efficiently access these order 0 coefficients without the indirections that would normally be necessary given the layout of FourierBSDFTable::a. We therefore keep an additional copy of these values in an array of size monospace n monospace upper M monospace u times monospace n monospace upper M monospace u in FourierBSDFTable::a0. This array is initialized by copying the corresponding entries from FourierBSDFTable::a at scene loading time in the FourierBSDFTable::Read() method.

<<FourierBSDFTable Public Data>>+=  
Float *a0;

With a marginal distribution at hand, we are now able to split the sampling operation into two lower dimensional steps: first, we use the a 0 coefficients to sample mu Subscript normal i given mu Subscript normal o . Second, with left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline right-parenthesis known, we can interpolate the Fourier coefficients that specify the BSDF’s dependence on the remaining azimuth difference angle parameter phi and sample from their distribution. These operations are all implemented as smooth mappings that preserve the properties of structured point sets, such as Sobol prime or Halton sequences. Given these angles, the last step is to compute the corresponding direction and value of the BSDF.

<<BxDF Method Definitions>>+=  
Spectrum FourierBSDF::Sample_f(const Vector3f &wo, Vector3f *wi, const Point2f &u, Float *pdf, BxDFType *sampledType) const { <<Sample zenith angle component for FourierBSDF>> 
Float muO = CosTheta(wo); Float pdfMu; Float muI = SampleCatmullRom2D(bsdfTable.nMu, bsdfTable.nMu, bsdfTable.mu, bsdfTable.mu, bsdfTable.a0, bsdfTable.cdf, muO, u[1], nullptr, &pdfMu);
<<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]; }
} }
<<Importance sample the luminance Fourier expansion>> 
Float phi, pdfPhi; Float Y = SampleFourier(ak, bsdfTable.recip, mMax, u[0], &pdfPhi, &phi); *pdf = std::max((Float)0, pdfPhi * pdfMu);
<<Compute the scattered direction for FourierBSDF>> 
Float sin2ThetaI = std::max((Float)0, 1 - muI * muI); Float norm = std::sqrt(sin2ThetaI / Sin2Theta(wo)); Float sinPhi = std::sin(phi), cosPhi = std::cos(phi); *wi = -Vector3f(norm * (cosPhi * wo.x - sinPhi * wo.y), norm * (sinPhi * wo.x + cosPhi * wo.y), muI);
<<Evaluate remaining Fourier expansions for angle phi >> 
Float scale = muI != 0 ? (1/std::abs(muI)) : (Float) 0; 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); 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();
}

Sampling the zenith angle is implemented using SampleCatmullRom2D(), which will be defined in a few pages. This helper function operates in stages: after first mapping a uniform variate to one of the spline segments, it then samples a specific position within the segment. To select a segment, the function requires an array of precomputed CDFs

upper I Subscript normal i comma normal o Baseline equals integral Subscript negative 1 Superscript mu Subscript normal i Baseline Baseline a 0 left-parenthesis mu prime comma mu Subscript normal o Baseline right-parenthesis normal d mu Superscript prime Baseline comma

where 0 less-than-or-equal-to normal i comma normal o less-than n . Each column of this matrix stores a discrete CDF over mu Subscript normal i for a different (fixed) value of mu Subscript normal o . The above integral is computed directly from the spline interpolant and can thus be used to efficiently select spline segments proportional to their definite integral.

<<FourierBSDFTable Public Data>>+=  
Float *cdf;

In the case of the FourierBSDF, this cdf array is already part of the input file, and we need not be concerned with its generation.

<<Sample zenith angle component for FourierBSDF>>= 
Float muO = CosTheta(wo); Float pdfMu; Float muI = SampleCatmullRom2D(bsdfTable.nMu, bsdfTable.nMu, bsdfTable.mu, bsdfTable.mu, bsdfTable.a0, bsdfTable.cdf, muO, u[1], nullptr, &pdfMu);

After SampleCatmullRom2D() returns, muI records the sampled incident zenith angle cosine, and pdfMu contains the associated probability density in the same domain.

We can now interpolate the nearby Fourier coefficients, reusing the fragment <<Compute Fourier coefficients a Subscript k for left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline right-parenthesis >> from FourierBSDF::f() in Section 8.6. Given the coefficients ak, sampling of the azimuth difference angle is also implemented as a separate function SampleFourier(), also to be defined in a few pages. This function returns the sampled difference angle phi, the value Y of the luminance Fourier expansion evaluated at phi, and the sample probability pdfPhi per radian. The final sample probability per unit solid angle is the product of the azimuthal and zenith angle cosine PDFs. (As with values computed via Fourier series in Section 8.6, negative values must be clamped to 0.)

<<Importance sample the luminance Fourier expansion>>= 
Float phi, pdfPhi; Float Y = SampleFourier(ak, bsdfTable.recip, mMax, u[0], &pdfPhi, &phi); *pdf = std::max((Float)0, pdfPhi * pdfMu);

SampleFourier() takes an additional input array recip containing precomputed integer reciprocals 1 slash i for all mMax Fourier orders. These reciprocals are frequently accessed within the function—precomputing them is an optimization to avoid costly arithmetic to generate them over and over again, causing pipeline stalls due to the high latency of division operations on current processor architectures. This reciprocal array is initialized in FourierBSDFTable::Read().

<<FourierBSDFTable Public Data>>+= 
Float *recip;

We now have the values mu Subscript normal i Baseline equals cosine theta Subscript normal i and phi . The sampled incident direction’s x y coordinates are given by rotating the x y components of omega Subscript normal o by an angle phi about the surface normal, and its z component is given by mu Subscript normal i (using spherical coordinates).

There are two details to note in the computation of the direction omega Subscript normal i . First, the x y components are scaled by a factor sine theta Subscript normal i slash sine theta Subscript normal o , which ensures that the resulting vector is normalized. Second, the computed direction is negated before being assigned to *wi; this follows the coordinate system convention for the FourierBSDF that was described in Section 8.6.

<<Compute the scattered direction for FourierBSDF>>= 
Float sin2ThetaI = std::max((Float)0, 1 - muI * muI); Float norm = std::sqrt(sin2ThetaI / Sin2Theta(wo)); Float sinPhi = std::sin(phi), cosPhi = std::cos(phi); *wi = -Vector3f(norm * (cosPhi * wo.x - sinPhi * wo.y), norm * (sinPhi * wo.x + cosPhi * wo.y), muI);

The fragment <<Evaluate remaining Fourier expansions for angle phi >> is identical to <<Evaluate Fourier expansion for angle phi >> defined in Section 8.6 except that doesn’t evaluate the luminance channel, which was already done by SampleFourier() above.

The FourierBSDF::Pdf() method returns the solid angle density for the preceding sampling method. Since this method produces samples that are exactly distributed according to the product f left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline comma phi right-parenthesis StartAbsoluteValue mu Subscript normal i Baseline EndAbsoluteValue , we could simply copy the implementation of FourierBSDF::f() except for the division that cancels StartAbsoluteValue mu Subscript normal i Baseline EndAbsoluteValue . However, doing so would underestimate the probability when the BSDF doesn’t reflect all of the incident illumination.

To correct for this, we scale the unnormalized f left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline comma phi right-parenthesis StartAbsoluteValue mu Subscript normal i Baseline EndAbsoluteValue by a suitable normalization factor rho Superscript negative 1 to ensure that the product integrates to 1:

integral Subscript 0 Superscript 2 pi Baseline integral Subscript 0 Superscript pi Baseline StartFraction 1 Over rho EndFraction f left-parenthesis cosine theta Subscript normal i Baseline comma cosine theta Subscript o Baseline comma phi right-parenthesis StartAbsoluteValue cosine theta Subscript normal i Baseline EndAbsoluteValue sine theta Subscript normal i Baseline normal d theta Subscript normal i Baseline normal d phi Subscript Baseline equals 1 period

Note that the outgoing zenith angle cosine cosine theta Subscript normal o Baseline equals mu Subscript normal o was left unspecified in the above equation. In general, the normalization factor rho is not constant and, instead, it depends on the current value of mu Subscript normal o . rho left-parenthesis mu Subscript normal o Baseline right-parenthesis has a simple interpretation: it is the hemispherical–directional reflectance of a surface that is illuminated from the zenith angle cosine Superscript negative 1 Baseline mu Subscript normal o .

Suppose briefly that mu Subscript normal o happens to be part of the discretized set of zenith angle cosines mu 0 comma ellipsis comma mu Subscript n minus 1 Baseline stored in the array FourierBSDFTable::mu. Then

StartLayout 1st Row 1st Column rho left-parenthesis mu Subscript normal o Baseline right-parenthesis 2nd Column equals integral Subscript 0 Superscript 2 pi Baseline integral Subscript 0 Superscript pi Baseline f left-parenthesis cosine theta Subscript normal i Baseline comma cosine theta Subscript normal o Baseline comma phi right-parenthesis StartAbsoluteValue cosine theta Subscript normal i Baseline EndAbsoluteValue sine theta Subscript normal i Baseline normal d theta Subscript normal i Baseline normal d phi Subscript Baseline 2nd Row 1st Column Blank 2nd Column equals integral Subscript 0 Superscript 2 pi Baseline integral Subscript negative 1 Superscript 1 Baseline f left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline comma phi right-parenthesis StartAbsoluteValue mu Subscript normal i Baseline EndAbsoluteValue normal d mu Subscript normal i Baseline normal d phi Subscript Baseline 3rd Row 1st Column Blank 2nd Column equals 2 pi integral Subscript negative 1 Superscript 1 Baseline left-bracket StartFraction 1 Over pi EndFraction integral Subscript 0 Superscript pi Baseline f left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline comma phi right-parenthesis StartAbsoluteValue mu Subscript normal i Baseline EndAbsoluteValue normal d phi Subscript Baseline right-bracket normal d mu Subscript normal i Baseline 4th Row 1st Column Blank 2nd Column equals 2 pi integral Subscript negative 1 Superscript 1 Baseline a 0 left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline right-parenthesis normal d mu Subscript normal i Baseline 5th Row 1st Column Blank 2nd Column equals 2 pi upper I Subscript n minus 1 comma normal o Baseline comma EndLayout

where upper I was defined in Equation (14.3). In other words, the needed normalization factor is readily available in the FourierBSDFTable::cdf array. For intermediate values of mu Subscript normal o , we can simply interpolate the neighboring four entries of upper I Subscript n minus 1 comma i using the usual spline interpolation scheme—the linearity of this interpolation coupled with the linearity of the analytic integrals in (14.4) ensures a result that is consistent with FourierBSDF::f().

<<BxDF Method Definitions>>+=  
Float FourierBSDF::Pdf(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 luminance Fourier coefficients a Subscript k for left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline right-parenthesis >> 
int offsetI, offsetO; Float weightsI[4], weightsO[4]; if (!bsdfTable.GetWeightsAndOffset(muI, &offsetI, weightsI) || !bsdfTable.GetWeightsAndOffset(muO, &offsetO, weightsO)) return 0; Float *ak = ALLOCA(Float, bsdfTable.mMax * bsdfTable.nChannels); memset(ak, 0, bsdfTable.mMax * bsdfTable.nChannels * sizeof(Float)); int mMax = 0; for (int o = 0; o < 4; ++o) { for (int i = 0; i < 4; ++i) { Float weight = weightsI[i] * weightsO[o]; if (weight == 0) continue; int order; const Float *coeffs = bsdfTable.GetAk(offsetI + i, offsetO + o, &order); mMax = std::max(mMax, order); for (int k = 0; k < order; ++k) ak[k] += *coeffs++ * weight; } }
<<Evaluate probability of sampling wi>> 
Float rho = 0; for (int o = 0; o < 4; ++o) { if (weightsO[o] == 0) continue; rho += weightsO[o] * bsdfTable.cdf[(offsetO + o) * bsdfTable.nMu + bsdfTable.nMu - 1] * (2 * Pi); } Float Y = Fourier(ak, mMax, cosPhi); return (rho > 0 && Y > 0) ? (Y / rho) : 0;
}

We won’t include the second fragment here—it is almost identical to <<Compute Fourier coefficients a Subscript k for left-parenthesis mu Subscript normal i Baseline comma mu Subscript normal o Baseline right-parenthesis >>, the only difference being that it only interpolates the luminance coefficients (samples are generated proportional to luminance; hence the other two channels are not relevant here).

The last fragment interpolates the directional albedo from Equation (14.4) and uses it to correct the result of Fourier() for absorption.

<<Evaluate probability of sampling wi>>= 
Float rho = 0; for (int o = 0; o < 4; ++o) { if (weightsO[o] == 0) continue; rho += weightsO[o] * bsdfTable.cdf[(offsetO + o) * bsdfTable.nMu + bsdfTable.nMu - 1] * (2 * Pi); } Float Y = Fourier(ak, mMax, cosPhi); return (rho > 0 && Y > 0) ? (Y / rho) : 0;

Sampling 1D Spline Interpolants

Before defining the SampleCatmullRom2D() function used in the previously discussed FourierBSDF::Sample_f() method, we’ll first focus on a simpler 1D case: suppose that a function f was evaluated at n positions x 0 comma ellipsis comma x Subscript n minus 1 Baseline , resulting in a piecewise cubic Catmull-Rom spline interpolant ModifyingAbove f With caret with n minus 1 spline segments ModifyingAbove f With caret Subscript i Baseline left-parenthesis x right-parenthesis left-parenthesis i equals 0 comma ellipsis comma n minus 2 right-parenthesis . Given a precomputed discrete CDF over these spline segments defined as

upper F Subscript i Baseline equals StartLayout Enlarged left-brace 1st Row 1st Column 0 comma 2nd Column if i equals 0 comma 2nd Row 1st Column sigma-summation Underscript k equals 0 Overscript i minus 1 Endscripts integral Subscript x Subscript k Baseline Superscript x Subscript k plus 1 Baseline Baseline StartFraction 1 Over c EndFraction ModifyingAbove f With caret Subscript k Baseline left-parenthesis x Superscript prime Baseline right-parenthesis normal d x Superscript prime Baseline comma 2nd Column if i greater-than 0 comma EndLayout

where c is the normalization term,

c equals integral Subscript x 0 Superscript x Subscript n minus 1 Baseline Baseline ModifyingAbove f With caret left-parenthesis x right-parenthesis normal d x comma

a straightforward way of importance sampling ModifyingAbove f With caret in two steps using the inversion method entails first finding the interval i such that

upper F Subscript i Baseline less-than-or-equal-to xi 1 less-than upper F Subscript i plus 1 Baseline comma

where xi 1 is a random variate on the interval left-bracket 0 comma 1 right-parenthesis , and then sampling an x prime value in the i th interval. Since the values upper F Subscript i are monotonically increasing, the interval can be found using an efficient binary search.

In the following, we won’t actually normalize the upper F Subscript i values, effectively setting c equals 1 . We can equivalently sample i by first multiplying the random variable xi 1 by the last upper F Subscript i entry, upper F Subscript n minus 1 , which is the total integral over all spline segments and is thus equal to c . Thus, the binary search looks for

upper F Subscript i Baseline less-than-or-equal-to xi 1 upper F Subscript n minus 1 Baseline less-than-or-equal-to upper F Subscript i plus 1 Baseline comma

Having selected a segment i , we can offset and re-scale xi 1 to obtain a second uniform variate in left-bracket 0 comma 1 right-parenthesis :

xi 2 equals StartFraction xi 1 upper F Subscript n minus 1 Baseline minus upper F Subscript i Baseline Over upper F Subscript i plus 1 Baseline minus upper F Subscript i Baseline EndFraction period

We then use xi 2 to sample a position x within the interval left-bracket x Subscript i Baseline comma x Subscript i plus 1 Baseline right-bracket using that segment’s integral,

ModifyingAbove upper F With caret Subscript i Baseline left-parenthesis x right-parenthesis equals integral Subscript x Subscript i Baseline Superscript x Baseline ModifyingAbove f With caret Subscript i Baseline left-parenthesis x Superscript prime Baseline right-parenthesis normal d x Superscript prime Baseline comma

where again we won’t compute a properly normalized CDF but will instead multiply xi 2 by ModifyingAbove upper F With caret Subscript i Baseline left-parenthesis x Subscript i plus 1 Baseline right-parenthesis rather than normalizing ModifyingAbove upper F With caret Subscript i . We must then compute

StartLayout 1st Row 1st Column x 2nd Column equals ModifyingAbove upper F With caret Subscript i Superscript negative 1 Baseline left-parenthesis ModifyingAbove upper F With caret Subscript i Baseline left-parenthesis x Subscript i plus 1 Baseline right-parenthesis xi 2 right-parenthesis equals ModifyingAbove upper F With caret Subscript i Superscript negative 1 Baseline left-parenthesis left-parenthesis upper F Subscript i italic plus italic 1 Baseline minus upper F Subscript i Baseline right-parenthesis StartFraction xi 1 upper F Subscript n minus 1 Baseline minus upper F Subscript i Baseline Over upper F Subscript i plus 1 Baseline minus upper F Subscript i Baseline EndFraction right-parenthesis 2nd Row 1st Column Blank 2nd Column equals ModifyingAbove upper F With caret Subscript i Superscript negative 1 Baseline left-parenthesis xi 1 upper F Subscript n minus 1 Baseline minus upper F Subscript i Baseline right-parenthesis period EndLayout

This approach is implemented in SampleCatmullRom(), which takes the following inputs: the number of samples n; x contains the locations x 0 comma ellipsis comma x Subscript n minus 1 Baseline where the original function f was evaluated; f stores the value of the function at each point x Subscript i ; u is used to pass the uniform variate xi ; and integrated upper F Subscript i values must be provided via the F parameter—these values can be precomputed with IntegrateCatmullRom() when necessary. fval and pdf are used to return the function value and associated PDF value.

<<Spline Interpolation Definitions>>+=  
Float SampleCatmullRom(int n, const Float *x, const Float *f, const Float *F, Float u, Float *fval, Float *pdf) { <<Map u to a spline interval by inverting F>> 
u *= F[n - 1]; int i = FindInterval(n, [&](int i) { return F[i] <= u; });
<<Look up x Subscript i and function values of spline segment i>> 
Float x0 = x[i], x1 = x[i + 1]; Float f0 = f[i], f1 = f[i + 1]; Float width = x1 - x0;
<<Approximate derivatives using finite differences>> 
Float d0, d1; if (i > 0) d0 = width * (f1 - f[i - 1]) / (x1 - x[i - 1]); else d0 = f1 - f0; if (i + 2 < n) d1 = width * (f[i + 2] - f0) / (x[i + 2] - x0); else d1 = f1 - f0;
<<Re-scale u for continous spline sampling step>> 
u = (u - F[i]) / width;
<<Invert definite integral over spline segment and return solution>> 
<<Set initial guess for t by importance sampling a linear interpolant>> 
Float t; if (f0 != f1) t = (f0 - std::sqrt( std::max((Float)0, f0 * f0 + 2 * u * (f1 - f0)))) / (f0 - f1); else t = u / f0;
Float a = 0, b = 1, Fhat, fhat; while (true) { <<Fall back to a bisection step when t is out of bounds>> 
if (!(t >= a && t <= b)) t = 0.5f * (a + b);
<<Evaluate target function and its derivative in Horner form>> 
Fhat = t * (f0 + t * (.5f * d0 + t * ((1.f/3.f) * (-2 * d0 - d1) + f1 - f0 + t * (.25f * (d0 + d1) + .5f * (f0 - f1))))); fhat = f0 + t * (d0 + t * (-2 * d0 - d1 + 3 * (f1 - f0) + t * (d0 + d1 + 2 * (f0 - f1))));
<<Stop the iteration if converged>> 
if (std::abs(Fhat - u) < 1e-6f || b - a < 1e-6f) break;
<<Update bisection bounds using updated t>> 
if (Fhat - u < 0) a = t; else b = t;
<<Perform a Newton step>> 
t -= (Fhat - u) / fhat;
} <<Return the sample position and function value>> 
if (fval) *fval = fhat; if (pdf) *pdf = fhat / F[n - 1]; return x0 + width * t;
}

The function begins by scaling the uniform variate u by the last entry of F following Equation (14.6). Following this, u is mapped to a spline interval via the FindInterval() helper function, which returns the last interval satisfying F[i] <= u while clamping to the array bounds in case of rounding errors.

<<Map u to a spline interval by inverting F>>= 
u *= F[n - 1]; int i = FindInterval(n, [&](int i) { return F[i] <= u; });

The next fragment fetches the associated function values and node positions from f and x; the variable width contains the segment length.

<<Look up x Subscript i and function values of spline segment i>>= 
Float x0 = x[i], x1 = x[i + 1]; Float f0 = f[i], f1 = f[i + 1]; Float width = x1 - x0;

Recall that Catmull-Rom splines require an approximation of the first derivative of the function (Section 8.6.1) at the segment endpoints. Depending on i, this derivative is computed using forward, backward, or central finite difference approximations.

<<Approximate derivatives using finite differences>>= 
Float d0, d1; if (i > 0) d0 = width * (f1 - f[i - 1]) / (x1 - x[i - 1]); else d0 = f1 - f0; if (i + 2 < n) d1 = width * (f[i + 2] - f0) / (x[i + 2] - x0); else d1 = f1 - f0;

The remainder of the function then has to find the inverse of the continuous cumulative distribution function from Equation (14.8):

upper F Subscript i Superscript negative 1 Baseline left-parenthesis xi 1 upper F Subscript n minus 1 Baseline minus upper F Subscript i Baseline right-parenthesis period

Since the scaling by upper F Subscript n minus 1 was already applied in the first fragment, we need only subtract  upper F Subscript i .

The actual inversion is done in <<Invert definite integral over spline segment and return solution>>, whose discussion we postpone for the following discussion of the 2D cases. The internals of this inversion operate on a scaled and shifted spline segment defined on the interval left-bracket 0 comma 1 right-bracket , which requires an additional scaling by the associated change of variable factor equal to the reciprocal of width.

<<Re-scale u for continous spline sampling step>>= 
u = (u - F[i]) / width;

Sampling 2D Spline Interpolants

The main use cases of spline interpolants in pbrt actually importance sample 2D functions f left-parenthesis alpha comma x right-parenthesis , where alpha is considered a fixed parameter for the purpose of sampling (e.g., the albedo of the underlying material or the outgoing zenith angle cosine mu Subscript normal o in the case of the FourierBSDF). This case is handled by SampleCatmullRom2D().

<<Spline Interpolation Definitions>>+=  
Float SampleCatmullRom2D(int size1, int size2, const Float *nodes1, const Float *nodes2, const Float *values, const Float *cdf, Float alpha, Float u, Float *fval, Float *pdf) { <<Determine offset and coefficients for the alpha parameter>> 
int offset; Float weights[4]; if (!CatmullRomWeights(size1, nodes1, alpha, &offset, weights)) return 0;
<<Define a lambda function to interpolate table entries>> 
auto interpolate = [&](const Float *array, int idx) { Float value = 0; for (int i = 0; i < 4; ++i) if (weights[i] != 0) value += array[(offset + i) * size2 + idx] * weights[i]; return value; };
<<Map u to a spline interval by inverting the interpolated cdf>> 
Float maximum = interpolate(cdf, size2-1); u *= maximum; int idx = FindInterval(size2, [&](int i) { return interpolate(cdf, i) <= u; });
<<Look up node positions and interpolated function values>> 
Float f0 = interpolate(values, idx), f1 = interpolate(values, idx+1); Float x0 = nodes2[idx], x1 = nodes2[idx+1]; Float width = x1 - x0; Float d0, d1;
<<Re-scale u using the interpolated cdf>> 
u = (u - interpolate(cdf, idx)) / width;
<<Approximate derivatives using finite differences of the interpolant>> 
if (idx > 0) d0 = width * (f1 - interpolate(values, idx - 1)) / (x1 - nodes2[idx - 1]); else d0 = f1 - f0; if (idx + 2 < size2) d1 = width * (interpolate(values, idx+2) - f0) / (nodes2[idx + 2] - x0); else d1 = f1 - f0;
<<Invert definite integral over spline segment and return solution>> 
<<Set initial guess for t by importance sampling a linear interpolant>> 
Float t; if (f0 != f1) t = (f0 - std::sqrt( std::max((Float)0, f0 * f0 + 2 * u * (f1 - f0)))) / (f0 - f1); else t = u / f0;
Float a = 0, b = 1, Fhat, fhat; while (true) { <<Fall back to a bisection step when t is out of bounds>> 
if (!(t >= a && t <= b)) t = 0.5f * (a + b);
<<Evaluate target function and its derivative in Horner form>> 
Fhat = t * (f0 + t * (.5f * d0 + t * ((1.f/3.f) * (-2 * d0 - d1) + f1 - f0 + t * (.25f * (d0 + d1) + .5f * (f0 - f1))))); fhat = f0 + t * (d0 + t * (-2 * d0 - d1 + 3 * (f1 - f0) + t * (d0 + d1 + 2 * (f0 - f1))));
<<Stop the iteration if converged>> 
if (std::abs(Fhat - u) < 1e-6f || b - a < 1e-6f) break;
<<Update bisection bounds using updated t>> 
if (Fhat - u < 0) a = t; else b = t;
<<Perform a Newton step>> 
t -= (Fhat - u) / fhat;
} <<Return the sample position and function value>> 
if (fval) *fval = fhat; if (pdf) *pdf = fhat / F[n - 1]; return x0 + width * t;
}

The parameters size1, size2, nodes1, and nodes2 specify separate discretizations for each dimension. The values argument supplies a matrix of function values in row-major order, with rows corresponding to sets of samples that have the same position along the first dimension. The function uses the parameter alpha to choose a slice in the first dimension that is then importance sampled along the second dimension. The parameter cdf supplies a matrix of discrete CDFs, where each row was obtained by running IntegrateCatmullRom() on the corresponding row of values.

The first fragment of SampleCatmullRom2D() calls CatmullRomWeights() to select four adjacent rows of the values array along with interpolation weights.

<<Determine offset and coefficients for the alpha parameter>>= 
int offset; Float weights[4]; if (!CatmullRomWeights(size1, nodes1, alpha, &offset, weights)) return 0;

To proceed, we could now simply interpolate the selected rows of values and cdf and finish by calling the 1D sampling function SampleCatmullRom(). However, only a few entries of values and cdf are truly needed to generate a sample in practice, making such an approach unnecessarily slow. Instead, we define a C++11 lambda function that interpolates entries of these arrays on demand:

<<Define a lambda function to interpolate table entries>>= 
auto interpolate = [&](const Float *array, int idx) { Float value = 0; for (int i = 0; i < 4; ++i) if (weights[i] != 0) value += array[(offset + i) * size2 + idx] * weights[i]; return value; };

The rest of the function is identical to SampleCatmullRom() except that every access to x[i] is replaced by interpolate(values, i) and similarly for cdf. For brevity, this code is omitted in the book.

We now return to the inversion of the integral in Equation (14.8), which we glossed over. Recall that ModifyingAbove upper F With caret Subscript i was defined as an integral over the cubic spline segment ModifyingAbove f With caret Subscript i , making it a quartic polynomial. It is possible but burdensome to invert this function analytically. We prefer a numerical approach that is facilitated by a useful pair of properties:

  1. The function ModifyingAbove upper F With caret Subscript i is the definite integral over the (assumed nonnegative) interpolant ModifyingAbove f With caret Subscript i , and so it increases monotonically.
  2. The interval left-bracket x Subscript i Baseline comma x Subscript i plus 1 Baseline right-bracket selected by the function FindInterval() contains exactly one solution to Equation (14.8).

In this case, the interval left-bracket x Subscript i Baseline comma x Subscript i plus 1 Baseline right-bracket is known as a bracketing interval. The existence of such an interval allows using bisection search, a simple iterative root-finding technique that is guaranteed to converge to the solution. In each iteration, bisection search splits the interval into two parts and discards the subinterval that does not bracket the solution—in this way, it can be interpreted as a continuous extension of binary search. The method’s robustness is clearly desirable, but its relatively slow (linear) convergence can still be improved. We use Newton-Bisection, which is a combination of the quadratically converging but potentially unsafe Newton’s method with the safety of bisection search as a fallback.

As mentioned earlier, all of the following steps assume that the spline segment under consideration is defined on the interval left-bracket 0 comma 1 right-bracket with endpoint values f0 and f1 and derivative estimates d0 and d1. We will use the variable t to denote positions in this shifted and scaled interval and the values a and b store the current interval extent. Fhat stores the value of ModifyingAbove upper F With caret left-parenthesis t right-parenthesis and fhat stores ModifyingAbove f With caret left-parenthesis t right-parenthesis .

<<Invert definite integral over spline segment and return solution>>= 
<<Set initial guess for t by importance sampling a linear interpolant>> 
Float t; if (f0 != f1) t = (f0 - std::sqrt( std::max((Float)0, f0 * f0 + 2 * u * (f1 - f0)))) / (f0 - f1); else t = u / f0;
Float a = 0, b = 1, Fhat, fhat; while (true) { <<Fall back to a bisection step when t is out of bounds>> 
if (!(t >= a && t <= b)) t = 0.5f * (a + b);
<<Evaluate target function and its derivative in Horner form>> 
Fhat = t * (f0 + t * (.5f * d0 + t * ((1.f/3.f) * (-2 * d0 - d1) + f1 - f0 + t * (.25f * (d0 + d1) + .5f * (f0 - f1))))); fhat = f0 + t * (d0 + t * (-2 * d0 - d1 + 3 * (f1 - f0) + t * (d0 + d1 + 2 * (f0 - f1))));
<<Stop the iteration if converged>> 
if (std::abs(Fhat - u) < 1e-6f || b - a < 1e-6f) break;
<<Update bisection bounds using updated t>> 
if (Fhat - u < 0) a = t; else b = t;
<<Perform a Newton step>> 
t -= (Fhat - u) / fhat;
} <<Return the sample position and function value>> 
if (fval) *fval = fhat; if (pdf) *pdf = fhat / F[n - 1]; return x0 + width * t;

The number of required Newton-Bisection iterations can be reduced by starting the algorithm with a good initial guess. We use a heuristic that assumes that the spline segment is linear, i.e.,

ModifyingAbove f With caret left-parenthesis t right-parenthesis equals left-parenthesis 1 minus t right-parenthesis f left-parenthesis 0 right-parenthesis plus t f left-parenthesis 1 right-parenthesis period

Then the definite integral

ModifyingAbove upper F With caret left-parenthesis t right-parenthesis equals integral Subscript 0 Superscript t Baseline ModifyingAbove f With caret left-parenthesis t prime right-parenthesis normal d t Superscript prime Baseline equals StartFraction t Over 2 EndFraction left-parenthesis t f left-parenthesis 1 right-parenthesis minus left-parenthesis t minus 2 right-parenthesis f left-parenthesis 0 right-parenthesis right-parenthesis

has the inverse

ModifyingAbove upper F With caret Superscript negative 1 Baseline left-parenthesis xi right-parenthesis equals StartLayout Enlarged left-brace 1st Row 1st Column StartFraction f left-parenthesis 0 right-parenthesis plus-or-minus StartRoot f left-parenthesis 0 right-parenthesis squared minus 2 f left-parenthesis 0 right-parenthesis xi plus 2 f left-parenthesis 1 right-parenthesis xi EndRoot Over f left-parenthesis 0 right-parenthesis minus f left-parenthesis 1 right-parenthesis EndFraction 2nd Column f left-parenthesis 0 right-parenthesis not-equals f left-parenthesis 1 right-parenthesis 2nd Row 1st Column StartFraction xi Over f left-parenthesis 0 right-parenthesis EndFraction 2nd Column otherwise comma EndLayout

of which only one of the quadratic roots is relevant (the other one yields values outside of left-bracket 0 comma 1 right-bracket ).

<<Set initial guess for t by importance sampling a linear interpolant>>= 
Float t; if (f0 != f1) t = (f0 - std::sqrt( std::max((Float)0, f0 * f0 + 2 * u * (f1 - f0)))) / (f0 - f1); else t = u / f0;

The first fragment in the inner loop checks if the current iterate is inside the bracketing interval left-bracket monospace a comma monospace b right-bracket . Otherwise it is reset to the interval center, resulting in a standard bisection step (Figure 14.6).

Figure 14.6: The Robustness of Newton-Bisection. (left) this function increases monotonically and contains a single root on the shown interval, but a naive application of Newton’s method diverges. (right) the bisection feature of the robust root-finder enables recovery from the third Newton step, which jumps far away from the root (the bisection interval is highlighted). The method converges a few iterations later.

<<Fall back to a bisection step when t is out of bounds>>= 
if (!(t >= a && t <= b)) t = 0.5f * (a + b);

Next, F is initialized by evaluating the quartic ModifyingAbove upper F With caret left-parenthesis t right-parenthesis from Equation (14.7). For Newton’s method, we also require the derivative of this function, which is simply the original cubic spline—thus, f is set to the spline evaluated at t. The following expressions result after converting both functions to Horner form:

<<Evaluate target function and its derivative in Horner form>>= 
Fhat = t * (f0 + t * (.5f * d0 + t * ((1.f/3.f) * (-2 * d0 - d1) + f1 - f0 + t * (.25f * (d0 + d1) + .5f * (f0 - f1))))); fhat = f0 + t * (d0 + t * (-2 * d0 - d1 + 3 * (f1 - f0) + t * (d0 + d1 + 2 * (f0 - f1))));

The iteration should stop either if Fhat - u is close to 0, or if the bracketing interval has become sufficiently small.

<<Stop the iteration if converged>>= 
if (std::abs(Fhat - u) < 1e-6f || b - a < 1e-6f) break;

If ModifyingAbove upper F With caret left-parenthesis monospace t right-parenthesis minus monospace u less-than 0 , then the monotonicity of ModifyingAbove upper F With caret implies that the interval left-bracket monospace a comma monospace t right-bracket cannot possibly contain the solution (and a similar statement holds for b ). The next fragment uses this information to update the bracketing interval.

<<Update bisection bounds using updated t>>= 
if (Fhat - u < 0) a = t; else b = t;

Finally, the function and derivative values are used in a Newton step.

<<Perform a Newton step>>= 
t -= (Fhat - u) / fhat;

Once converged, the last fragment maps t back to the original interval. The function optionally returns the spline value and probability density at this position.

<<Return the sample position and function value>>= 
if (fval) *fval = fhat; if (pdf) *pdf = fhat / F[n - 1]; return x0 + width * t;

Sampling Fourier Expansions

Next, we’ll discuss sample generation for the Fourier series, Equation (8.21), using a Newton-Bisection-type method that is very similar to what was used in SampleCatmullRom(). We’d like to sample from a distribution that matches

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

for given Fourier coefficients a Subscript k . Integrating gives a simple analytic expression:

upper F left-parenthesis phi right-parenthesis equals integral Subscript 0 Superscript phi Baseline f left-parenthesis phi Superscript prime Baseline right-parenthesis normal d phi Superscript prime Baseline equals a 0 phi plus sigma-summation Underscript k equals 1 Overscript m minus 1 Endscripts StartFraction a Subscript k Baseline Over k EndFraction sine left-parenthesis k phi right-parenthesis comma

though note that this isn’t necessarily a normalized CDF: upper F left-parenthesis 2 pi right-parenthesis equals 2 pi a 0 , since the sine left-parenthesis k phi right-parenthesis terms are all zero at phi equals 2 pi .

The function SampleFourier() numerically inverts upper F left-parenthesis phi right-parenthesis to sample azimuths using the inversion method. It takes an array ak of Fourier coefficients of length m as input. The u parameter is used to pass a uniform variate, and recip should be a pointer to an array of m integer reciprocals. SampleFourier() returns the value of the Fourier expansion at the sampled position, which is stored in *phiPtr along with a probability density in *pdf.

<<Fourier Interpolation Definitions>>+= 
Float SampleFourier(const Float *ak, const Float *recip, int m, Float u, Float *pdf, Float *phiPtr) { <<Pick a side and declare bisection variables>> 
bool flip = (u >= 0.5); if (flip) u = 1 - 2 * (u - .5f); else u *= 2; double a = 0, b = Pi, phi = 0.5 * Pi; double F, f;
while (true) { <<Evaluate upper F left-parenthesis phi right-parenthesis and its derivative f left-parenthesis phi right-parenthesis >> 
<<Initialize sine and cosine iterates>> 
double cosPhi = std::cos(phi); double sinPhi = std::sqrt(1 - cosPhi * cosPhi); double cosPhiPrev = cosPhi, cosPhiCur = 1; double sinPhiPrev = -sinPhi, sinPhiCur = 0;
<<Initialize F and f with the first series term>> 
F = ak[0] * phi; f = ak[0];
for (int k = 1; k < m; ++k) { <<Compute next sine and cosine iterates>> 
double sinPhiNext = 2 * cosPhi * sinPhiCur - sinPhiPrev; double cosPhiNext = 2 * cosPhi * cosPhiCur - cosPhiPrev; sinPhiPrev = sinPhiCur; sinPhiCur = sinPhiNext; cosPhiPrev = cosPhiCur; cosPhiCur = cosPhiNext;
<<Add the next series term to F and f>> 
F += ak[k] * recip[k] * sinPhiNext; f += ak[k] * cosPhiNext;
} F -= u * ak[0] * Pi;
<<Update bisection bounds using updated phi >> 
if (F > 0) b = phi; else a = phi;
<<Stop the Fourier bisection iteration if converged>> 
if (std::abs(F) < 1e-6f || b - a < 1e-6f) break;
<<Perform a Newton step given f left-parenthesis phi right-parenthesis and upper F left-parenthesis phi right-parenthesis >> 
phi -= F / f;
<<Fall back to a bisection step when phi is out of bounds>> 
if (!(phi >= a && phi <= b)) phi = 0.5f * (a + b);
} <<Potentially flip phi and return the result>> 
if (flip) phi = 2 * Pi - phi; *pdf = (Float)f / (2 * Pi * ak[0]); *phiPtr = (Float)phi; return f;
}

Since SampleFourier() operates on even functions that are periodic on the interval left-bracket 0 comma 2 pi right-bracket , the probability of generating a sample in each of the two subintervals left-bracket 0 comma pi right-bracket and left-bracket pi comma 2 pi right-bracket