## 11.4 The BSSRDF

The bidirectional scattering-surface reflectance distribution function (BSSRDF) was introduced in Section 5.6.2; it gives exitant radiance at a point on a surface given incident differential irradiance at another point : . Accurately rendering translucent surfaces with subsurface scattering requires integrating over both area—points on the surface of the object being rendered—and incident direction, evaluating the BSSRDF and computing reflection with the subsurface scattering equation

Subsurface light transport is described by the volumetric scattering processes introduced in Sections 11.1 and 11.2 as well as the volume light transport equation that will be introduced in Section 15.1. The BSSRDF is a summarized representation modeling the outcome of these scattering processes between a given pair of points and directions on the boundary.

A variety of BSSRDF models have been developed to model subsurface reflection;
they generally include some simplifications to the underlying scattering
processes to make them tractable. One such model will be introduced in
Section 15.5. Here, we will begin by specifying a fairly abstract
interface analogous to `BSDF`.
All of the code related to BSSRDFs is in the files
`core/bssrdf.h` and
`core/bssrdf.cpp`.

`BSSRDF` implementations must pass the current (outgoing) surface
interaction as well as the index of refraction of the scattering medium to
the base class constructor. There is thus an implicit assumption that the
index of refraction is constant throughout the medium, which is a
widely used assumption in BSSRDF models.

The key method that `BSSRDF` implementations must provide is one that
evaluates the
eight-dimensional distribution function `S()`, which quantifies the ratio
of differential radiance at point in direction to the
incident differential flux at from direction
(Section 5.6.2). Since the and arguments are already
available via the `BSSRDF::po` and `Interaction::wo` fields, they
aren’t included in the method signature.

Like the `BSDF`, the `BSSRDF` interface also defines functions to
sample the distribution and to evaluate the probability density of the
implemented sampling scheme. The specifics of this part of the interface are
discussed in Section 15.4.

During the shading process, the current `Material`’s
`ComputeScatteringFunctions()` method initializes the
`SurfaceInteraction::bssrdf` member variable with an appropriate
`BSSRDF` if the material exhibits subsurface scattering.
(Section 11.4.3 will define two materials for
subsurface scattering.)

### 11.4.1 Separable BSSRDFs

One issue with the `BSSRDF` interface as defined above is its extreme
generality. Finding solutions to the subsurface light transport even in simple
planar or spherical geometries is already a fairly challenging problem, and the
fact that `BSSRDF` implementations can be attached to arbitrary and
considerably more complex `Shape`s leads to an impracticably difficult
context. To retain the ability to support general `Shape`s, we’ll
introduce a simpler BSSRDF representation in `SeparableBSSRDF`.

The constructor of `SeparableBSSRDF` initializes a local coordinate frame
defined by `ss`, `ts`, and `ns`, records the current light
transport mode `mode`, and keeps a pointer to the underlying `Material`.
The need for these values will be clarified in Section 15.4.

The simplified `SeparableBSSRDF` interface casts the BSSRDF into a
separable form with three independent components (one spatial and two
directional):

The Fresnel term at the beginning models the fraction of the light that is transmitted into direction after exiting the material. A second Fresnel term contained inside accounts for the influence of the boundary on the directional distribution of light entering the object from direction . The profile term is a spatial distribution characterizing how far light travels after entering the material.

For high-albedo media, the scattered radiance distribution is generally fairly isotropic and the Fresnel transmittance is the most important factor for defining the final directional distribution. However, directional variation can be meaningful for low-albedo media; in that case, this approximation is less accurate.

Given the separable expression in Equation (11.6), the integral for determining the outgoing illumination due to subsurface scattering (Section 15.5) simplifies to

The normalization factor is chosen so that integrates to one over the cosine-weighted hemisphere:

In other words,

This integral is referred to as the *first moment*
of the Fresnel reflectance function. Other moments involving
higher powers of the cosine function also exist and frequently occur in
subsurface scattering-related computations; the general definition of the
th Fresnel moment is

`pbrt` provides two functions `FresnelMoment1()` and
`FresnelMoment2()` that evaluate the corresponding moments based on
polynomial fits to these functions. (We won’t include their implementations
in the book here.) One subtlety is that these functions follow a
convention that is slightly different from the above definition—they are
actually called with the reciprocal of . This is due to their main
application in Section 15.5, where they account for the
effect light that is reflected back into the material due to a reflection
at the internal boundary with a relative index of refraction of .

Using `FresnelMoment1()`, the definition of `SeparableBSSRDF::Sw()`
based on Equation (11.7) can be easily
implemented.

