## 15.4 Sampling Subsurface Reflection Functions

We’ll now implement techniques to sample the subsurface scattering equation
introduced in Section 5.6.2, building on the `BSSRDF`
interface introduced in Section 11.4. Our task is to
estimate

Figure 15.6 suggests the complexity of evaluating the integral. To compute the standard Monte Carlo estimate of this equation given a point at which to compute outgoing radiance, we need a technique to sample points on the surface and to compute the incident radiance at these points, as well as an efficient way to compute the specific value of the BSSRDF for each sampled point and incident direction.

The `VolPathIntegrator` could be used to evaluate the BSSRDF:
given a pair of points on the surface and a pair of directions, the
integrator can be used to compute the fraction of incident light from
direction at the point that exits the object at the point
in direction by following light-carrying paths through the
multiple scattering events in the medium. Beyond standard path-tracing or
bidirectional path-tracing techniques, many other light transport
algorithms are applicable to this task.

However, many translucent objects are characterized by having very high albedos, which are not efficiently handled by classic approaches. For example, Jensen et al. (2001b) measured the scattering properties of skim milk and found an albedo of . When essentially all of the light is scattered at each interaction in the medium and almost none of it is absorbed, light easily travels far from where it first enters the medium. Hundreds or even thousands of scattering events must be considered to compute an accurate result; given the high albedo of milk, after 100 scattering events, of the incident light is still carried by a path, after 500 scattering events, and still after 1000.

`BSSRDF` class implementations represent the aggregate scattering behavior
of these sorts of media,
making it possible to render them fairly efficiently.
Figure 15.7 shows an
example of the dragon model rendered with a BSSRDF.
The main sampling operation that must be provided by implementations of the
`BSSRDF` interface, `BSSRDF::Sample_S()`, determines the surface
position where a ray re-emerges following internal scattering.

The value of the BSSRDF for the two points and directions is returned
directly, and the associated surface interaction record and probability
density are returned via the `si` and `pdf` parameters. Two
samples must be provided: a 1D sample for discrete sampling decisions
(e.g., choosing a specific spectral channel of the profile) and a 2D sample
that is mapped onto `si`. As we will see shortly, it’s useful for
BSSRDF implementations to be able to trace rays against the scene geometry
to find `si`, so the scene is also provided as an argument.

### 15.4.1 Sampling the SeparableBSSRDF

Recall the simplifying assumption introduced in Section 11.4.1, which factored the BSSRDF into spatial and directional components that can be sampled independently from one another. Specifically, Equation (11.6) defined as a product of a single spatial term and a pair of directional terms related to the incident and outgoing directions.

The spatial term was further simplified to a radial profile function :

We will now explain how each of these factors is handled by the
`SeparableBSSRDF`’s sampling routines. This class implements an
abstract sampling interface
that works for any radial profile function . The `TabulatedBSSRDF`
class, discussed in Section 15.4.2, derives from
`SeparableBSSRDF` and provides a specific tabulated representation of this
profile with support for efficient evaluation and exact importance sampling.

Returning to Equation (15.14), if we assume that the BSSRDF is only sampled for rays that are transmitted through the surface boundary, where transmission is selected with probability , then nothing needs to be done for the part here. (This is the case for the fragment <<Account for subsurface scattering, if applicable>>.) This is a reasonable expectation to place on calling code, as this approach gives good Monte Carlo efficiency.

This leaves the and terms—the former is handled by a call to
`SeparableBSSRDF::Sample_Sp()` (to be discussed shortly), which returns the
position `si`.

If sampling a position is successful, the method initializes
`si->bsdf` with an instance of the class `SeparableBSSRDFAdapter`,
which represents the directional term as a `BxDF`. Although
this `BxDF` does not truly depend on the outgoing direction
`si->wo`, we still need to initialize it with a dummy direction.

The `SeparableBSSRDFAdapter` class is a thin wrapper around
`SeparableBSSRDF::Sw()`. Recall that from
Equation (11.7) was defined as a
diffuse-like term scaled by the normalized Fresnel transmission. For this reason, the
`SeparableBSSRDFAdapter` classifies itself as `BSDF_DIFFUSE` and
just uses the default cosine-weighted sampling routine provided by
`BxDF::Sample_f()`.

