## 12.6 Light Sampling

Due to the linearity assumption in radiometry, illumination at a point in
a scene with multiple light sources can be computed by summing the independent
contributions of each light. As we have seen before, however, correctness
alone is not always sufficient—if it were, we might have sampled
`ImageInfiniteLight`s uniformly with the suggestion that one take
thousands of samples per pixel until error has been reduced
sufficiently.
Especially in scenes with thousands or more independent light sources,
considering all of them carries too much of a performance cost.

Fortunately, here, too, is a setting where stochastic sampling can be applied. An unbiased estimator for a sum of terms is given by

where the probability mass function (PMF) for any term where is nonzero and where . This is the discrete analog to the integral Monte Carlo estimator, Equation (2.7). Therefore, we can replace any sum over all the scene’s light sources with a sum over just one or a few of them, where the contributions are weighted by one over the probability of sampling the ones selected.

Figure 12.23 is a rendered image of a scene with 8,878 light sources. A few observations motivate some of the light sampling algorithms to come. At any given point in the scene, some lights are facing away and others are occluded. Ideally, such lights would be given a zero sampling probability. Furthermore, often many lights are both far away from a given point and have relatively low power; such lights should have a low probability of being sampled. (Consider, for example, the small yellow lights inset in the machinery.) Of course, even a small and dim light is important to points close to it. Therefore, the most effective light sampling probabilities will vary across the scene depending on position, surface normal, the BSDF, and so forth.

*Zero Day*Scene, with 8,878 Light Sources. It is infeasible to consider every light when computing reflected radiance at a point on a surface, and therefore light sampling methods from this section are necessary to render this scene efficiently.

*(Scene courtesy of Beeple.)*

The `LightSampler` class defines the
`LightSampler` interface for
sampling light sources.
It is defined in the file
`base/lightsampler.h`. `LightSampler` implementations can
be found in `lightsamplers.h` and
`lightsamplers.cpp`.

The key `LightSampler` method is `Sample()`, which takes a uniform
1D sample and information about a reference point in the form of a
`LightSampleContext`. When sampling is successful, a
`SampledLight` is returned. Otherwise, the optional value is left
unset, as may happen if the light sampler is able to determine that no
lights illuminate the provided point.

`SampledLight` just wraps up a light and the discrete probability
for it having been sampled.

In order to compute the MIS weighting term when a ray happens to hit a
light source, it is necessary to be able to find the value of the
probability mass function for sampling a particular light. This task is
handled by `PMF()` method implementations.

`LightSampler`s must also provide methods to sample a light and
return the corresponding probability independent
of a specific point being illuminated. These
methods are useful for light transport algorithms like
bidirectional path tracing that start paths at the light sources.

### 12.6.1 Uniform Light Sampling

`UniformLightSampler` is the simplest possible light sampler: it
samples all lights with uniform probability. In practice, more
sophisticated sampling algorithms are usually much more effective, but this
one is easy to implement and provides a useful baseline for comparing light
sampling techniques.

As with all light samplers, an array of all the lights in the scene is
provided to the constructor; `UniformLightSampler` makes a copy of
them in a member variable.

Since the light sampling probabilities do not depend on the lookup point,
we will only include the variants of `Sample()` and `PMF()`
that do not take a `LightSampleContext` here. The versions of these methods
that do take a context just call these variants. For sampling, a light is chosen by
scaling the provided uniform sample by the array size and returning the
corresponding light.

Given uniform sampling probabilities, the value of the PMF is always one over the number of lights.

### 12.6.2 Power Light Sampler

`PowerLightSampler` sets the probability for sampling each light
according to its power. Doing so generally gives better results than
sampling uniformly, but the lack of spatial variation in sampling
probabilities limits its effectiveness. (We will return to this topic at the
end of this section where some comparisons between the two techniques are
presented.)

Its constructor also makes a copy of the provided lights but initializes some additional data structures as well.

`lightToIndex`hash table>>

To efficiently return the value of the PMF for a given light, it is necessary
to be able to find the index in the
`lights` array of a given light. Therefore, the constructor also initializes a hash
table that maps from `Light`s to indices.

