11.4 The BSSRDF

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

upper L Subscript normal o Superscript Baseline left-parenthesis normal p Subscript normal o Baseline comma omega Subscript normal o Baseline right-parenthesis equals integral Underscript upper A Endscripts integral Underscript script upper H squared left-parenthesis bold n Subscript Baseline right-parenthesis Endscripts upper S left-parenthesis normal p Subscript normal o Baseline comma omega Subscript normal o Baseline comma normal p Subscript normal i Baseline comma omega Subscript normal i Baseline right-parenthesis upper L Subscript normal i Superscript Baseline left-parenthesis normal p Subscript normal i Baseline comma omega Subscript normal i Baseline right-parenthesis StartAbsoluteValue cosine theta Subscript normal i Baseline EndAbsoluteValue normal d omega Subscript normal i Baseline normal d upper A Subscript Baseline period

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

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

<<BSSRDF Declarations>>= 
class BSSRDF { public: <<BSSRDF Public Methods>> 
BSSRDF(const SurfaceInteraction &po, Float eta) : po(po), eta(eta) { }
<<BSSRDF Interface>> 
virtual Spectrum S(const SurfaceInteraction &pi, const Vector3f &wi) = 0; virtual Spectrum Sample_S(const Scene &scene, Float u1, const Point2f &u2, MemoryArena &arena, SurfaceInteraction *si, Float *pdf) const = 0;
protected: <<BSSRDF Protected Data>> 
const SurfaceInteraction &po; Float eta;
};

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

<<BSSRDF Public Methods>>= 
BSSRDF(const SurfaceInteraction &po, Float eta) : po(po), eta(eta) { }

<<BSSRDF Protected Data>>= 
const SurfaceInteraction &po; Float eta;

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

<<BSSRDF Interface>>= 
virtual Spectrum S(const SurfaceInteraction &pi, const Vector3f &wi) = 0;

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

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

11.4.1 Separable BSSRDFs

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

<<BSSRDF Declarations>>+=  
class SeparableBSSRDF : public BSSRDF { public: <<SeparableBSSRDF Public Methods>> 
SeparableBSSRDF(const SurfaceInteraction &po, Float eta, const Material *material, TransportMode mode) : BSSRDF(po, eta), ns(po.shading.n), ss(Normalize(po.shading.dpdu)), ts(Cross(ns, ss)), material(material), mode(mode) { } Spectrum S(const SurfaceInteraction &pi, const Vector3f &wi) { Float Ft = 1 - FrDielectric(Dot(po.wo, po.shading.n), 1, eta); return Ft * Sp(pi) * Sw(wi); } Spectrum Sw(const Vector3f &w) const { Float c = 1 - 2 * FresnelMoment1(1 / eta); return (1 - FrDielectric(CosTheta(w), 1, eta)) / (c * Pi); } Spectrum Sp(const SurfaceInteraction &pi) const { return Sr(Distance(po.p, pi.p)); } Spectrum Sample_S(const Scene &scene, Float u1, const Point2f &u2, MemoryArena &arena, SurfaceInteraction *si, Float *pdf) const; Spectrum Sample_Sp(const Scene &scene, Float u1, const Point2f &u2, MemoryArena &arena, SurfaceInteraction *si, Float *pdf) const; Float Pdf_Sp(const SurfaceInteraction &si) const;
<<SeparableBSSRDF Interface>> 
virtual Spectrum Sr(Float d) const = 0; virtual Float Sample_Sr(int ch, Float u) const = 0; virtual Float Pdf_Sr(int ch, Float r) const = 0;
private: <<SeparableBSSRDF Private Data>> 
const Normal3f ns; const Vector3f ss, ts; const Material *material; const TransportMode mode;
};

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

<<SeparableBSSRDF Public Methods>>= 
SeparableBSSRDF(const SurfaceInteraction &po, Float eta, const Material *material, TransportMode mode) : BSSRDF(po, eta), ns(po.shading.n), ss(Normalize(po.shading.dpdu)), ts(Cross(ns, ss)), material(material), mode(mode) { }

