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.
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.
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.)
The constructor makes its own copy of the function values and computes the function’s CDF. Note that the constructor allocates n+1 Floats 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.)
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].)
Mapping u to an interval matching the above criterion is carried out using the efficient binary search implemented in FindInterval().
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.