15.5 Subsurface Scattering Using the Diffusion Equation

Our last task to complete the subsurface scattering implementation is to be able to initialize the TabulatedBSSRDF with a radial profile function that accurately describes subsurface scattering for given properties of the scattering medium (, , the phase function asymmetry parameter , and the relative index of refraction ). The technique we’ll discuss in this section is based on the photon beam diffusion (PBD) technique by Habel et al. (2013). The resulting profile takes all orders of scattering into account, effectively accounting for all of the light transport that occurs within the surface.

Beam diffusion makes several significant assumptions and approximations: first, the distribution of light in the translucent medium is modeled with the diffusion approximation, which describes the equilibrium distribution of illumination in highly scattering optically thick participating media. Second, it assumes homogeneous scattering properties throughout the medium, and it implicitly assumes that the medium is semi-infinite (it continues infinitely beneath a planar surface of infinite lateral extent). Finally, PBD builds upon the separable BSSRDF approximation of Equation (11.6), which imposes a simple multiplicative relationship between the spatial and directional scattering distribution. When these approximations are satisfied, the solutions computed by PBD are in close agreement with ground truth simulations performed using the equation of transfer, Equation (15.2).

Of course, many of these assumptions won’t be valid when the profile is applied to an arbitrary shape, potentially also with spatially varying material properties. The appeal of diffusion-type methods in the context of computer graphics is that they degrade in a graceful manner, producing visually reasonable results even in cases where some or all of their fundamental assumptions are violated. See the “Further Reading” section and exercises at the end of this chapter for references to improvements to this approach that generalize it to handle a wider range of settings more accurately.

While PBD can compute the profile for any radius and material parameters, profile evaluations tend to be fairly expensive, as they involve a numerical integration step. Furthermore, we need to be able to invert the CDF of in polar coordinates to importance sample the model, but this inverse is not available in closed form. The TabulatedBSSRDF from Section 11.4.2 nicely addresses these issues, providing efficient evaluation and sampling operations. Our approach will thus be to precompute for a range of radii and albedo values and use the results to populate the BSSRDFTable of a TabulatedBSSRDF. This precomputation is performed during scene construction.

In the following, we discuss the key ingredients of the PBD method, starting with the principle of similarity and diffusion theory.

15.5.1 Principle of Similarity

A number of important ideas are used in the process of transforming the fully general equation of transfer to the diffusion equation, which can be approximately solved for subsurface scattering. The first is the principle of similarity, which says that for an anisotropically scattering medium with a high albedo, the medium can instead be modeled as having an isotropic phase function with appropriately modified scattering and attenuation coefficients. Light transport solutions computed based on the modified coefficients correspond well to those with the original coefficients and phase function, while allowing simplifications due to the assumption of isotropic scattering.

The principle of similarity is based on the observation that after many scattering events the distribution of light in media with high albedos becomes more and more uniformly directionally distributed regardless of the original illumination distribution or the anisotropy of the phase function. One way to see how this happens is to consider an expression derived by Yanovitskij (1997); it describes isotropization due to multiple scattering events from the Henyey–Greenstein phase function. He showed that the distribution of light that has been scattered times is given by

As grows large, this converges to the isotropic phase function, . Figure 15.14 plots this function for a few values of . When coupled with the observation from Section 15.4.1 about how much of the light energy remains after tens or even hundreds of scattering events in high-albedo materials, one can see intuitively why it is reasonable to work with an isotropic phase function approximation for high-albedo media.

If the principle of similarity is applied to treat the phase function as isotropic, modified versions of various scattering properties should be used. The reduced scattering coefficient is defined as , and the reduced attenuation coefficient is . The albedo also changes to a reduced albedo defined as , which is generally different from . These new coefficients account for the effect of using the isotropic phase function approximation.