`lightToIndex`hash table>>=

The `PowerLightSampler` uses an `AliasTable` for sampling. It is
initialized here with weights based on the emitted power returned by each
light’s `Phi()` method.
Note that if the light’s emission distribution is spiky (e.g., as with many
fluorescent lights), there is a risk of underestimating its power if a spike
is missed. We have not found this to be a problem in practice, however.

Given the alias table, sampling is easy.

To evaluate the PMF, the hash table gives the mapping to an index in the array of lights. In turn, the PMF returned by the alias table for the corresponding entry is the probability of sampling the light.

As with the `UniformLightSampler`, the `Sample()` and `PMF()`
methods that do take a `LightSampleContext` just call the corresponding
methods that do not take one.

Sampling lights based on their power usually works well.
Figure 12.24
compares both sampling methods using the *Zero Day* scene. For this scene, noise
is visibly reduced when sampling according to power, and mean squared error
(MSE) is improved by a factor of 12.4.

*Zero Day*Scene. (a) Rendered with uniform light sampling. (b) Rendered with lights sampled according to power. Both images are rendered with 16 samples per pixel and rendering time is nearly the same. Sampling lights according to power reduces MSE by a factor of 12.4 for this scene.

*(Scene courtesy of Beeple.)*

Although sampling according to power generally works well, it is not optimal. Like uniform sampling, it is hindered by not taking the geometry of emitters and the relationship between emitters and the receiving point into account. Relatively dim light sources may make the greatest visual contribution in a scene, especially if the bright ones are far away, mostly occluded, or not visible at all.

As an extreme example of this problem with sampling according to power, consider a large triangular light source that emits a small amount of radiance. The triangle’s emitted power can be made arbitrarily large by scaling it to increase its total area. However, at any point in the scene the triangle can do no more than subtend a hemisphere, which limits its effect on the total incident radiance at a point. Sampling by power can devote far too many samples to such lights.

### 12.6.3 BVH Light Sampling

Varying the light sampling probabilities based on the point being shaded can be an effective light sampling strategy, though if there are more than a handful of lights, some sort of data structure is necessary to do this without having to consider every light at each point being shaded. One widely used approach is to construct a hierarchy over the light sources with the effect of multiple lights aggregated into the higher nodes of the hierarchy. This representation makes it possible to traverse the hierarchy to find the important lights near a given point.

Given a good light hierarchy, it is possible to render scenes with hundreds
of thousands or even millions of light sources nearly as efficiently as a
scene with just one light. In this section, we will describe the
implementation of the `BVHLightSampler`, which applies bounding volume
hierarchies to this task.

#### Bounding Lights

When bounding volume hierarchies (BVHs) were used
for intersection acceleration structures in
Section 7.3, it was necessary to abstract away the details of the
various types of primitives and underlying shapes that they stored so that
the `BVHAggregate` did not have to be explicitly aware of each of them.
There, the primitives’ rendering-space bounding boxes were used for
building the hierarchy. Although there were cases where the quality of the
acceleration structure might have been improved using shape-specific information
(e.g., if the acceleration structure was aware of skinny diagonal triangles
with large bounding boxes with
respect to the triangle’s area), the `BVHAggregate`’s implementation
was substantially simplified with that approach.

We would like to have a similar decoupling for the `BVHLightSampler`,
though it is less obvious what the right abstraction should be. For
example, we might want to know that a spotlight only emits light in a
particular cone, so that the sampler does not choose it for points outside
the cone. Similarly, we might want to know that a one-sided area light
only shines light on one side of a particular plane. For all sorts of
lights, knowing their total power would be helpful so that brighter lights
can be sampled preferentially to dimmer ones. Of course, power does not
tell the whole story, as the light’s spatial extent and relationship to a
receiving point affect how much power is potentially received at that point.

The `LightBounds` structure provides the abstraction used by `pbrt` for these purposes. It stores a variety of values that make it possible to
represent the emission distribution of a variety of types of light.

