## 15.2 Sampling Volume Scattering

Before proceeding to algorithms that model the effect of light scattering in participating media, we’ll first define some building-block functionality for sampling from distributions related to participating media and for computing the beam transmittance for spatially varying media.

The `Medium` interface defines a `Sample()` method, which takes a
world space ray and possibly samples a medium scattering
interaction along it. The input ray will generally have been intersected
against the scene geometry; thus, implementations of this method shouldn’t ever
sample a medium interaction at a point on the ray beyond its
value. Without loss of generality, the following discussion assumes that there
is always a surface at some distance .

The objective of this method is to sample the integral form of the equation of transfer, Equation (15.2), which consists of a surface and a medium-related term:

where is the point on the surface. We will neglect the effect of medium emission and assume directionally constant medium properties, in which case the source term is given by

Two cases can occur: if `Sample()` doesn’t sample an
interaction on the given ray interval , then
the surface-related term
should be estimated. If
it does sample an interaction, the
second integral term is to be estimated, and the provided `MediumInteraction`
should be initialized accordingly.

Suppose that denotes the probability per unit distance of generating an interaction at position . Due to the possibility of not sampling a medium interaction, this function generally doesn’t integrate to 1, and we define as the associated discrete probability of sampling the surface term:

With these definitions, we can now specify the semantics of `Sample()`, which
differs from previously encountered techniques for scattering functions
like `BSDF::Sample_f()` in that it does not provide the caller with
separate information about the function value and PDF at the
sampled position. This information is not generally needed, and some medium
models (specifically, the heterogeneous medium) admit more efficient sampling
schemes when it is possible to compute *ratios* of these quantities
instead.

When the surface term is selected, the method should return a weight equal to

which corresponds to sampling the first summand. Note that the value of the outgoing radiance is not included in ; it is the responsibility of the caller to account for this term. In the medium case, the method returns

which corresponds to sampling all medium-related terms except for the integral over in-scattered light in Equation (15.6), which must be handled separately.

The scattering coefficient and transmittance allow for spectral variation,
hence this method returns a `Spectrum`-valued weighting factor to
update the path throughput weight up to the surface or
medium scattering event.

As is generally the case for Monte Carlo integration, estimators like and admit a variety of sampling techniques that all produce the desired distribution. The implementation of the heterogeneous medium will make use of this fact to provide an implementation that is considerably more efficient than the canonical sampling approach based on the inversion method.

So that calling code can easily determine whether the provided
`MediumInteraction` was initialized by `Sample()`,
`MediumInteraction` provides an `IsValid()` method that takes
advantage of the fact that any time a medium scattering event has been
sampled, the phase function pointer will be set.

### 15.2.1 Homogeneous Medium

The `HomogeneousMedium` implementation of this method is fairly
straightforward; the only complexities come from needing to handle
attenuation coefficients that vary by wavelength.

In Section 13.3.1 we derived the sampling method for an exponential distribution defined over . For , it is

with PDF

However, the attenuation coefficient in general varies by wavelength. It is not desirable to sample multiple points in the medium, so a uniform sample is first used to select a spectral channel ; the corresponding scalar value is then used to sample a distance along the distribution

using the technique from Equation (15.9). The resulting sampling density is the average of the individual strategies :

The (discrete) probability of sampling a surface interaction at is the complement of generating a medium scattering event between and . This works out to a probability equal to the average transmittance over all spectral channels:

The implementation draws a sample according to Equation (15.11);
if the sampled distance is before the ray–primitive intersection (if any),
then a medium scattering event is recorded by initializing the
`MediumInteraction`. Otherwise, the sampled point in the medium is
ignored, and corresponding surface interaction should be used as the next
path vertex by the integrator. This sampling
approach is naturally efficient: the probability of generating a medium
interaction instead of the surface interaction is exactly equal to 1
minus the beam transmittance for the selected wavelength. Thus, given
optically thin media (or a short ray extent), the surface interaction is
more often used, and for thick media (or longer rays), a
medium interaction is more likely to be sampled.

In either case, the beam transmittance `Tr` is easily computed using
Beer’s law, Equation (11.3), just as in the
`HomogeneousMedium::Tr()` method.

Finally, the method computes the sample density using
Equations (15.11)
or (15.12) and returns the resulting sampling weight
or , depending on the value of
`sampledMedium`.

### 15.2.2 Heterogeneous Medium

In the case of the `GridDensityMedium`, extra effort is necessary to
deal with the medium’s heterogeneous nature. When the spatial variation can be
decomposed into uniform regions (e.g., piecewise constant voxels), a technique
known as *regular tracking* applies standard homogeneous medium techniques
to the voxels individually; a disadvantage of this approach is that it becomes
costly when there are many voxels. Since the `GridDensityMedium`
relies on linear interpolation, this approach cannot be used.

Other techniques build on a straightforward generalization of the homogeneous sampling PDF from Equation (15.10) with a spatially varying attenuation coefficient:

where evaluates the attenuation at distance
along the ray. The most commonly used method for importance sampling
Equation (15.13), is known as *ray marching*. This method inverts an
approximate cumulative distribution by partitioning the range
into a number of subintervals, numerically approximating
the integral in each interval, and finally inverting this discrete
representation. Unfortunately discretizing the problem in this way introduces
systemic statistical bias, which means that an `Integrator` using ray
marching generally won’t converge to the right result (even when an
infinite number of samples per pixel is used). Furthermore,
this bias can manifest itself in the
form of distracting visual artifacts.

