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 Lights 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 SpotLights, ProjectionLights, GonioPhotometricLights, and DistantLights 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.

Shapes 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.

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.

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.

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 Disks, 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.

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.

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.

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.

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:

(See Figure 14.9(c).) Solving for , we have

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.

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 Shapes, 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 InfiniteAreaLights 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 InfiniteAreaLights 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:

1. 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.
2. We apply the sampling method from Section 13.6.7 to transform 2D samples to samples drawn from the piecewise-constant distribution.
3. 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.

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. ^†

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.

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.