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 upper S Subscript normal r that accurately describes subsurface scattering for given properties of the scattering medium ( sigma Subscript normal a , sigma Subscript normal s , the phase function asymmetry parameter g , and the relative index of refraction eta ). 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 upper S Subscript normal r 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 upper S Subscript normal r 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 upper S Subscript normal r 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 n times is given by

p Subscript n Baseline left-parenthesis omega Subscript Baseline right-arrow omega prime Subscript Baseline right-parenthesis equals StartFraction 1 minus g Superscript 2 n Baseline Over 4 pi left-parenthesis 1 plus g Superscript 2 n Baseline minus 2 g StartAbsoluteValue g Superscript n minus 1 Baseline EndAbsoluteValue left-parenthesis minus omega Subscript Baseline dot omega prime Subscript Baseline right-parenthesis right-parenthesis Superscript 3 slash 2 Baseline EndFraction period

As n grows large, this converges to the isotropic phase function, 1 slash left-parenthesis 4 pi right-parenthesis . Figure 15.14 plots this function for a few values of n . 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.

Figure 15.14: Light Distribution after Many Scattering Events. (a) Directional distribution of a single incident ray of light after 10 scattering events in a highly anisotropic medium with g equals 0.9 , (b) 20 scattering events, (c) 200 scattering events. The distribution becomes increasingly isotropic, even though it was initially very anisotropic.

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 sigma prime Subscript normal s Baseline equals left-parenthesis 1 minus g right-parenthesis sigma Subscript normal s , and the reduced attenuation coefficient is sigma prime Subscript normal t Baseline equals sigma Subscript normal a Baseline plus sigma prime Subscript normal s . The albedo also changes to a reduced albedo defined as rho prime equals sigma prime Subscript normal s Baseline slash sigma prime Subscript normal t , which is generally different from rho . 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 g right-arrow 1 . 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 sigma prime Subscript normal s Baseline equals left-parenthesis 1 minus g right-parenthesis sigma Subscript normal s is much smaller than sigma Subscript normal s , 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 g right-arrow negative 1 . 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.

Figure 15.15: Representative Light Paths for Highly Anisotropic Scattering Media. (a) Forward-scattering medium, with g equals 0.9 . Light generally scatters in the same direction it was originally traveling. (b) Backward-scattering medium, with g equals negative 0.9 . Light frequently bounces back and forth, making relatively little forward progress with respect to its original direction.

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

StartLayout 1st Row 1st Column StartFraction partial-differential Over partial-differential t EndFraction upper L Subscript normal o Superscript Baseline left-parenthesis normal p Subscript Baseline plus t omega comma omega Subscript Baseline right-parenthesis equals 2nd Column minus sigma Subscript normal t Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis upper L Subscript normal i Superscript Baseline left-parenthesis normal p Subscript Baseline comma minus omega Subscript Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column plus sigma Subscript normal s Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis integral Underscript script upper S squared Endscripts p left-parenthesis normal p Subscript Baseline comma minus omega prime Subscript Baseline comma omega Subscript Baseline right-parenthesis upper L Subscript normal i Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega prime Subscript Baseline right-parenthesis normal d omega prime Subscript plus upper L Subscript normal e Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis comma EndLayout

we assume spatially uniform material parameters and switch to an isotropic phase function p equals 1 slash 4 pi , making a corresponding change to the scattering and attenuation coefficient using similarity theory. We also replace upper L Subscript normal o Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis equals upper L Subscript normal i Baseline left-parenthesis normal p Subscript Baseline comma minus omega Subscript Baseline right-parenthesis with a single function upper L left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis .

StartFraction partial-differential Over partial-differential t EndFraction upper L Subscript Superscript Baseline left-parenthesis normal p Subscript Baseline plus t omega Subscript Baseline comma omega Subscript Baseline right-parenthesis equals minus sigma prime Subscript normal t Baseline upper L Subscript Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis plus StartFraction sigma prime Subscript normal s Baseline Over 4 pi EndFraction integral Underscript script upper S squared Endscripts upper L Subscript Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega prime Subscript Baseline right-parenthesis normal d omega Subscript Baseline Superscript prime Baseline plus upper L Subscript normal e Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis period

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 f colon script upper S squared right-arrow bold upper R Superscript , the n -th moment on the unit sphere is defined as

left-parenthesis mu Subscript n Baseline left-bracket f right-bracket right-parenthesis Subscript i comma j comma k comma ellipsis Baseline equals integral Underscript upper S squared Endscripts ModifyingBelow omega Subscript normal i Baseline omega Subscript j Baseline omega Subscript k Baseline midline-horizontal-ellipsis With bottom-brace Underscript n factors Endscripts f left-parenthesis omega right-parenthesis normal d omega Subscript period Baseline

In other words, to get the i comma j comma k comma ellipsis entry of the n -tensor mu Subscript n Baseline left-bracket f right-bracket , we integrate the product of f and the i comma j comma k comma ellipsis components of the direction omega written in Cartesian coordinates. Note that there is some notational overlap with the angle cosines mu Subscript k 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 3 times 3 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 n less-than-or-equal-to 1 .

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

phi left-parenthesis normal p Subscript Baseline right-parenthesis equals mu 0 left-bracket upper L left-parenthesis normal p Subscript Baseline comma dot right-parenthesis right-bracket equals integral Underscript script upper S squared Endscripts upper L Subscript Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis normal d omega Subscript Baseline period

Note that this expression differs from the fluence function upper H left-parenthesis normal p Subscript Baseline right-parenthesis 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:

bold upper E left-parenthesis normal p Subscript Baseline right-parenthesis equals mu 1 left-bracket upper L left-parenthesis normal p Subscript Baseline comma dot right-parenthesis right-bracket equals integral Underscript script upper S squared Endscripts omega Subscript Baseline upper L Subscript Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis normal d omega Subscript Baseline period

The two-term expansion mentioned before is defined as

upper L Subscript normal d Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis equals StartFraction 1 Over 4 pi EndFraction phi left-parenthesis normal p Subscript Baseline right-parenthesis plus StartFraction 3 Over 4 pi EndFraction omega Subscript Baseline dot bold upper E left-parenthesis normal p Subscript Baseline right-parenthesis comma

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

mu 0 left-bracket upper L Subscript normal d Superscript Baseline left-parenthesis normal p Subscript Baseline comma dot right-parenthesis right-bracket equals phi left-parenthesis normal p Subscript Baseline right-parenthesis and mu 1 left-bracket upper L Subscript normal d Superscript Baseline left-parenthesis normal p Subscript Baseline comma dot right-parenthesis right-bracket equals bold upper E left-parenthesis normal p Subscript Baseline right-parenthesis period