It is evident that the spatial bounds of the light and its emitted power
will be useful quantities, so those are included in `LightBounds`.
However, this representation excludes light sources at infinity such as the
`DistantLight` and the various infinite lights. That limitation is
fine, however, since it is unclear how such light sources would be stored in a
BVH anyway. (The `BVHLightSampler` therefore handles these types of lights
separately.)

As suggested earlier, bounding a light’s directional emission distribution is important for sampling lights effectively. The representation used here is based on a unit vector that specifies a principal direction for the emitter’s surface normal and two angles that specify its variation. First, specifies the maximum deviation of the emitter’s surface normal from the principal normal direction . Second, specifies the angle beyond up to which there may be emission (see Figure 12.25). Thus, directions that make an angle up to with may receive illumination from a light and those that make a greater angle certainly do not.

While this representation may seem overly specialized for emissive shapes
alone, it works well for all of `pbrt`’s (noninfinite) light types. For
example, a point light can be represented with an arbitrary average normal
and an angle of for . A spotlight can use
the direction it is facing for , its central cone angle for
, and the angular width of its falloff region for
.

Our implementation stores the cosine of these angles rather than the angles themselves; this representation will make it possible to avoid the expense of evaluating a number of trigonometric functions in the following.

The last part of the emission bounds for a light is a `twoSided` flag,
which indicates whether the direction should be negated to specify a
second cone that uses the same pair of angles.

The `LightBounds` constructor takes corresponding parameters and
initializes the member variables. The implementation is straightforward,
and so is not included here.

To cluster lights into a hierarchy, we will need to be able to
find the bounds that encompass two specified `LightBounds` objects.
This capability is provided by the `Union()` function.

`LightBounds`has zero power, return the other>> <<Find average direction and updated angles for

`LightBounds`>>

`LightBounds`union>> }

It is worthwhile to start out by checking for the easy case of one or the
other specified `LightBounds` having zero power. In this case, the
other can be returned immediately.

`LightBounds`has zero power, return the other>>=

Otherwise, a new average normal direction and updated angles
and must be computed. Most of the
work involved is handled by the `DirectionCone`’s
`Union()` method, which takes a
pair of cones of directions and returns one that bounds the two of them.
The cosine of the new angle is then given by the cosine
of the spread angle of that cone.

The value of should be the maximum of the
values for the two cones. However, because
`LightBounds` stores the cosines of the angles and because the cosine
function is monotonically decreasing over the range of possible
values, , we take the minimum cosine of the two angles.

`LightBounds`>>=

The remainder of the parameter values for the `LightBounds`
constructor are easily computed from the two `LightBounds` that were provided.

`LightBounds`union>>=

A utility method returns the centroid of the spatial bounds; this value will be handy when building the light BVH.

The `Importance()` method provides the key `LightBounds`
functionality: it returns a scalar value that estimates the contribution of
the light or lights represented in the `LightBounds` at a given point.
If the provided normal is nondegenerate, it is assumed to be the surface
normal at the receiving point. Scattering points in participating media
pass a zero-valued `Normal3f`.

`w`, >>

It is necessary to make a number of assumptions in order to estimate the
amount of light arriving at a point given a `LightBounds`. While it
will be possible to make use of principles such as the received power
falling off with the squared distance from the emitter or the incident
irradiance at a surface varying according to Lambert’s law, some
approximations are inevitable, given the loss of specifics that comes with
adopting the `LightBounds` representation.

`w`, >>

Even computing the squared distance for the falloff of received power is
challenging if `bounds` is not degenerate: to which point in
`bounds` should the distance be computed? It may seem that finding
the minimum distance from the point `p` to the bounds would be a safe
choice, though this would imply a very small distance for a point close to
the bounds and a zero distance for a point inside it. Either of these
would lead to a very large factor and potentially high error due to
giving too much preference to such a `LightBounds`. Further, choosing
between the two child `LightBounds` of a node when a point is inside
both would be impossible, given infinite weights for each.

Therefore, our first fudge is to compute the distance from the center of the bounding box but further to ensure that the squared distance is not too small with respect to the length of the diagonal of the bounding box. Thus, for larger bounding boxes with corresponding uncertainty about what the actual spatial distribution of emission is, the inverse squared distance factor cannot become too large.

In the following computations, we will need to produce a series of values of the form and . Given the sine and cosine of and , it is possible to do so without evaluating any trigonometric functions. To see how, consider the cosine: implies that and that therefore . We thus start by checking that case and returning if it is true. We are otherwise left with , which can be expressed in terms of the sines and cosines of the two angles using a trigonometric identity, . The case for sine follows analogously.

Two simple lambda functions provide these capabilities. (Only the one for
cosine is included in the text, as `sinSubClamped` follows a similar
form.)

There are a number of angles involved in the importance computation. In addition to the ones that specify the directional emission bounds, and , we will start by computing the sine and cosine of , the angle between the principal normal direction and the vector from the center of the light bounds to the reference point (Figure 12.26(a)).

`w`, >>=

To bound the variation of various angles across the extent of the
bounding box, we will also make use of the angle that the bounding box
subtends with respect to the reference point (see
Figure 12.26(b)). We
will denote this angle . The preexisting
`DirectionCone::BoundSubtendedDirections()` function takes care of
computing its cosine. Its sine follows directly.

The last angle we will use is the minimum angle between the emitter’s normal and the vector to the reference point. We will denote it by , and it is given by

see Figure 12.27. As with the other angles, we only need its sine and cosine, which can be computed one subtraction at a time.

If this angle is greater than (or, here, if its cosine is less than ), then it is certain that the lights represented by the bounds do not illuminate the reference point and an importance value of 0 can be returned immediately.

The importance value can now be computed. It starts with the product of the power of the lights, the factor that accounts for the cosine at the emitter, and the inverse squared distance.

At a surface, the importance also accounts for a conservative estimate of
the incident cosine factor there. We have `wi`, the unit vector from
the reference point to the center of the `LightBounds`, but would like
to conservatively set the importance based on the maximum value of the
incident cosine over the entire bounding box. The corresponding minimum
angle with the surface normal is given by (see Figure 12.28).

Our implementation of this computation uses the `cosSubClamped()`
lambda function introduced earlier to compute the cosine of the angle
directly using the sines and cosines of the two
contributing angles.

#### Bounds for Light Implementations

Given the definition of `LightBounds`, we will add another method to
the `Light` interface to allow lights to return bounds on their
emission.

Lights at infinity return an unset `optional` value. Here, for
example, is the implementation of this method for `ImageInfiniteLight`. The
other infinite lights and the `DistantLight` do likewise.

The `PointLight`’s implementation is just a few lines of code, as
befitting the simplicity of that type of light source. The spatial bounds
are given by the light’s rendering space position and the total emitted power
is easily computed following the approach in `PointLight::Phi()`.
Because this light shines in all directions, the average normal direction
is arbitrary and is , corresponding to the full
sphere of directions.

The `SpotLight`’s bounding method is a bit more interesting: now the average
normal vector is relevant; it is set here to be the light’s direction. The
range is set to be the angular width of the inner cone
of the light and corresponds to the width of its falloff
at the edges. While this falloff does not exactly match the cosine-weighted
falloff used in the `LightBounds::Importance()` method, it is
close enough for these purposes.

There is a subtlety in the computation of `phi` for this light: it is
computed as if the light source was an isotropic point source and is not
affected by the spotlight’s cone angle, like the computation in
`SpotLight::Phi()` is. To understand the reason for this, consider two
spotlights with the same radiant intensity, one with a very narrow cone
and one with a wide cone, both illuminating some point in the scene. The
total power emitted by the former is much less than the latter, though for
a point inside both of their cones, both should be sampled with equal
probability—effectively, the cone is accounted for in the light
importance function and so should not be included in the `phi` value
supplied here.

We will skip past the implementations of the `ProjectionLight` and
`GoniometricLight` `Bounds()` methods, which are along similar lines.