For this reason, we prefer an alternative unbiased approach proposed by
Woodcock et al. (1965) that was originally developed to simulate volumetric
scattering of neutrons in atomic reactors. This technique is known as
*delta tracking* and is easiest to realize when the attenuation
coefficient is monochromatic. Our implementation includes an
assertion test (not shown here) to verify that this is indeed the case. Note
that the scattering and absorption coefficients are still permitted to vary
with respect to wavelength—however, their sum must
be uniform.

Figure 15.3 compares regular tracking, ray marching, and delta tracking. Delta tracking can be interpreted as filling the medium with additional (virtual) particles until its attenuation coefficient is constant everywhere. Sampling the resulting homogeneous medium is then easily accomplished using the basic exponential scheme from Equation (15.9). However, whenever an interaction with a particle occurs, it is still necessary to determine if it involved a “real” or a “virtual” particle (in which case the interaction is disregarded). The elegant insight of Woodcock et al. was that this decision can be made randomly based on the local fraction of “real” particles, which leads to a distribution of samples matching Equation (15.13).

The following fragment is part of the
`GridDensityMedium::GridDensityMedium()` constructor; its purpose is to
precompute the inverse of the maximum density scale factor over the entire
medium, which will be a useful quantity in the delta tracking implementation
discussed next.

`GridDensityMedium`>>=

The `Sample()` method begins by transforming the ray into the medium
coordinate system and normalizing the ray direction; `ray.tMax` is scaled
appropriately to account for the normalization.

`ray`’s overlap with medium bounds>>

`mi`with medium interaction information and return>>

Next, the implementation computes the parametric range of the ray’s overlap with the medium’s bounds, which are the unit cube . This step is technically not required for correct operation but is generally a good idea: reducing the length of the considered ray segment translates into a correspondingly smaller number of delta tracking iterations.

`ray`’s overlap with medium bounds>>=

Assuming that the maximum extinction value throughout the medium is given by , each delta-tracking iteration performs a standard exponential step through the uniform medium:

where . These steps are repeated until one of two stopping criteria is
satisfied: first, if then we have left the medium without an
interaction and `Medium::Sample()` hasn’t sampled a scattering event.
Alternatively, the loop may be terminated
at each iteration with probability , the local fraction of “real” particles.
This random decision consumes , the
second of two uniform samples per iteration .

`mi`with medium interaction information and return>>

The probability of not sampling a medium interaction is equal to the
transmittance of the ray segment ; hence
`1.0` is returned for the sampling weight
according to Equation (15.7). The medium interaction case
resembles the fragment <<Sample a channel and distance along the
ray>>.

`mi`with medium interaction information and return>>=

Finally, we must also provide an implementation of the `Tr()` method to
compute the transmittance along a ray segment. Consider the pseudocode of the
following simplistic implementation that performs a call to `Sample()` and
returns `1.0` if the ray passed through the segment and
`0.0` when a medium interaction occurred along the way. This effectively
turns the transmittance function into a binary random variable.

Since the probability of passing through the medium is equal to the
transmittance, this random variable has the correct mean and could be used in the
context of unbiased Monte Carlo integration. Calling `Tr()` many times and
averaging the result would produce an increasingly accurate estimate of the
transmittance, though this will generally be too costly to do in practice.
On the other hand, using the naive binary implementation leads to a high
amount of variance.

Novák et al. (2014) observed that this binary-valued function can be
interpreted as an instance of Russian roulette. However, instead of randomly
terminating the algorithm with a value of zero in each iteration, we could also
remove the Russian roulette logic and simply multiply the transmittance by the
probability of continuation. The resulting estimator has the same mean with a
considerably lower variance.
We will use this approach in the
implementation of `GridDensityMedium::Tr()`.

`ray`’s overlap with medium bounds>>

The beginning of the `Tr()` method matches `Sample()`. The loop body
is also identical except for the last line, which multiplies a running product
by the ratio of real particles to hypothetical particles.
(Novák referred to this scheme as *ratio tracking*).

### 15.2.3 Sampling Phase Functions

It is also useful to be able to draw samples from the distribution described by
phase functions—applications include applying multiple importance
sampling to computing direct lighting in participating media as well as for
sampling scattered directions for indirect lighting samples in
participating media. For these applications, `PhaseFunction`
implementations must implement the `Sample_p()` method, which samples
an incident direction given the outgoing direction and
a sample value in .

Note that, unlike the `BxDF` sampling methods, `Sample_p()`
doesn’t return both the phase function’s value and its PDF. Rather, `pbrt` assumes that phase functions are sampled with PDFs that perfectly match
their distributions. In conjunction with the requirement that phase
functions themselves be normalized
(Equation (11.4)), a single return value
encodes both values. When the value of the PDF alone is needed, a call to
`PhaseFunction::p()` suffices.

The PDF for the Henyey–Greenstein phase function is separable into and components, with as usual. The main task is to sample .

`wi`for Henyey–Greenstein sample>>

For Henyey–Greenstein, the distribution for is

if ; otherwise, gives a uniform sampling over the sphere of directions.

Given , what should now be a familiar approach converts them to the direction .

`wi`for Henyey–Greenstein sample>>=