(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 upper L Subscript normal d Superscript 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

StartLayout 1st Row 1st Column mu Subscript i Baseline left-bracket StartFraction partial-differential Over partial-differential t EndFraction upper L Subscript normal d Superscript Baseline left-parenthesis normal p Subscript Baseline plus t omega Subscript Baseline comma omega Subscript Baseline right-parenthesis right-bracket equals mu Subscript i Baseline left-bracket 2nd Column minus sigma prime Subscript normal t Baseline upper L Subscript normal d Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column plus StartFraction sigma prime Subscript normal s Baseline Over 4 pi EndFraction integral Underscript script upper S squared Endscripts upper L Subscript normal d Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega prime Subscript Baseline right-parenthesis normal d omega prime Subscript plus upper L Subscript normal e Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis right-bracket EndLayout

for i equals 0 and i equals 1 . 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:

div of bold upper E left-parenthesis normal p Subscript Baseline right-parenthesis equals left-parenthesis minus sigma prime Subscript normal t plus sigma prime Subscript normal s right-parenthesis phi left-parenthesis normal p Subscript Baseline right-parenthesis equals minus sigma Subscript normal a Baseline phi left-parenthesis normal p Subscript Baseline right-parenthesis plus upper Q 0 left-parenthesis normal p Subscript Baseline right-parenthesis comma

where div of bold upper E left-parenthesis normal p Subscript Baseline right-parenthesis equals StartFraction partial-differential Over partial-differential x EndFraction bold upper E left-parenthesis normal p Subscript Baseline right-parenthesis plus StartFraction partial-differential Over partial-differential y EndFraction bold upper E left-parenthesis normal p Subscript Baseline right-parenthesis plus StartFraction partial-differential Over partial-differential z EndFraction bold upper E left-parenthesis normal p Subscript Baseline right-parenthesis is the divergence of bold upper E left-parenthesis normal p Subscript Baseline right-parenthesis and

upper Q Subscript i Baseline left-parenthesis normal p Subscript Baseline right-parenthesis equals mu Subscript i Baseline left-bracket upper L Subscript normal e Superscript Baseline left-parenthesis normal p Subscript Baseline comma dot right-parenthesis right-bracket

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

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

one-third nabla phi left-parenthesis normal p Subscript Baseline right-parenthesis equals minus sigma prime Subscript normal t Baseline bold upper E left-parenthesis normal p Subscript Baseline right-parenthesis plus upper Q 1 left-parenthesis normal p Subscript Baseline right-parenthesis comma

A reasonable simplification at this point is to assume light sources in the medium emit light uniformly in all directions, in which case upper Q 1 left-parenthesis normal p Subscript Baseline right-parenthesis equals 0 .

The next step of the traditional derivation is to solve the above equation for bold upper E and substitute it into the equation relating zeroth moments. The substitution removes bold upper E left-parenthesis normal p Subscript Baseline right-parenthesis and yields the diffusion equation, which now only involves the fluence rate phi left-parenthesis normal p Subscript Baseline right-parenthesis :

StartFraction 1 Over 3 sigma prime Subscript normal t EndFraction div of nabla phi left-parenthesis normal p Subscript Baseline right-parenthesis equals sigma Subscript normal a Baseline phi left-parenthesis normal p Subscript Baseline right-parenthesis minus upper Q 0 left-parenthesis normal p Subscript Baseline right-parenthesis plus StartFraction 1 Over sigma prime Subscript normal t EndFraction nabla dot upper Q 1 left-parenthesis normal p Subscript Baseline right-parenthesis period

Assuming that upper Q 1 left-parenthesis normal p Subscript Baseline right-parenthesis equals 0 , the diffusion equation can be written more compactly as

upper D nabla squared phi left-parenthesis normal p Subscript Baseline right-parenthesis minus sigma Subscript normal a Baseline phi left-parenthesis normal p Subscript Baseline right-parenthesis equals minus upper Q 0 left-parenthesis normal p Subscript Baseline right-parenthesis comma

where upper D equals 1 slash left-parenthesis 3 sigma prime Subscript normal t right-parenthesis is the classical diffusion coefficient and nabla squared is a shorter notation for div of nabla , 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 upper L Subscript normal e Superscript for this setup is given by

upper L Subscript normal e Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis equals StartFraction 1 Over 4 pi EndFraction delta left-parenthesis normal p Subscript Baseline right-parenthesis comma

and the corresponding moments are equal to

upper Q Subscript i Baseline left-parenthesis normal p Subscript Baseline right-parenthesis equals StartLayout Enlarged left-brace 1st Row 1st Column delta left-parenthesis normal p Subscript Baseline right-parenthesis comma 2nd Column i equals 0 comma 2nd Row 1st Column 0 comma 2nd Column i equals 1 period EndLayout

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

phi Subscript normal upper M Baseline left-parenthesis r right-parenthesis equals StartFraction 1 Over 4 pi upper D EndFraction StartFraction normal e Superscript minus sigma Super Subscript normal t normal r Superscript r Baseline Over r EndFraction comma

where r is the distance from the light source. The constant sigma Subscript normal t normal r Baseline equals StartRoot sigma Subscript normal a Baseline slash upper D EndRoot is called the effective transport coefficient; it occurs in an exponential falloff term that accounts for absorption in the medium. Observe that sigma Subscript normal t normal r Baseline not-equals sigma prime Subscript normal t : 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 nabla squared of a function f that only depends on the radius in spherical coordinates is given by

nabla squared f equals r Superscript negative 2 Baseline StartFraction partial-differential Over partial-differential r EndFraction left-parenthesis r squared StartFraction partial-differential Over partial-differential r EndFraction f left-parenthesis r right-parenthesis right-parenthesis period

Taking the inner derivative of phi Subscript normal upper M Baseline left-parenthesis r right-parenthesis and multiplying by r squared yields

r squared StartFraction partial-differential Over partial-differential r EndFraction phi Subscript normal upper M Baseline left-parenthesis r right-parenthesis equals minus r phi Subscript normal upper M Baseline left-parenthesis r right-parenthesis left-parenthesis 1 plus sigma Subscript normal t normal r Baseline r right-parenthesis period

Taking the outer derivative, multiplying by r Superscript negative 2 , and simplifying gives

r Superscript negative 2 Baseline StartFraction partial-differential Over partial-differential r EndFraction left-parenthesis minus r phi Subscript normal upper M Baseline left-parenthesis r right-parenthesis left-parenthesis 1 plus sigma Subscript normal t normal r Baseline r right-parenthesis right-parenthesis equals sigma Subscript normal t normal r Superscript 2 Baseline phi Subscript normal upper M Baseline left-parenthesis r right-parenthesis equals StartFraction sigma Subscript normal a Baseline Over upper D EndFraction phi Subscript normal upper M Baseline left-parenthesis r right-parenthesis comma

which is exactly the ratio predicted by Equation (15.19); i.e., upper D nabla squared phi Subscript normal upper M Baseline minus sigma Subscript normal a Baseline phi Subscript normal upper M Baseline equals 0 .

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

StartLayout 1st Row 1st Column bold upper E Subscript normal upper M Baseline left-parenthesis normal p Subscript Baseline right-parenthesis 2nd Column equals minus upper D nabla phi Subscript normal upper M Baseline left-parenthesis normal p Subscript Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column equals left-bracket minus upper D StartFraction partial-differential Over partial-differential r EndFraction phi Subscript normal upper M Baseline left-parenthesis r right-parenthesis right-bracket ModifyingAbove bold r With caret 3rd Row 1st Column Blank 2nd Column equals StartFraction 1 plus r sigma Subscript normal t normal r Baseline Over 4 pi r squared EndFraction normal e Superscript minus sigma Super Subscript normal t normal r Superscript r Baseline ModifyingAbove bold r With caret comma EndLayout

where ModifyingAbove bold r With caret is a unit vector pointing away from the source.

15.5.4 Non-classical Diffusion

As we have seen, the monopole fluence rate phi Subscript normal upper M 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:

phi Subscript normal upper G Baseline left-parenthesis r right-parenthesis equals StartFraction normal e Superscript minus sigma prime Super Subscript normal t r Baseline Over 4 pi r squared EndFraction plus ModifyingAbove phi Subscript normal upper M Baseline With tilde left-parenthesis r right-parenthesis period

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 rho prime accounts for the energy reduction due to the extra scattering event:

ModifyingAbove phi Subscript normal upper M Baseline With tilde left-parenthesis r right-parenthesis equals rho prime phi Subscript normal upper M Baseline left-parenthesis r right-parenthesis period

This expression uses the previous monopole solution phi Subscript normal upper M , though its diffusion coefficient upper D must be replaced with a non-classical version given by

upper D Subscript normal upper G Baseline equals StartFraction 2 sigma Subscript normal a Baseline plus sigma prime Subscript normal s Over 3 left-parenthesis sigma Subscript normal a Baseline plus sigma prime Subscript normal s right-parenthesis squared EndFraction period

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

Figure 15.16: Comparison of Classical and Non-Classical Diffusion Solutions. The two plots compare the fluence rate due to a point source for (a) a low albedo ( rho equals 1 slash 3 ) and (b) a high albedo ( rho equals 0.9 ). In both cases, Grosjean’s non-classical monopole (blue) is a considerably better match to the exact solution (black) compared to the classical diffusion solution (red).

In the following sections, we will just focus on the non-classical diffusive part phi Subscript normal upper M Baseline overTilde 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 z equals 0 with a normal of bold n Subscript Baseline equals left-parenthesis 0 comma 0 comma negative 1 right-parenthesis so that positive values on the z axis correspond to points inside the medium. Let eta 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 z -axis at position left-parenthesis 0 comma 0 comma z Subscript normal r Baseline right-parenthesis . Let’s assume that z Subscript normal r Baseline greater-than 0 , 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 z equals 0 ; 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 left-parenthesis 0 comma 0 comma z Subscript normal v Baseline right-parenthesis with z Subscript normal v Baseline less-than 0 . 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 z Subscript normal v 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).