<<SeparableBSSRDF Private Data>>= 
const Normal3f ns; const Vector3f ss, ts; const Material *material; const TransportMode mode;

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

upper S left-parenthesis normal p Subscript normal o Baseline comma omega Subscript normal o Baseline comma normal p Subscript normal i Baseline comma omega Subscript normal i Baseline right-parenthesis almost-equals left-parenthesis 1 minus upper F Subscript normal r Baseline left-parenthesis cosine theta Subscript normal o Baseline right-parenthesis right-parenthesis upper S Subscript normal p Baseline left-parenthesis normal p Subscript normal o Baseline comma normal p Subscript normal i Baseline right-parenthesis upper S Subscript omega Baseline left-parenthesis omega Subscript normal i Baseline right-parenthesis period

The Fresnel term at the beginning models the fraction of the light that is transmitted into direction omega Subscript normal o after exiting the material. A second Fresnel term contained inside upper S Subscript omega Baseline left-parenthesis omega Subscript normal i Baseline right-parenthesis accounts for the influence of the boundary on the directional distribution of light entering the object from direction omega Subscript normal i . The profile term upper S Subscript normal p is a spatial distribution characterizing how far light travels after entering the material.

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

<<SeparableBSSRDF Public Methods>>+=  
Spectrum S(const SurfaceInteraction &pi, const Vector3f &wi) { Float Ft = 1 - FrDielectric(Dot(po.wo, po.shading.n), 1, eta); return Ft * Sp(pi) * Sw(wi); }

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

StartLayout 1st Row 1st Column upper L Subscript normal o Superscript Baseline left-parenthesis normal p Subscript normal o Baseline comma omega Subscript normal o Baseline right-parenthesis 2nd Column equals integral Underscript upper A Endscripts integral Underscript script upper H squared left-parenthesis bold n Subscript Baseline right-parenthesis Endscripts upper S left-parenthesis normal p Subscript normal o Baseline comma omega Subscript normal o Baseline comma normal p Subscript normal i Baseline comma omega Subscript normal i Baseline right-parenthesis upper L Subscript normal i Superscript Baseline left-parenthesis normal p Subscript normal i Baseline comma omega Subscript normal i Baseline right-parenthesis StartAbsoluteValue cosine theta Subscript normal i Baseline EndAbsoluteValue normal d omega Subscript normal i Baseline normal d upper A Subscript Baseline left-parenthesis normal p Subscript normal i Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column equals left-parenthesis 1 minus upper F Subscript normal r Baseline left-parenthesis cosine theta Subscript o Baseline right-parenthesis right-parenthesis integral Underscript upper A Endscripts upper S Subscript normal p Baseline left-parenthesis normal p Subscript normal o Baseline comma normal p Subscript normal i Baseline right-parenthesis integral Underscript script upper H squared left-parenthesis bold n Subscript Baseline right-parenthesis Endscripts upper S Subscript omega Baseline left-parenthesis omega Subscript normal i Baseline right-parenthesis upper L Subscript normal i Superscript Baseline left-parenthesis normal p Subscript normal i Baseline comma omega Subscript normal i Baseline right-parenthesis StartAbsoluteValue cosine theta Subscript normal i Baseline EndAbsoluteValue normal d omega Subscript normal i Baseline normal d upper A Subscript Baseline left-parenthesis normal p Subscript normal i Baseline right-parenthesis period EndLayout

We define the directional term upper S Subscript omega Baseline left-parenthesis omega Subscript normal i Baseline right-parenthesis as a scaled version of the Fresnel transmittance (Section 8.2):

upper S Subscript omega Baseline left-parenthesis omega Subscript normal i Baseline right-parenthesis equals StartFraction 1 minus upper F Subscript normal r Baseline left-parenthesis cosine theta Subscript normal i Baseline right-parenthesis Over c pi EndFraction period

The normalization factor c is chosen so that upper S Subscript omega integrates to one over the cosine-weighted hemisphere:

integral Underscript script upper H squared Endscripts upper S Subscript omega Baseline left-parenthesis omega right-parenthesis cosine theta normal d omega Subscript Baseline equals 1 period

In other words,

