## A.4 Sampling 1D Functions

Throughout the implementation of `pbrt` we have found it useful to draw
samples from a wide variety of functions. This section therefore presents
the implementations of additional functions for sampling in 1D to augment
the ones in Section 2.3.2. All are based
on the inversion method and most introduce useful tricks for
sampling that are helpful to know when deriving new sampling algorithms.

### A.4.1 Sampling the Tent Function

`SampleTent()` uses `SampleLinear()` to sample the “tent”
function with radius ,

The sampling algorithm first uses the provided uniform
sample `u` with the `SampleDiscrete()` function to choose whether to sample a value greater than or less
than zero, with each possibility having equal probability. Note the use of
`SampleDiscrete()`’s capability of returning a new uniform random
sample here, overwriting `u`’s original value. In turn, one of the two
linear functions is sampled, with the result scaled so that the interval
is sampled.

One thing to note in this function is that the cases and expressions have
been carefully crafted so that `u==0` maps to `-r` and then
as `u` increases, the sampled value increases monotonically until
`u==1` maps to `r`, without
any jumps or reversals. This property is helpful for preserving well-distributed
sample points (e.g., if they have low discrepancy).

The tent function is easily normalized to find its PDF.

The inversion function is based on `InvertLinearSample()`.

### A.4.2 Sampling Exponential Distributions

Sampling the transmittance function when rendering images with participating media often requires samples from a distribution . As before, the first step is to find a constant that normalizes this distribution so that it integrates to one. In this case, we will assume for now that the range of values we’d like the generated samples to cover is rather than , so

Thus, and our PDF is .

We can integrate to find :

which gives a function that is easy to invert:

Therefore, we can draw samples using

It may be tempting to simplify the log term from to , under the theory that because , these are effectively the same and a subtraction can thus be saved. The problem with this idea is that may have the value 0 but never has the value 1. With the simplification, it is possible that we would try to take the logarithm of 0, which is undefined; this danger is avoided with the first formulation. While a value of 0 may seem very unlikely, it is possible, especially in the world of floating-point arithmetic and not the real numbers. Sample generation algorithms based on the radical inverse function (Section 8.6.1) are particularly prone to generating the value 0.

As before, the inverse sampling function is given by evaluating .

### A.4.3 Sampling the Gaussian

The Gaussian function is parameterized by its center and standard deviation :

The probability distribution it defines is called the *normal
distribution*. The Gaussian is already normalized, so the PDF follows
directly.

However, the Gaussian’s CDF cannot be expressed with elementary functions. It is

where is the error function. If we equate and solve, we find that:

The inverse error function can be well approximated with a polynomial, which in turn gives a sampling technique.

`InvertNormalSample()`, not included
here, evaluates .

The *Box-Muller transform* is an alternative sampling technique for
the normal distribution; it takes a pair of random samples and returns a
pair of normally distributed samples. It makes use of the fact that if two
normally distributed variables are considered as a 2D point and transformed to 2D polar coordinates
, then and .

### A.4.4 Sampling the Logistic Function

The *logistic function* is shaped similarly to the Gaussian, but can
be sampled directly. It is therefore useful in cases where a distribution
similar to the Gaussian is useful but an exact Gaussian is not needed. (It
is used, for example, in the implementation of `pbrt`’s scattering model for
hair.) The logistic function centered at the origin is

where is a parameter that controls its rate of falloff similar to in the Gaussian. Figure A.4 shows a plot of the logistic and Gaussian functions with parameter values that lead to curves with similar shapes.

The logistic function is normalized by design, so the PDF evaluation function follows directly.

Its CDF,

is easily found, and can be inverted to derive a sampling routine. The
result is implemented in `SampleLogistic()`.

As usual in 1D, the sample inversion method is performed by evaluating the CDF.

### A.4.5 Sampling a Function over an Interval

It is sometimes useful to sample from a function’s distribution over a specified interval . It turns out that this is easy to do if we are able to evaluate the function’s CDF. We will use the logistic function as an example here, though the underlying technique applies more generally.

First consider the task of finding the PDF of the function limited to the interval, : we need to renormalize it. Doing so requires being able to integrate , which is otherwise known as finding its CDF:

The function to evaluate the PDF follows directly. Here we have wrapped a
call to `InvertLogisticSample()` in a simple lambda expression in
order to make the relationship to Equation (A.17)
more clear.

Next, consider sampling using the inversion method. Following the definition of , we can see that the CDF associated with is

Setting and solving for , we have

Thus, if we compute a new value (that, in a slight abuse of notation, is not between 0 and 1) by using to linearly interpolate between and and then apply the original sampling algorithm, we will generate a sample from the distribution over the interval .

The inversion routine follows directly from Equation (A.17).

### A.4.6 Sampling Non-Invertible CDFs

It was not possible to invert the normal distribution’s CDF to derive a
sampling technique, so there we used a polynomial approximation of the
inverse CDF. In cases like that, another option is to use numerical
root–finding techniques. We will demonstrate that approach using the
*smoothstep* function as an example.