To understand the idea these coefficients embody, consider a strongly forward-scattering phase function, where . With the original phase function, light will mostly continue in the same direction once it scatters. In this example, the value of the reduced scattering coefficient is much smaller than , which means that light travels a larger distance in the medium before scattering; the medium is approximated as being thinner, allowing light to travel farther. The change thus has the same effect as a highly forward-scattering phase function.

Conversely, consider the case of . In this case, at a scattering event the light will tend to scatter back in the direction it came from. But then the next time it scatters after that, it will generally reverse course again; it bounces back and forth without making very much forward progress. In this case, the reduced scattering coefficient is larger than the original scattering coefficient, indicating greater probability of a scattering interaction. In other words, the medium is treated as being thicker than it actually is, which approximates the effect of light having relatively more trouble making forward progress. Figure 15.15 illustrates these ideas, showing representative paths of scattering interactions in highly forward-scattering and highly backward-scattering media.

15.5.2 Diffusion Theory

Diffusion theory provides a way of transitioning from the equation of transfer to a simpler diffusion equation, which provides a solution to the equation of transfer for the case of homogeneous, optically thick, highly scattering materials (i.e., those with relatively large albedos). For the application to subsurface scattering in pbrt, it can be derived by writing the equation of transfer using the reduced scattering and attenuation coefficients and an isotropic phase function.

Starting with the integro-differential form of the equation of transfer from Equation (15.1),

we assume spatially uniform material parameters and switch to an isotropic phase function , making a corresponding change to the scattering and attenuation coefficient using similarity theory. We also replace with a single function .

The key assumption of diffusion theory is that because each scattering event effectively blurs the incident illumination, high frequencies disappear from the angular radiance distribution as light propagates farther into the medium; in dense and isotropically scattering media, all directionality is eventually lost. Motivated by this observation, the radiance function is restricted to a simple two-term expansion based on spherical moments. Formally, for a function , the -th moment on the unit sphere is defined as ^†

In other words, to get the entry of the -tensor , we integrate the product of and the components of the direction written in Cartesian coordinates. Note that there is some notational overlap with the angle cosines from Section 8.6; the remainder of this chapter will exclusively refer to the definition above.

The zeroth moment, for instance, gives the function’s integral over the sphere, the first moment can be interpreted as a “center of mass” 3-vector, and the second moment is a positive definite 33 matrix. Higher order moments have many symmetries: for instance, exchanging any pair of indices leaves the value unchanged. High-order moments are useful to derive extended versions of diffusion theory that allow for more pronounced directional behavior; here, we will just focus on degrees .

The moments of the radiance function have special designations: the zeroth moment is referred to as the fluence rate:

Note that this expression differs from the fluence function in Equation (6.6), which was defined as the fluence rate on a surface boundary integrated over time.

The first moment is the vector irradiance:

The two-term expansion mentioned before is defined as

so that the moments can be exactly recovered, i.e.,

(Here, the “d” subscript denotes the diffusion approximation, not the direct lighting term as it did in Section 14.3.)

To derive the diffusion equation from the equation of transfer, we simply substitute the two-term radiance function into Equation (15.15). The resulting expression is unfortunately not guaranteed to have a solution, but this issue can be addressed with a simple trick: by only enforcing equality of its moments, i.e., by requiring that

for and . Computing these moments is a fairly lengthy and mechanical exercise in trigonometric calculus that we skip here. ^† The end result is an equation equating the zeroth moments:

where is the divergence of and

is the -th moment of the medium emission. This equation states that the divergence of the irradiance vector field is negative in the presence of absorption (i.e., light is being removed) and positive when light is being added by .

Another similar equation for the first moments states that the irradiance vector field , which represents the overall flow of energy, points from regions with a higher fluence rate to regions with a lower rate.

A reasonable simplification at this point is to assume light sources in the medium emit light uniformly in all directions, in which case .

The next step of the traditional derivation is to solve the above equation for and substitute it into the equation relating zeroth moments. The substitution removes and yields the diffusion equation, which now only involves the fluence rate :