StartLayout 1st Row 1st Column c 2nd Column equals integral Subscript 0 Superscript 2 pi Baseline integral Subscript 0 Superscript StartFraction pi Over 2 EndFraction Baseline StartFraction 1 minus upper F Subscript normal r Baseline left-parenthesis eta comma cosine theta right-parenthesis Over pi EndFraction sine theta cosine theta normal d theta Subscript Baseline normal d phi Subscript Baseline 2nd Row 1st Column Blank 2nd Column equals 1 minus 2 integral Subscript 0 Superscript StartFraction pi Over 2 EndFraction Baseline upper F Subscript normal r Baseline left-parenthesis eta comma cosine theta right-parenthesis sine theta cosine theta normal d theta Subscript Baseline period EndLayout

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

ModifyingAbove upper F With bar Subscript normal r comma i Baseline left-parenthesis eta right-parenthesis equals integral Subscript 0 Superscript StartFraction pi Over 2 EndFraction Baseline upper F Subscript normal r Baseline left-parenthesis eta comma cosine theta right-parenthesis sine theta cosine Superscript i Baseline theta normal d theta Subscript Baseline

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

<<BSSRDF Utility Declarations>>= 
Float FresnelMoment1(Float invEta); Float FresnelMoment2(Float invEta);

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

<<SeparableBSSRDF Public Methods>>+=  
Spectrum Sw(const Vector3f &w) const { Float c = 1 - 2 * FresnelMoment1(1 / eta); return (1 - FrDielectric(CosTheta(w), 1, eta)) / (c * Pi); }

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

upper S Subscript normal p Baseline left-parenthesis normal p Subscript normal o Baseline comma normal p Subscript normal i Baseline right-parenthesis almost-equals upper S Subscript normal r Baseline left-parenthesis double-vertical-bar normal p Subscript normal o Baseline minus normal p Subscript normal i Baseline double-vertical-bar right-parenthesis period

As before, the actual implementation of the spatial term upper S Subscript normal p doesn’t take normal p Subscript normal o as an argument, since it is already available in BSSRDF::po.

<<SeparableBSSRDF Public Methods>>+= 
Spectrum Sp(const SurfaceInteraction &pi) const { return Sr(Distance(po.p, pi.p)); }

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

<<SeparableBSSRDF Interface>>= 
virtual Spectrum Sr(Float d) const = 0;

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

11.4.2 Tabulated BSSRDF

<<BSSRDF Declarations>>+=  
class TabulatedBSSRDF : public SeparableBSSRDF { public: <<TabulatedBSSRDF Public Methods>> 
TabulatedBSSRDF(const SurfaceInteraction &po, const Material *material, TransportMode mode, Float eta, const Spectrum &sigma_a, const Spectrum &sigma_s, const BSSRDFTable &table) : SeparableBSSRDF(po, eta, material, mode), table(table) { sigma_t = sigma_a + sigma_s; for (int c = 0; c < Spectrum::nSamples; ++c) rho[c] = sigma_t[c] != 0 ? (sigma_s[c] / sigma_t[c]) : 0; } Spectrum Sr(Float distance) const; Float Pdf_Sr(int ch, Float distance) const; Float Sample_Sr(int ch, Float sample) const;
private: <<TabulatedBSSRDF Private Data>> 
const BSSRDFTable &table; Spectrum sigma_t, rho;
};

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

Figure 11.15: Objects rendered using the TabulatedBSSRDF using a variety of measured BSSRDFs. From left to right: cola, apple, skin, and ketchup.

It is important to note that the radial profile upper S Subscript normal r is only a 1D function when all BSSRDF material properties are fixed. More generally, it depends on four additional parameters: the index of refraction eta , the scattering anisotropy g , the albedo rho , and the extinction coefficient sigma Subscript normal t , leading to a complete function signature upper S Subscript normal r Baseline left-parenthesis eta comma g comma rho comma sigma Subscript normal t Baseline comma r right-parenthesis that is unfortunately too high-dimensional for discretization. We must thus either remove or fix some of the parameters.

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