The `DiffuseAreaLight`’s `Bounds()` implementation is different
than the previous ones. The utility of the
`Shape::NormalBounds()` method may now be better appreciated; the
cone of directions that it returns gives the average normal direction
and its spread angle .
For area lights, , since illumination is emitted in
the entire hemisphere around each surface normal.

`phi`for diffuse area light bounds>>

`DiffuseAreaLight`image channel value>>

The `phi` value is found by integrating over the light’s area. For
lights that use an `Image` for spatially varying emission, the
<<Compute average `DiffuseAreaLight` image channel value>>
fragment, not included here, computes its average value. Because
`LightBounds` accounts for whether the emitter is one- or two-sided, it
is important not to double the shape’s area if it is two-sided; that factor
is already included in the importance computation. (This subtlety is
similar to the one for the `SpotLight`’s `phi` value.)
See Figure 12.29
for an illustration of how this detail makes a difference.

`DiffuseAreaLight::Bounds()`method includes an additional factor of 2 for the two-sided light’s importance, then it receives more samples than it should. (b) Without this factor, the light importance values are more accurate, which in turn gives a visible reduction in error. The MSE improvement is a factor of 1.42.

`phi`for diffuse area light bounds>>=

`DiffuseAreaLight`image channel value>>

#### Compactly Bounding Lights

The `LightBounds` class uses 52 bytes of storage. This is not a
problem as far as the total amount of memory consumed for the lights in the
scene, but it does affect
performance from the amount of space it uses in the cache. For scenes with
thousands of lights, multiple instances of the `LightBounds` will be
accessed when traversing the light BVH, and so minimizing its storage
requirements improves cache performance and thus overall performance.
(This is especially the case on the GPU, since many threads run
concurrently on each processor and each will generally follow a different
path through the light BVH and thus access different `LightBounds`
instances.)

Therefore, we have also implemented a `CompactLightBounds` class, which
applies a number of techniques to reduce storage requirements for the
`LightBounds` information. It uses just 24 bytes of storage. We use
both classes in `pbrt`: the uncompressed `LightBounds` is convenient for
lights to return from their `Bounds()` methods and is also a good
representation to use when building the light BVH.
`CompactLightBounds` is used solely in the in-memory representation of
light BVH nodes, where minimizing size is beneficial to performance.

`qb`>>

`w`, >>

Its constructor takes both a `LightBounds` instance and a bounding box
`allb` that must completely bound `LightBounds::bounds`.
This bounding box is used to compute quantized bounding box vertex
positions to reduce their storage requirements.

`qb`>>

The `OctahedralVector` class from Section 3.8.3
stores a unit vector in 4 bytes, saving 8 from the `Vector3` used in
`LightBounds`. Then, the two cosines and the `twoSided` flag are
packed into another 4 bytes using a bitfield, saving another 8. We have
left `phi` alone, since the various compactions already implemented
are sufficient for `pbrt`’s current requirements.

`QuantizeCos()` maps the provided value (which is expected to be the
cosine of an angle and thus between and 1) to a 15-bit unsigned
integer. After being remapped to be in , multiplying by the largest
representable 15-bit unsigned integer, , gives a value
that spans the valid range.

Note the use of `pstd::floor()` to round the quantized cosine value
down before returning it: this is preferable to, say, rounding to the
nearest integer, since it ensures that any quantization error serves to
slightly increase the corresponding angle rather than decreasing it.
(Decreasing it could lead to inadvertently determining that the light did
not illuminate a point that it actually did.)

The bounding box corners are also quantized. Each coordinate of each
corner gets 16 bits, all of them stored in the `qb` member variable.
This brings the storage for the bounds down to 12 bytes, from 24 before.
Here the quantization is also conservative, rounding down at the lower end
of the extent and rounding up at the upper end.

`qb`>>=

`QuantizeBounds()` remaps a coordinate value `c` between
`min` and `max` to the range , the range of values
that an unsigned 16-bit integer can store.

A few convenience methods make the values of various member variables
available. For the two quantized cosines, the inverse computation of
`QuantizeCos()` is performed.

The `Bounds()` method returns the `Bounds3f` for the
`CompactLightBounds`. It must be passed the same `Bounds3f` as was
originally passed to its constructor for the correct result to be returned.