Figure 15.17: Basic Setting for the Dipole Approximation to the Solution to the Diffusion Equation. A light source with positive flux is placed at a position z equals z Subscript normal r inside the medium, below the point where incident illumination arrives, and a second light with an equal amount of negative flux is placed at z equals z Subscript normal v above the medium. These sources are placed so that the flux for the two of them cancels out at a height z Subscript normal e above the boundary, fulfilling the linearized boundary condition. The fluence rate that results from subtracting the closed-form solution for each of them, Equation (15.25), at the medium’s boundary z equals 0 is a reasonable approximation of the fluence rate at the boundary due to subsurface scattering.

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 left-parenthesis r comma 0 comma 0 right-parenthesis on the boundary is simply the sum of a positive and negative phi Subscript normal upper M Baseline overTilde term:

phi Subscript normal upper D Baseline left-parenthesis r right-parenthesis equals ModifyingAbove phi Subscript normal upper M Baseline With tilde left-parenthesis d Subscript normal r Baseline right-parenthesis minus ModifyingAbove phi Subscript normal upper M Baseline With tilde left-parenthesis d Subscript normal v Baseline right-parenthesis

where d Subscript normal r Baseline equals StartRoot r squared plus z Subscript normal r Superscript 2 Baseline EndRoot and d Subscript normal v Baseline equals StartRoot r squared plus z Subscript normal v Superscript 2 Baseline EndRoot 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):

bold upper E Subscript normal upper D Baseline left-parenthesis r right-parenthesis equals ModifyingAbove bold upper E Subscript normal upper M Baseline With tilde left-parenthesis d Subscript normal r Baseline right-parenthesis minus ModifyingAbove bold upper E Subscript normal upper M Baseline With tilde left-parenthesis d Subscript normal v Baseline right-parenthesis

The tilde in bold upper E Subscript normal upper M Baseline overTilde above indicates the use of Grosjean’s modified diffusion coefficient. We will later require the z component of this expression, which is given by

minus bold n Subscript Baseline dot bold upper E Subscript normal upper D Baseline left-parenthesis r right-parenthesis equals StartFraction 1 Over 4 pi EndFraction left-bracket StartFraction z Subscript normal r Baseline left-parenthesis 1 plus d Subscript normal r Baseline sigma Subscript normal t normal r Baseline right-parenthesis Over d Subscript normal r Superscript 3 Baseline EndFraction normal e Superscript minus sigma Super Subscript normal t normal r Superscript d Super Subscript normal r Superscript Baseline minus StartFraction z Subscript normal v Baseline left-parenthesis 1 plus d Subscript normal v Baseline sigma Subscript normal t normal r Baseline right-parenthesis Over d Subscript normal v Superscript 3 Baseline EndFraction normal e Superscript minus sigma Super Subscript normal t normal r Superscript d Super Subscript normal v Superscript Baseline right-bracket period

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 z less-than 0 ). We’ll assume that the real light source depth z Subscript normal r is specified, and that z Subscript normal v 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 z -axis using a first-order Taylor expansion of the boundary conditions at z equals 0 , then this linear function eventually reaches 0 (and then negative values) at a linear extrapolation depth of z Subscript normal e Baseline less-than 0 .

The idea of this approach then is to mirror the real light source across the mirror plane z equals z Subscript normal e defined by this extrapolation depth to obtain the virtual light source depth z Subscript normal v Baseline equals 2 z Subscript normal e Baseline minus z Subscript normal r ; this ensures that the fluence rate of the dipole is exactly equal to 0 at z equals z Subscript normal e 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

z Subscript normal e Baseline equals minus 2 upper D Subscript normal upper G Baseline StartFraction 1 plus 3 ModifyingAbove upper F With bar Subscript normal r comma 2 Baseline left-parenthesis eta right-parenthesis Over 1 minus 2 ModifyingAbove upper F With bar Subscript normal r comma 1 Baseline left-parenthesis eta right-parenthesis EndFraction comma

where upper F overbar Subscript normal r comma 1 and upper F overbar Subscript normal r comma 2 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

upper L Subscript normal d Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis equals StartFraction 1 Over 4 pi EndFraction phi Subscript normal upper D Baseline left-parenthesis double-vertical-bar normal p Subscript Baseline double-vertical-bar right-parenthesis plus StartFraction 3 Over 4 pi EndFraction omega Subscript Baseline dot bold upper E Subscript normal upper D Baseline left-parenthesis double-vertical-bar normal p Subscript Baseline double-vertical-bar right-parenthesis comma

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 upper L Subscript normal d Superscript against the Fresnel transmittance and a cosine factor due to the normal d upper A Subscript Superscript up-tack 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

StartLayout 1st Row 1st Column upper E Subscript normal d Baseline left-parenthesis normal p Subscript Baseline right-parenthesis 2nd Column equals integral Underscript script upper H squared left-parenthesis bold n Subscript Baseline right-parenthesis Endscripts left-parenthesis 1 minus upper F Subscript normal r Baseline left-parenthesis eta Superscript negative 1 Baseline comma cosine theta right-parenthesis right-parenthesis upper L Subscript normal d Baseline left-parenthesis normal p Subscript Baseline comma omega right-parenthesis cosine theta normal d omega Subscript Baseline 2nd Row 1st Column Blank 2nd Column equals upper E Subscript normal d comma phi Sub Subscript normal upper D Subscript Baseline left-parenthesis double-vertical-bar normal p Subscript Baseline double-vertical-bar right-parenthesis plus upper E Subscript normal d comma bold upper E Sub Subscript normal upper D Subscript Baseline left-parenthesis double-vertical-bar normal p Subscript Baseline double-vertical-bar right-parenthesis period EndLayout