Assuming that , the diffusion equation can be written more compactly as

where is the classical diffusion coefficient and is a shorter notation for , which is known as the Laplace operator.

With the diffusion equation at hand, we’ll proceed as follows: starting from a solution for a point source that is only correct in a space where the medium infinitely extends in all directions, we will consider ways of improving the solution’s accuracy in more challenging cases and introduce an approximation that can account for the effect of a refractive boundary.

We’ll initially focus on a point light source that is placed below the surface to approximate the effect of incident illumination striking the surface. Later, we switch to a more accurate light source approximation and derive the beam diffusion solution to the multiple scattering component as well as a single scattering correction that is based on the classical equation of transfer.

15.5.3 Monopole Solution

Consider an infinite homogeneous medium with a point light source of unit power (a monopole) located at its origin. The emitted radiance function for this setup is given by

and the corresponding moments are equal to

The fluence rate due to this type of source has a simple analytic expression:

where is the distance from the light source. The constant is called the effective transport coefficient; it occurs in an exponential falloff term that accounts for absorption in the medium. Observe that : instead, this modified attenuation coefficient additionally depends on the albedo to model the effect of multiple scattering inside the medium, hence the term effective.

Let’s verify that this expression satisfies the diffusion equation away from the origin, where the fluence has a pole singularity. The Laplace operator of a function that only depends on the radius in spherical coordinates is given by

Taking the inner derivative of and multiplying by yields

Taking the outer derivative, multiplying by , and simplifying gives

which is exactly the ratio predicted by Equation (15.19); i.e., .

Using the identity from Equation (15.18), we can also find the irradiance vector field induced by , which will be useful later on:

where is a unit vector pointing away from the source.

15.5.4 Non-classical Diffusion

As we have seen, the monopole fluence rate in Equation (15.20) exactly solves the diffusion equation. Despite this, it turns out that the solution can still have significant errors compared to the original equation of transfer when the assumptions of the underlying diffusion approximation are violated.

There are two important cases: the first occurs when absorption prevents the radiance from reaching an isotropic equilibrium distribution. The second occurs close to the source, where the true radiance function is dominated by a (highly non-isotropic) spherical Dirac delta function.

A number of modified diffusion theories have been proposed that improve the accuracy in a number of different cases. An effective one is to switch to the modified monopole solution developed by Grosjean (1956) in the field of neutron transport:

The first term in this solution separates out an attenuated source term modeled using standard radiative transport—this effectively removes the portion that cannot easily be handled by diffusion. The remainder is a scaled diffusive term that accounts for light which has been scattered at least once, where the reduced albedo accounts for the energy reduction due to the extra scattering event:

This expression uses the previous monopole solution , though its diffusion coefficient must be replaced with a non-classical version given by

Figure 15.16 illustrates the superior accuracy of Grosjean’s solution in both absorption-dominated and scattering-dominated media.

In the following sections, we will just focus on the non-classical diffusive part and ignore the attenuated source term in Equation (15.22), as it is simple to handle separately later on.

15.5.5 Dipole Solution

To apply these results to subsurface scattering for rendering, it is clear that the solution must account for the presence of a surface. We will now switch to the simplest kind of geometric setup that satisfies this requirement: a semi-infinite half space, i.e., a medium that fills all space below a planar surface of infinite lateral extent. The region above is modeled as a dielectric without any kind of scattering (i.e., a vacuum).

For simplicity, we assume that the boundary is located at with a normal of so that positive values on the axis correspond to points inside the medium. Let denote the relative index of refraction over the boundary. We are still interested in the fluence rate due to a point light source that we (arbitrarily) place on the -axis at position . Let’s assume that , i.e., the light source is located inside the medium (more on this later). Due to the newly added boundary, a portion of the light traveling upward can escape from the layer and undergo no further scattering. Another portion is specularly reflected at ; the monopole solution from Equation (15.20) accounts for neither effect and thus no longer yields accurate results.

