15.4 Sampling Subsurface Reflection Functions

We’ll now implement techniques to sample the subsurface scattering equation introduced in Section 5.6.2, building on the BSSRDF interface introduced in Section 11.4. Our task is to estimate

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

Figure 15.6 suggests the complexity of evaluating the integral. To compute the standard Monte Carlo estimate of this equation given a point at which to compute outgoing radiance, we need a technique to sample points normal p Subscript normal i on the surface and to compute the incident radiance at these points, as well as an efficient way to compute the specific value of the BSSRDF 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 for each sampled point normal p Subscript normal i and incident direction.

Figure 15.6: Computing Subsurface Reflection. When a surface is translucent, in order to compute outgoing radiance from a point normal p Subscript normal o in direction omega Subscript normal o , it’s necessary to integrate the illumination arriving from directions omega Subscript normal i at nearby points normal p Subscript normal i weighted by the BSSRDF 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 . The BSSRDF can be difficult to evaluate efficiently, since it represents all scattering within the volume for light that enters at one point and exits at the other.

The VolPathIntegrator could be used to evaluate the BSSRDF: given a pair of points on the surface and a pair of directions, the integrator can be used to compute the fraction of incident light from direction omega Subscript normal i at the point normal p Subscript normal i that exits the object at the point normal p Subscript normal o in direction omega Subscript normal o by following light-carrying paths through the multiple scattering events in the medium. Beyond standard path-tracing or bidirectional path-tracing techniques, many other light transport algorithms are applicable to this task.

However, many translucent objects are characterized by having very high albedos, which are not efficiently handled by classic approaches. For example, Jensen et al. (2001b) measured the scattering properties of skim milk and found an albedo of 0.9987 . When essentially all of the light is scattered at each interaction in the medium and almost none of it is absorbed, light easily travels far from where it first enters the medium. Hundreds or even thousands of scattering events must be considered to compute an accurate result; given the high albedo of milk, after 100 scattering events, 87.5 percent-sign of the incident light is still carried by a path, 51 percent-sign after 500 scattering events, and still 26 percent-sign after 1000.

BSSRDF class implementations represent the aggregate scattering behavior of these sorts of media, making it possible to render them fairly efficiently. Figure 15.7 shows an example of the dragon model rendered with a BSSRDF. The main sampling operation that must be provided by implementations of the BSSRDF interface, BSSRDF::Sample_S(), determines the surface position where a ray re-emerges following internal scattering.

Figure 15.7: Subsurface Scattering from the Dragon Model Rendered Using Different Material Densities. (a) Although incident illumination is arriving from behind the model, the front of the model has light exiting from it due to subsurface light transport. (b)  sigma Subscript normal s and sigma Subscript normal a scaled by a factor of 5. (c) scaled by 25. Note how the dragon becomes increasingly opaque as the scattering coefficients increase.

<<BSSRDF Interface>>+= 
virtual Spectrum Sample_S(const Scene &scene, Float u1, const Point2f &u2, MemoryArena &arena, SurfaceInteraction *si, Float *pdf) const = 0;

The value of the BSSRDF for the two points and directions is returned directly, and the associated surface intersection record and probability density are returned via the si and pdf parameters. Two samples must be provided: a 1D sample for discrete sampling decisions (e.g., choosing a specific spectral channel of the profile) and a 2D sample that is mapped onto si. As we will see shortly, it’s useful for BSSRDF implementations to be able to trace rays against the scene geometry to find si, so the scene is also provided as an argument.

15.4.1 Sampling the SeparableBSSRDF

Recall the simplifying assumption introduced in Section 11.4.1, which factored the BSSRDF into spatial and directional components that can be sampled independently from one another. Specifically, Equation (11.6) defined upper S as a product of a single spatial term and a pair of directional terms related to the incident and outgoing directions.

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 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 spatial term upper S Subscript normal p was further simplified to a radial profile function upper S Subscript normal r :

upper S Subscript normal p Baseline left-parenthesis normal p Subscript normal o Baseline comma normal p Subscript normal i Baseline right-parenthesis 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

We will now explain how each of these factors is handled by the SeparableBSSRDF’s sampling routines. This class implements an abstract sampling interface that works for any radial profile function upper S Subscript normal r . The TabulatedBSSRDF class, discussed in Section 15.4.2, derives from SeparableBSSRDF and provides a specific tabulated representation of this profile with support for efficient evaluation and exact importance sampling.

