13.6 2D Sampling with Multidimensional Transformations

Suppose we have a 2D joint density function that we wish to draw samples from. Sometimes multidimensional densities are separable and can be expressed as the product of 1D densities—for example,

for some and . In this case, random variables can be found by independently sampling from and from . Many useful densities aren’t separable, however, so we will introduce the theory of how to sample from multidimensional distributions in the general case.

Given a 2D density function, the marginal density function is obtained by “integrating out” one of the dimensions:

This can be thought of as the density function for alone. More precisely, it is the average density for a particular over all possible values.

The conditional density function is the density function for given that some particular has been chosen (it is read “ of given ”):

The basic idea for 2D sampling from joint distributions is to first compute the marginal density to isolate one particular variable and draw a sample from that density using standard 1D techniques. Once that sample is drawn, one can then compute the conditional density function given that value and draw a sample from that distribution, again using standard 1D sampling techniques.

13.6.1 Uniformly Sampling a Hemisphere

As an example, consider the task of choosing a direction on the hemisphere uniformly with respect to solid angle. Remember that a uniform distribution means that the density function is a constant, so we know that . In conjunction with the fact that the density function must integrate to one over its domain, we have

This tells us that , or (using a result from the previous example about spherical coordinates). Note that this density function is separable. Nevertheless, we will use the marginal and conditional densities to illustrate the multidimensional sampling technique.

Consider sampling first. To do so, we need ’s marginal density function :

Now, compute the conditional density for :

Notice that the density function for is itself uniform; this should make intuitive sense given the symmetry of the hemisphere. Now, we use the 1D inversion technique to sample each of these PDFs in turn:

Inverting these functions is straightforward, and here we can safely replace with , giving

Converting these back to Cartesian coordinates, we get the final sampling formulae:

This sampling strategy is implemented in the following code. Two uniform random numbers are provided in u, and a vector on the hemisphere is returned.

For each sampling routine like this in pbrt, there is a corresponding function that returns the value of the PDF for a particular sample. For such functions, it is important to be clear which PDF is being evaluated—for example, for a direction on the hemisphere, we have already seen these densities expressed differently in terms of solid angle and in terms of . For hemispheres (and all other directional sampling), these functions return values with respect to solid angle. For the hemisphere, the solid angle PDF is a constant .

Sampling the full sphere uniformly over its area follows almost exactly the same derivation, which we omit here. The end result is

13.6.2 Sampling a Unit Disk

Although the disk seems a simpler shape to sample than the hemisphere, it can be trickier to sample uniformly because it has an incorrect intuitive solution. The wrong approach is the seemingly obvious one: . Although the resulting point is both random and inside the circle, it is not uniformly distributed; it actually clumps samples near the center of the circle. Figure 13.10(a) shows a plot of samples on the unit disk when this mapping was used for a set of uniform random samples . Figure 13.10(b) shows uniformly distributed samples resulting from the following correct approach.

Since we’re going to sample uniformly with respect to area, the PDF must be a constant. By the normalization constraint, . If we transform into polar coordinates (see the example in Section 13.5.2), we have . Now we compute the marginal and conditional densities as before:

As with the hemisphere case, the fact that is a constant should make sense because of the symmetry of the circle. Integrating and inverting to find , , , and , we can find that the correct solution to generate uniformly distributed samples on a disk is

Taking the square root of effectively pushes the samples back toward the edge of the disk, counteracting the clumping referred to earlier.

Although this mapping solves the problem at hand, it distorts areas on the disk; areas on the unit square are elongated and/or compressed when mapped to the disk (Figure 13.11). (Section 13.8.3 will discuss in more detail why this distortion is a disadvantage.)

A better approach is a “concentric” mapping from the unit square to the unit circle that avoids this problem. The concentric mapping takes points in the square to the unit disk by uniformly mapping concentric squares to concentric circles (Figure 13.12).

The mapping turns wedges of the square into slices of the disk. For example, points in the shaded area in Figure 13.12 are mapped to  by

See Figure 13.13. The other seven wedges are handled analogously.

13.6.3 Cosine-Weighted Hemisphere Sampling

As we will see in Section 13.10, it is often useful to sample from a distribution that has a shape similar to that of the integrand being estimated. For example, because the scattering equation weights the product of the BSDF and the incident radiance with a cosine term, it is useful to have a method that generates directions that are more likely to be close to the top of the hemisphere, where the cosine term has a large value, than the bottom, where the cosine term is small.

Mathematically, this means that we would like to sample directions from a PDF

Normalizing as usual,

so

We could compute the marginal and conditional densities as before, but instead we can use a technique known as Malley’s method to generate these cosine-weighted points. The idea behind Malley’s method is that if we choose points uniformly from the unit disk and then generate directions by projecting the points on the disk up to the hemisphere above it, the resulting distribution of directions will have a cosine distribution (Figure 13.14).

Why does this work? Let be the polar coordinates of the point chosen on the disk (note that we’re using instead of the usual here). From our calculations before, we know that the joint density gives the density of a point sampled on the disk.

Now, we map this to the hemisphere. The vertical projection gives , which is easily seen from Figure 13.14. To complete the transformation, we need the determinant of the Jacobian

Therefore,

which is exactly what we wanted! We have used the transformation method to prove that Malley’s method generates directions with a cosine-weighted distribution. Note that this technique works regardless of the method used to sample points from the circle, so we can use the earlier concentric mapping just as well as the simpler method.

Remember that all of the directional PDF evaluation routines in pbrt are defined with respect to solid angle, not spherical coordinates, so the PDF function returns a weight of .

13.6.4 Sampling a Cone

For both area light sources based on Spheres as well as for the SpotLight, it’s useful to be able to uniformly sample rays in a cone of directions. Such distributions are separable in , with and so we therefore need to derive a method to sample a direction uniformly over the cone of directions around a central direction up to the maximum angle of the beam, . Incorporating the term from the measure on the unit sphere from Equation (5.5),

So and .

The PDF can be integrated to find the CDF, and the sampling technique,

follows. There are two UniformSampleCone() functions that implement this sampling technique: the first samples about the axis, and the second (not shown here) takes three basis vectors for the coordinate system to be used where samples taken are with respect to the axis of the given coordinate system.

13.6.5 Sampling a Triangle

Although uniformly sampling a triangle may seem like a simple task, it turns out to be more complex than the ones we’ve seen so far. ^† To simplify the problem, we will assume we are sampling an isosceles right triangle of area . The output of the sampling routine that we will derive will be barycentric coordinates, however, so the technique will actually work for any triangle despite this simplification. Figure 13.15 shows the shape to be sampled.

We will denote the two barycentric coordinates here by . Since we are sampling with respect to area, we know that the PDF must be a constant equal to the reciprocal of the shape’s area, , so .

First, we find the marginal density :

and the conditional density :

The CDFs are, as always, found by integration:

Inverting these functions and assigning them to uniform random variables gives the final sampling strategy:

Notice that the two variables in this case are not independent!

We won’t provide a PDF evaluation function for this sampling strategy since the PDF’s value is just one over the triangle’s area.

13.6.6 Sampling Cameras

Section 6.4.7 described the radiometry of the measurement that is performed by film in a realistic camera. For reference, the equation that gave joules arriving over an area of the sensor, Equation (6.7), is

where the areas integrated over are the area of a pixel on the film, , and an area over the rear lens element that bounds the lens’s exit pupil, . The integration over the pixel in the film is handled by the Integrator and the Sampler, so here we’ll consider the estimate

for a fixed point on the film plane .

Given a PDF for sampling a time and a PDF for sampling a point on the exit pupil, , we have the estimator

For a uniformly sampled time value, . The weighting factor that should be applied to the incident radiance function to compute an estimate of Equation (6.7) is then

Recall from Section 6.4.5 that points on the exit pupil were uniformly sampled within a 2D bounding box. Therefore, to compute , we just need one over the area of this bounding box. Conveniently, RealisticCamera::SampleExitPupil() returns this value, which GenerateRay() stores in exitPupilBoundsArea.

The computation of this weight is implemented in <<Return weighting for RealisticCamera ray>>; we can now see that if the simpleWeighting member variable is true, then the RealisticCamera computes a modified version of this term that only accounts for the term. While the value it computes isn’t a useful physical quantity, this option gives images with pixel values that are closely related to the scene radiance (with vignetting), which is often convenient.

13.6.7 Piecewise-Constant 2D Distributions

Our final example will show how to sample from discrete 2D distributions. We will consider the case of a 2D function defined over by a 2D array of sample values. This case is particularly useful for generating samples from distributions defined by texture maps and environment maps.

Consider a 2D function defined by a set of values where , , and give the constant value of over the range . Given continuous values , we will use to denote the corresponding discrete indices, with and so that .

Integrals of are simple sums of values, so that, for example, the integral of over the domain is

Using the definition of the PDF and the integral of , we can find ’s PDF,

Recalling Equation (13.11), the marginal density can be computed as a sum of values

Because this function only depends on , it is thus itself a piecewise constant 1D function, , defined by values. The 1D sampling machinery from Section 13.3.1 can be applied to sampling from its distribution.

Given a sample, the conditional density is then

Note that, given a particular value of , is a piecewise-constant 1D function of , that can be sampled with the usual 1D approach. There are such distinct 1D conditional densities, one for each possible value of .

Putting this all together, the Distribution2D structure provides functionality similar to Distribution1D, except that it generates samples from piecewise-constant 2D distributions.

The constructor has two tasks. First, it computes a 1D conditional sampling density for each of the individual values using Equation (13.14). It then computes the marginal sampling density with Equation (13.13).

Distribution1D can directly compute the distributions from a pointer to each of the rows of function values, since they’re laid out linearly in memory. The and terms from Equation (13.14) don’t need to be included in the values passed to Distribution1D, since they have the same value for all of the values and are thus just a constant scale that doesn’t affect the normalized distribution that Distribution1D computes.

Given the conditional densities for each value, we can find the 1D marginal density for sampling each value, . The Distribution1D class stores the integral of the piecewise-constant function it represents in its funcInt member variable, so it’s just necessary to copy these values to the marginalFunc buffer so they’re stored linearly in memory for the Distribution1D constructor.

As described previously, in order to sample from the 2D distribution, first a sample is drawn from the marginal distribution in order to find the coordinate of the sample. The offset of the sampled function value gives the integer value that determines which of the precomputed conditional distributions should be used for sampling the value. Figure 13.16 illustrates this idea using a low-resolution image as an example.

The value of the PDF for a given sample value is computed as the product of the conditional and marginal PDFs for sampling it from the distribution.