When the fluence-related part is written in spherical coordinates, the integral over azimuth turns into a multiplication by 2 pi as no terms depend on it, and phi Subscript normal upper D Baseline left-parenthesis r right-parenthesis can be moved out of the integral. The remaining expression reduces to a constant plus a scaled Fresnel moment (Equation (11.8)).

StartLayout 1st Row 1st Column upper E Subscript normal d comma phi Sub Subscript normal upper D Subscript Baseline left-parenthesis r right-parenthesis 2nd Column equals integral Subscript 0 Superscript 2 pi Baseline integral Subscript 0 Superscript StartFraction pi Over 2 EndFraction Baseline left-parenthesis 1 minus upper F Subscript normal r Baseline left-parenthesis eta Superscript negative 1 Baseline comma cosine theta right-parenthesis right-parenthesis StartFraction 1 Over 4 pi EndFraction phi Subscript normal upper D Baseline left-parenthesis r right-parenthesis cosine theta sine theta normal d theta Subscript Baseline normal d phi Subscript Baseline 2nd Row 1st Column Blank 2nd Column equals one-half phi Subscript normal upper D Baseline left-parenthesis r right-parenthesis integral Subscript 0 Superscript StartFraction pi Over 2 EndFraction Baseline left-parenthesis 1 minus upper F Subscript normal r Baseline left-parenthesis eta Superscript negative 1 Baseline comma cosine theta right-parenthesis right-parenthesis cosine theta sine theta normal d theta Subscript Baseline 3rd Row 1st Column Blank 2nd Column equals phi Subscript normal upper D Baseline left-parenthesis r right-parenthesis left-parenthesis one-fourth minus one-half upper F overbar Subscript normal r comma 1 Baseline right-parenthesis period EndLayout

For the part depending on bold upper E Subscript normal upper D , we obtain

StartLayout 1st Row 1st Column upper E Subscript normal d comma bold upper E Sub Subscript normal upper D Subscript Baseline left-parenthesis r right-parenthesis 2nd Column equals integral Subscript 0 Superscript 2 pi Baseline integral Subscript 0 Superscript StartFraction pi Over 2 EndFraction Baseline left-parenthesis 1 minus upper F Subscript normal r Baseline left-parenthesis eta Superscript negative 1 Baseline comma cosine theta right-parenthesis right-parenthesis left-parenthesis StartFraction 3 Over 4 pi EndFraction omega Subscript Baseline dot bold upper E Subscript normal upper D Baseline left-parenthesis r right-parenthesis right-parenthesis cosine theta sine theta normal d theta Subscript Baseline normal d phi Subscript Baseline 2nd Row 1st Column Blank 2nd Column equals integral Subscript 0 Superscript StartFraction pi Over 2 EndFraction Baseline left-parenthesis 1 minus upper F Subscript normal r Baseline left-parenthesis eta Superscript negative 1 Baseline comma cosine theta right-parenthesis right-parenthesis left-parenthesis StartFraction 3 cosine theta Over 2 EndFraction bold n Subscript Baseline dot bold upper E Subscript normal upper D Baseline left-parenthesis r right-parenthesis right-parenthesis cosine theta sine theta normal d theta Subscript Baseline 3rd Row 1st Column Blank 2nd Column equals bold n Subscript Baseline dot bold upper E Subscript normal upper D Baseline left-parenthesis r right-parenthesis left-parenthesis one-half minus three-halves upper F overbar Subscript normal r comma 2 Baseline right-parenthesis period EndLayout

In summary: after putting together all the pieces, we have a method of evaluating the radiant exitance at a position normal p Subscript on the boundary due to an internal source at depth z Subscript normal r . 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 upper E Subscript normal d Baseline left-parenthesis normal p Subscript Baseline comma z Subscript normal r Baseline right-parenthesis .

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 z Subscript normal r 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., z Subscript normal r Baseline equals 1 slash sigma prime Subscript normal t —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 z Subscript normal r Baseline element-of left-bracket 0 comma normal infinity right-parenthesis 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

upper E Subscript normal d Baseline left-parenthesis normal p Subscript Baseline right-parenthesis equals integral Subscript 0 Superscript normal infinity Baseline sigma prime Subscript normal s Baseline normal e Superscript minus sigma prime Super Subscript normal t Superscript z Super Subscript normal r Baseline upper E Subscript normal d Baseline left-parenthesis normal p Subscript Baseline comma z Subscript normal r Baseline right-parenthesis normal d z Subscript normal r Baseline

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 sigma Subscript normal s Baseline comma sigma Subscript normal a Baseline comma g comma eta , and a radius r and returns an average of 100 samples of the integrand.

<<BSSRDF Utility Functions>>= 
Float BeamDiffusionMS(Float sigma_s, Float sigma_a, Float g, Float eta, Float r) { const int nSamples = 100; Float Ed = 0; <<Precompute information for dipole integrand>>  for (int i = 0; i < nSamples; ++i) { <<Sample real point source depth z Subscript normal r >> 
Float zr = -std::log(1 - (i + .5f) / nSamples) / sigmap_t;
<<Evaluate dipole integrand upper E Subscript normal d at z Subscript normal r and add to Ed>> 
Float zv = -zr + 2 * ze; Float dr = std::sqrt(r * r + zr * zr), dv = std::sqrt(r * r + zv * zv); <<Compute dipole fluence rate phi Subscript normal upper D Baseline left-parenthesis r right-parenthesis using Equation left-parenthesis 15.25 right-parenthesis >> 
Float phiD = Inv4Pi / D_g * (std::exp(-sigma_tr * dr) / dr - std::exp(-sigma_tr * dv) / dv);
<<Compute dipole vector irradiance minus bold n Subscript dot bold upper E Subscript normal upper D Baseline left-parenthesis r right-parenthesis using Equation left-parenthesis 15.27 right-parenthesis >> 
Float EDn = Inv4Pi * (zr * (1 + sigma_tr * dr) * std::exp(-sigma_tr * dr) / (dr*dr*dr) - zv * (1 + sigma_tr * dv) * std::exp(-sigma_tr * dv) / (dv*dv*dv));
<<Add contribution from dipole for depth z Subscript normal r to Ed>> 
Float E = phiD * cPhi + EDn * cE; Float kappa = 1 - std::exp(-2 * sigmap_t * (dr + zr)); Ed += kappa * rhop * rhop * E;
} return Ed / nSamples; }

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

<<Precompute information for dipole integrand>>= 

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

<<Compute reduced scattering coefficients sigma prime Subscript normal s Baseline comma sigma prime Subscript normal t and albedo rho prime >>= 
Float sigmap_s = sigma_s * (1 - g); Float sigmap_t = sigma_a + sigmap_s; Float rhop = sigmap_s / sigmap_t;

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

<<Compute non-classical diffusion coefficient upper D Subscript normal upper G using Equation left-parenthesis 15.24 right-parenthesis >>= 
Float D_g = (2 * sigma_a + sigmap_s) / (3 * sigmap_t * sigmap_t);

<<Compute effective transport coefficient sigma Subscript normal t normal r based on upper D Subscript normal upper G >>= 
Float sigma_tr = std::sqrt(sigma_a / D_g);

Neither the linear extrapolation depth z Subscript normal e nor the scale factors from the radiant exitance computation depend on z Subscript normal r , so they can also be computed. The FresnelMoment1() and FresnelMoment2() functions previously encountered in Section 11.4.1 evaluate upper F overbar Subscript normal r comma 1 and upper F overbar Subscript normal r comma 2 .