Returning to Equation (15.14), if we assume that the BSSRDF is only sampled for rays that are transmitted through the surface boundary, where transmission is selected with probability left-parenthesis 1 minus upper F Subscript normal r Baseline left-parenthesis cosine theta Subscript normal o Baseline right-parenthesis right-parenthesis , then nothing needs to be done for the left-parenthesis 1 minus upper F Subscript normal r Baseline left-parenthesis cosine theta Subscript normal o Baseline right-parenthesis right-parenthesis part here. (This is the case for the fragment <<Account for subsurface scattering, if applicable>>.) This is a reasonable expectation to place on calling code, as this approach gives good Monte Carlo efficiency.

This leaves the upper S Subscript normal p and upper S Subscript omega terms—the former is handled by a call to SeparableBSSRDF::Sample_Sp() (to be discussed shortly), which returns the position si.

<<BSSRDF Method Definitions>>+=  
Spectrum SeparableBSSRDF::Sample_S(const Scene &scene, Float u1, const Point2f &u2, MemoryArena &arena, SurfaceInteraction *si, Float *pdf) const { Spectrum Sp = Sample_Sp(scene, u1, u2, arena, si, pdf); if (!Sp.IsBlack()) { <<Initialize material model at sampled surface interaction>> 
si->bsdf = ARENA_ALLOC(arena, BSDF)(*si); si->bsdf->Add(ARENA_ALLOC(arena, SeparableBSSRDFAdapter)(this)); si->wo = Vector3f(si->shading.n);
} return Sp; }

If sampling a position is successful, the method initializes si->bsdf with an instance of the class SeparableBSSRDFAdapter, which represents the directional term upper S Subscript omega Baseline left-parenthesis omega Subscript normal i Baseline right-parenthesis as a BxDF. Although this BxDF does not truly depend on the outgoing direction si->wo, we still need to initialize it with a dummy direction.

<<Initialize material model at sampled surface interaction>>= 
si->bsdf = ARENA_ALLOC(arena, BSDF)(*si); si->bsdf->Add(ARENA_ALLOC(arena, SeparableBSSRDFAdapter)(this)); si->wo = Vector3f(si->shading.n);

The SeparableBSSRDFAdapter class is a thin wrapper around SeparableBSSRDF::Sw(). Recall that upper S Subscript omega from Equation (11.7) was defined as a diffuse-like term scaled by the normalized Fresnel transmission. For this reason, the SeparableBSSRDFAdapter classifies itself as BSDF_DIFFUSE and just uses the default cosine-weighted sampling routine provided by BxDF::Sample_f().

<<BSSRDF Declarations>>+= 
class SeparableBSSRDFAdapter : public BxDF { public: <<SeparableBSSRDFAdapter Public Methods>> 
SeparableBSSRDFAdapter(const SeparableBSSRDF *bssrdf) : BxDF(BxDFType(BSDF_REFLECTION | BSDF_DIFFUSE)), bssrdf(bssrdf) { } Spectrum f(const Vector3f &wo, const Vector3f &wi) const { Spectrum f = bssrdf->Sw(wi); <<Update BSSRDF transmission term to account for adjoint light transport>> 
if (bssrdf->mode == TransportMode::Radiance) f *= bssrdf->eta * bssrdf->eta;
return f; }
private: const SeparableBSSRDF *bssrdf; };

<<SeparableBSSRDFAdapter Public Methods>>= 
SeparableBSSRDFAdapter(const SeparableBSSRDF *bssrdf) : BxDF(BxDFType(BSDF_REFLECTION | BSDF_DIFFUSE)), bssrdf(bssrdf) { }

Similar to refractive BSDFs, a scaling factor related to the light transport mode must be applied to the value the f() method returns for the upper S Subscript omega term. This issue is discussed in more detail in Section 16.1, and the fragment that applies this scaling, <<Update BSSRDF transmission term to account for adjoint light transport>>, is defined there.

<<SeparableBSSRDFAdapter Public Methods>>+= 
Spectrum f(const Vector3f &wo, const Vector3f &wi) const { Spectrum f = bssrdf->Sw(wi); <<Update BSSRDF transmission term to account for adjoint light transport>> 
if (bssrdf->mode == TransportMode::Radiance) f *= bssrdf->eta * bssrdf->eta;
return f; }

To sample the spatial component upper S Subscript normal p , we need a way of mapping a 2D distribution function onto an arbitrary surface using a parameterization of the surface in the neighborhood of the outgoing position. A conceptually straightforward way to obtain such a parameterization is by means of geodesics, but finding and evaluating them is non-trivial and requires significant implementation effort for each shape that is supported. We use a much simpler approach that uses ray tracing to map the radial profile upper S Subscript normal r onto the scene geometry.