When a lookup for a given extinction coefficient sigma Subscript normal t and radius r occurs at run time we find the corresponding unitless optical radius r Subscript normal o normal p normal t normal i normal c normal a normal l Baseline equals sigma Subscript normal t Baseline r and evaluate the lower-dimensional tabulation as follows:

upper S Subscript normal r Baseline left-parenthesis eta comma g comma rho comma sigma Subscript normal t Baseline comma r right-parenthesis equals sigma Subscript normal t Baseline Superscript 2 Baseline upper S Subscript normal r Baseline left-parenthesis eta comma g comma rho comma 1 comma r Subscript normal o normal p normal t normal i normal c normal a normal l Baseline right-parenthesis

Since upper S Subscript normal r is a 2D density function in polar coordinates left-parenthesis r comma phi right-parenthesis , a corresponding scale factor of sigma Subscript normal t Baseline Superscript 2 is needed to account for this change of variables (see also Section 13.5.2).

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

The TabulatedBSSRDF constructor takes all parameters of the BSSRDF constructor in addition to spectrally varying absorption and scattering coefficients sigma Subscript normal a and sigma Subscript normal s . It precomputes the extinction coefficient sigma Subscript normal t Baseline equals sigma Subscript normal a Baseline plus sigma Subscript normal s and albedo rho equals sigma Subscript normal s Baseline slash sigma Subscript normal t for the current surface location, avoiding issues with a division by zero when there is no extinction for a given spectral channel.

<<TabulatedBSSRDF Public Methods>>= 
TabulatedBSSRDF(const SurfaceInteraction &po, const Material *material, TransportMode mode, Float eta, const Spectrum &sigma_a, const Spectrum &sigma_s, const BSSRDFTable &table) : SeparableBSSRDF(po, eta, material, mode), table(table) { sigma_t = sigma_a + sigma_s; for (int c = 0; c < Spectrum::nSamples; ++c) rho[c] = sigma_t[c] != 0 ? (sigma_s[c] / sigma_t[c]) : 0; }

<<TabulatedBSSRDF Private Data>>= 
const BSSRDFTable &table; Spectrum sigma_t, rho;

Detailed information about the scattering profile upper S Subscript normal r is supplied via the table parameter, which is an instance of the BSSRDFTable data structure:

<<BSSRDF Declarations>>+=  
struct BSSRDFTable { <<BSSRDFTable Public Data>> 
const int nRhoSamples, nRadiusSamples; std::unique_ptr<Float[]> rhoSamples, radiusSamples; std::unique_ptr<Float[]> profile; std::unique_ptr<Float[]> rhoEff; std::unique_ptr<Float[]> profileCDF;
<<BSSRDFTable Public Methods>> 
BSSRDFTable(int nRhoSamples, int nRadiusSamples); inline Float EvalProfile(int rhoIndex, int radiusIndex) const { return profile[rhoIndex * nRadiusSamples + radiusIndex]; }
};

Instances of BSSRDFTable record samples of the function upper S Subscript normal r taken at a set of single scattering albedos left-parenthesis rho 1 comma rho 2 comma ellipsis comma rho Subscript n Baseline right-parenthesis and radii left-parenthesis r 1 comma r 2 comma ellipsis comma r Subscript m Baseline right-parenthesis . The spacing between radius and albedo samples will generally be non-uniform in order to more accurately represent the underlying function. Section 15.5.8 will show how to initialize the BSSRDFTable for a specific BSSRDF model.

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

<<BSSRDFTable Public Methods>>= 
BSSRDFTable(int nRhoSamples, int nRadiusSamples);

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

<<BSSRDFTable Public Data>>= 
const int nRhoSamples, nRadiusSamples; std::unique_ptr<Float[]> rhoSamples, radiusSamples;

A sample value is stored in the profile member variable for each of the m times n pairs left-parenthesis rho Subscript i Baseline comma r Subscript j Baseline right-parenthesis .

<<BSSRDFTable Public Data>>+=  
std::unique_ptr<Float[]> profile;

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

We define the effective albedo as the following integral of the profile upper S Subscript normal r in polar coordinates.