<<Determine linear extrapolation distance z Subscript normal e using Equation left-parenthesis 15.28 right-parenthesis >>= 
Float fm1 = FresnelMoment1(eta), fm2 = FresnelMoment2(eta); Float ze = -2 * D_g * (1 + 3 * fm2) / (1 - 2 * fm1);

<<Determine exitance scale factors using Equations left-parenthesis 15.31 right-parenthesis and left-parenthesis 15.32 right-parenthesis >>= 
Float cPhi = .25f * (1 - 2 * fm1), cE = .5f * (1 - 3 * fm2);

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

z Subscript normal r Baseline equals minus StartFraction ln left-parenthesis 1 minus xi Subscript i Baseline right-parenthesis Over sigma prime Subscript normal t EndFraction

for xi Subscript i Baseline element-of left-bracket 0 comma 1 right-parenthesis . Because we are integrating a smooth 1D function, Monte Carlo isn’t required. We therefore use equal-spaced positions left-parenthesis i plus one-half right-parenthesis slash upper N , where 0 less-than-or-equal-to i less-than upper N .

<<Sample real point source depth z Subscript normal r >>= 
Float zr = -std::log(1 - (i + .5f) / nSamples) / sigmap_t;

Given z Subscript normal r , 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 z equals z Subscript normal e .

<<Evaluate dipole integrand upper E Subscript normal d at z Subscript normal r and add to Ed>>= 
Float zv = -zr + 2 * ze; Float dr = std::sqrt(r * r + zr * zr), dv = std::sqrt(r * r + zv * zv); <<Compute dipole fluence rate phi Subscript normal upper D Baseline left-parenthesis r right-parenthesis using Equation left-parenthesis 15.25 right-parenthesis >> 
Float phiD = Inv4Pi / D_g * (std::exp(-sigma_tr * dr) / dr - std::exp(-sigma_tr * dv) / dv);
<<Compute dipole vector irradiance minus bold n Subscript dot bold upper E Subscript normal upper D Baseline left-parenthesis r right-parenthesis using Equation left-parenthesis 15.27 right-parenthesis >> 
Float EDn = Inv4Pi * (zr * (1 + sigma_tr * dr) * std::exp(-sigma_tr * dr) / (dr*dr*dr) - zv * (1 + sigma_tr * dv) * std::exp(-sigma_tr * dv) / (dv*dv*dv));
<<Add contribution from dipole for depth z Subscript normal r to Ed>> 
Float E = phiD * cPhi + EDn * cE; Float kappa = 1 - std::exp(-2 * sigmap_t * (dr + zr)); Ed += kappa * rhop * rhop * E;

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

<<Compute dipole fluence rate phi Subscript normal upper D Baseline left-parenthesis r right-parenthesis using Equation left-parenthesis 15.25 right-parenthesis >>= 
Float phiD = Inv4Pi / D_g * (std::exp(-sigma_tr * dr) / dr - std::exp(-sigma_tr * dv) / dv);

<<Compute dipole vector irradiance minus bold n Subscript dot bold upper E Subscript normal upper D Baseline left-parenthesis r right-parenthesis using Equation left-parenthesis 15.27 right-parenthesis >>= 
Float EDn = Inv4Pi * (zr * (1 + sigma_tr * dr) * std::exp(-sigma_tr * dr) / (dr*dr*dr) - zv * (1 + sigma_tr * dv) * std::exp(-sigma_tr * dv) / (dv*dv*dv));

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 sigma prime Subscript normal s 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

kappa equals 1 minus normal e Superscript minus 2 sigma prime Super Subscript normal t Superscript left-parenthesis d Super Subscript normal r Superscript plus z Super Subscript normal r Superscript right-parenthesis

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

<<Add contribution from dipole for depth z Subscript normal r to Ed>>= 
Float E = phiD * cPhi + EDn * cE; Float kappa = 1 - std::exp(-2 * sigmap_t * (dr + zr)); Ed += kappa * rhop * rhop * E;

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 r equals 0 . 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 normal p 0 after scattering precisely once in the volume at position normal p 1 is given by

upper L Subscript normal s normal s Baseline left-parenthesis normal p 1 right-arrow normal p 0 right-parenthesis equals ModifyingAbove upper P With caret left-parenthesis normal p Subscript Baseline overbar Subscript 2 Baseline right-parenthesis equals integral Underscript script upper P 1 Endscripts upper L Subscript normal e Baseline left-parenthesis normal p Subscript Baseline Subscript 2 Baseline right-arrow normal p Subscript Baseline Subscript 1 Baseline right-parenthesis ModifyingAbove upper T With caret left-parenthesis normal p 0 comma normal p 1 comma normal p 2 right-parenthesis normal d mu 1 left-parenthesis normal p 2 right-parenthesis comma

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

upper L Subscript normal e Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis equals delta left-parenthesis normal p Subscript Baseline minus normal p Subscript normal i Baseline right-parenthesis delta left-parenthesis omega Subscript Baseline plus bold n Subscript Baseline right-parenthesis period

Similar to the previous section, we are interested in the outgoing irradiance at a point normal p Subscript normal o 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 z equals 0 to simplify the derivation but will account for it later on in the final result.

Figure 15.18: Computing the Single Scattering Portion of the BSSRDF Profile. Illumination at perpendicular incidence is scattered once inside the material. Scattering sites before the depth t Subscript normal c normal r normal i normal t undergo internal reflection and cannot directly contribute to single scattering. We therefore integrate the source function over the remaining depths t element-of left-bracket t Subscript normal c normal r normal i normal t Baseline comma normal infinity right-parenthesis to account for the full effect of single scattering.

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

StartLayout 1st Row 1st Column upper E Subscript normal s normal s Baseline left-parenthesis normal p Subscript normal o Baseline right-parenthesis 2nd Column equals integral Underscript script upper P 1 Endscripts upper L Subscript normal s normal s Baseline left-parenthesis normal p 1 right-arrow normal p Subscript normal o Baseline right-parenthesis ModifyingAbove upper G With caret left-parenthesis normal p 1 left-right-arrow normal p Subscript normal o Baseline right-parenthesis normal d mu 1 left-parenthesis normal p 1 right-parenthesis 2nd Row 1st Column Blank 2nd Column equals integral Underscript script upper P 2 Endscripts upper L Subscript normal e Baseline left-parenthesis normal p Subscript Baseline Subscript 2 Baseline right-arrow normal p Subscript Baseline Subscript 1 Baseline right-parenthesis ModifyingAbove upper T With caret left-parenthesis normal p Subscript normal o Baseline comma normal p 1 comma normal p 2 right-parenthesis ModifyingAbove upper G With caret left-parenthesis normal p 1 left-right-arrow normal p Subscript normal o Baseline right-parenthesis normal d mu 2 left-parenthesis normal p 1 comma normal p 2 right-parenthesis period EndLayout

After expanding upper L Subscript normal e Superscript , the spatial delta function from Equation (15.36) removes the integration over normal p 2 , and the directional term reduces Equation (15.37) to a 1D integral along the ray left-parenthesis normal p Subscript normal i Baseline comma minus bold n Subscript Baseline right-parenthesis :