Decoupling the spatial and directional arguments considerably reduces the
dimension of but does not address the fundamental difficulty with regard to
supporting general `Shape` implementations. We introduce a second
approximation, which assumes that the surface is not only locally planar but
that it is the distance between the points rather than their actual locations
that affects the value of the BSSRDF. This reduces to a
function that only involves distance of the two points
and :

As before, the actual implementation of the spatial term doesn’t
take as an argument, since it is already available in `BSSRDF::po`.

The `SeparableBSSRDF::Sr()` method remains virtual—it is overridden
in subclasses that implement specific 1D subsurface scattering
profiles. Note that the dependence on the distance introduces an implicit
assumption that the scattering medium is relatively homogeneous and does
not strongly vary as a function of position—any variation should be
larger than the mean free path length.

BSSRDF models are usually models of the function `Sr()` derived through a
careful analysis of the light transport within a homogeneous slab. This means
models such as the `SeparableBSSRDF` will yield good approximations in the planar
setting, but with increasing error as the underlying geometry deviates from
this assumption.

### 11.4.2 Tabulated BSSRDF

The single current implementation of the `SeparableBSSRDF` interface in
`pbrt` is the
`TabulatedBSSRDF` class. It provides access to a tabulated BSSRDF
representation that can handle a wide range of scattering profiles including
measured real-world BSSRDFs. `TabulatedBSSRDF` uses the same type of
adaptive spline-based interpolation method that was also used by the
`FourierBSDF` reflectance model (Section 8.6); in this
case, we are interpolating the distance-dependent scattering profile function
from Equation (11.9). The `FourierBSDF`’s
second stage of capturing directional variation using Fourier series is not
needed. Figure 11.15 shows
spheres rendered using the `TabulatedBSSRDF`.

`TabulatedBSSRDF`using a variety of measured BSSRDFs. From left to right: cola, apple, skin, and ketchup.

It is important to note that the radial profile is only a 1D function when all BSSRDF material properties are fixed. More generally, it depends on four additional parameters: the index of refraction , the scattering anisotropy , the albedo , and the extinction coefficient , leading to a complete function signature that is unfortunately too high-dimensional for discretization. We must thus either remove or fix some of the parameters.

Consider that the only parameter with physical units (apart from ) is
. This parameter quantifies the rate of scattering or absorption
interactions per unit distance.
The effect of is simple: it only controls the spatial scale of the
BSSRDF profile. To reduce the dimension of the necessary tables, we thus fix
and tabulate a *unitless* version of the BSSRDF profile.

When a lookup for a given extinction coefficient and radius
occurs at run time we find the corresponding unitless *optical radius*
and evaluate the lower-dimensional
tabulation as follows:

Since is a 2D density function in polar coordinates , a corresponding scale factor of is needed to account for this change of variables (see also Section 13.5.2).

We will also fix the index of refraction and scattering anisotropy
parameter —in practice, this means that these parameters cannot be
textured over an object that has a material that uses a
`TabulatedBSSRDF`. These
simplifications leave us with a fairly manageable 2D function that is only
discretized over albedos and optical radii .

The `TabulatedBSSRDF` constructor takes all parameters of the
`BSSRDF` constructor in addition to spectrally varying absorption and
scattering coefficients and . It
precomputes the extinction coefficient
and albedo
for the current
surface location, avoiding issues with a division by zero when there is no
extinction for a given spectral channel.

Detailed information about the scattering profile is supplied via the
`table` parameter, which is an instance of the `BSSRDFTable` data
structure:

Instances of `BSSRDFTable` record samples of the function taken
at a set of single scattering albedos
and radii . The spacing between radius and albedo
samples will generally be non-uniform in order to more accurately represent
the underlying function. Section 15.5.8 will
show how to initialize the `BSSRDFTable` for a specific BSSRDF model.

We omit the constructor implementation, which takes the desired resolution and allocates memory for the representation.

The sample locations and counts are exposed as public member variables.

A sample value is stored in the `profile` member variable for each of
the pairs .

Note that the `TabulatedBSSRDF::rho` member variable gives the
reduction in energy after a single scattering event; this is different
from the material’s overall albedo, which takes all orders of scattering into
account. To stress
the difference, we will refer to these different types of albedos as the
*single scattering albedo* and the *effective albedo*
.

We define the effective albedo as the following integral of the profile in polar coordinates.

The value of will frequently be accessed both by the
profile sampling code and by the `KdSubsurfaceMaterial`. We introduce an
array `rhoEff` of length `BSSRDFTable::nRhoSamples`, which maps every
albedo sample to its corresponding effective albedo.

Computation of will be discussed in Section 15.5. For now, we only note that it is a nonlinear and strictly monotonically increasing function of the single scattering albedo .