Similar to refractive BSDFs, a scaling factor related to the light
transport mode must be applied to the value the `f()` method returns
for the term. This issue is discussed in more detail in
Section 16.1, and the fragment that applies this
scaling, <<Update BSSRDF transmission term to account for adjoint
light transport>>, is defined there.

To sample the spatial component , we need a way of mapping a 2D distribution function onto an arbitrary surface using a parameterization of the surface in the neighborhood of the outgoing position. A conceptually straightforward way to obtain such a parameterization is by means of geodesics, but finding and evaluating them is non-trivial and requires significant implementation effort for each shape that is supported. We use a much simpler approach that uses ray tracing to map the radial profile onto the scene geometry.

Figure 15.8 illustrates the basic idea: the
position and associated normal define a planar approximation
to the surface. Using 2D polar coordinates, we first sample an azimuth and a
radius value centered around and then map this position onto the
actual surface by intersecting an offset perpendicular ray with the
primitive, producing the position .
The `SeparableBSSRDF` class only supports radially symmetric profile
functions; hence is drawn from a uniform distribution on , and
is distributed according to the radial profile function .

There are still several difficulties with this basic approach:

- The radial profile is not necessarily uniform across wavelengths—in practice, the mean free path can differ by orders of magnitude between different spectral channels.
- If the surface geometry is poorly approximated by a plane and , where is the surface normal at , the probe rays will hit the surface at a grazing angle so that positions with comparatively high values of may be sampled with too low a probability. The result is high variance in renderings (Figure 15.9).
- Finally, the probe ray may intersect multiple surface locations along its length, all of which may contribute to reflected radiance.

The first two problems can be addressed with a familiar approach; namely, by introducing additional tailored sampling distributions and combining them using multiple importance sampling. The third will be addressed shortly.

We use a different sampling technique per wavelength to deal with
spectral variation, and each technique is additionally replicated three
times with different projection axes given by the basis vectors of a local
frame, resulting in a total of `3 * Spectrum::nSamples`
sampling techniques. This ensures that every point where takes on non-negligible
values is intersected with a reasonable probability. This combination of
techniques is implemented in `SeparableBSSRDF::Sample_Sp()`.

`ss`>>

`ts`>>

`IntersectionChain`and linked list>>

`IntersectionChain`>>

We begin by choosing a projection axis.
Note that when the surface is close to planar, projecting
along the normal `SeparableBSSRDF::ns` is clearly the best sampling
strategy, as probe rays along the other two axes are likely to miss the
surface. We therefore allocate a fairly large portion (50%) of the sample
budget to perpendicular rays. The other half is equally shared between
tangential projections along `SeparableBSSRDF::ss` and
`SeparableBSSRDF::ts`. The three axes of the chosen coordinate system
are stored in `vx`, `vy`, and `vz`, and follow our usual
convention of measuring angles in spherical coordinates with
respect to the axis.

After this discrete sampling operation, we scale and offset `u1`
so that additional sampling operations can reuse it as a uniform variate.

`ss`>>

`ts`>>

The fragments for the other two axes are similar and therefore not included here.

Next, we uniformly choose a spectral channel and re-scale `u1` once more.

The 2D profile sampling operation is then carried out in polar coordinates
using the `SeparableBSSRDF::Sample_Sr()` method. This method returns a
negative radius to indicate a failure (e.g., when there is no scattering
from channel `ch`); the implementation here returns a BSSRDF value of 0
in this case.

Both the radius sampling method `SeparableBSSRDF::Sample_Sr()` and its
associated density function `SeparableBSSRDF::Pdf_Sr()` are declared
as pure virtual functions; an implementation for
`TabulatedBSSRDF` is presented in the next section.

Because the profile falls off fairly quickly, we are not interested in
positions that are too far
away from . In order to
reduce the computational expense of the ray-tracing step, the probe ray is
clamped to a sphere of radius around . Another call
to `SeparableBSSRDF::Sample_Sr()` is used to determine
. Assuming that this function implements a perfect
importance sampling scheme based on the inversion method
(Section 13.3.1), `Sample_Sr()` maps a sample value
to the radius of a sphere containing a fraction of the
scattered energy.

