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.

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.

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.

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 .

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

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().

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 .

For Henyey–Greenstein with pbrt’s convention for the orientation of direction vectors, 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 .