Given radius value `r` and single scattering albedo, the function
`Sr()` implements a spline-interpolated lookup into the tabulated profile.
The albedo parameter is taken to be the `TabulatedBSSRDF::rho` value.
There is thus an
implicit assumption that the albedo does not vary within the support of the
BSSRDF profile centered around .

Since the albedo is of type `Spectrum`, the function performs a separate
profile lookup for each spectral channel and returns a `Spectrum` as a
result. The return value must be clamped in case the interpolation produces
slightly negative values.

`ch`>>

`Sr[ch]`using tensor spline interpolation>>

The first line in the loop applies the scaling identity from
Equation (11.10) to obtain an optical radius for
the current channel `ch`.

Given the adjusted radius `rOptical` and the albedo
`TabulatedBSSRDF::rho``[ch]` at location `BSSRDF::po`, we next
call `CatmullRomWeights()` to obtain offsets and cubic spline weights to
interpolate the profile values. This step is identical to the `FourierBSDF`
interpolation in Section 8.6.

`ch`>>=

We can now sum over the product of the spline weights and the profile values.

`Sr[ch]`using tensor spline interpolation>>=

A convenience method in `BSSRDFTable` helps with finding profile
values.

It’s necessary to cancel a multiplicative factor of
that came from the entries of
`BSSRDFTable::profile` related to
Equation (11.11). This factor is present in the
tabularized values to facilitate importance sampling (more about this in
Section 15.4).
Since it is not part of the definition of the
BSSRDF, the term must be removed here.

Finally, we apply the change of variables factor from
Equation (11.10) to convert the interpolated
unitless BSSRDF value in `Sr` into world space units.

### 11.4.3 Subsurface Scattering Materials

There are two `Material`s for translucent
objects: `SubsurfaceMaterial`, which is defined in
`materials/subsurface.h` and
`materials/subsurface.cpp`,
and `KdSubsurfaceMaterial`, defined in
`materials/kdsubsurface.h` and
`materials/kdsubsurface.cpp`.
The only difference between these two materials is how the scattering
properties of the medium are specified.

`SubsurfaceMaterial` stores textures that allow the scattering properties
to vary as a function of the position on the surface. Note that this isn’t the
same as scattering properties that vary as a function of 3D inside the
scattering medium, but it can give a reasonable approximation to heterogeneous
media in some cases. (Note, however, that if used with spatially varying
textures, this feature destroys reciprocity of the BSSRDF, since these textures
are evaluated at just one of the two scattering points, and so interchanging
them will generally result in different values from the texture.
)

In addition to the volumetric scattering properties, a number of textures allow the user to specify coefficients for a BSDF that represents perfect or glossy specular reflection and transmission at the surface.

The `ComputeScatteringFunctions()` method uses the textures to
compute the values of the scattering properties at the point. The
absorption and scattering coefficients are then scaled by
the `scale` member variable, which provides an easy way to change the
units of the scattering properties. (Recall that they’re expected to be
specified in terms of inverse meters.)
Finally, the method creates a `TabulatedBSSRDF` with these parameters.

`bumpMap`, if present>> <<Initialize BSDF for

`SubsurfaceMaterial`>>

`bsdf`for smooth or rough dielectric>>

The fragment <<Initialize BSDF for `SubsurfaceMaterial`>> won’t
be included here; it follows the now familiar approach of allocating
appropriate `BxDF` components for the `BSDF`
according to which textures have nonzero SPDs.

Directly setting the absorption and scattering coefficients to achieve a
desired visual look is difficult. The parameters have a nonlinear and
unintuitive effect on the result. The
`KdSubsurfaceMaterial` allows the
user to specify the subsurface scattering properties in terms of the
diffuse reflectance of the surface and the mean free path
. It then uses the `SubsurfaceFromDiffuse()`
utility function, which will be defined in Section 15.5, to
compute the corresponding intrinsic scattering properties.

Being able to specify translucent materials in this manner is particularly useful in that it makes it possible to use standard texture maps that might otherwise be used for diffuse reflection to define scattering properties (again with the caveat that varying properties on the surface don’t properly correspond to varying properties in the medium).

We won’t include the definition of `KdSubsurfaceMaterial` here since
its implementation just evaluates `Texture`s to compute the diffuse
reflection and mean free path values and
calls `SubsurfaceFromDiffuse()` to compute the scattering properties
needed by the `BSSRDF`.

Finally, `GetMediumScatteringProperties()` is a utility function that
has a small library of measured scattering data for translucent materials;
it returns the corresponding scattering properties if it has an entry for
the given name. (For a list of the valid names, see the implementation
in `core/media.cpp`.) The data provided by this function is from
papers by Jensen et al. (2001b) and Narasimhan et al. (2006).