Here, we set so that the sphere from Figure 15.8 contains 99.9% of the scattered energy. When lies outside of , sampling fails—this helps to keep the probe rays short, which significantly improves the run-time performance. Given and , the length of the intersection of the probe ray with the sphere of radius is

(See Figure 15.10.)

Given the sampled polar coordinate value, we can compute the world space
origin of a ray that lies on the boundary of the sphere and a target point,
`pTarget`, where it exits the sphere.

In practice, there could be more than just one intersection along the probe ray, and we want to collect all of them here. We’ll create a linked list of all of the found interactions.

`IntersectionChain`and linked list>>

`IntersectionChain`>>

`IntersectionChain` lets us maintain this list. Once again, the
`MemoryArena` makes it possible to efficiently perform allocations, here
for the list nodes.

`IntersectionChain`and linked list>>=

We now start by finding intersections along the segment within the sphere.
The list’s tail node’s `SurfaceInteraction` is initialized with each
intersection’s information, and `base` `Interaction` is updated so that the next ray
can be spawned on the other side of the intersected surface. (See
Figure 15.11.)

`SeparableBSSRDF::Sample_Sp()`method finds all of the intersections of a ray with the surface of the primitive, where ray extents are limited to a sphere around the intersection point (red dot). At each intersection (blue dots), the corresponding

`SurfaceInteraction`is stored in a linked list before a new ray leaving the other side of the intersected surface is generated.

`IntersectionChain`>>

When tracing rays to sample nearby points on the surface of the primitive,
it’s important to ignore any intersections on other primitives in the
scene. (There is an implicit assumption that scattering between primitives
will be handled by the integrator, and the BSSRDF should be limited to
account for single primitives’ scattering.) The implementation here
uses equality of `Material` pointers as a proxy to determine if an
intersection is on the same primitive. Valid intersections are appended to
the chain, and the variable `nFound` records their total count when
the loop terminates.

`IntersectionChain`>>=

With the set of intersections at hand, we must now choose one of them, as
`Sample_Sp()` can only return a single position . The following
fragment uses the variable `u1` one last time to pick one of the list
entries with uniform probability.

Finally, we can call `SeparableBSSRDF::Pdf_Sp()` (to be defined
shortly) to evaluate the combined PDF that takes all of the sampling
strategies into account. The probability it returns is divided by
`nFound` to account for the discrete probability of selecting
`pi` from the `IntersectionChain`. Finally, the value of
is returned.

`SeparableBSSRDF::Pdf_Sp()` returns the probability per unit area of
sampling the position `pi` with the total of `3 *
Spectrum::nSamples` sampling techniques available to
`SeparableBSSRDF::Sample_Sp()`.

First, `nLocal` is initialized with the surface normal at and
`dLocal` with the difference vector , both expressed using
local coordinates at .

To determine the combined PDF, we must query the probability of sampling a
radial profile radius matching the pair for each
technique. This radius is measured in 2D and thus depends on the chosen
projection axis (Figure 15.12). The
`rProj` variable records radii for projections perpendicular to
`ss`, `ts`, and `ns`.

The remainder of the implementation simply loops over all combinations of spectral channels and projection axes and sums up the product of the probability of selecting each technique and its area density under projection onto the surface at .

The alert reader may have noticed a slight inconsistency in the above
definitions: the probability of choosing one of several (`nFound`)
intersections in `SeparableBSSRDF::Sample_Sp()` should really have
been part of the density function computed in the
`SeparableBSSRDF::Pdf_Sp()` method rather
than the ad hoc division that occurs in the fragment <<Compute
sample PDF and return the spatial BSSRDF term >>. In practice, the
number of detected intersections varies with respect of the projection axis
and spectral channel; correctly accounting for this in the PDF
computation requires counting the number of intersections along each of a
total of `3 * Spectrum::nSamples` probe rays for every sample!
We neglect this issue, trading a more efficient implementation for a small
amount of bias.