Finally, `CompactLightBounds()` also provides an `Importance()`
method. Its implementation also requires that the original `Bounds3f`
be provided so that the `Bounds()` method can be called. Given the
unquantized bounds and cosines made available in appropriately named local
variables, the remainder of the implementation can share the same fragments
as were used in the implementation of `LightBounds::Importance()`.

`w`, >>

#### Light Bounding Volume Hierarchies

Given a way of bounding lights as well as a compact representation of these
bounds, we can turn to the implementation of the `BVHLightSampler`.
This light sampler is the default for most of the integrators in `pbrt`.
Not only is it effective at efficiently sampling among large collections of
lights, it even reduces error in simple scenes with just a few lights.
Figures 12.30
and 12.31 show two examples.

`PowerLightSampler`. (b) Rendered with 1 sample per pixel using the

`BVHLightSampler`. Even with a small number of lights, error is visibly lower with a sampler that uses spatially varying sampling probabilities due to being able to choose nearby bright lights with higher probability. In this case, MSE is improved by a factor of 2.72.

*Zero Day*Scene, with 8,878 Area Lights. (a) Rendered with the

`PowerLightSampler`. (b) Rendered with the

`BVHLightSampler`. Both images are rendered with 16 samples per pixel. For a scene of this complexity, an effective light sampling algorithm is crucial. The

`BVHLightSampler`gives an MSE improvement of 2.37 with only a 5.8% increase in rendering time. Monte Carlo efficiency is improved by a factor of 2.25.

*(Scene courtesy of Beeple.)*

`pInfinite`>>

`light`PMF computation>>

`pInfinite`>>

`bitTrail`to find next node index and update its value>>

`LightBounds`>>

`LightBounds`>>

Its constructor starts by making a copy of the provided array of lights before proceeding to initialize the BVH and additional data structures.

`infiniteLights`array and light BVH>>

`infiniteLights`or

`bvhLights`>>

Because the BVH cannot store lights at infinity, the first step is to partition the lights into those that can be stored in the BVH and those that cannot. This is handled by a loop over all the provided lights after which the BVH is constructed.

`infiniteLights`array and light BVH>>=

`infiniteLights`or

`bvhLights`>>

Lights that are not able to provide a `LightBounds` are added to the
`infiniteLights` array and are sampled independently of the lights
stored in the BVH. As long as they have nonzero emitted power, the rest
are added to the `bvhLights` array, which is used during BVH
construction. Along the way, a bounding box that encompasses all the
BVH lights’ bounding boxes is maintained in `allLightBounds`; this is
the bounding box that will be passed to the `CompactLightBounds` for
quantizing the spatial bounds of individual lights.

`infiniteLights`or

`bvhLights`>>=

The light BVH is represented using an instance of the `LightBVHNode`
structure for each tree node, both interior and leaf. It uses a total of
28 bytes of storage, adding just 4 to the 24 used by
`CompactLightBounds`, though its declaration specifies 32-byte
alignment, ensuring that 2 of them fit neatly into a typical 64-byte cache
line on a CPU, and 4 of them fit into a 128-byte GPU cache line.

Naturally, each `LightBVHNode` stores the `CompactLightBounds` for
either a single light or a collection of them. Like the
`BVHAggregate`’s BVH, the light BVH is laid out in memory so that the
first child of each interior node is immediately after it. Therefore, it
is only necessary to store information about the second child’s location in
the `LightBVHNode`. The `BVHLightSampler` stores all nodes in a
contiguous array, so an index suffices; a full pointer is unnecessary.

Two object-creation methods return a `LightBVHNode` of the specified
type.

The `buildBVH()` method constructs the BVH by recursively partitioning
the lights until it reaches a single light, at which point a leaf node is
constructed. Its implementation closely follows the approach implemented
in the `BVHAggregate::buildRecursive()` method: along each dimension,
the light bounds are assigned to a fixed number of buckets according to
their centroids. Next, a cost model is evaluated for splitting the
lights at each bucket boundary. The minimum cost split is chosen and
the lights are partitioned into two sets, each one passed to a recursive
invocation of `buildBVH()`.