rho Subscript normal e normal f normal f Baseline equals integral Subscript 0 Superscript 2 pi Baseline integral Subscript 0 Superscript normal infinity Baseline r upper S Subscript normal r Baseline left-parenthesis r right-parenthesis normal d r Subscript Baseline normal d phi Subscript Baseline equals 2 pi integral Subscript 0 Superscript normal infinity Baseline r upper S Subscript normal r Baseline left-parenthesis r right-parenthesis normal d r Subscript Baseline period

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

<<BSSRDFTable Public Data>>+=  
std::unique_ptr<Float[]> rhoEff;

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

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

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

<<BSSRDF Method Definitions>>= 
Spectrum TabulatedBSSRDF::Sr(Float r) const { Spectrum Sr(0.f); for (int ch = 0; ch < Spectrum::nSamples; ++ch) { <<Convert r into unitless optical radius r Subscript normal o normal p normal t normal i normal c normal a normal l >> 
Float rOptical = r * sigma_t[ch];
<<Compute spline weights to interpolate BSSRDF on channel ch>> 
int rhoOffset, radiusOffset; Float rhoWeights[4], radiusWeights[4]; if (!CatmullRomWeights(table.nRhoSamples, table.rhoSamples.get(), rho[ch], &rhoOffset, rhoWeights) || !CatmullRomWeights(table.nRadiusSamples, table.radiusSamples.get(), rOptical, &radiusOffset, radiusWeights)) continue;
<<Set BSSRDF value Sr[ch] using tensor spline interpolation>> 
Float sr = 0; for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) { Float weight = rhoWeights[i] * radiusWeights[j]; if (weight != 0) sr += weight * table.EvalProfile(rhoOffset + i, radiusOffset + j); } } <<Cancel marginal PDF factor from tabulated BSSRDF profile>> 
if (rOptical != 0) sr /= 2 * Pi * rOptical;
Sr[ch] = sr;
} <<Transform BSSRDF value into world space units>> 
Sr *= sigma_t * sigma_t;
return Sr.Clamp(); }

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

<<Convert r into unitless optical radius r Subscript normal o normal p normal t normal i normal c normal a normal l >>= 
Float rOptical = r * sigma_t[ch];

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

<<Compute spline weights to interpolate BSSRDF on channel ch>>= 
int rhoOffset, radiusOffset; Float rhoWeights[4], radiusWeights[4]; if (!CatmullRomWeights(table.nRhoSamples, table.rhoSamples.get(), rho[ch], &rhoOffset, rhoWeights) || !CatmullRomWeights(table.nRadiusSamples, table.radiusSamples.get(), rOptical, &radiusOffset, radiusWeights)) continue;

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

<<Set BSSRDF value Sr[ch] using tensor spline interpolation>>= 
Float sr = 0; for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) { Float weight = rhoWeights[i] * radiusWeights[j]; if (weight != 0) sr += weight * table.EvalProfile(rhoOffset + i, radiusOffset + j); } } <<Cancel marginal PDF factor from tabulated BSSRDF profile>> 
if (rOptical != 0) sr /= 2 * Pi * rOptical;
Sr[ch] = sr;

A convenience method in BSSRDFTable helps with finding profile values.

<<BSSRDFTable Public Methods>>+= 
inline Float EvalProfile(int rhoIndex, int radiusIndex) const { return profile[rhoIndex * nRadiusSamples + radiusIndex]; }

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

<<Cancel marginal PDF factor from tabulated BSSRDF profile>>= 
if (rOptical != 0) sr /= 2 * Pi * rOptical;

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

<<Transform BSSRDF value into world space units>>= 
Sr *= sigma_t * sigma_t;

11.4.3 Subsurface Scattering Materials

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

