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 Shapes leads to an impracticably difficult context. To retain the ability to support general Shapes, 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

We define the directional term as a scaled version of the Fresnel transmittance (Section 8.2):

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

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.

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.

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

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

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 Textures 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/medium.cpp.) The data provided by this function is from papers by Jensen et al. (2001b) and Narasimhan et al. (2006).