2.3 Sampling Using the Inversion Method

To evaluate the Monte Carlo estimator in Equation (2.7), it is necessary to be able to draw random samples from a chosen probability distribution. There are a variety of techniques for doing so, but one of the most important for rendering is the inversion method, which maps uniform samples from left-bracket 0 comma 1 right-parenthesis to a given 1D probability distribution by inverting the distribution’s CDF. (In Section 2.4.2 we will see how this approach can be applied to higher-dimensional functions by considering a sequence of 1D distributions.) When used with well-distributed samples such as those generated by the samplers that are defined in Chapter 8, the inversion method can be particularly effective. Throughout the remainder of the book, we will see the application of the inversion method to generate samples from the distributions defined by BSDFs, light sources, cameras, and scattering media.

2.3.1 Discrete Case

Equation (2.2) leads to an algorithm for sampling from a set of discrete probabilities using a uniform random variable. Suppose we have a process with four possible outcomes where the probabilities of each of the four outcomes are given by p 1 , p 2 , p 3 , and p 4 , with sigma-summation Underscript i Endscripts p Subscript i Baseline equals 1 . The corresponding PMF is shown in Figure 2.3.

Figure 2.3: A PMF for Four Events, Each with a Probability p Subscript i . The sum of their probabilities sigma-summation Underscript i Endscripts p Subscript i is necessarily 1.

There is a direct connection between the sums in Equation (2.2) and the definition of the CDF. The discrete CDF is given by

upper P Subscript i Baseline equals sigma-summation Underscript j equals 1 Overscript i Endscripts p Subscript j Baseline comma

which can be interpreted graphically by stacking the bars of the PMF on top of each other, starting at the left. This idea is shown in Figure 2.4.

Figure 2.4: A Discrete CDF, Corresponding to the PMF in Figure 2.3. Each column’s height is given by the PMF for the event that it represents plus the sum of the PMFs for the previous events, upper P Subscript i Baseline equals sigma-summation Underscript j equals 1 Overscript i Endscripts p Subscript j .

The sampling operation of Equation (2.2) can be expressed as finding i such that

upper P Subscript i minus 1 Baseline less-than-or-equal-to xi Subscript Baseline less-than upper P Subscript i Baseline comma
(2.19)

which can be interpreted as inverting the CDF upper P , and thus, the name of the technique. Continuing the graphical interpretation, this sampling operation can be considered in terms of projecting the events’ probabilities onto the vertical axis where they cover the range left-bracket 0 comma 1 right-bracket and using a random variable xi Subscript to select among them (see Figure 2.5). It should be clear that this draws from the correct distribution—the probability of the uniform sample hitting any particular bar is exactly equal to the height of that bar.

The SampleDiscrete() function implements this algorithm. It takes a not-necessarily normalized set of nonnegative weights, a uniform random sample u, and returns the index of one of the weights with probability proportional to its weight. The sampling operation it performs corresponds to finding i such that

sigma-summation Underscript j equals 1 Overscript i minus 1 Endscripts w Subscript j Baseline less-than-or-equal-to xi sigma-summation w Subscript i Baseline less-than sigma-summation Underscript j equals 1 Overscript i Endscripts w Subscript j Baseline comma
(2.20)

which corresponds to multiplying Equation (2.19) by sigma-summation w Subscript i . (Not requiring a normalized PMF is a convenience for calling code and not much more work in the function’s implementation.) Two optional parameters are provided to return the value of the PMF for the sample as well as a new uniform random sample that is derived from u.

This function is designed for the case where only a single sample needs to be generated from the weights’ distribution; if multiple samples are required, the AliasTable, which will be introduced in Section A.1, should generally be used instead: it generates samples in upper O left-parenthesis 1 right-parenthesis time after an upper O left-parenthesis n right-parenthesis preprocessing step, whereas SampleDiscrete() requires upper O left-parenthesis n right-parenthesis time for each sample generated.

<<Sampling Inline Functions>>+=  
int SampleDiscrete(pstd::span<const Float> weights, Float u, Float *pmf, Float *uRemapped) { <<Handle empty weights for discrete sampling>> 
if (weights.empty()) { if (pmf) *pmf = 0; return -1; }
<<Compute sum of weights>> 
Float sumWeights = 0; for (Float w : weights) sumWeights += w;
<<Compute rescaled u prime sample>> 
Float up = u * sumWeights; if (up == sumWeights) up = NextFloatDown(up);
<<Find offset in weights corresponding to u prime >> 
int offset = 0; Float sum = 0; while (sum + weights[offset] <= up) sum += weights[offset++];
<<Compute PMF and remapped u value, if necessary>> 
if (pmf) *pmf = weights[offset] / sumWeights; if (uRemapped) *uRemapped = std::min((up - sum) / weights[offset], OneMinusEpsilon);
return offset; }

The case of weights being empty is handled first so that subsequent code can assume that there is at least one weight.

<<Handle empty weights for discrete sampling>>= 
if (weights.empty()) { if (pmf) *pmf = 0; return -1; }

The discrete probability of sampling the ith element is given by weights[i] divided by the sum of all weight values. Therefore, the function computes that sum next.

<<Compute sum of weights>>= 
Float sumWeights = 0; for (Float w : weights) sumWeights += w;

Following Equation (2.20), the uniform sample u is scaled by the sum of the weights to get a value u prime that will be used to sample from them. Even though the provided u value should be in the range left-bracket 0 comma 1 right-parenthesis , it is possible that u * sumWeights will be equal to sumWeights due to floating-point round-off. In that rare case, up is bumped down to the next lower floating-point value so that subsequent code can assume that up < sumWeights.

<<Compute rescaled u prime sample>>= 
Float up = u * sumWeights; if (up == sumWeights) up = NextFloatDown(up);