upper E Subscript normal s normal s Baseline left-parenthesis normal p Subscript normal o Baseline right-parenthesis equals integral Subscript 0 Superscript normal infinity Baseline t squared ModifyingAbove upper T With caret left-parenthesis normal p Subscript normal o Baseline comma normal p Subscript normal i Baseline minus t bold n Subscript Baseline comma normal p Subscript normal i Baseline right-parenthesis ModifyingAbove upper G With caret left-parenthesis left-parenthesis normal p Subscript normal i Baseline minus t bold n Subscript Baseline right-parenthesis left-right-arrow normal p Subscript normal o Baseline right-parenthesis normal d t period

The extra t squared term is a change of variables factor resulting from an intermediate step (not shown here), where the volumetric integral over normal p 1 is expressed in terms of spherical coordinates left-parenthesis t comma phi comma theta right-parenthesis . The phi and theta variables subsequently drop out due to the directional delta function delta left-parenthesis omega Subscript Baseline plus bold n Subscript Baseline right-parenthesis in upper L Subscript normal e Superscript .

Expanding the generalized throughput ModifyingAbove upper T With caret using Equation (15.3) yields

StartLayout 1st Row 1st Column upper E Subscript normal s normal s Baseline left-parenthesis normal p Subscript normal o Baseline right-parenthesis equals integral Subscript 0 Superscript normal infinity Baseline t squared 2nd Column ModifyingAbove upper G With caret left-parenthesis normal p Subscript normal i Baseline left-right-arrow left-parenthesis normal p Subscript normal i Baseline minus t bold n Subscript Baseline right-parenthesis right-parenthesis ModifyingAbove f With caret left-parenthesis normal p Subscript normal i Baseline right-arrow left-parenthesis normal p Subscript normal i Baseline minus t bold n Subscript Baseline right-parenthesis right-arrow normal p Subscript Baseline Subscript 0 Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column times ModifyingAbove upper G With caret left-parenthesis left-parenthesis normal p Subscript normal i Baseline minus t bold n Subscript Baseline right-parenthesis left-right-arrow normal p Subscript normal o Baseline right-parenthesis normal d t EndLayout

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 cosine theta Subscript normal s imply that ModifyingAbove f With caret can be written as

ModifyingAbove f With caret left-parenthesis normal p Subscript normal i Baseline right-arrow left-parenthesis normal p Subscript normal i Baseline minus t bold n Subscript Baseline right-parenthesis right-arrow normal p Subscript Baseline Subscript normal o Baseline right-parenthesis equals sigma Subscript s Baseline p left-parenthesis minus cosine theta Subscript normal s Baseline right-parenthesis period

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

cosine theta Subscript normal s Baseline equals StartFraction t Over d EndFraction period

The hypotenuse d is the distance from the scattering location normal p Subscript normal i Baseline minus t bold n Subscript to the exit point normal p Subscript normal o and is given by the Pythagorean theorem:

d equals StartRoot r squared plus t squared EndRoot comma

where r equals double-vertical-bar normal p Subscript normal i Baseline minus normal p Subscript normal o Baseline double-vertical-bar .

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

ModifyingAbove upper G With caret left-parenthesis normal p Subscript normal i Baseline left-right-arrow left-parenthesis normal p Subscript normal i Baseline minus t bold n Subscript Baseline right-parenthesis right-parenthesis equals StartFraction normal e Superscript minus sigma Super Subscript normal t Superscript t Baseline Over t squared EndFraction period

Note that this equation uses the original attenuation coefficient sigma Subscript normal t rather than the reduced version sigma prime Subscript normal t 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 bold n Subscript when intersecting the boundary at normal p Subscript normal o :

ModifyingAbove upper G With caret left-parenthesis left-parenthesis normal p Subscript normal i Baseline minus t bold n Subscript Baseline right-parenthesis left-right-arrow normal p Subscript normal o Baseline right-parenthesis equals StartFraction normal e Superscript minus sigma Super Subscript normal t Superscript d Baseline Over d squared EndFraction StartAbsoluteValue cosine theta Subscript normal o Baseline EndAbsoluteValue period

The two cosine terms cosine theta Subscript normal o and cosine theta Subscript normal s 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:

upper E Subscript normal s normal s Baseline left-parenthesis normal p Subscript normal o Baseline right-parenthesis equals integral Subscript 0 Superscript normal infinity Baseline StartFraction sigma Subscript normal s Baseline normal e Superscript minus sigma Super Subscript normal t Superscript left-parenthesis t plus d right-parenthesis Baseline Over d squared EndFraction p left-parenthesis cosine theta Subscript normal s Baseline right-parenthesis StartAbsoluteValue cosine theta Subscript normal o Baseline EndAbsoluteValue normal d t period

At this point, we’ll reintroduce the effect of the refractive boundary by adding a Fresnel transmission factor left-parenthesis 1 minus upper F Subscript normal r Baseline left-parenthesis eta comma cosine theta Subscript normal o Baseline right-parenthesis right-parenthesis . 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 normal p Subscript normal i is delegated to the Material’s BSDF and thus there is just a single Fresnel transmission term,

upper E Subscript normal s normal s comma upper F Sub Subscript normal r Subscript Baseline left-parenthesis normal p Subscript normal o Baseline right-parenthesis equals integral Subscript 0 Superscript normal infinity Baseline StartFraction sigma Subscript normal s Baseline normal e Superscript minus sigma Super Subscript normal t Superscript left-parenthesis t plus d right-parenthesis Baseline Over d squared EndFraction p left-parenthesis cosine theta Subscript normal s Baseline right-parenthesis left-parenthesis 1 minus upper F Subscript normal r Baseline left-parenthesis eta comma cosine theta Subscript normal o Baseline right-parenthesis right-parenthesis StartAbsoluteValue cosine theta Subscript normal o Baseline EndAbsoluteValue normal d t period

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

theta Subscript normal c normal r normal i normal t Baseline equals sine Superscript negative 1 Baseline StartFraction 1 Over eta EndFraction period

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

cosine theta Subscript normal o Baseline less-than minus cosine theta Subscript normal c normal r normal i normal t Baseline equals minus StartRoot 1 minus StartFraction 1 Over eta squared EndFraction EndRoot period

Combining this with Equation (15.41) and solving for t yields

t greater-than t Subscript normal c normal r normal i normal t Baseline equals r StartRoot eta squared minus 1 EndRoot period

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

t Subscript i Baseline equals t Subscript normal c normal r normal i normal t Baseline minus StartFraction ln left-parenthesis 1 minus xi Subscript i Baseline right-parenthesis Over sigma Subscript normal t Baseline EndFraction

for evenly spaced xi Subscript i Baseline element-of left-bracket 0 comma 1 right-parenthesis . The associated PDF is

p Subscript t Baseline left-parenthesis t right-parenthesis equals sigma Subscript normal t Baseline normal e Superscript minus sigma Super Subscript normal t Superscript left-parenthesis t minus t Super Subscript normal c normal r normal i normal t Superscript right-parenthesis Baseline comma

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

beta equals StartFraction rho normal e Superscript minus sigma Super Subscript normal t Superscript left-parenthesis t Super Subscript normal c normal r normal i normal t Superscript plus d right-parenthesis Baseline Over d squared EndFraction p left-parenthesis minus cosine theta Subscript normal o Baseline right-parenthesis left-parenthesis 1 minus upper F Subscript normal r Baseline left-parenthesis eta comma minus cosine theta Subscript normal o Baseline right-parenthesis right-parenthesis StartAbsoluteValue cosine theta Subscript normal o Baseline EndAbsoluteValue comma

where rho was previously defined as sigma Subscript normal s Baseline slash sigma Subscript normal t .

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

