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

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

`Ed`>>

`Ed`>>

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 .

`Ed`>>=

`Ed`>>

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

`Ed`>>=

### 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: