## 14.2 Volume Scattering Integrators

The path space expression of the null-scattering equation of transfer allows a variety of sampling techniques to be applied to the light transport problem. This section defines two integrators that are based on path tracing starting from the camera.

First is the `SimpleVolPathIntegrator`, which uses simple sampling
techniques, giving an implementation that is short and easily verified.
This integrator is particularly useful for computing ground-truth results
when debugging more sophisticated volumetric sampling and integration algorithms.

The `VolPathIntegrator` is defined next. This integrator is fairly
complex, but it applies state-of-the-art sampling techniques to volume
light transport while handling surface scattering similarly to the
`PathIntegrator`. It is `pbrt`’s default integrator and is also the
template for the wavefront integrator in Chapter 15.

### 14.2.1 A Simple Volumetric Integrator

The `SimpleVolPathIntegrator` implements a basic volumetric path
tracer, following the sampling approach described in
Section 14.1.2.
Its `Li()` method is under 100 lines of
code, none of them too tricky. However, with this simplicity comes a
number of limitations. First, like the `RandomWalkIntegrator`, it does
not perform any explicit light sampling, so it requires that rays are able
to randomly intersect the lights in the scene. Second, it does
not handle scattering from surfaces. An error message is therefore issued
if it is used with a
scene that contains delta distribution light sources or has surfaces with
nonzero-valued BSDFs. (These defects are all addressed in the
`VolPathIntegrator` discussed in
Section 14.2.2.)
Nevertheless, this integrator is capable of
rendering complex volumetric effects; see Figure 14.6.

`SimpleVolPathIntegrator`. With 256 samples per pixel, this integrator gives a reasonably accurate rendering of the volumetric model, though there are variance spikes in some pixels (especially visible toward the bottom of the volume) due to error from the integrator not directly sampling the scene’s light sources. The

`VolPathIntegrator`, which uses more sophisticated sampling strategies, renders this scene with 1,288 times lower MSE; it is discussed in Section 14.2.2.

*(Scene courtesy of Jim Price.)*

This integrator’s only parameter is the maximum path length, which is set via a value passed to the constructor (not included here).

The general form of the `Li()` method follows that of the
`PathIntegrator`.

`RNG`for sampling the majorant transmittance>>

A few familiar variables track the path state, including `L` to
accumulate the radiance estimate for the path. For this integrator,
`beta`, which tracks the path throughput weight, is just a single
`Float` value, since the product of ratios of phase function values and
sampling PDFs from Equation (14.19) is a
scalar value.

Media with scattering properties that vary according to wavelength
introduce a number of complexities in sampling and evaluating Monte Carlo
estimators. We will defer addressing them until we cover the
`VolPathIntegrator`. The `SimpleVolPathIntegrator` instead
estimates radiance at a single wavelength by terminating all but the first wavelength
sample.

Here is a case where we have chosen simplicity over efficiency for this integrator’s implementation: we might instead have accounted for all wavelengths until the point that spectrally varying scattering properties were encountered, enjoying the variance reduction benefits of estimating all of them for scenes where doing so is possible. However, doing this would have led to a more complex integrator implementation.

The first step in the loop is to find the ray’s intersection with the scene geometry, if any. This gives the parametric distance beyond which no samples should be taken for the current ray, as the intersection either represents a transition to a different medium or a surface that occludes farther-away points.

The `scattered` and `terminated` variables declared here will allow the
lambda function that is passed to
`SampleT_maj()` to report
back the state of the path after sampling terminates.

`RNG`for sampling the majorant transmittance>>

An `RNG` is required for the call to the `SampleT_maj()`
function. We derive seeds for it based on two random values from the
sampler, hashing them to convert `Float`s into integers.

`RNG`for sampling the majorant transmittance>>=

With that, a call to `SampleT_maj()` starts the generation of samples
according to
. The `Sampler` is used to generate the first uniform
sample `u` that is passed to the method; recall from
Section 14.1.3 that subsequent ones will be
generated using the provided `RNG`. In a similar fashion, the
`Sampler` is used for the initial value of `uMode` here. It will be
used to choose among the three types of scattering event at the first
sampled point. For `uMode` as well, the `RNG` will provide
subsequent values.

In this case, the transmittance that `SampleT_maj()` returns for the
final segment is unneeded, so it is ignored.

