## 14.2 Sampling Light Sources

Being able to take a point and sample the directions around it where
direct illumination may be incident is another important sampling operation for
rendering. Consider a diffuse surface illuminated by a small spherical
area light source (Figure 14.7): sampling
directions using the BSDF’s sampling distribution is likely to be very
inefficient because the light is only visible along a small cone of
directions from the point. A much better approach is to instead use a
sampling distribution that is based on the light source. For example, the
sampling routine could choose from among only those directions where the
sphere is potentially visible. This section introduces the
`Light::Sample_Li()` method, which allows this operation to be
implemented for the lights in `pbrt`.

There are two sampling methods that `Light`s must implement. The
first, `Sample_Li()`, samples an incident direction at a point in the
scene along which illumination from the light may be arriving. The second,
`Light::Sample_Le()`, will be defined in
Section 16.1.2; it returns a
light-carrying ray leaving the light source. Both have corresponding
methods that return the PDF for an incident direction or a ray, respectively.

The `Light::Sample_Li()` method was first introduced in
Section 12.2. For reference, here is its declaration:

We can now understand the meaning of its `u` and `pdf` parameters:
`u` provides a 2D sample value for sampling the
light source, and the PDF for sampling the chosen direction is returned in
`*pdf`.

The `Light`’s `Pdf_Li()` method returns the probability density with
respect to solid angle for the light’s `Sample_Li()` method to sample
the direction `wi` from the reference point `ref`.

### 14.2.1 Lights with Singularities

Just as with perfect specular reflection and transmission, light sources
that are defined in terms of delta distributions fit naturally into this
sampling framework, although they require care on the part of the routines
that call their sampling methods, since there are implicit delta
distributions in the radiance and PDF values that they return. For the
most part, these delta distributions naturally cancel out when estimators
are evaluated, although multiple importance sampling code must be aware of
this case, just as was the case with BSDFs. (The `IsDeltaLight()`
utility function can be used to see if a light source is described by a
delta distribution.)

Point lights are described by a delta distribution such that they only
illuminate a receiving point from a single direction. Thus, the sampling
problem is trivial. The `PointLight::Sample_Li()` method was
already implemented in Section 12.3, where we glossed
over the nuance that the method performs Monte Carlo sampling of a delta
distribution and thus always returns a single direction and doesn’t need
the random sample values.

Due to the delta distribution, the `PointLight::Pdf_Li()` method
returns 0. This value reflects the fact that there is no chance for some other
sampling process to randomly generate a direction that would intersect an
infinitesimal light source.

The `Sample_Li()` methods for `SpotLight`s,
`ProjectionLight`s, `GonioPhotometricLight`s, and
`DistantLight`s were also previously implemented in
Sections 12.3.1, 12.3.2,
12.3.3, and 12.4, respectively. All
also return 0 from their `Pdf_Li()` methods.

### 14.2.2 Sampling Shapes

In `pbrt`, area lights are defined by attaching an emission profile to a
`Shape`. Therefore, in order to sample incident illumination from
such light sources, it is useful to be able to generate samples over the
surface of shapes. To make this possible, we will add sampling methods to
the `Shape` class that sample points on their surfaces. The
`AreaLight` sampling methods to be defined shortly will in turn call
these methods.

There are two shape sampling methods, both named `Shape::Sample()`.
The first chooses points on the surface of the shape using a sampling distribution
with respect to surface area and returns the local geometric
information about the sampled point in an `Interaction`. In addition
to initializing the position `p` and normal `n` of the sampled
point, the `Sample()` method should set `Interaction::pError` with
bounds on the floating-point rounding error in the computed `p` value.
`pError` is used to compute robust origins for rays leaving the
surface of the light—see Section 3.9.5.

`Shape`s almost always sample uniformly by area on their surface.
Therefore, we will provide a default implementation of the
`Shape::Pdf()` method corresponding to this sampling approach that
returns the corresponding PDF: 1 over the surface area.

The second shape sampling method takes the point from which the surface of the shape is being integrated over as a parameter. This method is particularly useful for lighting, since the caller can pass in the point to be lit and allow shape implementations to ensure that they only sample the portion of the shape that is potentially visible from that point. The default implementation ignores the additional point and calls the earlier sampling method.

Unlike the first `Shape` sampling method, which generates points on the
shape according to a probability density with respect to surface area on
the shape, the second one uses a density with respect to solid angle from
the reference point `ref`. This difference stems from the fact that
the area light sampling routines evaluate the direct lighting integral as
an integral over directions from the reference point—expressing these
sampling densities with respect to solid angle at the point is more
convenient. Therefore, the standard implementation of the second
`Pdf()` method here transforms the density from one defined over area
to one defined over solid angle.