<<BSSRDF Utility Functions>>+=  
Float BeamDiffusionSS(Float sigma_s, Float sigma_a, Float g, Float eta, Float r) { <<Compute material parameters and minimum t below the critical angle>> 
Float sigma_t = sigma_a + sigma_s, rho = sigma_s / sigma_t; Float tCrit = r * std::sqrt(eta * eta - 1);
Float Ess = 0; const int nSamples = 100; for (int i = 0; i < nSamples; ++i) { <<Evaluate single-scattering integrand and add to Ess>> 
Float ti = tCrit - std::log(1 - (i + .5f) / nSamples) / sigma_t; <<Determine length d of connecting segment and cosine theta Subscript normal o >> 
Float d = std::sqrt(r * r + ti * ti); Float cosThetaO = ti / d;
<<Add contribution of single scattering at depth t >> 
Ess += rho * std::exp(-sigma_t * (d + tCrit)) / (d * d) * PhaseHG(cosThetaO, g) * (1 - FrDielectric(-cosThetaO, 1, eta)) * std::abs(cosThetaO);
} return Ess / nSamples; }

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

<<Compute material parameters and minimum t below the critical angle>>= 
Float sigma_t = sigma_a + sigma_s, rho = sigma_s / sigma_t; Float tCrit = r * std::sqrt(eta * eta - 1);

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

<<Evaluate single-scattering integrand and add to Ess>>= 
Float ti = tCrit - std::log(1 - (i + .5f) / nSamples) / sigma_t; <<Determine length d of connecting segment and cosine theta Subscript normal o >> 
Float d = std::sqrt(r * r + ti * ti); Float cosThetaO = ti / d;
<<Add contribution of single scattering at depth t >> 
Ess += rho * std::exp(-sigma_t * (d + tCrit)) / (d * d) * PhaseHG(cosThetaO, g) * (1 - FrDielectric(-cosThetaO, 1, eta)) * std::abs(cosThetaO);

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

<<Determine length d of connecting segment and cosine theta Subscript normal o >>= 
Float d = std::sqrt(r * r + ti * ti); Float cosThetaO = ti / d;

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

<<Add contribution of single scattering at depth t >>= 
Ess += rho * std::exp(-sigma_t * (d + tCrit)) / (d * d) * PhaseHG(cosThetaO, g) * (1 - FrDielectric(-cosThetaO, 1, eta)) * std::abs(cosThetaO);

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 g and relative index of refraction eta 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.

<<BSSRDF Utility Functions>>+=  
void ComputeBeamDiffusionBSSRDF(Float g, Float eta, BSSRDFTable *t) { <<Choose radius values of the diffusion profile discretization>> 
t->radiusSamples[0] = 0; t->radiusSamples[1] = 2.5e-3f; for (int i = 2; i < t->nRadiusSamples; ++i) t->radiusSamples[i] = t->radiusSamples[i - 1] * 1.2f;
<<Choose albedo values of the diffusion profile discretization>> 
for (int i = 0; i < t->nRhoSamples; ++i) t->rhoSamples[i] = (1 - std::exp(-8 * i / (Float)(t->nRhoSamples - 1))) / (1 - std::exp(-8));
ParallelFor( [&](int i) { <<Compute the diffusion profile for the ith albedo sample>> 
<<Compute scattering profile for chosen albedo rho >> 
for (int j = 0; j < t->nRadiusSamples; ++j) { Float rho = t->rhoSamples[i], r = t->radiusSamples[j]; t->profile[i * t->nRadiusSamples + j] = 2 * Pi * r * (BeamDiffusionSS(rho, 1 - rho, g, eta, r) + BeamDiffusionMS(rho, 1 - rho, g, eta, r)); }
<<Compute effective albedo rho Subscript normal e normal f normal f and CDF for importance sampling>> 
}, t->nRhoSamples); }

Both single- and multiple-scattering components of the profile function upper S Subscript normal r Baseline left-parenthesis r right-parenthesis are characterized by an exponential decay with increasing radius r . 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 sigma Subscript normal t . These unitless optical radii were previously introduced in Section 11.4.2.)

<<Choose radius values of the diffusion profile discretization>>= 
t->radiusSamples[0] = 0; t->radiusSamples[1] = 2.5e-3f; for (int i = 2; i < t->nRadiusSamples; ++i) t->radiusSamples[i] = t->radiusSamples[i - 1] * 1.2f;

Next, we need to decide on the locations of upper N albedo samples rho Subscript i on the interval left-bracket 0 comma 1 right-bracket . Some precautions must be taken here as well, since the material’s scattering behavior has a highly nonlinear dependence on the rho parameter.

Consider the example of a medium with a single scattering albedo of rho equals 0.8 : this produces a surprisingly absorptive BSSRDF with an effective albedo rho Subscript normal e normal f normal f of less than 0.15 ! The reason for this behavior is that most incident illumination is scattered many times before it finally leaves the half space; each scattering event incurs an energy reduction by rho , leading to this striking nonlinearity. The left side of Figure 15.19 plots rho Subscript normal e normal f normal f as a function of rho , which shows that most of the interesting behavior is crammed into a small region near rho almost-equals 1 .

Figure 15.19: (a) The relationship between the single scattering albedo rho and the effective albedo rho Subscript normal e normal f normal f is highly nonlinear. (b) We use the reparameterization in Equation (15.48) to achieve a more effective sample placement with an approximately uniform sample spacing in effective albedo space.

Clearly, uniform sample placement rho Subscript i Baseline equals i slash left-parenthesis upper N minus 1 right-parenthesis will not capture this behavior in a satisfactory way. Instead, we use a heuristic that approximately inverts the nonlinear mapping between rho and rho Subscript normal e normal f normal f :

rho Subscript i Baseline equals StartFraction 1 minus normal e Superscript minus 8 i slash left-parenthesis upper N minus 1 right-parenthesis Baseline Over 1 minus normal e Superscript negative 8 Baseline EndFraction comma

<<Choose albedo values of the diffusion profile discretization>>= 
for (int i = 0; i < t->nRhoSamples; ++i) t->rhoSamples[i] = (1 - std::exp(-8 * i / (Float)(t->nRhoSamples - 1))) / (1 - std::exp(-8));

This results in a more perceptually uniform placement of the albedo samples, which can be seen on the right side of Figure 15.19: while not a perfect straight line, the effective albedo associated with increasing sample indices i is now reasonably close to linear.

The loop body then iterates over albedo samples and computes a diffusion profile and associated effective albedo for each one.

<<Compute the diffusion profile for the ith albedo sample>>= 
<<Compute scattering profile for chosen albedo rho >> 
for (int j = 0; j < t->nRadiusSamples; ++j) { Float rho = t->rhoSamples[i], r = t->radiusSamples[j]; t->profile[i * t->nRadiusSamples + j] = 2 * Pi * r * (BeamDiffusionSS(rho, 1 - rho, g, eta, r) + BeamDiffusionMS(rho, 1 - rho, g, eta, r)); }
<<Compute effective albedo rho Subscript normal e normal f normal f and CDF for importance sampling>> 

The first fragment interprets the scattering profile upper S Subscript normal r Baseline left-parenthesis r right-parenthesis as a distribution in polar coordinates left-parenthesis r comma phi right-parenthesis . The marginal distribution over the radius parameter is then given by 2 pi r upper S Subscript normal r Baseline left-parenthesis r right-parenthesis , where upper S Subscript normal r Baseline left-parenthesis r right-parenthesis equals upper S Subscript normal r comma normal s normal s Baseline left-parenthesis r right-parenthesis plus upper S Subscript normal r comma normal m normal s Baseline left-parenthesis r right-parenthesis is the sum of single and multiple scattering profiles.

