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

`UniformSampleDisk()`distorts areas substantially. Each section of the disk here has equal area and represents of the unit square of uniform random samples in each direction. In general, we’d prefer a mapping that did a better job at mapping nearby values to nearby points on the disk.

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

`RealisticCamera`ray>>=

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

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.

*(Grace Cathedral environment map courtesy of Paul Debevec.)*

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.