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 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 in *pdf. The PDF value returned should be measured with respect to solid angle, and both the outgoing direction and the sampled incident direction 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: for all 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.)

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 in the reflection coordinate system. If is in the opposite hemisphere, then must be flipped to lie in the same hemisphere as . 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 direction.

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 ), it is first necessary to check that and lie on the same side of the surface; if not, the sampling probability is 0. Otherwise, the method computes . One potential pitfall with this method is that the order of the and arguments is significant. For the cosine-weighted distribution, in general. Code that calls this method must be careful to use the correct argument ordering.

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 where each microfacet exhibited perfect specular reflection and/or transmission. Because the 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.

The classic approach to sampling a microfacet orientation is to sample from 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.

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

The sampling method for the Beckmann–Spizzichino distribution’s full distribution of normals returns and the angle in spherical coordinates, which in turn are converted to a normalized direction vector 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 , and so the PDF for this distribution, , is separable into and .

is constant with a value of , and thus values can be sampled by

For , we have

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

To find the sampling technique, we need to solve

for in terms of . In this case, suffices to find the microfacet orientation and is more efficient to compute, so we will compute

The sampling code follows directly.

The algorithm to sample and 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 , we can compute using the identities and . follows, and we have enough information to compute the microfacet orientation using the spherical coordinates formula. Because pbrt transforms the normal to in the reflection coordinate system, we can almost use the computed direction from spherical coordinates directly. The last detail to handle is that if is in the opposite hemisphere than the normal, then the half-angle vector needs to be flipped to be in that hemisphere as well.

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

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 :

Here, the term accounts for microfacet self-shadowing, and the 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 have a higher probability of being visible than they do in the overall distribution .

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

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

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

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 , then its direction is under the surface and no light is reflected. ^†

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 to one in terms of using the techniques introduced in Section 13.5. Doing so requires applying the adjustment for a change of variables .

A simple geometric construction gives the relationship between the two distributions. Consider the spherical coordinate system oriented about (Figure 14.4). Measured with respect to , the differential solid angles and are and , respectively; thus,

Because is computed by reflecting about , . Furthermore, because , we can find the desired conversion factor:

Therefore, the PDF after transformation is

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

The approach for transmission is analogous: given a sampled 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.

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.

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

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 is less than or greater than . In both cases, it remaps to cover the range after using it to make this decision. (Otherwise, values of used for the cosine-weighted sampling would always be less than , for example.) Using the sample for two purposes in this manner slightly reduces the quality of the stratification of the values that are actually used for sampling directions. ^†

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

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

and

Thus, it is a probability density function, where the PDF has a value of 0 for all . 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 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 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

where is the hemispherical–directional reflectance and is the direction for perfect specular reflection or transmission.

Because the PDF has a delta term as well, , the two delta distributions cancel out, and the estimator is

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.

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

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

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

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.

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

where the function was discretized over the incident and outgoing zenith angle cosines with endpoints and . An even Fourier expansion with real coefficients was used to model the dependence on the azimuth angle difference parameter .

The task now is to first sample given and then sample the angle relative to . 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 satisfying these criteria can be found using an inner product between and the individual basis functions analogous to orthogonal projections of vectors on Euclidean vector spaces. Here, we are dealing with continuous functions on , where a suitable inner product can be defined as

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

This quantity turns out to be very helpful in constructing a method for importance sampling the BSDF: disregarding normalization factors, this average over can be interpreted as the marginal distribution of the cosine-weighted BSDF with respect to pairs of 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 in FourierBSDFTable::a0. This array is initialized by copying the corresponding entries from FourierBSDFTable::a at scene loading time in the FourierBSDFTable::Read() method.

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

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

where . Each column of this matrix stores a discrete CDF over for a different (fixed) value of . 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.

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.

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

SampleFourier() takes an additional input array recip containing precomputed integer reciprocals 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().

We now have the values and . The sampled incident direction’s coordinates are given by rotating the components of by an angle about the surface normal, and its component is given by (using spherical coordinates).

There are two details to note in the computation of the direction . First, the components are scaled by a factor , 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.

The fragment <<Evaluate remaining Fourier expansions for angle >> is identical to <<Evaluate Fourier expansion for angle >> 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 , we could simply copy the implementation of FourierBSDF::f() except for the division that cancels . 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 by a suitable normalization factor to ensure that the product integrates to 1:

Note that the outgoing zenith angle cosine was left unspecified in the above equation. In general, the normalization factor is not constant and, instead, it depends on the current value of . has a simple interpretation: it is the hemispherical–directional reflectance of a surface that is illuminated from the zenith angle .

Suppose briefly that happens to be part of the discretized set of zenith angle cosines stored in the array FourierBSDFTable::mu. Then

where 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 , we can simply interpolate the neighboring four entries of 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().

We won’t include the second fragment here—it is almost identical to <<Compute Fourier coefficients for >>, 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.

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 was evaluated at positions , resulting in a piecewise cubic Catmull-Rom spline interpolant with spline segments . Given a precomputed discrete CDF over these spline segments defined as

where is the normalization term,

a straightforward way of importance sampling in two steps using the inversion method entails first finding the interval such that

where is a random variate on the interval , and then sampling an value in the th interval. Since the values are monotonically increasing, the interval can be found using an efficient binary search.

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

Having selected a segment , we can offset and re-scale to obtain a second uniform variate in :

We then use to sample a position within the interval using that segment’s integral,

where again we won’t compute a properly normalized CDF but will instead multiply by rather than normalizing . We must then compute

This approach is implemented in SampleCatmullRom(), which takes the following inputs: the number of samples n; x contains the locations where the original function was evaluated; f stores the value of the function at each point ; u is used to pass the uniform variate ; and integrated 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.

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.

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

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.

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

Since the scaling by was already applied in the first fragment, we need only subtract .

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 , which requires an additional scaling by the associated change of variable factor equal to the reciprocal of width.

Sampling 2D Spline Interpolants

The main use cases of spline interpolants in pbrt actually importance sample 2D functions , where is considered a fixed parameter for the purpose of sampling (e.g., the albedo of the underlying material or the outgoing zenith angle cosine in the case of the FourierBSDF). This case is handled by SampleCatmullRom2D().

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.

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:

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 was defined as an integral over the cubic spline segment , 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 is the definite integral over the (assumed nonnegative) interpolant , and so it increases monotonically.
2. The interval selected by the function FindInterval() contains exactly one solution to Equation (14.8).

In this case, the interval 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 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 and fhat stores .

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

Then the definite integral

has the inverse

of which only one of the quadratic roots is relevant (the other one yields values outside of ).

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

Next, F is initialized by evaluating the quartic 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:

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

If , then the monotonicity of implies that the interval cannot possibly contain the solution (and a similar statement holds for ). The next fragment uses this information to update the bracketing interval.

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

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.

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

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

though note that this isn’t necessarily a normalized CDF: , since the terms are all zero at .

The function SampleFourier() numerically inverts 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.

Since SampleFourier() operates on even functions that are periodic on the interval , the probability of generating a sample in each of the two subintervals and