Figure 15.8 illustrates the basic idea: the position normal p Subscript normal o and associated normal bold n Subscript bold o define a planar approximation to the surface. Using 2D polar coordinates, we first sample an azimuth phi and a radius value r centered around normal p Subscript normal o and then map this position onto the actual surface by intersecting an offset perpendicular ray with the primitive, producing the position normal p Subscript normal i . The SeparableBSSRDF class only supports radially symmetric profile functions; hence phi is drawn from a uniform distribution on left-bracket 0 comma 2 pi right-parenthesis , and r is distributed according to the radial profile function upper S Subscript normal r .

Figure 15.8: Sampling the Spatial Component of the Separable BSSRDF. To compute the outgoing radiance at a point normal p Subscript normal o on a translucent surface, we sample a radius r from a radial scattering profile and map it onto the surface by tracing a probe ray back toward the surface, in the opposite direction of the surface normal bold n Subscript bold o . To improve efficiency, the probe ray is clamped to a sphere of radius r Subscript normal m normal a normal x , after which the value of the BSSRDF becomes negligible.

There are still several difficulties with this basic approach:

  • The radial profile upper S Subscript normal r is not necessarily uniform across wavelengths—in practice, the mean free path can differ by orders of magnitude between different spectral channels.
  • If the surface geometry is poorly approximated by a plane and bold n Subscript bold o Baseline dot bold n Subscript bold i Baseline almost-equals 0 , where bold n Subscript bold i is the surface normal at normal p Subscript normal i , the probe rays will hit the surface at a grazing angle so that positions normal p Subscript normal i with comparatively high values of upper S left-parenthesis normal p Subscript normal o Baseline comma omega Subscript normal o Baseline comma normal p Subscript normal i Baseline comma dot right-parenthesis may be sampled with too low a probability. The result is high variance in renderings (Figure 15.9).
  • Finally, the probe ray may intersect multiple surface locations along its length, all of which may contribute to reflected radiance.

The first two problems can be addressed with a familiar approach; namely, by introducing additional tailored sampling distributions and combining them using multiple importance sampling. The third will be addressed shortly.

Figure 15.9: Comparison of Scattering Profile Projection. (a) Projecting the scattering profile perpendicularly along the current normal direction usually works well, but occasionally there are surface regions that are sampled with a very small probability despite a large corresponding BSSRDF value upper S , producing high variance in renderings. (b) Projecting along multiple axes and combining the resulting sampling techniques using multiple importance sampling greatly reduces the maximum variance at the cost of an increase in the overall amount of variance in well-converged regions. As more samples are taken, the latter approach shows better overall convergence.

We use a different sampling technique per wavelength to deal with spectral variation, and each technique is additionally replicated three times with different projection axes given by the basis vectors of a local frame, resulting in a total of 3 * Spectrum::nSamples sampling techniques. This ensures that every point where upper S takes on non-negligible values is intersected with a reasonable probability. This combination of techniques is implemented in SeparableBSSRDF::Sample_Sp().

<<BSSRDF Method Definitions>>+=  
Spectrum SeparableBSSRDF::Sample_Sp(const Scene &scene, Float u1, const Point2f &u2, MemoryArena &arena, SurfaceInteraction *pi, Float *pdf) const { <<Choose projection axis for BSSRDF sampling>> 
Vector3f vx, vy, vz; if (u1 < .5f) { vx = ss; vy = ts; vz = Vector3f(ns); u1 *= 2; } else if (u1 < .75f) { <<Prepare for sampling rays with respect to ss>> 
vx = ts; vy = Vector3f(ns); vz = ss; u1 = (u1 - .5f) * 4;
} else { <<Prepare for sampling rays with respect to ts>> 
vx = Vector3f(ns); vy = ss; vz = ts; u1 = (u1 - .75f) * 4;
}
<<Choose spectral channel for BSSRDF sampling>> 
int ch = Clamp((int)(u1 * Spectrum::nSamples), 0, Spectrum::nSamples - 1); u1 = u1 * Spectrum::nSamples - ch;
<<Sample BSSRDF profile in polar coordinates>> 
Float r = Sample_Sr(ch, u2[0]); if (r < 0) return Spectrum(0.f); Float phi = 2 * Pi * u2[1];
<<Compute BSSRDF profile bounds and intersection height>> 
Float rMax = Sample_Sr(ch, 0.999f); if (r > rMax) return Spectrum(0.f); Float l = 2 * std::sqrt(rMax * rMax - r * r);
<<Compute BSSRDF sampling ray segment>> 
Interaction base; base.p = po.p + r * (vx * std::cos(phi) + vy * std::sin(phi)) - l * vz * 0.5f; base.time = po.time; Point3f pTarget = base.p + l * vz;
<<Intersect BSSRDF sampling ray against the scene geometry>> 
<<Declare IntersectionChain and linked list>> 
struct IntersectionChain { SurfaceInteraction si; IntersectionChain *next = nullptr; }; IntersectionChain *chain = ARENA_ALLOC(arena, IntersectionChain)();
<<Accumulate chain of intersections along ray>> 
IntersectionChain *ptr = chain; int nFound = 0; while (scene.Intersect(base.SpawnRayTo(pTarget), &ptr->si)) { base = ptr->si; <<Append admissible intersection to IntersectionChain>> 
if (ptr->si.primitive->GetMaterial() == material) { IntersectionChain *next = ARENA_ALLOC(arena, IntersectionChain)(); ptr->next = next; ptr = next; nFound++; }
}
<<Randomly choose one of several intersections during BSSRDF sampling>> 
if (nFound == 0) return Spectrum(0.0f); int selected = Clamp((int)(u1 * nFound), 0, nFound - 1); while (selected-- > 0) chain = chain->next; *pi = chain->si;
<<Compute sample PDF and return the spatial BSSRDF term upper S Subscript normal p >> 
*pdf = Pdf_Sp(*pi) / nFound; return Sp(*pi);
}