The influence of the boundary can be approximated by the method of images, where a negative source is placed on the vacuum side of the boundary at position with . The negative contribution of this “virtual” light source subtracts a portion of the “real” light source to account for the combined effect of internal reflections and illumination that has left the medium and can no longer scatter. The choice of the virtual depth is crucial to ensure that this indeed works out—we discuss this step shortly.

This arrangement of a positive and a negative light source is known as a dipole (Figure 15.17).

Due to the linearity of the diffusion equation, (15.19), superpositions of solutions still solve the diffusion equation; hence the dipole fluence rate at a point on the boundary is simply the sum of a positive and negative term:

where and are the straight-line distances from the evaluation point to the real and virtual light source. Once more, we can find a matching vector irradiance value using the monopole solution from Equation (15.21):

The tilde in above indicates the use of Grosjean’s modified diffusion coefficient. We will later require the component of this expression, which is given by

Boundary Conditions

Having defined the diffusion dipole, we must still specify how the two light sources should be placed in relation to each other so that they fulfill the appropriate boundary conditions (internal reflections, no scattering for ). We’ll assume that the real light source depth is specified, and that must be set to correct for the boundary.

Moulton (1990) showed that a very simple approximation is enough to get a reasonable answer. As we move from the boundary into the region filled by vacuum, it is intuitive that the fluence rate should decrease fairly quickly due to the inverse squared falloff in a region of space that does not scatter. If we model the fluence rate along the -axis using a first-order Taylor expansion of the boundary conditions at , then this linear function eventually reaches 0 (and then negative values) at a linear extrapolation depth of .

The idea of this approach then is to mirror the real light source across the mirror plane defined by this extrapolation depth to obtain the virtual light source depth ; this ensures that the fluence rate of the dipole is exactly equal to 0 at above the half-space. Note that we would generally expect a point outside the medium to have a nonzero positive fluence rate due to radiance leaving the medium; here we are only interested in computing a good solution at the boundary, where the somewhat non-physical negative light source does not pose problems.

An improved variational approximation of boundary conditions for interfaces with internal Fresnel reflection was derived by Pomraning and Ganapol (1995). With their approach, the linear extrapolation depth is given by

where and are the Fresnel moments first encountered in Equation (11.8).

Radiant Exitance

At this point, we have all ingredients to compute the fluence rate due to a point source inside the surface, with corrections to account for the half-space geometry and internal reflections. To use this solution in a light transport simulation, we’ll need to know how much light actually leaves the surface.

Recall Equation (15.16), which related radiance to the fluence rate and vector irradiance. Plugging the dipole solutions into this expression yields

but this is only valid inside the surface. To find the amount of diffusely scattered light leaving at the boundary, we can integrate the internal radiance distribution against the Fresnel transmittance and a cosine factor due to the term in the definition of radiance (Section 5.4).

Using linearity to split the integral into two parts that are related to the fluence rate and vector irradiance, we have

When the fluence-related part is written in spherical coordinates, the integral over azimuth turns into a multiplication by as no terms depend on it, and can be moved out of the integral. The remaining expression reduces to a constant plus a scaled Fresnel moment (Equation (11.8)).

For the part depending on , we obtain

In summary: after putting together all the pieces, we have a method of evaluating the radiant exitance at a position on the boundary due to an internal source at depth . So far, this depth was assumed to be a fixed parameter, but now we’ll promote it to an argument of the radiant exitance we just derived, which is now written .

15.5.6 Beam Solution

At this point, we are almost ready to implement the model, though one missing piece that must still be addressed is the depth of the positive point light source—this choice will clearly have an effect on the accuracy of the final solution.

The first dipole-based BSSRDF model for computer graphics proposed by Jensen et al. (2001b) placed the source at a depth of one mean free path—i.e., —inside the medium, which is the expected distance that light will travel after entering the surface. This is a reasonable approximation, though it leads to significant errors close to the source.