The reason for tabulating the marginal distribution rather than the profile itself is to facilitate sample generation: in this way, BSSRDFTable::profile can be treated as parameterized sequence of 1D distributions that can be sampled using existing tools like SampleCatmullRom2D() (as is done in TabulatedBSSRDF::Sample_Sr()). However, we must be careful to cancel the extra factor of 2 pi r during normal profile evaluations—this was implemented in the fragment <<Cancel marginal PDF factor from tabulated BSSRDF profile>> in Section 11.4.2.

<<Compute scattering profile for chosen albedo rho >>= 
for (int j = 0; j < t->nRadiusSamples; ++j) { Float rho = t->rhoSamples[i], r = t->radiusSamples[j]; t->profile[i * t->nRadiusSamples + j] = 2 * Pi * r * (BeamDiffusionSS(rho, 1 - rho, g, eta, r) + BeamDiffusionMS(rho, 1 - rho, g, eta, r)); }

Integrating the (unnormalized) marginal PDF over the interval r element-of left-bracket 0 comma normal infinity right-parenthesis yields the effective albedo, which was defined in Equation (11.11). In practice, we limit the integral to a finite interval left-bracket 0 comma r Subscript normal m normal a normal x Baseline right-bracket —this is not an issue, as upper S Subscript normal r Baseline left-parenthesis r right-parenthesis is negligible for r greater-than r Subscript normal m normal a normal x Baseline period The IntegrateCatmullRom() function computes this integral in addition to an auxiliary CDF for importance sampling.

<<Compute effective albedo rho Subscript normal e normal f normal f and CDF for importance sampling>>= 

pbrt uses the inversion method to importance sample piecewise cubic spline functions. This entails first choosing a spline segment according to a discrete probability mass function and then picking a position inside the segment. The IntegrateCatmullRom() function therefore computes a cumulative distribution function that is useful for implementing this sampling operation in an efficient way.

<<Spline Interpolation Definitions>>+= 
Float IntegrateCatmullRom(int n, const Float *x, const Float *values, Float *cdf) { Float sum = 0; cdf[0] = 0; for (int i = 0; i < n - 1; ++i) { <<Look up x Subscript i and function values of spline segment i>> 
Float x0 = x[i], x1 = x[i + 1]; Float f0 = f[i], f1 = f[i + 1]; Float width = x1 - x0;
<<Approximate derivatives using finite differences>> 
Float d0, d1; if (i > 0) d0 = width * (f1 - f[i - 1]) / (x1 - x[i - 1]); else d0 = f1 - f0; if (i + 2 < n) d1 = width * (f[i + 2] - f0) / (x[i + 2] - x0); else d1 = f1 - f0;
<<Keep a running sum and build a cumulative distribution function>> 
sum += ((d0 - d1) * (1.f / 12.f) + (f0 + f1) * .5f) * width; cdf[i + 1] = sum;
} return sum; }

The loop iterates over each spline segment and computes the definite integral of the associated cubic spline interpolant, Equation (8.26), over the associated interval. This definite integral has a simple analytic solution:

integral Subscript x Subscript i Baseline Superscript x Subscript i plus 1 Baseline Baseline p Subscript i Baseline left-parenthesis x right-parenthesis normal d x equals left-parenthesis x Subscript i plus 1 Baseline minus x Subscript i Baseline right-parenthesis left-bracket StartFraction f left-parenthesis x Subscript i Baseline right-parenthesis plus f left-parenthesis x Subscript i plus 1 Baseline right-parenthesis Over 2 EndFraction plus StartFraction f prime left-parenthesis x Subscript i Baseline right-parenthesis minus f prime left-parenthesis x Subscript i plus 1 Baseline right-parenthesis Over 12 EndFraction right-bracket comma

where p Subscript i is the spline interpolant defined on the interval left-bracket x Subscript i Baseline comma x Subscript i plus 1 Baseline right-bracket , the values f left-parenthesis x Subscript i Baseline right-parenthesis and f left-parenthesis x Subscript i plus 1 Baseline right-parenthesis are function evaluations at the endpoints, and f prime left-parenthesis x Subscript i Baseline right-parenthesis and f prime left-parenthesis x Subscript i plus 1 Baseline right-parenthesis are derivative estimates.

The first two fragments in the loop body were previously discussed—they look up the function values (denoted by f0 and f1) and derivative estimates (denoted by d0 and d1) in the current interval and initialize width with the interval length.

We can now compute the definite integral and add it to a running sum. Partial sums are written into the cdf array.

<<Keep a running sum and build a cumulative distribution function>>= 
sum += ((d0 - d1) * (1.f / 12.f) + (f0 + f1) * .5f) * width; cdf[i + 1] = sum;

15.5.9 Setting Scattering Properties

It is remarkably unintuitive to set values of the absorption and scattering coefficients sigma Subscript normal a and sigma Subscript normal s to achieve a desired visual result. If measured values of these parameters aren’t available (e.g., from values from the GetMediumScatteringProperties() utility function of Section 11.4.3), then the task for an artist trying to render subsurface scattering can be difficult.

In this section, we define a convenience function SubsurfaceFromDiffuse() used by the KdSubsurfaceMaterial, which solves an inverse problem using information stored in the BSSRDFTable to derive the medium’s scattering properties. The function takes considerably more intuitive parameters as input: in addition to the BSSRDFTable, it requires an effective albedo and the average distance light travels in the medium before scattering (the mean free path length).

The medium can be made more transparent by increasing the mean free path length or denser by decreasing it. The amount of multiple scattering can be controlled by how close the effective albedo is to 1. The remaining medium properties (index of refraction and scattering anisotropy) are assumed to be fixed.

We perform the inversion separately for each wavelength using InvertCatmullRom() to map from effective albedo to a single scattering albedo rho . With rho known, the desired coefficients are given by sigma Subscript normal s Baseline equals rho sigma Subscript normal t and sigma Subscript normal a Baseline equals left-parenthesis 1 minus rho right-parenthesis sigma Subscript normal t , where sigma Subscript normal t is the reciprocal of the mean free path length.

<<BSSRDF Utility Functions>>+= 
void SubsurfaceFromDiffuse(const BSSRDFTable &t, const Spectrum &rhoEff, const Spectrum &mfp, Spectrum *sigma_a, Spectrum *sigma_s) { for (int c = 0; c < Spectrum::nSamples; ++c) { Float rho = InvertCatmullRom(t.nRhoSamples, t.rhoSamples.get(), t.rhoEff.get(), rhoEff[c]); (*sigma_s)[c] = rho / mfp[c]; (*sigma_a)[c] = (1 - rho) / mfp[c]; } }

InvertCatmullRom() is very similar to the already defined SampleCatmullRom() function except that it directly inverts the spline function and not its definite integral—it otherwise applies the same Newton-Bisection algorithm (and reuses a number of code fragments). We therefore won’t include its implementation here. Note that using this approach requires that the underlying function be either monotonically increasing or monotonically decreasing; in the case of SubsurfaceFromDiffuse(), the function is monotonically increasing.

<<Spline Interpolation Declarations>>= 
Float InvertCatmullRom(int n, const Float *x, const Float *values, Float u);