We begin by choosing a projection axis. Note that when the surface is close to planar, projecting along the normal SeparableBSSRDF::ns is clearly the best sampling strategy, as probe rays along the other two axes are likely to miss the surface. We therefore allocate a fairly large portion (50%) of the sample budget to perpendicular rays. The other half is equally shared between tangential projections along SeparableBSSRDF::ss and SeparableBSSRDF::ts. The three axes of the chosen coordinate system are stored in vx, vy, and vz, and follow our usual convention of measuring angles theta in spherical coordinates with respect to the z axis.

After this discrete sampling operation, we scale and offset u1 so that additional sampling operations can reuse it as a uniform variate.

<<Choose projection axis for BSSRDF sampling>>= 
Vector3f vx, vy, vz; if (u1 < .5f) { vx = ss; vy = ts; vz = Vector3f(ns); u1 *= 2; } else if (u1 < .75f) { <<Prepare for sampling rays with respect to ss>> 
vx = ts; vy = Vector3f(ns); vz = ss; u1 = (u1 - .5f) * 4;
} else { <<Prepare for sampling rays with respect to ts>> 
vx = Vector3f(ns); vy = ss; vz = ts; u1 = (u1 - .75f) * 4;
}

The fragments for the other two axes are similar and therefore not included here.

Next, we uniformly choose a spectral channel and re-scale u1 once more.

<<Choose spectral channel for BSSRDF sampling>>= 
int ch = Clamp((int)(u1 * Spectrum::nSamples), 0, Spectrum::nSamples - 1); u1 = u1 * Spectrum::nSamples - ch;

The 2D profile sampling operation is then carried out in polar coordinates using the SeparableBSSRDF::Sample_Sr() method. This method returns a negative radius to indicate a failure (e.g., when there is no scattering from channel ch); the implementation here returns a BSSRDF value of 0 in this case.

<<Sample BSSRDF profile in polar coordinates>>= 
Float r = Sample_Sr(ch, u2[0]); if (r < 0) return Spectrum(0.f); Float phi = 2 * Pi * u2[1];

Both the radius sampling method SeparableBSSRDF::Sample_Sr() and its associated density function SeparableBSSRDF::Pdf_Sr() are declared as pure virtual functions; an implementation for TabulatedBSSRDF is presented in the next section.

<<SeparableBSSRDF Interface>>+= 
virtual Float Sample_Sr(int ch, Float u) const = 0; virtual Float Pdf_Sr(int ch, Float r) const = 0;

Because the profile falls off fairly quickly, we are not interested in positions normal p Subscript normal i that are too far away from normal p Subscript normal o . In order to reduce the computational expense of the ray-tracing step, the probe ray is clamped to a sphere of radius r Subscript normal m normal a normal x around normal p Subscript normal o . Another call to SeparableBSSRDF::Sample_Sr() is used to determine r Subscript normal m normal a normal x . Assuming that this function implements a perfect importance sampling scheme based on the inversion method (Section 13.3.1), Sample_Sr() maps a sample value x to the radius of a sphere containing a fraction x of the scattered energy.

Here, we set monospace r monospace upper M monospace a monospace x so that the sphere from Figure 15.8 contains 99.9% of the scattered energy. When r lies outside of r Subscript normal m normal a normal x , sampling fails—this helps to keep the probe rays short, which significantly improves the run-time performance. Given r and r Subscript normal m normal a normal x , the length of the intersection of the probe ray with the sphere of radius r Subscript normal m normal a normal x is

l equals 2 StartRoot r Subscript normal m normal a normal x Superscript 2 Baseline minus r squared EndRoot period