For each sample returned by `SampleT_maj()`, it is necessary
to select which type of scattering it represents. The first step is to
compute the probability of each possibility. Because we have specified
such that it is nonnegative and , the null-scattering probability can be found as one minus the
other two probabilities. A call to `std::max()` ensures that any
slightly negative values due to floating-point round-off error are clamped at
zero.

A call to `SampleDiscrete()` then selects one of the three terms of
using the specified probabilities.

If absorption is chosen, the path terminates. Any emission is added to the
radiance estimate, and evaluation of
Equation (14.19) is complete. The
fragment therefore sets `terminated` to indicate that the path is
finished and returns `false` from the lambda function so that no further
samples are generated along the ray.

For a scattering event, `beta` is updated according to the ratio of phase
function and its directional sampling probability from
Equation (14.19).

The counter for the number of scattering events is only incremented for real-scattering events; we do not want the number of null-scattering events to affect path termination. If this scattering event causes the limit to be reached, the path is terminated.

If the path is not terminated, then a new direction is sampled from the phase function’s distribution.

Given a sampled direction, the `beta` factor must be updated.
Volumetric path-tracing
implementations often assume that the phase function sampling
distribution matches the phase function’s actual distribution and dispense
with `beta` entirely since it is always equal to 1. This variation is worth pausing to consider: in
that case, emitted radiance at the end of the path is always returned,
unscaled. All of the effect of transmittance, phase functions, and so forth
is entirely encapsulated in the distribution of how often
various terms are evaluated and in the distribution of scattered ray
directions.
`pbrt` does not impose the requirement on phase functions that their
importance sampling technique be perfect, though this is the case for the
Henyey–Greenstein phase function in `pbrt`.

Be it with `beta` or without, there is no need to do any further work
along the current ray after a scattering event, so after
the following code updates the path state to account for
scattering, it too returns `false` to direct that no further
samples should be taken along the ray.

Null-scattering events are ignored, so there is nothing to do but to return
`true` to indicate that additional samples along the
current ray should be taken. Similar to the real-scattering case, this can be
interpreted as starting a recursive evaluation of
Equation (14.3) from the current sampled position without
incurring the
overhead of actually doing so. Since this is the only case that may lead
to another invocation of the lambda function, `uMode` must be
refreshed with a new uniform sample value in case another sample is
generated.

If the path was terminated due to absorption,
then there is no more work to do in the `Li()`
method; the final
radiance value can be returned. Further, if the ray was scattered, then
there is nothing more to do but to restart the `while` loop and start
sampling the scattered ray. Otherwise, the ray either underwent no
scattering events or only underwent null scattering.

If the ray is unscattered and unabsorbed, then any emitters it interacts with contribute radiance to the path. Either surface emission or emission from infinite light sources is accounted for, depending on whether an intersection with a surface was found. Further, if the ray did not intersect a surface, then the path is finished and the radiance estimate can be returned.

It is still necessary to consider surface intersections, even if scattering from them is not handled by this integrator. There are three cases to consider:

- If the surface has no BSDF, it represents a transition between different
types of participating media. A call to
`SkipIntersection()`moves the ray past the intersection and updates its medium appropriately. - If there is a valid
`BSDF`and that`BDSF`also returns a valid sample from`Sample_f()`, then we have a BSDF that scatters; an error is issued and rendering stops. - A valid but zero-valued BSDF is allowed; such a BSDF should be assigned to area light sources in scenes to be rendered using this integrator.

### 14.2.2 Improving the Sampling Techniques

The `VolPathIntegrator`
adds three significant improvements to the approach
implemented in `SimpleVolPathIntegrator`: it supports scattering from surfaces
as well as from volumes; it
handles spectrally varying medium scattering properties without falling
back to sampling a single wavelength;
and it samples lights directly, using multiple
importance sampling to reduce variance when doing so.
The first improvement—including surface scattering—is mostly a matter
of applying the ideas of
Equation (14.7), sampling distances
in volumes but then choosing surface scattering if the sampled distance is
past the closest intersection. For the other two, we will here discuss
the underlying foundations before turning to their implementation.

#### Chromatic Media

We have thus far glossed over some of the implications of spectrally varying medium
properties. Because `pbrt` uses point-sampled spectra, they introduce no
complications in terms of evaluating things
like the modified path throughput or the path throughput
weight : given a set of path
vertices, such quantities can be evaluated for all the wavelength
samples simultaneously using the `SampledSpectrum` class.