<<SubsurfaceMaterial Declarations>>= 
class SubsurfaceMaterial : public Material { public: <<SubsurfaceMaterial Public Methods>> 
SubsurfaceMaterial(Float scale, const std::shared_ptr<Texture<Spectrum>> &Kr, const std::shared_ptr<Texture<Spectrum>> &Kt, const std::shared_ptr<Texture<Spectrum>> &sigma_a, const std::shared_ptr<Texture<Spectrum>> &sigma_s, Float g, Float eta, const std::shared_ptr<Texture<Float>> &uRoughness, const std::shared_ptr<Texture<Float>> &vRoughness, const std::shared_ptr<Texture<Float>> &bumpMap, bool remapRoughness) : scale(scale), Kr(Kr), Kt(Kt), sigma_a(sigma_a), sigma_s(sigma_s), uRoughness(uRoughness), vRoughness(vRoughness), bumpMap(bumpMap), eta(eta), remapRoughness(remapRoughness), table(100, 64) { ComputeBeamDiffusionBSSRDF(g, eta, &table); } void ComputeScatteringFunctions(SurfaceInteraction *si, MemoryArena &arena, TransportMode mode, bool allowMultipleLobes) const;
private: <<SubsurfaceMaterial Private Data>> 
const Float scale; std::shared_ptr<Texture<Spectrum>> Kr, Kt, sigma_a, sigma_s; std::shared_ptr<Texture<Float>> uRoughness, vRoughness; std::shared_ptr<Texture<Float>> bumpMap; const Float eta; const bool remapRoughness; BSSRDFTable table;
};

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

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

<<SubsurfaceMaterial Private Data>>= 
const Float scale; std::shared_ptr<Texture<Spectrum>> Kr, Kt, sigma_a, sigma_s; std::shared_ptr<Texture<Float>> uRoughness, vRoughness; std::shared_ptr<Texture<Float>> bumpMap; const Float eta; const bool remapRoughness; BSSRDFTable table;

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

<<SubsurfaceMaterial Method Definitions>>= 
void SubsurfaceMaterial::ComputeScatteringFunctions( SurfaceInteraction *si, MemoryArena &arena, TransportMode mode, bool allowMultipleLobes) const { <<Perform bump mapping with bumpMap, if present>> 
if (bumpMap) Bump(bumpMap, si);
<<Initialize BSDF for SubsurfaceMaterial>> 
Spectrum R = Kr->Evaluate(*si).Clamp(); Spectrum T = Kt->Evaluate(*si).Clamp(); Float urough = uRoughness->Evaluate(*si); Float vrough = vRoughness->Evaluate(*si); <<Initialize bsdf for smooth or rough dielectric>> 
si->bsdf = ARENA_ALLOC(arena, BSDF)(*si, eta); if (R.IsBlack() && T.IsBlack()) return; bool isSpecular = urough == 0 && vrough == 0; if (isSpecular && allowMultipleLobes) { si->bsdf->Add(ARENA_ALLOC(arena, FresnelSpecular)(R, T, 1.f, eta, mode)); } else { if (remapRoughness) { urough = TrowbridgeReitzDistribution::RoughnessToAlpha(urough); vrough = TrowbridgeReitzDistribution::RoughnessToAlpha(vrough); } MicrofacetDistribution *distrib = isSpecular ? nullptr : ARENA_ALLOC(arena, TrowbridgeReitzDistribution)(urough, vrough); if (!R.IsBlack()) { Fresnel *fresnel = ARENA_ALLOC(arena, FresnelDielectric)(1.f, eta); if (isSpecular) si->bsdf->Add(ARENA_ALLOC(arena, SpecularReflection)(R, fresnel)); else si->bsdf->Add(ARENA_ALLOC(arena, MicrofacetReflection)(R, distrib, fresnel)); } if (!T.IsBlack()) { if (isSpecular) si->bsdf->Add(ARENA_ALLOC(arena, SpecularTransmission)(T, 1.f, eta, mode)); else si->bsdf->Add(ARENA_ALLOC(arena, MicrofacetTransmission)(T, distrib, 1.f, eta, mode)); } }
Spectrum sig_a = scale * sigma_a->Evaluate(*si).Clamp(); Spectrum sig_s = scale * sigma_s->Evaluate(*si).Clamp(); si->bssrdf = ARENA_ALLOC(arena, TabulatedBSSRDF)( *si, this, mode, eta, sig_a, sig_s, table); }

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

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

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

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

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

<<Media Declarations>>+= 
bool GetMediumScatteringProperties(const std::string &name, Spectrum *sigma_a, Spectrum *sigma_s);