(See Figure 15.10.)

Figure 15.10: Given a sampled radius r that is less than the maximum radius r Subscript normal m normal a normal x , the length of the segment l , which is the total length of the ray in the sphere, can be found using the Pythagorean theorem.

<<Compute BSSRDF profile bounds and intersection height>>= 
Float rMax = Sample_Sr(ch, 0.999f); if (r > rMax) return Spectrum(0.f); Float l = 2 * std::sqrt(rMax * rMax - r * r);

Given the sampled polar coordinate value, we can compute the world space origin of a ray that lies on the boundary of the sphere and a target point, pTarget, where it exits the sphere.

<<Compute BSSRDF sampling ray segment>>= 
Interaction base; base.p = po.p + r * (vx * std::cos(phi) + vy * std::sin(phi)) - l * vz * 0.5f; base.time = po.time; Point3f pTarget = base.p + l * vz;

In practice, there could be more than just one intersection along the probe ray, and we want to collect all of them here. We’ll create a linked list of all of the found interactions.

<<Intersect BSSRDF sampling ray against the scene geometry>>= 
<<Declare IntersectionChain and linked list>> 
struct IntersectionChain { SurfaceInteraction si; IntersectionChain *next = nullptr; }; IntersectionChain *chain = ARENA_ALLOC(arena, IntersectionChain)();
<<Accumulate chain of intersections along ray>> 
IntersectionChain *ptr = chain; int nFound = 0; while (scene.Intersect(base.SpawnRayTo(pTarget), &ptr->si)) { base = ptr->si; <<Append admissible intersection to IntersectionChain>> 
if (ptr->si.primitive->GetMaterial() == material) { IntersectionChain *next = ARENA_ALLOC(arena, IntersectionChain)(); ptr->next = next; ptr = next; nFound++; }
}

IntersectionChain lets us maintain this list. Once again, the MemoryArena makes it possible to efficiently perform allocations, here for the list nodes.

<<Declare IntersectionChain and linked list>>= 
struct IntersectionChain { SurfaceInteraction si; IntersectionChain *next = nullptr; }; IntersectionChain *chain = ARENA_ALLOC(arena, IntersectionChain)();

We now start by finding intersections along the segment within the sphere. The list’s tail node’s SurfaceInteraction is initialized with each intersection’s information, and base Interaction is updated so that the next ray can be spawned on the other side of the intersected surface. (See Figure 15.11.)

Figure 15.11: Accumulating Surface Intersections along a Sample Ray. The SeparableBSSRDF::Sample_Sp() method finds all of the intersections of a ray with the surface of the primitive, where ray extents are limited to a sphere around the intersection point (red dot). At each intersection (blue dots), the corresponding SurfaceInteraction is stored in a linked list before a new ray leaving the other side of the intersected surface is generated.

<<Accumulate chain of intersections along ray>>= 
IntersectionChain *ptr = chain; int nFound = 0; while (scene.Intersect(base.SpawnRayTo(pTarget), &ptr->si)) { base = ptr->si; <<Append admissible intersection to IntersectionChain>> 
if (ptr->si.primitive->GetMaterial() == material) { IntersectionChain *next = ARENA_ALLOC(arena, IntersectionChain)(); ptr->next = next; ptr = next; nFound++; }
}

When tracing rays to sample nearby points on the surface of the primitive, it’s important to ignore any intersections on other primitives in the scene. (There is an implicit assumption that scattering between primitives will be handled by the integrator, and the BSSRDF should be limited to account for single primitives’ scattering.) The implementation here uses equality of Material pointers as a proxy to determine if an intersection is on the same primitive. Valid intersections are appended to the chain, and the variable nFound records their total count when the loop terminates.

<<Append admissible intersection to IntersectionChain>>= 
if (ptr->si.primitive->GetMaterial() == material) { IntersectionChain *next = ARENA_ALLOC(arena, IntersectionChain)(); ptr->next = next; ptr = next; nFound++; }

With the set of intersections at hand, we must now choose one of them, as Sample_Sp() can only return a single position normal p Subscript normal i . The following fragment uses the variable u1 one last time to pick one of the list entries with uniform probability.

<<Randomly choose one of several intersections during BSSRDF sampling>>= 
if (nFound == 0) return Spectrum(0.0f); int selected = Clamp((int)(u1 * nFound), 0, nFound - 1); while (selected-- > 0) chain = chain->next; *pi = chain->si;

Finally, we can call SeparableBSSRDF::Pdf_Sp() (to be defined shortly) to evaluate the combined PDF that takes all of the sampling strategies into account. The probability it returns is divided by nFound to account for the discrete probability of selecting pi from the IntersectionChain. Finally, the value of upper S Subscript normal p Baseline left-parenthesis normal p Subscript normal i Baseline right-parenthesis is returned.

