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

`tan2Theta`and

`phi`for anisotropic Beckmann distribution>>

`wh`>>

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

`tan2Theta`and

`phi`for anisotropic Beckmann distribution>>

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

`tan2Theta`and

`phi`for anisotropic Beckmann distribution>>

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.

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

`wi`for microfacet reflection>> 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
, 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

`wi`for microfacet reflection>>=

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.

`dwh_dwi`for microfacet transmission>>

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

`FresnelSpecular`>>

`FresnelSpecular`>>

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

`FresnelSpecular`>>=

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

`FresnelSpecular`>>=

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

`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

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.

`FourierBSDF`>>

`ak`coefficients>>

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

`FourierBSDF`>>

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.

`FourierBSDF`>>=

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.

`FourierBSDF`>>=

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

`wi`>> }

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.

`wi`>>=

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

`u`to a spline interval by inverting

`F`>>

`i`>>

`u`for continous spline sampling step>>

`t`is out of bounds>>

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

`u`to a spline interval by inverting

`F`>>=

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

`i`>>=

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

`u`for continous spline sampling step>>=

#### 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()`.

`alpha`parameter>>

`u`to a spline interval by inverting the interpolated

`cdf`>>

`u`using the interpolated

`cdf`>>

`t`is out of bounds>>

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

`alpha`parameter>>=

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 `values[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:

- The function is the definite integral over the (assumed nonnegative) interpolant , and so it increases monotonically.
- 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
.

`t`is out of bounds>>

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

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

`t`is out of bounds>>=

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.

`t`>>=

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

`F`and

`f`with the first series term>>

`F`and

`f`>>

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