### 15.4.2 Sampling the TabulatedBSSRDF

The previous section completed the discussion of BSSRDF sampling with the
exception of the
`Pdf_Sr()` and `Sample_Sr()` methods that were declared as pure
virtual functions in
the `SeparableBSSRDF` interface. The `TabulatedBSSRDF` subclass
implements this missing functionality.

The `TabulatedBSSRDF::Sample_Sr()` method samples radius values proportional
to the radial profile function . Recall from
Section 11.4.2 that the profile has an implicit dependence
on the albedo at the current surface position and that the
`TabulatedBSSRDF` provides interpolated evaluations of using
2D tensor product spline basis functions. `TabulatedBSSRDF::Sample_Sr()`
then determines the albedo for the given spectral channel `ch`
and draws samples proportional to the remaining 1D function .
Sample generation fails if there is neither scattering nor absorption on
channel `ch` (this case is indicated by returning a negative radius).

As in Section 11.4.2, there is a considerable amount of
overlap with the `FourierBSDF` implementation. The sampling
operation here actually reduces to a single call to `SampleCatmullRom2D()`, which was
previously used in `FourierBSDF::Sample_f()`.

Recall that this function depends on a precomputed CDF array, which is
initialized when the `BSSRDFTable` is created.

The `Pdf_Sr()` method returns the PDF of samples obtained via
`Sample_Sr()`. It evaluates the profile function divided by the
normalizing constant defined in
Equation (11.11).

The beginning is analogous to the spline evaluation code in
`TabulatedBSSRDF::Sr()`. The fragment <<Compute spline weights
to interpolate BSSRDF density on channel `ch`>> matches
<<Compute spline weights to interpolate BSSRDF on channel
`ch`>> in that method except that this method immediately returns zero
if the optical radius is outside the range represented by the spline.

`ch`>>

`ch`>>

The remainder of the implementation is very similar to fragment
<<Set BSSRDF value `Sr[ch]` using tensor spline interpolation>>
except that here, we also interpolate from the
tabulation and include it in the division at the end.

`ch`>>=

### 15.4.3 Subsurface Scattering in the Path Tracer

We now have the capability to apply Monte Carlo integration to the generalized scattering equation, (5.11), from Section 5.6.2. We will compute estimates of the form

where represents incident direct radiance and
is incident indirect radiance. The sample is generated in two steps.
First, given and , a call to
`BSSRDF::Sample_S()` returns a position whose distribution
is similar to the marginal distribution of with respect to .

Next, we sample the incident direction . Recall that the
`BSSRDF::Sample_S()` interface intentionally keeps these two steps
apart: instead of generating both and at the same time, it
returns a special `BSDF` instance via the `bsdf` field of
`si` that is used for the direction sampling step. No generality
is lost with such an approach: the returned `BSDF` can be completely
arbitrary and is explicitly allowed to depend on information computed
within `BSSRDF::Sample_S()`. The benefit is that we can re-use a
considerable amount of existing infrastructure for computing product
integrals of a BSDF and .

The `PathIntegrator`’s <<Find next path vertex and accumulate
contribution>> fragment invokes the following code to compute this estimate.

When the BSSRDF sampling case is triggered, the path tracer begins by calling
`BSSRDF::Sample_S()` to generate and incorporates the resulting
sampling weight into its throughput weight variable `beta` upon success.

Because `BSSRDF::Sample_S()` also initializes `pi`’s `bsdf`
with a BSDF that characterizes the dependence of on , we can reuse
the existing infrastructure for direct illumination computations.
Only a single line of code is necessary to compute the contribution
of direct lighting at `pi` to reflected radiance at `po`.

Similarly, accounting for indirect illumination at the newly sampled
incident point is almost the same as the BSDF indirect illumination
computation in the `PathIntegrator` except that `pi` is used for
the next path vertex instead of `isect`.

With this, the path tracer (and the volumetric path tracer) support subsurface scattering. See Figure 15.13 for an example.

`PathIntegrator`. These dragons both have BSSRDFs that describe subsurface scattering in their interiors. (

*Model courtesy of Christian Schüller.*)