<<Compute sample PDF and return the spatial BSSRDF term upper S Subscript normal p >>= 
*pdf = Pdf_Sp(*pi) / nFound; return Sp(*pi);

SeparableBSSRDF::Pdf_Sp() returns the probability per unit area of sampling the position pi with the total of 3 * Spectrum::nSamples sampling techniques available to SeparableBSSRDF::Sample_Sp().

<<BSSRDF Method Definitions>>+=  
Float SeparableBSSRDF::Pdf_Sp(const SurfaceInteraction &pi) const { <<Express normal p Subscript normal i Baseline minus normal p Subscript normal o and bold n Subscript i with respect to local coordinates at normal p Subscript normal o >> 
Vector3f d = po.p - pi.p; Vector3f dLocal(Dot(ss, d), Dot(ts, d), Dot(ns, d)); Normal3f nLocal(Dot(ss, pi.n), Dot(ts, pi.n), Dot(ns, pi.n));
<<Compute BSSRDF profile radius under projection along each axis>> 
Float rProj[3] = { std::sqrt(dLocal.y * dLocal.y + dLocal.z * dLocal.z), std::sqrt(dLocal.z * dLocal.z + dLocal.x * dLocal.x), std::sqrt(dLocal.x * dLocal.x + dLocal.y * dLocal.y) };
<<Return combined probability from all BSSRDF sampling strategies>> 
Float pdf = 0, axisProb[3] = { .25f, .25f, .5f }; Float chProb = 1 / (Float)Spectrum::nSamples; for (int axis = 0; axis < 3; ++axis) for (int ch = 0; ch < Spectrum::nSamples; ++ch) pdf += Pdf_Sr(ch, rProj[axis]) * std::abs(nLocal[axis]) * chProb * axisProb[axis]; return pdf;
}

First, nLocal is initialized with the surface normal at normal p Subscript normal i and dLocal with the difference vector normal p Subscript normal o Baseline minus normal p Subscript normal i , both expressed using local coordinates at normal p Subscript normal o .

<<Express normal p Subscript normal i Baseline minus normal p Subscript normal o and bold n Subscript i with respect to local coordinates at normal p Subscript normal o >>= 
Vector3f d = po.p - pi.p; Vector3f dLocal(Dot(ss, d), Dot(ts, d), Dot(ns, d)); Normal3f nLocal(Dot(ss, pi.n), Dot(ts, pi.n), Dot(ns, pi.n));

To determine the combined PDF, we must query the probability of sampling a radial profile radius matching the pair left-parenthesis normal p Subscript normal o Baseline comma normal p Subscript normal i Baseline right-parenthesis for each technique. This radius is measured in 2D and thus depends on the chosen projection axis (Figure 15.12). The rProj variable records radii for projections perpendicular to ss, ts, and ns.

Figure 15.12: Probability of Alternative Sampling Strategies. To determine the combined probability of a BSSRDF position sample normal p Subscript normal i (left), we must evaluate the radial PDF for each of the radii r prime corresponding to alternative projection axes (right).

<<Compute BSSRDF profile radius under projection along each axis>>= 
Float rProj[3] = { std::sqrt(dLocal.y * dLocal.y + dLocal.z * dLocal.z), std::sqrt(dLocal.z * dLocal.z + dLocal.x * dLocal.x), std::sqrt(dLocal.x * dLocal.x + dLocal.y * dLocal.y) };

The remainder of the implementation simply loops over all combinations of spectral channels and projection axes and sums up the product of the probability of selecting each technique and its area density under projection onto the surface at normal p Subscript normal o .

<<Return combined probability from all BSSRDF sampling strategies>>= 
Float pdf = 0, axisProb[3] = { .25f, .25f, .5f }; Float chProb = 1 / (Float)Spectrum::nSamples; for (int axis = 0; axis < 3; ++axis) for (int ch = 0; ch < Spectrum::nSamples; ++ch) pdf += Pdf_Sr(ch, rProj[axis]) * std::abs(nLocal[axis]) * chProb * axisProb[axis]; return pdf;

The alert reader may have noticed a slight inconsistency in the above definitions: the probability of choosing one of several (nFound) intersections in SeparableBSSRDF::Sample_Sp() should really have been part of the density function computed in the SeparableBSSRDF::Pdf_Sp() method rather than the ad hoc division that occurs in the fragment <<Compute sample PDF and return the spatial BSSRDF term upper S Subscript normal p >>. In practice, the number of detected intersections varies with respect of the projection axis and spectral channel; correctly accounting for this in the PDF computation requires counting the number of intersections along each of a total of 3 * Spectrum::nSamples probe rays for every sample! We neglect this issue, trading a more efficient implementation for a small amount of bias.