The problem with spectrally varying medium properties comes from sampling. Consider a wavelength-dependent function that we would like to integrate at wavelengths . If we draw samples from a wavelength-dependent PDF based on the first wavelength and then evaluate at all the wavelengths, we have the estimators

Even if the PDF that was used for sampling matches well, it may be a poor match for at the other wavelengths. It may not even be a valid PDF for them, if it is zero-valued where the function is nonzero. However, falling back to integrating a single wavelength at a time would be unfortunately inefficient, as shown in Section 4.5.4.

This problem of a single sampling PDF possibly mismatched with a wavelength-dependent function comes up repeatedly in volumetric path tracing. For example, sampling the majorant transmittance at one wavelength may be a poor approach for sampling it at others. That could be handled by selecting a majorant that bounds all wavelengths’ extinction coefficients, but such a majorant would lead to many null-scattering events at wavelengths that could have used a much lower majorant, which would harm performance.

The path tracer’s choice among absorption, real scattering, and null scattering at a sampled point cannot be sidestepped in a similar way: different wavelengths may have quite different probabilities for each of these types of medium interaction, yet with path tracing the integrator must choose only one of them. Splitting up the computation to let each wavelength choose individually would be nearly as inefficient as only considering a single wavelength at a time.

However, if a single type of interaction is chosen based on a single wavelength and we evaluate the modified path contribution function for all wavelengths, we could have arbitrarily high variance in the other wavelengths. To see why, note how all the factors that came from the factors in Equation (14.18) canceled out to give the delta-tracking estimator, Equation (14.19). In the spectral case, if, for example, real scattering is chosen based on a wavelength ’s scattering coefficient and if a wavelength has scattering coefficient , then the final estimator for will include a factor of that can be arbitrarily large.

The fact that `SampleT_maj()` nevertheless samples according to a single
wavelength’s majorant transmittance suggests that there is a solution to
this problem. That solution, yet again, is multiple importance sampling.
In this case, we are using a single sampling technique rather than
MIS-weighting multiple techniques, so we use the single-sample MIS
estimator from Equation (2.16), which here gives

where is the discrete probability of sampling using the wavelength , here uniform at with the number of spectral samples.

The balance heuristic is optimal for single-sample MIS. It gives the MIS weight

which gives the estimator

See Figure 14.7 for an example that shows the benefits of MIS for chromatic media.

`VolPathIntegrator`. For this scene, MSE is reduced by a factor of 149.

*(Scene courtesty of Jim Price.)*

#### Direct Lighting

Multiple importance sampling is also at the heart of how the
`VolPathIntegrator` samples direct
illumination. As with the `PathIntegrator`, we would like to combine
the strategies of sampling the light sources with sampling the BSDF or
phase function to find light-carrying paths and
then to weight the contributions of each sampling technique using MIS.
Doing so is more complex than it is in the absence of volumetric
scattering, however, because not only does the sampling distribution used
at the last path vertex differ (as before) but the `VolPathIntegrator`
also uses ratio
tracking to estimate the transmittance along the shadow ray. That is a
different distance sampling technique than the delta-tracking approach used when
sampling ray paths, and so it leads to a different path PDF.

In the following, we will say that the two path-sampling techniques used in
the `VolPathIntegrator` are *unidirectional path sampling* and
*light path sampling*; we will write their respective path PDFs as
and . The first corresponds to the sampling
approach from Section 14.1.5, with
delta tracking used to find real-scattering vertices and with the phase function
or BSDF sampled to find the new direction at each vertex. Light path
sampling follows the same approach up to the last real-scattering vertex
before the light vertex; there, the light selects the direction and then
ratio tracking gives the transmittance along the last path segment. (See
Figure 14.8.)
Given a path , both approaches share the same path
throughput weight up to the vertex and the same
path PDF up to that vertex, .

`VolPathIntegrator`uses ratio tracking to compute the transmittance along the ray by accumulating the product at sampled points along the ray (open circles). For the MIS weight, it is necessary to be able not only to compute the PDF for sampling the corresponding direction at the last path vertex but also to compute the probability of generating these samples using delta tracking, since that is how the path would be sampled with unidirectional path sampling.