With the PBD method, the point source solution is integrated over a semi-infinite interval that considers all positions where light traveling along a perpendicularly incident collimated beam could scatter. The impulse response of the medium to such a spatio-directional Dirac delta function provides a more faithful description of its reflection behavior. More advanced variants of this model also allow for rays with non-perpendicular incidence. Formally, this method computes

The exponential models how the incident beam’s power diminishes due to extinction by the medium. Though there are a variety of high-accuracy approximations to this integral that use only a few samples, for the implementation in pbrt, performance is less critical since the diffusion solution is only computed once. For this reason, we use a rudimentary but simpler importance sampling scheme of the exponential term in Equation (15.33).

The function BeamDiffusionMS() takes the medium properties , and a radius and returns an average of 100 samples of the integrand.

A number of coefficients that do not depend on can be precomputed outside the loop.

We begin by setting the reduced scattering and attenuation coefficients and single scattering albedo using the principle of similarity from Section 15.5.1.

Following this, we compute Grosjean’s non-classical diffusion coefficient, Equation (15.24), and the corresponding effective transport coefficient (Section 15.5.3).

Neither the linear extrapolation depth nor the scale factors from the radiant exitance computation depend on , so they can also be computed. The FresnelMoment1() and FresnelMoment2() functions previously encountered in Section 11.4.1 evaluate and .

This concludes the precomputation—all following fragments occur in the loop over samples. To select the point source depth , we importance sample the homogeneous attenuation term using the same approach as in Section 15.2, setting

for . Because we are integrating a smooth 1D function, Monte Carlo isn’t required. We therefore use equal-spaced positions , where .

Given , we next determine the straight-line distances to the real and virtual light source, whose position is found by mirroring the real source across the linear extrapolation depth at .

The dipole fluence rate and normal component of the irradiance were previously discussed in Equations (15.25) and (15.27).

The last fragment computes the diffuse radiant exitance E due to the dipole using Equation (15.30) and adds a scaled version of it to the running sum in Ed. In this computation, the first rhop factor in the scale is needed due to the ratio of the importance sampling weight of the sampling strategy in Equation (15.34) and the factor in Equation (15.33). The second rhop factor accounts for the additional scattering event in Grosjean’s non-classical monopole in Equation (15.23). Finally an empirical correction factor

corrects an overestimation that can occur when is small and the light source is close to the surface (see Habel et al. (2013) for details).

15.5.7 Single Scattering Term

Having accounted for the diffusive multiple scattering term, it’s also necessary to account for single scattering, which is not handled well by the diffusion approximation. Single scattering contributes a significant amount of energy to the scattering profile close to . In the absence of multiple scattering, it is fortunately simple enough to compute its contribution directly with the original equation of transfer.

Using the generalized path integral from Section 15.1.1, the radiance in the medium traveling toward position after scattering precisely once in the volume at position is given by

where is the generalized throughput function from Equation (15.3). We model illumination coming from a collimated beam at fixed position pointing into the negative normal direction, i.e.:

Similar to the previous section, we are interested in the outgoing irradiance at a point on the half-space boundary due to the decaying beam of light. Relevant distances and angles in this arrangement can be found using simple triangle identities (Figure 15.18). We’ll briefly ignore the effect of refraction at to simplify the derivation but will account for it later on in the final result.

The irradiance at can be found by integrating the single-scattered radiance and a generalized geometry term of Equation (15.5) over the volume. Expanding the definition of from Equation (15.35) produces an integral over paths of length 2:

After expanding , the spatial delta function from Equation (15.36) removes the integration over , and the directional term reduces Equation (15.37) to a 1D integral along the ray :

The extra term is a change of variables factor resulting from an intermediate step (not shown here), where the volumetric integral over is expressed in terms of spherical coordinates . The and variables subsequently drop out due to the directional delta function in .

Expanding the generalized throughput using Equation (15.3) yields