Given a reference point and direction , the `Pdf()` method
determines if the ray from the point in direction intersects the
shape. If the ray doesn’t intersect the shape at all, the probability that
the shape would have chosen the direction can be assumed to be 0
(an effective sampling algorithm wouldn’t generate such a sample, and in
any case the light will not contribute to such directions, so
using a zero probability density is fine).

Note that this ray intersection test is only between the ray and the single shape under consideration. The rest of the scene geometry is ignored, and thus the intersection test is fairly efficient.

To compute the value of the PDF with respect to solid angle from the reference point, the method starts by computing the PDF with respect to surface area. Conversion from a density with respect to area to a density with respect to solid angle requires division by the factor

where is the angle between the direction of the ray from the point on the light to the reference point and the light’s surface normal, and is the distance between the point on the light and the point being shaded (recall the discussion about transforming between area and directional integration domains in Section 5.5).

#### Sampling Disks

The `Disk` sampling method uses the concentric disk sampling function
to find a point on the unit disk and then scales and offsets this point to
lie on the disk of a given radius and height. Note that this method does not
account for partial disks due to `Disk::innerRadius` being nonzero or
`Disk::phiMax` being less than . Fixing this bug is left for
an exercise at the end of the chapter.

Because the object space value of the sampled point is equal to
`Disk::height`, zero-extent bounds can be used for the error bounds for
rays leaving the sampled point, just as with ray–disk intersections.
(These bounds may later be expanded by the object to world transformation,
however.)

#### Sampling Cylinders

Uniform sampling on cylinders is straightforward. The height and value are sampled uniformly. Intuitively, it can be understood that this approach works because a cylinder is just a rolled-up rectangle.

`pObj`to cylinder surface and compute

`pObjError`>>

If the system’s `std::sin()` and `std::cos()` functions compute
results that are as precise as possible—i.e., they always return the
closest floating-point value to the fully-precise result, then it can be
shown that the and components of `pObj` are within a factor of
of the actual surface of the cylinder. While many
implementations of those functions are that precise, not all are,
especially on GPUs. To be safe, the implementation here reprojects the
sampled point to lie on the cylinder. In this case, the error bounds are
the same as were derived for reprojected ray–cylinder intersection points
in Section 3.9.4.

`pObj`to cylinder surface and compute

`pObjError`>>=

#### Sampling Triangles

The `UniformSampleTriangle()` function, defined in the previous
chapter, returns the barycentric coordinates for a uniformly sampled point
on a triangle. The point on a particular triangle for those barycentrics
is easily computed.

`p0`,

`p1`, and

`p2`>>

Because the sampled point is computed with barycentric interpolation, it has the same error bounds as were computed in Section 3.9.4 for triangle intersection points.

#### Sampling Spheres

As with `Disk`s, the sampling method here does not handle partial
spheres; an exercise at the end of the chapter discusses this issue
further. For the sampling method that is not given an external point
that’s being lit, sampling a point on a sphere is extremely simple.
`Sphere::Sample()` just uses the `UniformSampleSphere()` function
to generate a point on the unit sphere and scales the point by the sphere’s
radius.

`pObj`to sphere surface and compute

`pObjError`>>

Because `UniformSampleSphere()` uses `std::sin()` and
`std::cos()`, the error bounds on the computed `pObj` value
depend on the accuracy of those functions. Therefore, as with cylinders,
the sampled point is reprojected to the sphere’s surface, so that the error
bounds derived earlier in Equation (3.14) can be used
without needing to worry about those functions’ accuracy.

`pObj`to sphere surface and compute

`pObjError`>>=

For the sphere sampling method that is given a point being illuminated, we can do much better than sampling over the sphere’s entire area. While uniform sampling over its surface is perfectly correct, a better approach is to not sample points on the sphere that are definitely not visible from the point (such as those on the back side of the sphere as seen from the point). The sampling routine here instead uniformly samples directions over the solid angle subtended by the sphere from the reference point and then computes the point on the sphere corresponding to the sampled direction.

`Interaction`for sampled point on sphere>>

`pObj`to sphere surface and compute

`pObjError`>>

This process is most easily done if we first compute a coordinate system to use for sampling the sphere, where the axis is the vector between the sphere’s center and the point being illuminated:

For points that lie inside the surface of the sphere, the entire sphere
should be sampled, since the whole sphere is visible from inside it.
Note that the reference point used in this determination, `pOrigin`,
is computed using the `OffsetRayOrigin()` function. Doing so ensures
that if the reference point came from a ray intersecting the sphere, the
point tested lies on the correct side of the sphere.

Otherwise sampling within the cone proceeds.

`Interaction`for sampled point on sphere>>

`pObj`to sphere surface and compute

`pObjError`>>

If the reference point is outside the sphere, then as seen from the point being shaded the sphere subtends an angle

where is the radius of the sphere and is its center
(Figure 14.8). The sampling method here computes the
cosine of the
subtended angle using Equation (14.11) and then
uniformly samples directions inside this cone of directions using the
approach that was derived earlier for the `UniformSampleCone()`
function, sampling an offset from the center vector and
then uniformly sampling a rotation angle around the vector.
That function isn’t used here, however, as we will need some of the
intermediate values in the following fragments.

Given a sample angle with respect to the sampling coordinate system computed earlier, we can directly compute the corresponding sample point on the sphere. The derivation of this approach follows three steps, illustrated in Figure 14.9.

First, if we denote the distance from the reference point to the center of the sphere by and form a right triangle with angle at the reference point, then we can see that the lengths of the other two sides of the triangle, as shown in Figure 14.9(a), are and .

Next, consider the right triangle shown in Figure 14.9(b), where the hypotenuse is the line segment with length equal to the sphere’s radius that goes from the center of the sphere to the point where the line from the sampled angle intersects the sphere. From the Pythagorean theorem, we can see that the length of the third side of that triangle is

Subtracting this length from gives us the length of the segment from the reference point to the sampled point on the sphere:

We can now compute the angle from the center of the sphere to the sampled point on the sphere using the law of cosines, which relates the squared lengths of two sides of the triangle and the angle opposite them:

The angle and give the spherical coordinates for the sampled point on the unit sphere, with respect to a coordinate system centered around the vector from the sphere center to the reference point. Since we earlier computed a coordinate system from the reference point to the center, we can use that one here with each vector flipped.

As with the other `Sphere::Sample()` method, the sampled point is
reprojected onto the surface of the sphere; in turn, we can use the same
error bounds for the computed point as were derived earlier.

`Interaction`for sampled point on sphere>>=

`pObj`to sphere surface and compute

`pObjError`>>

The value of the PDF for sampling a direction toward a sphere from a given point depends on which of the two sampling strategies would be used for the point.

If the reference point was inside the sphere, a uniform sampling strategy
was used, in which case, the implementation hands off to the `Pdf()`
method of the `Shape` class, which takes care of the solid angle
conversion.

In the general case, we recompute the angle subtended by the sphere and
call `UniformConePdf()`. Note that no conversion of sampling measures
is required here because `UniformConePdf()` already returns a value
with respect to the solid angle measure.

### 14.2.3 Area Lights

Given shape sampling methods, the `DiffuseAreaLight::Sample_Li()`
method is quite straightforward. Most of the hard work is done by the
`Shape`s, and the `DiffuseAreaLight` just needs to copy
appropriate values to output parameters and compute
the emitted radiance value.

`Pdf_Li()` calls the variant of `Shape::Pdf()` that returns a
density with respect to solid angle, so the value it returns can be
returned directly.

### 14.2.4 Infinite Area Lights

The `InfiniteAreaLight`, defined in Section 12.6,
can be considered to be an infinitely large sphere that surrounds the entire
scene, illuminating it from all directions. The environment maps used with
`InfiniteAreaLight`s often have substantial variation along different
directions: consider, for example, an environment map of the sky during
daytime, where the relatively small number of directions that the sun subtends
will be thousands of times brighter than the rest of the directions. Given
this substantial variation, implementing a sampling method for
`InfiniteAreaLight`s that matches the illumination distribution will
generally substantially reduce variance in images.

Figure 14.10 shows two images of a car model illuminated by the morning skylight environment map from Figure 12.20. The top image was rendered using a simple cosine-weighted sampling distribution for selecting incident illumination directions, while the bottom image was rendered using the sampling method implemented in this section. Both images used just 32 shadow samples per pixel. For the same number of samples taken and with negligible additional computational cost, this sampling method computes a much better result with much lower variance.

There are three main steps to the sampling approach implemented here:

- We define a piecewise-constant 2D probability distribution function in image coordinates that corresponds to the distribution of the radiance represented by the environment map.
- We apply the sampling method from Section 13.6.7 to transform 2D samples to samples drawn from the piecewise-constant distribution.
- We define a probability density function over directions on the unit sphere based on the probability density over .