15.4.2 Sampling the TabulatedBSSRDF

The previous section completed the discussion of BSSRDF sampling with the exception of the Pdf_Sr() and Sample_Sr() methods that were declared as pure virtual functions in the SeparableBSSRDF interface. The TabulatedBSSRDF subclass implements this missing functionality.

The TabulatedBSSRDF::Sample_Sr() method samples radius values proportional to the radial profile function upper S Subscript normal r . Recall from Section 11.4.2 that the profile has an implicit dependence on the albedo rho at the current surface position and that the TabulatedBSSRDF provides interpolated evaluations of upper S Subscript normal r Baseline left-parenthesis rho comma r right-parenthesis using 2D tensor product spline basis functions. TabulatedBSSRDF::Sample_Sr() then determines the albedo rho for the given spectral channel ch and draws samples proportional to the remaining 1D function upper S Subscript normal r Baseline left-parenthesis rho comma dot right-parenthesis . Sample generation fails if there is neither scattering nor absorption on channel ch (this case is indicated by returning a negative radius).

As in Section 11.4.2, there is a considerable amount of overlap with the FourierBSDF implementation. The sampling operation here actually reduces to a single call to SampleCatmullRom2D(), which was previously used in FourierBSDF::Sample_f().

<<BSSRDF Method Definitions>>+=  
Float TabulatedBSSRDF::Sample_Sr(int ch, Float u) const { if (sigma_t[ch] == 0) return -1; return SampleCatmullRom2D(table.nRhoSamples, table.nRadiusSamples, table.rhoSamples.get(), table.radiusSamples.get(), table.profile.get(), table.profileCDF.get(), rho[ch], u) / sigma_t[ch]; }

Recall that this function depends on a precomputed CDF array, which is initialized when the BSSRDFTable is created.

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

The Pdf_Sr() method returns the PDF of samples obtained via Sample_Sr(). It evaluates the profile function divided by the normalizing constant rho Subscript normal e normal f normal f defined in Equation (11.11).

The beginning is analogous to the spline evaluation code in TabulatedBSSRDF::Sr(). The fragment <<Compute spline weights to interpolate BSSRDF density on channel ch>> matches <<Compute spline weights to interpolate BSSRDF on channel ch>> in that method except that this method immediately returns zero if the optical radius is outside the range represented by the spline.

<<BSSRDF Method Definitions>>+= 
Float TabulatedBSSRDF::Pdf_Sr(int ch, Float r) const { <<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 density 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)) return 0.f;
<<Return BSSRDF profile density for channel ch>> 
Float sr = 0, rhoEff = 0; for (int i = 0; i < 4; ++i) { if (rhoWeights[i] == 0) continue; rhoEff += table.rhoEff[rhoOffset + i] * rhoWeights[i]; for (int j = 0; j < 4; ++j) { if (radiusWeights[j] == 0) continue; sr += table.EvalProfile(rhoOffset + i, radiusOffset + j) * rhoWeights[i] * radiusWeights[j]; } } <<Cancel marginal PDF factor from tabulated BSSRDF profile>> 
if (rOptical != 0) sr /= 2 * Pi * rOptical;
return std::max((Float)0, sr * sigma_t[ch] * sigma_t[ch] / rhoEff);
}

The remainder of the implementation is very similar to fragment <<Set BSSRDF value Sr[ch] using tensor spline interpolation>> except that here, we also interpolate rho Subscript normal e normal f normal f from the tabulation and include it in the division at the end.

<<Return BSSRDF profile density for channel ch>>= 
Float sr = 0, rhoEff = 0; for (int i = 0; i < 4; ++i) { if (rhoWeights[i] == 0) continue; rhoEff += table.rhoEff[rhoOffset + i] * rhoWeights[i]; for (int j = 0; j < 4; ++j) { if (radiusWeights[j] == 0) continue; sr += table.EvalProfile(rhoOffset + i, radiusOffset + j) * rhoWeights[i] * radiusWeights[j]; } } <<Cancel marginal PDF factor from tabulated BSSRDF profile>> 
if (rOptical != 0) sr /= 2 * Pi * rOptical;
return std::max((Float)0, sr * sigma_t[ch] * sigma_t[ch] / rhoEff);

15.4.3 Subsurface Scattering in the Path Tracer

We now have the capability to apply Monte Carlo integration to the generalized scattering equation, (5.11), from Section 5.6.2. We will compute estimates of the form