For the full PDF for unidirectional path sampling, at the last scattering vertex we have the probability of scattering, times the directional probability for sampling the new direction , which is given by the sampling strategy used for the BSDF or phase function. Then, for the path to find an emitter at the vertex , it must have only sampled null-scattering vertices between and ; absorption or a real-scattering vertex preclude making it to .

Using the results from Section 14.1.5, we can find that the path PDF between two points and with intermediate null-scattering vertices indexed by is given by the product of

for all null-scattering vertices. The factors cancel and the null-scattering path probability is

The full unidirectional path probability is then given by

For light sampling, we again have the discrete probability for scattering at but the directional PDF at the vertex is determined by the light’s sampling distribution, which we will denote by . The only missing piece is the PDF of the last segment (the shadow ray), where ratio tracking is used. In that case, points are sampled according to the majorant transmittance and so the PDF for a path sampled between points and with intermediate vertices is

and the full light sampling path PDF is given by

The `VolPathIntegrator` samples both types of paths according to the
first wavelength but evaluates these PDFs at all wavelengths so
that MIS over wavelengths can be used. Given the path sampled
using unidirectional path sampling and then the path sampled
using light path sampling, the two-sample MIS estimator is

Note that because the paths share the same vertices for all of , not only do the two factors share common factors, but and do as well, following Equations (14.21) and (14.23).

In this case, the MIS weights can account not only for the differences between unidirectional and light path sampling but also for the different per-wavelength probabilities for each sampling strategy. For example, with the balance heuristic, the MIS weight for the unidirectional strategy works out to be

with the number of spectral samples. The MIS weight for light sampling is equivalent, but with the function in the numerator replaced with .

### 14.2.3 Improved Volumetric Integrator

`VolPathIntegrator`.

*(Scene courtesy of Jim Price.)*

`VolPathIntegrator`.

*(Scene courtesy of Angelo Ferretti.)*

The `VolPathIntegrator` pulls together all of these ideas to robustly
handle both surface and volume transport. See
Figures 14.9 and 14.10 for images
rendered with this integrator that show off the visual complexity that
comes from volumetric emission, chromatic media, and multiple scattering in
participating media.

As with the other `Integrator` constructors that we have seen so far,
the `VolPathIntegrator` constructor does not perform any meaningful
computation, but just initializes member variables with provided values.
These three are all equivalent to their parallels in the
`PathIntegrator`.

`RNG`for sampling the majorant transmittance>>

`L`for medium emission>> }

`beta`and

`r_u`for real-scattering event>>

`visibleSurf`at first intersection>>

`ucRho`and

`uRho`for reflectance estimate>>

`beta`and rescaled path probabilities for BSDF scattering>>

`BSSRDFSample`>>

There is a common factor of in the denominator of the first term of the two-sample MIS estimator, Equation (14.24), and the numerator of the MIS weights, Equation (14.25). There is a corresponding factor in the second term of the estimator and in the weight. It is tempting to cancel these out; in that case, the path state to be tracked by the integrator would consist of and the wavelength-dependent probabilities and . Doing so is mathematically valid and would provide all the quantities necessary to evaluate Equation (14.24), but suffers from the problem that the quantities involved may overflow or underflow the range of representable floating-point values.

To understand the problem, consider a highly specular surface—the BSDF will have a large
value for directions around its peak, but the PDF for sampling those
directions will also be large. That causes no problems in the
`PathIntegrator`, since its `beta` variable tracks their ratio,
which ends up being close to 1. However, with
maintained independently, a series of specular bounces could lead to
overflow. (Many null-scattering events along a path can cause similar
problems.)

Therefore, the `VolPathIntegrator` tracks the path throughput weight for the
sampled path

which is numerically well behaved. Directly tracking the probabilities
and would also stress the
range of floating-point numbers, so instead it tracks the *rescaled
path probabilities*

where is the probability for sampling the current path. It is equal to the light path probability for paths that end with a shadow ray from light path sampling and the unidirectional path probability otherwise. (Later in the implementation, we will take advantage of the fact that these two probabilities are the same until the last scattering vertex, which in turn implies that whichever of them is chosen for does not affect the values of and .)

These rescaled path probabilities are all easily incrementally updated during path sampling. If , then MIS weights like those in Equation (14.25) can be found with

and similarly for when .