Because these two methods are so similar, here we will only include the
fragments where the `BVHLightSampler` substantially diverges—in how
nodes are initialized and in the cost model used to evaluate candidate splits.

When this method is called with a range corresponding to a single light, a
leaf node is initialized and the recursion terminates. A
`CompactLightBounds` is created using the bounding box of all lights’
bounds to initialize its quantized bounding box coordinates and the BVH
tree node can be added to the `nodes` array.

In order to implement the `PMF()` method, it is necessary to
follow a path through the BVH from the root down to the leaf node for the
given light. We encode these paths using *bit trails*, integers
where each bit encodes which of the two child nodes should be visited at
each level of the tree in order to reach the light’s leaf node.
The lowest bit indicates
which child should be visited after the root node, and so forth. These
values are computed while the BVH is built and stored in a hash table that
uses `Light`s as keys.

When there are multiple lights to consider, the `EvaluateCost()`
method is called to evaluate the cost model for the two `LightBounds`
for each split candidate. In addition to the `LightBounds` for which to
compute the cost, it takes the bounding box of all the lights passed to
the current invocation of `buildBVH()` as well as the dimension in
which the split is being performed.

`LightBounds`>>

`LightBounds`>>

The principal surface normal direction and the angles
and that are stored in `LightBounds` are
worthwhile to include in the light BVH cost function. Doing so can
encourage partitions of primitives into groups that emit light in different
directions, which can be helpful for culling groups of lights
that do not illuminate a given receiving point.
To compute these costs, `pbrt` uses a weighted measure of the solid angle of directions
that the direction bounds subtend. A weight of 1 is used for all
directions inside the center cone up to and then the
remainder of directions up to are cosine-weighted,
following the importance computation earlier. (See
Figure 12.32.)
Integrating over the relevant directions gives us the direction bounds
measure,

The first term integrates to and the second has
a simple analytic form that is evaluated in the second term of
`M_omega`’s initializer below.

`LightBounds`>>=

Three other factors go into the full cost estimate:

- The power estimate
`phi`: in general, the higher the power of the lights in a`LightBounds`, the more important it is to minimize factors like the spatial and direction bounds. - A regularization factor
`Kr`that discourages long and thin bounding boxes. - The surface area of the
`LightBounds`’s bounding box.

The use of surface area in the cost metric deserves note: with the
`BVHAggregate`, the surface area heuristic was grounded in geometric
probability, as the surface area of a convex object is proportional to the
probability of a random ray intersecting it. In this case, no rays are
being traced. Arguably, minimizing the volume of the bounds would be a
more appropriate approach in this case. In practice, the surface area
seems to be more effective—one reason is that it penalizes bounding boxes
that span a large extent in two dimensions but little or none in the third.
Such bounding boxes are undesirable as they may subtend large solid angles,
adding more uncertainty to importance estimates.

`LightBounds`>>=

Once the lights have been partitioned, two fragments take care of recursively
initializing the child nodes and then initializing the interior node. The
first step is to take a spot in the `nodes` array for the interior
node; this spot must be claimed before the children are initialized in
order to ensure that the first child is the successor of the interior node
in the array.
Two recursive calls to `buildBVH()` then initialize the children. The main
thing to note in them is the maintenance of the `bitTrail` value
passed down into each one. For the first child, the corresponding bit
should be set to zero. `bitTrail` is zero-initialized in the initial
call to `buildBVH()`, so it has
this value already and there is nothing to do. For the second call, the
bit for the current tree depth is set to 1.

`LightBVHNode`and recursively initialize children>>=

The interior node can be initialized after the children have been. Its
light bounds are given by the union of its children’s, which allows
initializing a `CompactLightBounds` and then the `LightBVHNode`
itself.

Given the BVH, we can now implement the `Sample()` method, which samples a
light given a reference point in a `LightSampleContext`.

`pInfinite`>>