We would now like to find the last offset in the weights array i where the random sample up is greater than the sum of weights up to i . Sampling is performed using a linear search from the start of the array, accumulating a sum of weights until the sum would be greater than u prime .

<<Find offset in weights corresponding to u prime >>= 
int offset = 0; Float sum = 0; while (sum + weights[offset] <= up) sum += weights[offset++];

After the while loop terminates, the randomness in the provided sample u has only been used to select an element of the array—a discrete choice. The offset of a sample between the CDF values that bracket it is itself a uniform random value that can easily be remapped to left-bracket 0 comma 1 right-parenthesis . This value is returned to the caller in uRemapped, if requested.

One might ask: why bother? It is not too difficult to generate uniform random variables, so the benefit of providing this option may seem marginal. However, for some of the high-quality sample generation algorithms in Chapter 8, it can be beneficial to reuse samples in this way rather than generating new ones—thus, this option is provided.

<<Compute PMF and remapped u value, if necessary>>= 
if (pmf) *pmf = weights[offset] / sumWeights; if (uRemapped) *uRemapped = std::min((up - sum) / weights[offset], OneMinusEpsilon);

Figure 2.5: To use the inversion method to draw a sample from the distribution described by the PMF in Figure 2.3, a canonical uniform random variable is plotted on the vertical axis. By construction, the horizontal extension of xi Subscript will intersect the box representing the i th outcome with probability p Subscript i . If the corresponding event is chosen for a set of random variables xi Subscript , then the resulting distribution of events will be distributed according to the PMF.

2.3.2 Continuous Case

In order to generalize this technique to continuous distributions, consider what happens as the number of discrete possibilities approaches infinity. The PMF from Figure 2.3 becomes a PDF, and the CDF from Figure 2.4 becomes its integral. The projection process is still the same, but it has a convenient mathematical interpretation—it represents inverting the CDF and evaluating the inverse at xi Subscript .

More precisely, we can draw a sample upper X Subscript i from a PDF p left-parenthesis x right-parenthesis with the following steps:

  1. Integrate the PDF to find the CDF upper P left-parenthesis x right-parenthesis equals integral Subscript 0 Superscript x Baseline p left-parenthesis x prime right-parenthesis normal d x Superscript prime Baseline .
  2. Obtain a uniformly distributed random number xi Subscript .
  3. Generate a sample by solving xi equals upper P left-parenthesis upper X right-parenthesis for upper X ; in other words, find upper X equals upper P Superscript negative 1 Baseline left-parenthesis xi Subscript Baseline right-parenthesis .

We will illustrate this algorithm with a simple example; see Section A.4 for its application to a number of additional functions.

Sampling a Linear Function

The function f left-parenthesis x right-parenthesis equals left-parenthesis 1 minus x right-parenthesis a plus x b defined over left-bracket 0 comma 1 right-bracket linearly interpolates between a at x equals 0 and b at x equals 1 . Here we will assume that a comma b greater-than-or-equal-to 0 ; an exercise at the end of the chapter discusses the more general case.

<<Math Inline Functions>>= 
Float Lerp(Float x, Float a, Float b) { return (1 - x) * a + x * b; }

The function’s integral is integral Subscript 0 Superscript 1 Baseline f left-parenthesis x right-parenthesis normal d x equals left-parenthesis a plus b right-parenthesis slash 2 , which gives the normalization constant 2 slash left-parenthesis a plus b right-parenthesis to define its PDF,

p left-parenthesis x right-parenthesis equals StartFraction 2 f left-parenthesis x right-parenthesis Over a plus b EndFraction period

<<Sampling Inline Functions>>+=  
Float LinearPDF(Float x, Float a, Float b) { if (x < 0 || x > 1) return 0; return 2 * Lerp(x, a, b) / (a + b); }

Integrating the PDF gives the CDF, which is the quadratic function

upper P left-parenthesis x right-parenthesis equals StartFraction x left-parenthesis a left-parenthesis 2 minus x right-parenthesis plus b x right-parenthesis Over a plus b EndFraction period

Inverting xi equals upper P left-parenthesis upper X right-parenthesis gives a sampling recipe

upper X equals StartFraction a minus StartRoot left-parenthesis 1 minus xi right-parenthesis a squared plus xi b squared EndRoot Over a minus b EndFraction comma

though note that in this form, the case a equals b gives an indeterminate result. The more stable formulation

upper X equals StartFraction xi left-parenthesis a plus b right-parenthesis Over a plus StartRoot left-parenthesis 1 minus xi right-parenthesis a squared plus xi b squared EndRoot EndFraction

computes the same result and is implemented here.

<<Sampling Inline Functions>>+=  
Float SampleLinear(Float u, Float a, Float b) { if (u == 0 && a == 0) return 0; Float x = u * (a + b) / (a + std::sqrt(Lerp(u, Sqr(a), Sqr(b)))); return std::min(x, OneMinusEpsilon); }

One detail to note is the std::min call in the return statement, which ensures that the returned value is within the range left-bracket 0 comma 1 right-parenthesis . Although the sampling algorithm generates values in that range given xi element-of left-bracket 0 comma 1 right-parenthesis , round-off error may cause the result to be equal to 1. Because some of the code that calls the sampling routines depends on the returned values being in the specified range, the sampling routines must ensure this is so.

In addition to providing functions that sample from a distribution and compute the PDF of a sample, pbrt usually also provides functions that invert sampling operations, returning the random sample xi Subscript that corresponds to a value x . In the 1D case, this is equivalent to evaluating the CDF.

<<Sampling Inline Functions>>+=  
Float InvertLinearSample(Float x, Float a, Float b) { return x * (a * (2 - x) + b * x) / (a + b); }