Smoothstep defines an s-shaped curve based on a third-degree polynomial that goes from zero to one starting at a point and ending at a point . It is zero for values and one for values . Otherwise, it is defined as

with . In `pbrt` the smoothstep function is used
to define the falloff at the edges of a spotlight.

We will consider the task of sampling the function within the range . First, it is easy to show that the PDF is

Integrating the PDF is also easy; the resulting CDF is

The challenge in sampling is evident: doing so requires solving a fourth-degree polynomial.

The sampling task can be expressed as a zero-finding problem: to apply the
inversion method, we would like to solve for . Doing so is
equivalent to finding the value such that . The
`SampleSmoothStep()` function below uses a Newton-bisection solver that
is defined in Section B.2.10 to do this. That function
takes a callback that returns the value of the function and its derivative
at a given point; these values are easily computed given the equations
derived so far.

Sample inversion can be performed following the same approach as was used earlier in Equation (A.17) for the logistic over an interval.

### A.4.7 Sampling Piecewise-Constant 1D Functions

The inversion method can also be applied to tabularized functions; in this
section, we will consider piecewise-constant functions defined over
. The algorithms described here will provide the foundation for
sampling piecewise-constant 2D functions, used in multiple parts
of `pbrt` to sample from distributions defined by images.

Assume that the 1D function’s domain is split into equal-sized pieces of size . These regions start and end at points , where ranges from 0 to , inclusive. Within each region, the value of the function is a constant (Figure A.5(a)).

The value of is then

The function need not always be positive, though its PDF must be. Therefore, the absolute value of the function is taken to define its PDF. The integral is

and so it is easy to construct the PDF for as . By direct application of the relevant formulae, the CDF is a piecewise-linear function defined at points by

Between two points and , the CDF is linearly increasing with slope .

Recall that in order to sample we need to invert the CDF to find the value such that

Because the CDF is monotonically increasing, the value of must be between the and such that . Given an array of CDF values, this pair of values can be efficiently found with a binary search.

The `PiecewiseConstant1D` class brings these ideas together to provide
methods for efficient sampling and PDF evaluation of this class of
functions.

`func`>>

`offset`>>

The `PiecewiseConstant1D` constructor takes `n` values of a
piecewise-constant function `f` defined over a range . (The generalization to a non- interval simply requires
remapping returned samples to the specified range and renormalizing the
PDF based on its extent.)

`func`>>

The constructor makes its own copy of the function values and computes the
function’s CDF. Note that the constructor allocates `n+1`
`Float`s for the `cdf` array because if has step
values, then there are values that define the CDF.
Storing the final CDF value of 1 is
redundant but simplifies the sampling code later.

Because the specified function may be negative, the absolute value of it is
taken here first. (There is no further need for the original function in
the `PiecewiseConstant1D` implementation.)

`func`>>=

Next, the integral of at each point is computed using
Equation (A.17), with the result
stored in the `cdf` array for now.

With the value of the integral stored in
`cdf[n]`, this value can be copied into `funcInt` and the CDF can
be normalized by dividing through all entries by this value.
The case of a zero-valued function is handled by defining a linear CDF,
which leads to uniform sampling. That case occurs more frequently than one
might expect due to the use of this class when sampling piecewise-constant
2D functions; when that is used with images, images
with zero-valued scanlines lead to zero-valued functions here.

The integral of the absolute value of the function is made available via a
method and the `size()` method returns the number of tabularized values.

The `PiecewiseConstant1D::Sample()` method uses the given random sample
`u` to sample from its distribution. It returns the corresponding
value and the value of the PDF . If the optional
`offset` parameter is not `nullptr`, it returns the offset into the
array of function values of the largest index where the CDF was less than
or equal to `u`. (In other words, `cdf[*offset] <= u < cdf[*offset+1]`.)

`offset`>>

Mapping `u` to an interval matching the above criterion is carried out
using the efficient binary search implemented in `FindInterval()`.

`offset`>>=

Given the pair of CDF values that straddle `u`, we can compute .
First, we determine how far `u` is between `cdf[o]` and
`cdf[o+1]`. We denote this value with `du`, where `du` is 0 if `u ==
cdf[o]` and goes up to 1 if `u == cdf[o+1]`. Because the
CDF is piecewise-linear, the sample value is the same offset between
and (Figure A.5(b)).

The PDF for this sample is easily computed since we have the
function’s integral in `funcInt`. (Note that the offset `o`
into the CDF array has been computed in a way so that `func[o]`
gives the value of the function in the CDF range that the sample landed in.)

Finally, the appropriate value of is computed and returned. Here is where the sampled value in is remapped to the user-specified range .

As with the other sampling routines so far, `PiecewiseConstant1D`
provides an inversion method that takes a point in the range
and returns the sample value that maps
to it. As before, this is a matter of evaluating the CDF at the
given position.

Because the CDF is tabularized at regular steps over , if we remap to lie within , scale by the
number of CDF values, and take the floor of that value, we have the offset
to the entry in the `cdf` array that precedes .

Given those two points, we linearly interpolate between their values to evaluate the CDF.