The remaining variables in the following fragment have the same function as
the variables of the same names in the `PathIntegrator`.

The `while` loop for each ray segment starts out similarly to the
corresponding loop in the `SimpleVolPathIntegrator`:
the integrator traces a ray to find a value at the closest surface
intersection before sampling the medium, if any, between the ray origin and
the intersection point.

`RNG`for sampling the majorant transmittance>>

`L`for medium emission>> }

`beta`and

`r_u`for real-scattering event>>

`visibleSurf`at first intersection>>

`ucRho`and

`uRho`for reflectance estimate>>

`beta`and rescaled path probabilities for BSDF scattering>>

`BSSRDFSample`>>

The form of the fragment for sampling the medium is similar as well:
`tMax` is set using the ray intersection , if available, and an
`RNG` is prepared before medium sampling proceeds. If the path is
terminated or undergoes real scattering in the medium, then no further work
is done to sample surface scattering at a ray intersection point.

`RNG`for sampling the majorant transmittance>>

`L`for medium emission>> }

`beta`and

`r_u`for real-scattering event>>

Given a sampled point in the medium, the lambda function’s task is to evaluate the source function, taking care of the second case of Equation (14.7).

`L`for medium emission>> }

`beta`and

`r_u`for real-scattering event>>

A small difference from the `SimpleVolPathIntegrator` is that
volumetric emission is added at every point that is sampled in the medium
rather than only when the absorption case is sampled. There is no reason
not to do so, since emission is already available via the
`MediumProperties` passed to the lambda function.

`L`for medium emission>> }

In the following, we will sometimes use the notation to denote the path with the vertex appended to it. Thus, for example, . The estimator that gives the contribution for volumetric emission at is then

`beta` holds , so we can incrementally compute
by

From Section 14.1.5, we know that . Because we are always sampling absorption (at least as far as including emission goes), is 1 here.

Even though this is the only sampling technique for volumetric emission,
different wavelengths may sample this vertex with different probabilities,
so it is worthwhile to apply MIS over the wavelengths’ probabilities. With
`r_u` storing the rescaled unidirectional probabilities up to the
previous path vertex, the rescaled path probabilities for sampling the
emissive vertex, `r_e`, can be found by multiplying `r_u` by
the per-wavelength probabilities and dividing by the
probability for the wavelength that was used for sampling , which is
already available in `pdf`. (Note that in monochromatic media, these
ratios are all 1.)

Here we have a single-sample MIS estimator with balance heuristic weights given by

The absorption coefficient and emitted radiance for evaluating
Equation (14.28) are available in
`MediumProperties` and the `SampledSpectrum::Average()` method
conveniently computes the average of rescaled probabilities in the
denominator of Equation (14.29).

`L`for medium emission>>=

Briefly returning to the initialization of `betap` and `r_e` in
the previous fragments: it may seem tempting to cancel out the
`T_maj` factors from them, but note how the final estimator does not
perform a component-wise division of these two quantities but instead
divides by the average of the rescaled probabilities when computing the MIS
weight. Thus, performing such cancellations would lead to incorrect
results.

After emission is handled, the next step is to determine which term of
to evaluate; this follows the same
approach as in the `SimpleVolPathIntegrator`.

`beta`and

`r_u`for real-scattering event>>

As before, the ray path is terminated in the event of absorption. Since any volumetric emission at the sampled point has already been added, there is nothing to do but handle the details associated with ending the path.

For a real-scattering event, a shadow ray is traced to a light to sample
direct lighting, and the path state is updated to account for the new ray.
A `false` value returned from the lambda function prevents further
sample generation along the current ray.

`beta`and

`r_u`for real-scattering event>>

The PDF for real scattering at this vertex is the product of the PDF for sampling its distance along the ray, , and the probability for sampling real scattering, . The values cancel.

Given the PDF value, `beta` can be updated to include
along the segment up to the new vertex divided by the PDF. The rescaled
probabilities are computed in the same way as the path sampling PDF before
being divided by it, following
Equation (14.26). The rescaled
light path probabilities will be set shortly, after a new ray direction is
sampled.

`beta`and

`r_u`for real-scattering event>>=

Direct lighting is handled by the `SampleLd()` method, which we will
defer until later in this section.

Sampling the phase function gives a new direction at the scattering event.

There is a bit of bookkeeping to take care of after a real-scattering
event. We can now incorporate the phase function value into `beta`,
which completes the contribution of from
Equation (14.12). Because both unidirectional path
sampling and light path sampling use the same set of sampling operations up
to a real-scattering vertex, an initial value for the rescaled light path
sampling probabilities `r_l` comes from the value of the rescaled
unidirectional probabilities before scattering. It is divided by the
directional PDF from for this vertex here. The
associated directional PDF for light sampling at this vertex will be
incorporated into `r_l` later. There is no need to update
`r_u` here, since the scattering direction’s probability is the same
for all wavelengths and so the update factor would always be 1.

At this point, the integrator also updates various variables that record the scattering history and updates the current ray.

If the ray intersects a light source, the
`LightSampleContext` from the previous path vertex will be needed to
compute MIS weights; `prevIntrContext` is updated to store it
after each scattering event, whether in a medium or on a surface.

If null scattering is selected, the updates to `beta` and the rescaled
path sampling probabilities follow the same form as we have seen
previously: the former is given by
Equation (14.11) and the latter with a
factor where, as with real scattering, the
cancels with a corresponding factor from the
probability (Section 14.1.5).

In this case, we also must update the rescaled path probabilities for sampling this path vertex via light path sampling, which samples path vertices according to .

This fragment concludes the implementation of the lambda function that is
passed to the `SampleT_maj()` function.

Returning to the `Li()` method immediately after the
`SampleT_maj()` call, if the path terminated due to absorption, it is
only here that we can break out and return the radiance estimate to the
caller of the `Li()` method. Further, it is only here that we can
jump back to the start of the `while` loop for rays that were
scattered in the medium.

With those cases taken care of, we are left with rays that either underwent no scattering events in the medium or only underwent null scattering. For those cases, both the path throughput weight and the rescaled path probabilities must be updated. takes a factor of to account for the transmittance from either the last null-scattering event or the ray’s origin to the ray’s position. The rescaled unidirectional and light sampling probabilities also take the same , which corresponds to the final factors outside of the parenthesis in the definitions of and .

There is much more to do for rays that have either escaped the scene or have intersected a surface without medium scattering or absorption. We will defer discussion of the first following fragment, <<Add emitted light at volume path vertex or from the environment>>, until later in the section when we discuss the direct lighting calculation. A few of the others are implemented reusing fragments from earlier integrators.

`visibleSurf`at first intersection>>

`ucRho`and

`uRho`for reflectance estimate>>

`beta`and rescaled path probabilities for BSDF scattering>>

`BSSRDFSample`>>

As with the `PathIntegrator`, path termination due to reaching the
maximum depth only occurs after accounting for illumination from any
emissive surfaces that are intersected.

Sampling the light source at a surface intersection is handled by the same
`SampleLd()` method that is called for real-scattering vertices in the
medium. As with medium scattering, the `LightSampleContext`
corresponding to this scattering event is recorded for possible later use
in MIS weight calculations.

The logic for sampling scattering at a surface is very similar to the
corresponding logic in the `PathIntegrator`.

`beta`and rescaled path probabilities for BSDF scattering>>

Given a BSDF sample, is first multiplied by the value of the BSDF, which takes care of from Equation (14.12). This is also a good time to incorporate the cosine factor from the factor of the generalized geometric term, Equation (14.13).

Updates to the rescaled path probabilities follow how they were done for
medium scattering: first, there is no need to update `r_u` since the
probabilities are the same over all wavelengths. The rescaled light path
sampling probabilities are also initialized from `r_u`, here also
with only the factor included. The other
factors in `r_l` will only be computed and included if the ray
intersects an emitter; otherwise `r_l` is unused.

One nit in updating `r_l` is that the BSDF and PDF value
returned in the `BSDFSample` may only be correct up to a (common) scale
factor. This case comes up with sampling techniques like the random walk
used by the `LayeredBxDF` that is described in
Section 14.3.2. In that case, a call to
`BSDF::PDF()` gives an independent value for the PDF that can be used.

`beta`and rescaled path probabilities for BSDF scattering>>=

A few additional state variables must be updated after surface scattering, as well.

Russian roulette follows the same general approach as before, though we
scale `beta` by
the accumulated effect of radiance scaling for transmission that is encoded in
`etaScale` and use the balance heuristic over wavelengths.
If the Russian roulette test passes, `beta` is updated
with a factor that accounts for the survival probability, `1 - q`.