upper L Subscript normal o Superscript Baseline left-parenthesis normal p Subscript normal o Baseline comma omega Subscript normal o Baseline right-parenthesis almost-equals StartFraction 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 left-parenthesis upper L Subscript normal d Baseline left-parenthesis normal p Subscript normal i Baseline comma omega Subscript normal i Baseline right-parenthesis plus upper L Subscript normal i Baseline left-parenthesis normal p Subscript normal i Baseline comma omega Subscript normal i Baseline right-parenthesis right-parenthesis StartAbsoluteValue cosine theta Subscript normal i Baseline EndAbsoluteValue Over p left-parenthesis normal p Subscript normal i Baseline right-parenthesis p left-parenthesis omega Subscript normal i Baseline right-parenthesis EndFraction comma

where upper L Subscript normal d represents incident direct radiance and upper L Subscript normal i is incident indirect radiance. The sample left-parenthesis normal p Subscript normal i Baseline comma omega Subscript normal i Baseline right-parenthesis is generated in two steps. First, given normal p Subscript normal o and omega Subscript normal o , a call to BSSRDF::Sample_S() returns a position normal p Subscript normal i whose distribution is similar to the marginal distribution of upper S with respect to normal p Subscript normal i .

Next, we sample the incident direction omega Subscript normal i . Recall that the BSSRDF::Sample_S() interface intentionally keeps these two steps apart: instead of generating both normal p Subscript normal i and omega Subscript normal i at the same time, it returns a special BSDF instance via the bsdf field of si that is used for the direction sampling step. No generality is lost with such an approach: the returned BSDF can be completely arbitrary and is explicitly allowed to depend on information computed within BSSRDF::Sample_S(). The benefit is that we can re-use a considerable amount of existing infrastructure for computing product integrals of a BSDF and upper L Subscript normal i Superscript .

The PathIntegrator’s <<Find next path vertex and accumulate contribution>> fragment invokes the following code to compute this estimate.

<<Account for subsurface scattering, if applicable>>= 
if (isect.bssrdf && (flags & BSDF_TRANSMISSION)) { <<Importance sample the BSSRDF>> 
SurfaceInteraction pi; Spectrum S = isect.bssrdf->Sample_S(scene, sampler.Get1D(), sampler.Get2D(), arena, &pi, &pdf); if (S.IsBlack() || pdf == 0) break; beta *= S / pdf;
<<Account for the direct subsurface scattering component>> 
L += beta * UniformSampleOneLight(pi, scene, arena, sampler);
<<Account for the indirect subsurface scattering component>> 
Spectrum f = pi.bsdf->Sample_f(pi.wo, &wi, sampler.Get2D(), &pdf, BSDF_ALL, &flags); if (f.IsBlack() || pdf == 0) break; beta *= f * AbsDot(wi, pi.shading.n) / pdf; specularBounce = (flags & BSDF_SPECULAR) != 0; ray = pi.SpawnRay(wi);
}

When the BSSRDF sampling case is triggered, the path tracer begins by calling BSSRDF::Sample_S() to generate normal p Subscript normal i and incorporates the resulting sampling weight into its throughput weight variable beta upon success.

<<Importance sample the BSSRDF>>= 
SurfaceInteraction pi; Spectrum S = isect.bssrdf->Sample_S(scene, sampler.Get1D(), sampler.Get2D(), arena, &pi, &pdf); if (S.IsBlack() || pdf == 0) break; beta *= S / pdf;

Because BSSRDF::Sample_S() also initializes pi’s bsdf with a BSDF that characterizes the dependence of upper S on omega Subscript normal i , we can reuse the existing infrastructure for direct illumination computations. Only a single line of code is necessary to compute the contribution of direct lighting at pi to reflected radiance at po.

<<Account for the direct subsurface scattering component>>= 
L += beta * UniformSampleOneLight(pi, scene, arena, sampler);

Similarly, accounting for indirect illumination at the newly sampled incident point is almost the same as the BSDF indirect illumination computation in the PathIntegrator except that pi is used for the next path vertex instead of isect.

<<Account for the indirect subsurface scattering component>>= 
Spectrum f = pi.bsdf->Sample_f(pi.wo, &wi, sampler.Get2D(), &pdf, BSDF_ALL, &flags); if (f.IsBlack() || pdf == 0) break; beta *= f * AbsDot(wi, pi.shading.n) / pdf; specularBounce = (flags & BSDF_SPECULAR) != 0; ray = pi.SpawnRay(wi);

With this, the path tracer (and the volumetric path tracer) support subsurface scattering. See Figure 15.13 for an example.

Figure 15.13: Subsurface Scattering with the PathIntegrator. These dragons both have BSSRDFs that describe subsurface scattering in their interiors. (Model courtesy of Christian Schüller.)