The combination of these three steps makes it possible to generate samples on the sphere of directions according to a distribution that matches the radiance function very closely, leading to substantial variance reduction.

We will start by defining the fragment <<Initialize sampling PDFs for
infinite area light>> at the end of the `InfiniteAreaLight` constructor.

`img`from environment map>>

The first step is to transform the continuously defined spectral radiance
function defined by the environment map’s texels to a piecewise-constant
scalar function by computing its luminance at a set of sample points using
the `Spectrum::y()` method. There are three things to note in the code
below that does this computation.

First, it computes values of the radiance function at the same number of points as there are texels in the original image map. It could use either more or fewer points, leading to a corresponding increase or decrease in memory use while still generating a valid sampling distribution, however. These values work well, though, as fewer points would lead to a sampling distribution that didn’t match the function as well while more would mostly waste memory with little incremental benefit.

The second thing of note in this code is that the piecewise constant
function values being stored here in `img` are found by slightly
blurring the radiance function with the `MIPMap::Lookup()` method
(rather than just copying the corresponding texel values). The motivation
for this is subtle; Figure 14.11 illustrates the
idea in 1D. Because the continuous radiance function used for rendering is
reconstructed by bilinearly interpolating between texels in the image, just
because some texel is completely black, for example, the radiance function
may be nonzero a tiny distance away from it due to a neighboring texel’s
contribution. Because we are using a piecewise-constant function for
sampling rather than a piecewise-linear one, it must account for this issue
in order to ensure greater-than-zero probability of sampling any point
where the radiance function is nonzero.

`InfiniteAreaLight`sampling routines, is to find the average value of the function over some range and use that to define the piecewise-constant function.

Finally, each image value in the `img` buffer is multiplied by the
value of corresponding to the value each row has
when the latitude–longitude image is mapped to the sphere. Note that this
multiplication has no effect on the sampling method’s correctness: because
the value of is always greater than 0 over the
range, we are just reshaping the sampling PDF. The motivation for
adjusting the PDF is to eliminate the effect of the distortion from mapping
the 2D image to the unit sphere in the sampling method here; the details
will be fully explained later in this section.

`img`from environment map>>=

Given this filtered and scaled image, the `Distribution2D` structure
handles computing and storing the 2D PDF.

Given this precomputed data, the task of the sampling method is relatively
straightforward. Given a sample `u`
over , it draws a sample from the
function’s distribution using the sampling algorithm described in
Section 13.6.7, which gives a value and the value
of the PDF for taking this sample, .

The sample is mapped to spherical coordinates by

and then the spherical coordinates formula gives the direction .

Recall that the probability density values returned by the light source sampling routines must be defined in terms of the solid angle measure on the unit sphere. Therefore, this routine must now compute the transformation between the sampling density used, which was the image function over , and the corresponding density after the image has been mapped to the unit sphere with the latitude–longitude mapping. (Recall that the latitude–longitude parameterization of an image is , , and .)

First, consider the function that maps from to ,

The absolute value of the determinant of the Jacobian is . Applying the multidimensional change of variables equation from Section 13.5.1, we can find the density in terms of spherical coordinates .

From the definition of spherical coordinates, it is easy to determine that the absolute value of the determinant of the Jacobian for the mapping from to is . Since we are interested in the unit sphere, , and again applying the multidimensional change of variables equation, we can find the final relationship between probability densities,

This is the key relationship for applying this technique: it lets us sample from the piecewise-constant distribution defined by the image map and transform the sample and its probability density to be in terms of directions on the unit sphere.

We can now see why the initialization routines multiplied the values
of the piecewise-constant sampling function by a term.
Consider, for example, a constant-valued environment map: with the
sampling technique, all values are equally likely to be
chosen. Due to the mapping to directions on the sphere, however, this
would lead to more directions being sampled near the poles of the sphere,
*not* a uniform sampling of directions on the sphere, which would be
a more desirable result. The term in the PDF corrects for
this non-uniform sampling of directions so that correct results are
computed in Monte Carlo estimates. Given this state of affairs, however,
it’s better to have modified the sampling distribution so that
it’s less likely to select directions near the poles in the first place.

The method can finally initialize the `VisibilityTester` with a light
sample point outside the scene’s bounds and return the radiance value for
the chosen direction.

The `InfiniteAreaLight::Pdf_Li()` method needs to convert the direction
to the corresponding coordinates in the sampling
distribution. Given these, the PDF is computed as the product of
the two 1D PDFs by the `Distribution2D::Pdf()` method, which is
adjusted here for mapping to the sphere as was done in the
`Sample_Li()` method.