#### Estimating Direct Illumination

All that remains in the `VolPathIntegrator`’s implementation is direct
illumination. We will start with the `SampleLd()` method, which is
called to estimate scattered radiance due to direct illumination by
sampling a light source, both at scattering points in media and on
surfaces. (It is responsible for computing the second term of
Equation (14.24).)
The purpose of most of its parameters should be evident. The last,
`r_p`, gives the rescaled path probabilities up to the vertex
`intr`. (A separate variable named `r_u` will be used in
the function’s implementation, so a new name is needed here.)

`intr`>>

`LightSampleContext`for volumetric light sampling>>

`lightSampler`>>

`f_hat`and

`scatterPDF`accounting for the BSDF>> } else { <<Update

`f_hat`and

`scatterPDF`accounting for the phase function>>

`T_ray`and PDFs using ratio-tracking estimator>>

The overall structure of this method’s implementation is similar to the
`PathIntegrator`’s
`SampleLd()` method: a light
source and a point on it are sampled, the vertex’s scattering function is
evaluated, and then the light’s visibility is determined. Here we have the
added complexity of needing to compute the transmittance between the
scattering point and the point on the light rather than finding a binary visibility
factor, as well as the need to compute spectral path sampling weights for
MIS.

`intr`>>=

`LightSampleContext`for volumetric light sampling>>

`lightSampler`>>

`f_hat`and

`scatterPDF`accounting for the BSDF>> } else { <<Update

`f_hat`and

`scatterPDF`accounting for the phase function>>

`T_ray`and PDFs using ratio-tracking estimator>>

Because it is called for both surface and volumetric scattering path
vertices, `SampleLd()` takes a plain `Interaction` to represent
the scattering point. Some extra care is therefore needed when
initializing the `LightSampleContext`: if scattering is from a surface,
it is important to interpret that interaction as the
`SurfaceInteraction` that it is so that the shading normal is included
in the `LightSampleContext`. This case also presents an opportunity,
as was done in the `PathIntegrator`, to
shift the light sampling point to avoid incorrectly sampling
self-illumination from area lights.

`LightSampleContext`for volumetric light sampling>>=

Sampling a point on the light follows in the usual way. Note that the
implementation is careful to consume the two sample dimensions from the
`Sampler` regardless of whether sampling a light was successful, in
order to keep the association of sampler dimensions with integration
dimensions fixed across pixel samples.

`lightSampler`>>=

The light samples a direction from the reference point in the usual manner.
The `true` value passed for the `allowIncompletePDF` parameter of
`Light::SampleLi()` indicates the use of MIS here.

As in `PathIntegrator::SampleLd()`, it is worthwhile to evaluate the BSDF
or phase function before tracing the shadow ray: if it turns out to be
zero-valued for the direction to the light source, then it is possible to
exit early and perform no further work.

`f_hat`and

`scatterPDF`accounting for the BSDF>> } else { <<Update

`f_hat`and

`scatterPDF`accounting for the phase function>>

The `f_hat` variable that holds the value of the scattering function
is slightly misnamed: it also includes the cosine factor for scattering
from surfaces and does not include the for scattering from
participating media, as that has already been included in the provided
value of `beta`.

`f_hat`and

`scatterPDF`accounting for the BSDF>>=

`f_hat`and

`scatterPDF`accounting for the phase function>>=

A handful of variables keep track of some useful quantities for the
ray-tracing and medium sampling operations that are performed to compute
transmittance. `T_ray` holds the transmittance along the ray
and `r_u` and `r_l` respectively hold
the rescaled path probabilities for unidirectional sampling and light sampling,
though only along the ray. Maintaining these values independently of the
full path contribution and PDFs facilitates the use of Russian roulette in
the transmittance computation.

`SampleLd()` successively intersects the shadow ray with the scene
geometry, returning zero contribution if an opaque surface is found and
otherwise sampling the medium to estimate the transmittance up to the
intersection. For intersections that represent transitions between
different media, this process repeats until the ray reaches the light
source.

For some scenes, it could be more efficient to instead first check that there are no intersections with opaque surfaces before sampling the media to compute the transmittance. With the current implementation, it is possible to do wasted work estimating transmittance before finding an opaque surface farther along the ray.

`T_ray`and PDFs using ratio-tracking estimator>>

If an intersection is found with a surface that has a non-`nullptr`
`Material`, the visibility term is zero and the method can return
immediately.

Otherwise, if participating media is present, `SampleT_maj()` is
called to sample it along the ray up to whichever is closer—the surface
intersection or the sampled point on the light.

`T_ray`and PDFs using ratio-tracking estimator>>

For each sampled point in the medium, the transmittance and rescaled path probabilities are updated before Russian roulette is considered.

`T_ray`and PDFs using ratio-tracking estimator>>

In the context of the equation of transfer, using ratio tracking to compute transmittance can be seen as sampling distances along the ray according to the majorant transmittance and then only including the null-scattering component of the source function to correct any underestimate of transmittance from . Because only null-scattering vertices are sampled along transmittance rays, the logic for updating the transmittance and rescaled path probabilities at each vertex exactly follows that in the <<Handle null scattering along ray path>> fragment.

`T_ray`and PDFs using ratio-tracking estimator>>=

Russian roulette is used to randomly terminate rays with low transmittance. A natural choice might seem to be setting the survival probability equal to the transmittance—along the lines of how Russian roulette is used for terminating ray paths from the camera according to . However, doing so would effectively transform ratio tracking to delta tracking, with the transmittance always equal to zero or one. The implementation therefore applies a less aggressive termination probability, only to highly attenuated rays.

In the computation of the transmittance value used for the Russian roulette test, note that an MIS weight that accounts for both the unidirectional and light sampling strategies is used, along the lines of Equation (14.27).

After the `SampleT_maj()` call returns, the transmittance and
rescaled path probabilities
all must be multiplied by the `T_maj` returned from
`SampleT_maj()` for the final ray segment. (See the discussion for
the earlier <<Handle terminated, scattered, and unscattered medium rays>>
fragment for why each is updated as it is.)

If the transmittance is zero (e.g., due to Russian roulette termination), it
is possible to return immediately. Furthermore, if there is no surface
intersection, then there is no further medium sampling to be done and we
can move on to computing the scattered radiance from the light.
Alternatively, if there is an intersection, it must be with a surface that
represents the boundary between two media; the `SpawnRayTo()` method
call returns the continuation ray on the other side of the surface, with
its `medium` member variable set appropriately.

After the `while` loop terminates, we can compute the final rescaled
path probabilities, compute MIS weights, and return the final estimate of the path
contribution function for the light sample.

The `r_p` variable passed in to `SampleLd()` stores the rescaled
path probabilities for unidirectional sampling of the path up to the vertex
where direct lighting is being computed—though here, `r_u` and
`r_l` have been rescaled using the light path sampling probability,
since that is how the vertices were sampled along the shadow ray. However,
recall from Equations (14.21)
and (14.23) that
for
the path up to the scattering vertex. Thus, `r_p` can be interpreted
as being rescaled using . This allows
multiplying `r_l` and `r_u` by `r_p` to compute final
rescaled path probabilities.

If the light source is described by a delta distribution, only the light sampling technique is applicable; there is no chance of intersecting such a light via sampling the BSDF or phase function. In that case, we still apply MIS using all the wavelengths’ individual path PDFs in order to reduce variance in chromatic media.

For area lights, we are able to use both light source and the scattering function samples, giving two primary sampling strategies, each of which has a separate weight for each wavelength.

With `SampleLd()` implemented, we will return to the fragments in the
`Li()` method that handle the cases where a ray escapes from the
scene and possibly finds illumination from infinite area lights, as well as where a
ray intersects an emissive surface. These handle the first term in
the direct lighting MIS estimator, Equation (14.24).

As with the `PathIntegrator`, if the previous
scattering event was due to a delta-distribution scattering function, then
sampling the light is not a useful strategy. In that case, the MIS weight
is only based on the per-wavelength PDFs for the unidirectional sampling
strategy.

Otherwise, the MIS weight should account for both sampling techniques. At
this point, `r_l` has everything but the probabilities for sampling
the light itself. (Recall that we deferred that when initializing
`r_l` at the real-scattering vertex earlier.) After incorporating
that factor, all that is left is to compute the final weight, accounting for
both sampling strategies.

The work done in the <<Add contribution of emission from intersected surface>> fragment is very similar to that done for infinite lights, so it is not included here.