The definition of the generalized scattering function from Equation (15.4) and the assumption that the phase function only depends on the cosine of the scattering angle imply that can be written as

We will use the Henyey–Greenstein phase function model for . The angle cosine can be found from the triangle edge lengths (Figure 15.18):

The hypotenuse is the distance from the scattering location to the exit point and is given by the Pythagorean theorem:

where .

Due to the assumptions of a homogeneous medium without internal occluders, the first geometric term in Equation (15.38) takes on a simple form:

Note that this equation uses the original attenuation coefficient rather than the reduced version from Section 15.5.1. In contrast to diffusion theory, the equation of transfer easily accounts for anisotropy; hence there is no longer a need for this approximation.

For the second geometry term, we must include a cosine factor that accounts for the angle that the connecting ray segment makes with the surface normal when intersecting the boundary at :

The two cosine terms and are equal except for a difference in sign, which can be seen from the triangle geometry in Figure 15.18.

Plugging Equations (15.42), (15.43), and (15.39) into Equation (15.38) yields the following expression for the irradiance due to single scattering:

At this point, we’ll reintroduce the effect of the refractive boundary by adding a Fresnel transmission factor . Light that is internally reflected by the boundary is already included in the diffusion solution and should be excluded from the single-scattering profile. As before in Section 15.5.6, responsibility for the first refraction at is delegated to the Material’s BSDF and thus there is just a single Fresnel transmission term,

Before we begin with the implementation of this integral, there is one more effect that should be considered: when the relative index of refraction over the boundary is greater than 1, no illumination can directly leave the material at angles below the critical angle due to total internal reflection, where

To avoid unnecessary computations for scattering locations that cannot possibly contribute, we restrict the integral to the range satisfying

Combining this with Equation (15.41) and solving for yields

As in Section 15.5.6, we compute this integral by uniformly sampling distances according to an exponential distribution (which now starts at depth ) and set

for evenly spaced . The associated PDF is

and dividing the integrand in Equation (15.44) by this PDF results in a throughput weight of

where was previously defined as .

The PBD single-scattering profile computation is implemented in BeamDiffusionSS(), which takes the medium scattering properties and a radius as input. One hundred samples are used for the integral estimate, as in BeamDiffusionMS().

The function begins by precomputing derived material parameters and setting to the minimal distance below the critical angle as specified by Equation (15.45).

The loop body generates distances according to the sampling scheme in Equation (15.46).

Next, the function computes the length of the connection segment using Equation (15.41). is given by Equation (15.40) except for a sign flip that is needed since the half-space normal points away from the medium.

The last loop statement accumulates the throughput weight from Equation (15.47) into the running sum Ess.

15.5.8 Filling the BSSRDFTable

With the definitions of BeamDiffusionMS() and BeamDiffusionSS() at hand, we’ll now implement the function ComputeBeamDiffusionBSSRDF() that uses these functions to fill the BSSRDFTable of a TabulatedBSSRDF with profile data.

ComputeBeamDiffusionBSSRDF() takes the medium’s anisotropy parameter and relative index of refraction as input and initializes a BSSRDFTable that stores these quantities as a function of radius and albedo. This function is in turn invoked by the initialization routines of the SubsurfaceMaterial and KdSubsurfaceMaterial using a default BSSRDFTable with 100 albedo samples and 64 radius samples.

Both single- and multiple-scattering components of the profile function are characterized by an exponential decay with increasing radius . By placing many samples in high-valued regions and comparably fewer in low-valued regions, we can obtain a better spline approximation of the true profile. The following fragment places the first radius sample at 0 and distributes the remaining samples at exponentially increasing distances. (Note that radius values are unitless and independent of the actual medium density . These unitless optical radii were previously introduced in Section 11.4.2.)

Next, we need to decide on the locations of albedo samples on the interval . Some precautions must be taken here as well, since the material’s scattering behavior has a highly nonlinear dependence on the parameter.