The first choice to make is whether an infinite light should be sampled or
whether the light BVH should be used to choose a noninfinite light.
The `BVHLightSampler` gives equal probability to sampling each infinite
light and to sampling the BVH, from which the probability of sampling an
infinite light follows directly.

`pInfinite`>>=

If an infinite light is to be sampled, then the random sample `u` is
rescaled to provide a new uniform random sample that is used to index
into the `infiniteLights` array.

Otherwise a light is sampled by traversing the BVH. At each interior node,
probabilities are found for sampling each of the two children using
importance values returned by the `LightBounds` for the reference point.
A child node is then randomly chosen according to these probabilities. In
the end, the probability of choosing a leaf node is equal to the product of
probabilities along the path from the root to the leaf (see
Figure 12.33).
With this traversal scheme, there is no need to maintain a stack of
nodes to be processed as the `BVHAggregate` does—a single path is
taken down the tree from the root to a leaf.

A few values will be handy in the following traversal of the BVH. Among
them are the uniform sample `u`, which is remapped to a new uniform
random sample in . `pmf` starts out with the probability of
sampling the BVH in the first place; each time a child node of the tree is
randomly sampled, it will be multiplied by the discrete probability of that
sampling choice so that in the end it stores the complete probability for
sampling the light.

At each interior node, a child node is randomly sampled. Given a leaf node, we have reached the sampled light.

The first step at an interior node is to compute the importance values for the two child nodes. It may be the case that both of them are 0, indicating that neither child illuminates the reference point. That we may end up in this situation may be surprising: in that case, why would we have chosen to visit this node in the first place? This state of affairs is a natural consequence of the accuracy of light bounds improving on the way down the tree, which makes it possible for them to better differentiate regions that their respective subtrees do and do not illuminate.

Given at least one nonzero importance value, `SampleDiscrete()` takes care
of choosing a child node. The sampling PMF it returns is
incorporated into the running `pmf` product. We further use its
capability for remapping the sample `u` to a new uniform sample in ,
which allows the reuse of the `u` variable in subsequent loop
iterations.

When a leaf node is reached, we have found a light. The light should only be returned if it has a nonzero importance value, however: if the importance is zero, then it is better to return no light than to return one and cause the caller to go through some additional work to sample it before finding that it has no contribution. Most of the time, we have already determined that the node’s light bounds have a nonzero importance value by dint of sampling the node while traversing the BVH in the first place. It is thus only in the case of a single-node BVH with a single light stored in it that this test must be performed here.

Computing the PMF for sampling a specified light follows a set of computations similar to those of the sampling method: if the light is an infinite light, the infinite light sampling probability is returned and otherwise the BVH is traversed to compute the light’s sampling probability. In this case, BVH traversal is not stochastic, but is specified by the bit trail for the given light, which encodes the path to the leaf node that stores it.

`light`PMF computation>>

`pInfinite`>>

`bitTrail`to find next node index and update its value>>

If the given light is not in the bit trail hash table, then it is not stored in the BVH and therefore must be an infinite light. The probability of sampling it is one over the total number of infinite lights plus one if there is a light BVH.

`light`PMF computation>>=

A number of values will be useful as the tree is traversed, including the bit trail that points the way to the correct leaf, the PMF of the path taken so far, and the index of the current node being visited, starting here at the root.

`pInfinite`>>

For a light that is stored in the BVH, the probability of sampling it is again computed as the product of each discrete probability of sampling the child node that leads to its leaf node.

`bitTrail`to find next node index and update its value>>

The lowest bit of `bitTrail` encodes which of the two children of the
node is visited on a path down to the light’s leaf node. In turn,
it is possible to compute the probability of sampling that node given
the two child nodes’ importance values.

The low-order bit of `bitTrail` also points us to which node to visit
next on the way down the tree. After `nodeIndex` is updated,
`bitTrail` is shifted right by one bit so that the low-order bit encodes the
choice to make at the next level of the tree.

`bitTrail`to find next node index and update its value>>=

The basic `Sample()` and `PMF()` methods for when a reference
point is not specified sample all the lights uniformly and so are not
included here, as they parallel the implementations in the
`UniformLightSampler`.