13.4 A Better Path Tracer

The PathIntegrator is based on the same path tracing approach as the SimplePathIntegrator but incorporates a number of improvements. They include these:

  • The direct lighting calculation is performed by sampling both the BSDF and the sampled light source and weighting both samples using multiple importance sampling. This approach can substantially reduce variance compared to sampling the light alone.
  • Any LightSampler can be used, which makes it possible to use effective light sampling algorithms like the one implemented in BVHLightSampler to choose lights.
  • It initializes the VisibleSurface when it is provided, giving geometric information about the first intersection point to Film implementations like GBufferFilm.
  • Russian roulette is used to terminate paths, which can significantly boost the integrator’s efficiency.
  • A technique known as path regularization can be applied in order to reduce variance from difficult-to-sample paths.

While these additions make its implementation more complex, they also substantially improve efficiency; see Figure 13.7 for a comparison of the two.

Figure 13.7: Comparison of the SimplePathIntegrator and the PathIntegrator. (a) Rendered using the SimplePathIntegrator with 64 samples per pixel. (b) The PathIntegrator, also with 64 samples per pixel. Once again, improving the underlying sampling algorithms leads to a substantial reduction in error. Not only is MSE improved by a factor of 1.97 , but execution time is 4.44 times faster, giving an overall efficiency improvement of 8.75 times . (Scene courtesy of Guillermo M. Leal Llaguno.)

The most important of these differences is how the direct lighting calculation is performed. In the SimplePathIntegrator, a light was chosen with uniform probability and then that light sampled a direction; the corresponding estimator was given by Equation (13.9). More generally, the path contribution estimator can be expressed in terms of an arbitrary directional probability distribution p , which gives

upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis almost-equals StartFraction upper L Subscript normal e Baseline left-parenthesis normal p Subscript i Baseline right-arrow normal p Subscript i minus 1 Baseline right-parenthesis f Subscript Baseline left-parenthesis normal p Subscript i Baseline right-arrow normal p Subscript i minus 1 Baseline right-arrow normal p Subscript i minus 2 Baseline right-parenthesis StartAbsoluteValue cosine theta Subscript i Baseline EndAbsoluteValue upper V left-parenthesis normal p Subscript i Baseline left-right-arrow normal p Subscript i minus 1 Baseline right-parenthesis Over p left-parenthesis omega Subscript i Baseline right-parenthesis EndFraction beta period

It may seem that using only a sampling PDF that matches the upper L Subscript normal e factor to sample these directions, as done by the SimplePathIntegrator, would be a good strategy; after all, the radiance upper L Subscript normal e can then be expected to be nonzero for the sampled direction. If we instead drew samples using the BSDF’s sampling distribution, we might choose directions that did not intersect a light source at all, finding no emitted radiance after incurring the expense of tracing a ray in the hope of intersecting a light.

However, there are cases where sampling the BSDF can be the more effective strategy. For a very smooth surface, the BSDF is nonzero for a small set of directions. Sampling the light source will be unlikely to find directions that have a significant effect on scattering from the surface, especially if the light source is large and close by. Even worse, when such a light sample happens to lie in the BSDF lobe, an estimate with large magnitude will be the result due to the combination of a high contribution from the numerator and a small value for the PDF in the denominator. The estimator has high variance.

Figure 13.8 shows a variety of cases where each of these sampling methods is much better than the other. In this scene, four rectangular surfaces ranging from very smooth (top) to very rough (bottom) are illuminated by spherical light sources of decreasing size. Figures 13.8(a) and (b) show the BSDF and light sampling strategies on their own. As the example illustrates, sampling the BSDF is much more effective when it takes on large values on a narrow set of directions that is much smaller than the set of directions that would be obtained by sampling the light sources. This case is most visible in the top left reflection of a large light source in a low-roughness surface. On the other hand, sampling the light sources can be considerably more effective in the opposite case—when the light source is small and the BSDF lobe is less concentrated (this case is most visible in the bottom right reflection).

Figure 13.8: Four surfaces ranging from very smooth (top) to very rough (bottom) illuminated by spherical light sources of decreasing size and rendered with different sampling techniques (modeled after a scene by Eric Veach). (a) BSDF sampling, (b) light sampling, and (c) both techniques combined using MIS. Sampling the BSDF is generally more effective for highly specular materials and large light sources, as illumination is coming from many directions, but the BSDF’s value is large for only a few of them (top left reflection). The converse is true for small sources and rough materials (bottom right reflection), where sampling the light source is more effective.

Taking a single sample with each sampling technique and averaging the estimators would be of limited benefit. The resulting estimator would still have high variance in cases where one of the sampling strategies was ineffective and that strategy happened to sample a direction with nonzero contribution.

This situation is therefore a natural for the application of multiple importance sampling—we have multiple sampling techniques, each of which is sometimes effective and sometimes not. That approach is used in the PathIntegrator with one light sample omega Subscript normal l Baseline tilde p Subscript normal l and one BSDF sample omega Subscript normal b Baseline tilde p Subscript normal b , giving the estimator

StartLayout 1st Row 1st Column upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis almost-equals 2nd Column w Subscript normal l Baseline left-parenthesis omega Subscript normal l Baseline right-parenthesis StartFraction upper L Subscript normal e Baseline left-parenthesis normal p Subscript normal l Baseline right-arrow normal p Subscript i minus 1 Baseline right-parenthesis f Subscript Baseline left-parenthesis normal p Subscript normal l Baseline right-arrow normal p Subscript i minus 1 Baseline right-arrow normal p Subscript i minus 2 Baseline right-parenthesis StartAbsoluteValue cosine theta Subscript normal l Baseline EndAbsoluteValue upper V left-parenthesis normal p Subscript normal l Baseline left-right-arrow normal p Subscript i minus 1 Baseline right-parenthesis Over p Subscript normal l Baseline left-parenthesis omega Subscript normal l Baseline right-parenthesis EndFraction beta plus 2nd Row 1st Column Blank 2nd Column w Subscript normal b Baseline left-parenthesis omega Subscript normal b Baseline right-parenthesis StartFraction upper L Subscript normal e Baseline left-parenthesis normal p Subscript normal b Baseline right-arrow normal p Subscript i minus 1 Baseline right-parenthesis f Subscript Baseline left-parenthesis normal p Subscript normal b Baseline right-arrow normal p Subscript i minus 1 Baseline right-arrow normal p Subscript i minus 2 Baseline right-parenthesis StartAbsoluteValue cosine theta Subscript normal b Baseline EndAbsoluteValue upper V left-parenthesis normal p Subscript normal b Baseline left-right-arrow normal p Subscript i minus 1 Baseline right-parenthesis Over p Subscript normal b Baseline left-parenthesis omega Subscript normal b Baseline right-parenthesis EndFraction beta comma EndLayout

where the surface intersection points corresponding to the two sampled directions are respectively denoted normal p Subscript normal l and normal p Subscript normal b and each term includes a corresponding multiple importance sampling (MIS) weight w Subscript normal l or w Subscript normal b that can be computed, for example, using the balance heuristic from Equation (2.14) or the power heuristic from Equation (2.15). Figure 13.8(c) shows the effectiveness of combining these two sampling techniques with multiple importance sampling.

With that context established, we can start the implementation of the PathIntegrator. It is another RayIntegrator.

<<PathIntegrator Definition>>= 
class PathIntegrator : public RayIntegrator { public: <<PathIntegrator Public Methods>> 
PathIntegrator(int maxDepth, Camera camera, Sampler sampler, Primitive aggregate, std::vector<Light> lights, const std::string &lightSampleStrategy = "bvh", bool regularize = false); SampledSpectrum Li(RayDifferential ray, SampledWavelengths &lambda, Sampler sampler, ScratchBuffer &scratchBuffer, VisibleSurface *visibleSurface) const; static std::unique_ptr<PathIntegrator> Create( const ParameterDictionary &parameters, Camera camera, Sampler sampler, Primitive aggregate, std::vector<Light> lights, const FileLoc *loc); std::string ToString() const;
private: <<PathIntegrator Private Methods>> 
SampledSpectrum SampleLd(const SurfaceInteraction &intr, const BSDF *bsdf, SampledWavelengths &lambda, Sampler sampler) const;
<<PathIntegrator Private Members>> 
int maxDepth; LightSampler lightSampler; bool regularize;
};

Three member variables affect the PathIntegrator’s operation: a maximum path depth; the lightSampler used to sample a light source; and regularize, which controls whether path regularization is used.

<<PathIntegrator Private Members>>= 
int maxDepth; LightSampler lightSampler; bool regularize;

The form of the Li() method is similar to SimplePathIntegrator::Li().

<<PathIntegrator Method Definitions>>= 
SampledSpectrum PathIntegrator::Li(RayDifferential ray, SampledWavelengths &lambda, Sampler sampler, ScratchBuffer &scratchBuffer, VisibleSurface *visibleSurf) const { <<Declare local variables for PathIntegrator::Li()>> 
SampledSpectrum L(0.f), beta(1.f); int depth = 0; Float p_b, etaScale = 1; bool specularBounce = false, anyNonSpecularBounces = false; LightSampleContext prevIntrCtx;
<<Sample path from camera and accumulate radiance estimate>> 
while (true) { <<Trace ray and find closest path vertex and its BSDF>> 
pstd::optional<ShapeIntersection> si = Intersect(ray); <<Add emitted light at intersection point or from the environment>> 
if (!si) { <<Incorporate emission from infinite lights for escaped ray>> 
for (const auto &light : infiniteLights) { SampledSpectrum Le = light.Le(ray, lambda); if (depth == 0 || specularBounce) L += beta * Le; else { <<Compute MIS weight for infinite light>> 
Float p_l = lightSampler.PMF(prevIntrCtx, light) * light.PDF_Li(prevIntrCtx, ray.d, true); Float w_b = PowerHeuristic(1, p_b, 1, p_l);
L += beta * w_b * Le; } }
break; } <<Incorporate emission from surface hit by ray>> 
SampledSpectrum Le = si->intr.Le(-ray.d, lambda); if (Le) { if (depth == 0 || specularBounce) L += beta * Le; else { <<Compute MIS weight for area light>> 
Light areaLight(si->intr.areaLight); Float lightPDF = lightSampler.PMF(prevIntrCtx, areaLight) * areaLight.PDF_Li(prevIntrCtx, ray.d, true); Float w_l = PowerHeuristic(1, bsdfPDF, 1, lightPDF);
L += beta * w_l * Le; } }
SurfaceInteraction &isect = si->intr; <<Get BSDF and skip over medium boundaries>> 
BSDF bsdf = isect.GetBSDF(ray, lambda, camera, scratchBuffer, sampler); if (!bsdf) { isect.SkipIntersection(&ray, si->tHit); continue; }
<<Initialize visibleSurf at first intersection>> 
if (depth == 0 && visibleSurf) { <<Estimate BSDF’s albedo>> 
<<Define sample arrays ucRho and uRho for reflectance estimate>> 
constexpr int nRhoSamples = 16; const Float ucRho[nRhoSamples] = { 0.75741637, 0.37870818, 0.7083487, 0.18935409, 0.9149363, 0.35417435, 0.5990858, 0.09467703, 0.8578725, 0.45746812, 0.686759, 0.17708716, 0.9674518, 0.2995429, 0.5083201, 0.047338516 }; const Point2f uRho[nRhoSamples] = { Point2f(0.855985, 0.570367), Point2f(0.381823, 0.851844), Point2f(0.285328, 0.764262), Point2f(0.733380, 0.114073), Point2f(0.542663, 0.344465), Point2f(0.127274, 0.414848), Point2f(0.964700, 0.947162), Point2f(0.594089, 0.643463), Point2f(0.095109, 0.170369), Point2f(0.825444, 0.263359), Point2f(0.429467, 0.454469), Point2f(0.244460, 0.816459), Point2f(0.756135, 0.731258), Point2f(0.516165, 0.152852), Point2f(0.180888, 0.214174), Point2f(0.898579, 0.503897) };
SampledSpectrum albedo = bsdf.rho(isect.wo, ucRho, uRho);
*visibleSurf = VisibleSurface(isect, albedo, lambda); }
<<Possibly regularize the BSDF>> 
if (regularize && anyNonSpecularBounces) bsdf.Regularize();
<<End path if maximum depth reached>> 
if (depth++ == maxDepth) break;
<<Sample direct illumination from the light sources>> 
if (IsNonSpecular(bsdf.Flags())) { SampledSpectrum Ld = SampleLd(isect, &bsdf, lambda, sampler); L += beta * Ld; }
<<Sample BSDF to get new path direction>> 
Vector3f wo = -ray.d; Float u = sampler.Get1D(); pstd::optional<BSDFSample> bs = bsdf.Sample_f(wo, u, sampler.Get2D()); if (!bs) break; <<Update path state variables after surface scattering>> 
beta *= bs->f * AbsDot(bs->wi, isect.shading.n) / bs->pdf; p_b = bs->pdfIsProportional ? bsdf.PDF(wo, bs->wi) : bs->pdf; specularBounce = bs->IsSpecular(); anyNonSpecularBounces |= !bs->IsSpecular(); if (bs->IsTransmission()) etaScale *= Sqr(bs->eta); prevIntrCtx = si->intr;
ray = isect.SpawnRay(ray, bsdf, bs->wi, bs->flags, bs->eta);
<<Possibly terminate the path with Russian roulette>> 
SampledSpectrum rrBeta = beta * etaScale; if (rrBeta.MaxComponentValue() < 1 && depth > 1) { Float q = std::max<Float>(0, 1 - rrBeta.MaxComponentValue()); if (sampler.Get1D() < q) break; beta /= 1 - q; }
} return L;
}

The L, beta, and depth variables play the same role as the corresponding variables did in the SimplePathIntegrator.

<<Declare local variables for PathIntegrator::Li()>>= 
SampledSpectrum L(0.f), beta(1.f); int depth = 0;

Also similarly, each iteration of the while loop traces a ray to find its closest intersection and its BSDF. Note that a number of code fragments from the SimplePathIntegrator are reused here and in what follows to define the body of the while loop. The loop continues until either the maximum path length is reached or the path is terminated via Russian roulette.

<<Sample path from camera and accumulate radiance estimate>>= 
while (true) { <<Trace ray and find closest path vertex and its BSDF>> 
pstd::optional<ShapeIntersection> si = Intersect(ray); <<Add emitted light at intersection point or from the environment>> 
if (!si) { <<Incorporate emission from infinite lights for escaped ray>> 
for (const auto &light : infiniteLights) { SampledSpectrum Le = light.Le(ray, lambda); if (depth == 0 || specularBounce) L += beta * Le; else { <<Compute MIS weight for infinite light>> 
Float p_l = lightSampler.PMF(prevIntrCtx, light) * light.PDF_Li(prevIntrCtx, ray.d, true); Float w_b = PowerHeuristic(1, p_b, 1, p_l);
L += beta * w_b * Le; } }
break; } <<Incorporate emission from surface hit by ray>> 
SampledSpectrum Le = si->intr.Le(-ray.d, lambda); if (Le) { if (depth == 0 || specularBounce) L += beta * Le; else { <<Compute MIS weight for area light>> 
Light areaLight(si->intr.areaLight); Float lightPDF = lightSampler.PMF(prevIntrCtx, areaLight) * areaLight.PDF_Li(prevIntrCtx, ray.d, true); Float w_l = PowerHeuristic(1, bsdfPDF, 1, lightPDF);
L += beta * w_l * Le; } }
SurfaceInteraction &isect = si->intr; <<Get BSDF and skip over medium boundaries>> 
BSDF bsdf = isect.GetBSDF(ray, lambda, camera, scratchBuffer, sampler); if (!bsdf) { isect.SkipIntersection(&ray, si->tHit); continue; }
<<Initialize visibleSurf at first intersection>> 
if (depth == 0 && visibleSurf) { <<Estimate BSDF’s albedo>> 
<<Define sample arrays ucRho and uRho for reflectance estimate>> 
constexpr int nRhoSamples = 16; const Float ucRho[nRhoSamples] = { 0.75741637, 0.37870818, 0.7083487, 0.18935409, 0.9149363, 0.35417435, 0.5990858, 0.09467703, 0.8578725, 0.45746812, 0.686759, 0.17708716, 0.9674518, 0.2995429, 0.5083201, 0.047338516 }; const Point2f uRho[nRhoSamples] = { Point2f(0.855985, 0.570367), Point2f(0.381823, 0.851844), Point2f(0.285328, 0.764262), Point2f(0.733380, 0.114073), Point2f(0.542663, 0.344465), Point2f(0.127274, 0.414848), Point2f(0.964700, 0.947162), Point2f(0.594089, 0.643463), Point2f(0.095109, 0.170369), Point2f(0.825444, 0.263359), Point2f(0.429467, 0.454469), Point2f(0.244460, 0.816459), Point2f(0.756135, 0.731258), Point2f(0.516165, 0.152852), Point2f(0.180888, 0.214174), Point2f(0.898579, 0.503897) };
SampledSpectrum albedo = bsdf.rho(isect.wo, ucRho, uRho);
*visibleSurf = VisibleSurface(isect, albedo, lambda); }
<<Possibly regularize the BSDF>> 
if (regularize && anyNonSpecularBounces) bsdf.Regularize();
<<End path if maximum depth reached>> 
if (depth++ == maxDepth) break;
<<Sample direct illumination from the light sources>> 
if (IsNonSpecular(bsdf.Flags())) { SampledSpectrum Ld = SampleLd(isect, &bsdf, lambda, sampler); L += beta * Ld; }
<<Sample BSDF to get new path direction>> 
Vector3f wo = -ray.d; Float u = sampler.Get1D(); pstd::optional<BSDFSample> bs = bsdf.Sample_f(wo, u, sampler.Get2D()); if (!bs) break; <<Update path state variables after surface scattering>> 
beta *= bs->f * AbsDot(bs->wi, isect.shading.n) / bs->pdf; p_b = bs->pdfIsProportional ? bsdf.PDF(wo, bs->wi) : bs->pdf; specularBounce = bs->IsSpecular(); anyNonSpecularBounces |= !bs->IsSpecular(); if (bs->IsTransmission()) etaScale *= Sqr(bs->eta); prevIntrCtx = si->intr;
ray = isect.SpawnRay(ray, bsdf, bs->wi, bs->flags, bs->eta);
<<Possibly terminate the path with Russian roulette>> 
SampledSpectrum rrBeta = beta * etaScale; if (rrBeta.MaxComponentValue() < 1 && depth > 1) { Float q = std::max<Float>(0, 1 - rrBeta.MaxComponentValue()); if (sampler.Get1D() < q) break; beta /= 1 - q; }
} return L;

We will defer discussing the implementation of the first fragment used below, <<Add emitted light at intersection point or from the environment>>, until later in this section after more details of the implementation of the MIS direct lighting calculation have been introduced.

<<Trace ray and find closest path vertex and its BSDF>>= 
pstd::optional<ShapeIntersection> si = Intersect(ray); <<Add emitted light at intersection point or from the environment>> 
if (!si) { <<Incorporate emission from infinite lights for escaped ray>> 
for (const auto &light : infiniteLights) { SampledSpectrum Le = light.Le(ray, lambda); if (depth == 0 || specularBounce) L += beta * Le; else { <<Compute MIS weight for infinite light>> 
Float p_l = lightSampler.PMF(prevIntrCtx, light) * light.PDF_Li(prevIntrCtx, ray.d, true); Float w_b = PowerHeuristic(1, p_b, 1, p_l);
L += beta * w_b * Le; } }
break; } <<Incorporate emission from surface hit by ray>> 
SampledSpectrum Le = si->intr.Le(-ray.d, lambda); if (Le) { if (depth == 0 || specularBounce) L += beta * Le; else { <<Compute MIS weight for area light>> 
Light areaLight(si->intr.areaLight); Float lightPDF = lightSampler.PMF(prevIntrCtx, areaLight) * areaLight.PDF_Li(prevIntrCtx, ray.d, true); Float w_l = PowerHeuristic(1, bsdfPDF, 1, lightPDF);
L += beta * w_l * Le; } }
SurfaceInteraction &isect = si->intr; <<Get BSDF and skip over medium boundaries>> 
BSDF bsdf = isect.GetBSDF(ray, lambda, camera, scratchBuffer, sampler); if (!bsdf) { isect.SkipIntersection(&ray, si->tHit); continue; }
<<Initialize visibleSurf at first intersection>> 
if (depth == 0 && visibleSurf) { <<Estimate BSDF’s albedo>> 
<<Define sample arrays ucRho and uRho for reflectance estimate>> 
constexpr int nRhoSamples = 16; const Float ucRho[nRhoSamples] = { 0.75741637, 0.37870818, 0.7083487, 0.18935409, 0.9149363, 0.35417435, 0.5990858, 0.09467703, 0.8578725, 0.45746812, 0.686759, 0.17708716, 0.9674518, 0.2995429, 0.5083201, 0.047338516 }; const Point2f uRho[nRhoSamples] = { Point2f(0.855985, 0.570367), Point2f(0.381823, 0.851844), Point2f(0.285328, 0.764262), Point2f(0.733380, 0.114073), Point2f(0.542663, 0.344465), Point2f(0.127274, 0.414848), Point2f(0.964700, 0.947162), Point2f(0.594089, 0.643463), Point2f(0.095109, 0.170369), Point2f(0.825444, 0.263359), Point2f(0.429467, 0.454469), Point2f(0.244460, 0.816459), Point2f(0.756135, 0.731258), Point2f(0.516165, 0.152852), Point2f(0.180888, 0.214174), Point2f(0.898579, 0.503897) };
SampledSpectrum albedo = bsdf.rho(isect.wo, ucRho, uRho);
*visibleSurf = VisibleSurface(isect, albedo, lambda); }
<<Possibly regularize the BSDF>> 
if (regularize && anyNonSpecularBounces) bsdf.Regularize();

If the Film being used takes a VisibleSurface, then a non-nullptr VisibleSurface * is passed to the Li() method. It is initialized at the first intersection.

<<Initialize visibleSurf at first intersection>>= 
if (depth == 0 && visibleSurf) { <<Estimate BSDF’s albedo>> 
<<Define sample arrays ucRho and uRho for reflectance estimate>> 
constexpr int nRhoSamples = 16; const Float ucRho[nRhoSamples] = { 0.75741637, 0.37870818, 0.7083487, 0.18935409, 0.9149363, 0.35417435, 0.5990858, 0.09467703, 0.8578725, 0.45746812, 0.686759, 0.17708716, 0.9674518, 0.2995429, 0.5083201, 0.047338516 }; const Point2f uRho[nRhoSamples] = { Point2f(0.855985, 0.570367), Point2f(0.381823, 0.851844), Point2f(0.285328, 0.764262), Point2f(0.733380, 0.114073), Point2f(0.542663, 0.344465), Point2f(0.127274, 0.414848), Point2f(0.964700, 0.947162), Point2f(0.594089, 0.643463), Point2f(0.095109, 0.170369), Point2f(0.825444, 0.263359), Point2f(0.429467, 0.454469), Point2f(0.244460, 0.816459), Point2f(0.756135, 0.731258), Point2f(0.516165, 0.152852), Point2f(0.180888, 0.214174), Point2f(0.898579, 0.503897) };
SampledSpectrum albedo = bsdf.rho(isect.wo, ucRho, uRho);
*visibleSurf = VisibleSurface(isect, albedo, lambda); }

The only quantity that is not immediately available from the SurfaceInteraction is the albedo of the surface, which is computed here as the hemispherical-directional reflectance, Equation (4.12). Recall that the BSDF::rho() method estimates this value using Monte Carlo integration. Here, a set of 16 precomputed Owen-scrambled Halton points in arrays ucRho and uRho, not included in the text, are used for the estimate.

The use of Monte Carlo with this many samples is somewhat unsatisfying. The computed albedo is most commonly used for image-space denoising algorithms after rendering; most of these start by dividing the final color at each pixel by the first visible surface’s albedo in order to approximate the incident illumination alone. It is therefore important that the albedo value itself not have very much error. However, the albedo can be computed analytically for some BSDFs (e.g., the ideal Lambertian BRDF). In those cases, executing both the BSDF sampling and evaluation algorithms repeatedly is wasteful. An exercise at the end of the chapter discusses this matter further.

<<Estimate BSDF’s albedo>>= 
<<Define sample arrays ucRho and uRho for reflectance estimate>> 
constexpr int nRhoSamples = 16; const Float ucRho[nRhoSamples] = { 0.75741637, 0.37870818, 0.7083487, 0.18935409, 0.9149363, 0.35417435, 0.5990858, 0.09467703, 0.8578725, 0.45746812, 0.686759, 0.17708716, 0.9674518, 0.2995429, 0.5083201, 0.047338516 }; const Point2f uRho[nRhoSamples] = { Point2f(0.855985, 0.570367), Point2f(0.381823, 0.851844), Point2f(0.285328, 0.764262), Point2f(0.733380, 0.114073), Point2f(0.542663, 0.344465), Point2f(0.127274, 0.414848), Point2f(0.964700, 0.947162), Point2f(0.594089, 0.643463), Point2f(0.095109, 0.170369), Point2f(0.825444, 0.263359), Point2f(0.429467, 0.454469), Point2f(0.244460, 0.816459), Point2f(0.756135, 0.731258), Point2f(0.516165, 0.152852), Point2f(0.180888, 0.214174), Point2f(0.898579, 0.503897) };
SampledSpectrum albedo = bsdf.rho(isect.wo, ucRho, uRho);

The next task is to sample a light source to find a direction omega Subscript normal i to use to estimate the first term of Equation (13.10). However, if the BSDF is purely specular, there is no reason to do this work, since the value of the BSDF for a sampled point on a light will certainly be zero.

<<Sample direct illumination from the light sources>>= 
if (IsNonSpecular(bsdf.Flags())) { SampledSpectrum Ld = SampleLd(isect, &bsdf, lambda, sampler); L += beta * Ld; }

Although SampleLd() is only called once and thus could be expanded inline in the Li() method, there are multiple points along the way where it may return early. We therefore prefer a function here, as it avoids deeply nested if statements that would be needed otherwise.

<<PathIntegrator Method Definitions>>+= 
SampledSpectrum PathIntegrator::SampleLd( const SurfaceInteraction &intr, const BSDF *bsdf, SampledWavelengths &lambda, Sampler sampler) const { <<Initialize LightSampleContext for light sampling>> 
LightSampleContext ctx(intr); <<Try to nudge the light sampling position to correct side of the surface>> 
BxDFFlags flags = bsdf->Flags(); if (IsReflective(flags) && !IsTransmissive(flags)) ctx.pi = intr.OffsetRayOrigin(intr.wo); else if (IsTransmissive(flags) && !IsReflective(flags)) ctx.pi = intr.OffsetRayOrigin(-intr.wo);
<<Choose a light source for the direct lighting calculation>> 
Float u = sampler.Get1D(); pstd::optional<SampledLight> sampledLight = lightSampler.Sample(ctx, u); Point2f uLight = sampler.Get2D(); if (!sampledLight) return {};
<<Sample a point on the light source for direct lighting>> 
Light light = sampledLight->light; pstd::optional<LightLiSample> ls = light.SampleLi(ctx, uLight, lambda, true); if (!ls || !ls->L || ls->pdf == 0) return {};
<<Evaluate BSDF for light sample and check light visibility>> 
Vector3f wo = intr.wo, wi = ls->wi; SampledSpectrum f = bsdf->f(wo, wi) * AbsDot(wi, intr.shading.n); if (!f || !Unoccluded(intr, ls->pLight)) return {};
<<Return light’s contribution to reflected radiance>> 
Float p_l = sampledLight->p * ls->pdf; if (IsDeltaLight(light.Type())) return ls->L * f / p_l; else { Float p_b = bsdf->PDF(wo, wi); Float w_l = PowerHeuristic(1, p_l, 1, p_b); return w_l * ls->L * f / p_l; }
}

A LightSampleContext is necessary both for choosing a specific light source and for sampling a point on it. One is initialized using the constructor that takes a SurfaceInteraction.

<<Initialize LightSampleContext for light sampling>>= 
LightSampleContext ctx(intr); <<Try to nudge the light sampling position to correct side of the surface>> 
BxDFFlags flags = bsdf->Flags(); if (IsReflective(flags) && !IsTransmissive(flags)) ctx.pi = intr.OffsetRayOrigin(intr.wo); else if (IsTransmissive(flags) && !IsReflective(flags)) ctx.pi = intr.OffsetRayOrigin(-intr.wo);

If the surface is purely reflective or purely transmissive, then the reference point used for sampling pi is shifted slightly so that it lies on the side of the surface from which the outgoing ray will leave the intersection point toward the light. Doing so helps avoid a subtle error that is the result of the combination of floating-point round-off error in the computed intersection point and a ray that intersects an emitter that does not have a completely absorbing BSDF. The problem is illustrated in Figure 13.9.

Figure 13.9: The LightSampleContext stores the error bounds around the computed intersection point pi. Typically, the center of these bounds (filled circle) is used as the reference point for sampling a point on the light source. If a ray intersects the non-emissive side of a one-sided light, the light’s BSDF is nonzero, and if the center of the pi bounds is on the emissive side of the light, then it may seem that the intersection point is illuminated by the light. The result is an occasional bright pixel on the back side of light sources. Offsetting the reference point to the side of the surface from which the outgoing ray will leave (open circle) works around this problem.

<<Try to nudge the light sampling position to correct side of the surface>>= 
BxDFFlags flags = bsdf->Flags(); if (IsReflective(flags) && !IsTransmissive(flags)) ctx.pi = intr.OffsetRayOrigin(intr.wo); else if (IsTransmissive(flags) && !IsReflective(flags)) ctx.pi = intr.OffsetRayOrigin(-intr.wo);

Next, the LightSampler selects a light. One thing to note in the implementation here is that two more dimensions are consumed from the Sampler even if the LightSampler does not return a valid light. This is done in order to keep the allocation of Sampler dimensions consistent across all the pixel samples. (Recall the discussion of this issue in Section 8.3.)

<<Choose a light source for the direct lighting calculation>>= 
Float u = sampler.Get1D(); pstd::optional<SampledLight> sampledLight = lightSampler.Sample(ctx, u); Point2f uLight = sampler.Get2D(); if (!sampledLight) return {};

Sampling a direction with the light proceeds using Light::SampleLi(), though here a true value is passed for its allowIncompletePDF parameter. Because we will use a second sampling technique, BSDF sampling, for the estimator in Equation (13.10), and that technique has nonzero probability of sampling all directions omega Subscript normal i where the integrand is nonzero, the light sampling distribution may not include directions where the light’s emission is relatively low. (The motivation for this was discussed in Section 2.2.3 in the context of MIS compensation.)

Given a light sample, it is worth checking for various cases that require no further processing here. As an example, consider a spotlight where the intersection point is outside of its emission cone; the LightLiSample will have a zero radiance value in that case. It is worthwhile to find that there is no incident radiance before incurring the cost of evaluating the BSDF.

<<Sample a point on the light source for direct lighting>>= 
Light light = sampledLight->light; pstd::optional<LightLiSample> ls = light.SampleLi(ctx, uLight, lambda, true); if (!ls || !ls->L || ls->pdf == 0) return {};

A shadow ray is only traced if the BSDF for the sampled direction is nonzero. It is not unusual for the BSDF to be zero here: for example, given a surface that is reflective but not transmissive, any sampled direction that is on the other side of the surface than the incident ray will have zero contribution.

<<Evaluate BSDF for light sample and check light visibility>>= 
Vector3f wo = intr.wo, wi = ls->wi; SampledSpectrum f = bsdf->f(wo, wi) * AbsDot(wi, intr.shading.n); if (!f || !Unoccluded(intr, ls->pLight)) return {};

Figure 13.10: Comparison of the Balance and Power Heuristics for Direct Lighting. A zoomed-in region of Figure 13.8 is shown here. (a) Rendered using the balance heuristic to weight BSDF and light samples in the direct lighting calculation. (b) Rendered using the power heuristic. The pixels behind the light source have a visible reduction in noise.

The light sample’s contribution can now be computed; recall that the returned value corresponds to the first term of Equation (13.10), save for the beta factor. The case of a light that is described by a delta distribution receives special treatment here; recall from Section 12.1 that in that case there is an implied delta distribution in the emitted radiance value returned from SampleLi() as well as the PDF and that they cancel out when the estimator is evaluated. Further, BSDF sampling is unable to generate a light sample and therefore we must not try to apply multiple importance sampling but should evaluate the standard estimator, Equation (13.9), instead. If we do not have a delta distribution light source, then the value of the BSDF’s PDF for sampling the direction omega Subscript normal i is found by calling BSDF::PDF() and the MIS weight is computed using the power heuristic. (See Figure 13.10 for a comparison between the balance heuristic and power heuristic for this computation.)

<<Return light’s contribution to reflected radiance>>= 
Float p_l = sampledLight->p * ls->pdf; if (IsDeltaLight(light.Type())) return ls->L * f / p_l; else { Float p_b = bsdf->PDF(wo, wi); Float w_l = PowerHeuristic(1, p_l, 1, p_b); return w_l * ls->L * f / p_l; }

Returning now to the Li() method implementation, the next step is to sample the BSDF at the intersection to get an outgoing direction for the next ray to trace. That ray will be used to sample indirect illumination as well as for the BSDF sample for the direct lighting estimator.

<<Sample BSDF to get new path direction>>= 
Vector3f wo = -ray.d; Float u = sampler.Get1D(); pstd::optional<BSDFSample> bs = bsdf.Sample_f(wo, u, sampler.Get2D()); if (!bs) break; <<Update path state variables after surface scattering>> 
beta *= bs->f * AbsDot(bs->wi, isect.shading.n) / bs->pdf; p_b = bs->pdfIsProportional ? bsdf.PDF(wo, bs->wi) : bs->pdf; specularBounce = bs->IsSpecular(); anyNonSpecularBounces |= !bs->IsSpecular(); if (bs->IsTransmission()) etaScale *= Sqr(bs->eta); prevIntrCtx = si->intr;
ray = isect.SpawnRay(ray, bsdf, bs->wi, bs->flags, bs->eta);

In addition to the path throughput weight beta, a number of additional values related to the path are maintained, as follows:

  • p_b is the PDF for sampling the direction bs->wi; this value is needed for the MIS-based direct lighting estimate. One nit comes from BSDFs like the LayeredBxDF that return a BSDFSample where the f and pdf are only proportional to their true values. In that case, an explicit call to BSDF::PDF() is required to get an estimate of the true PDF.
  • As in the SimplePathIntegrator, specularBounce tracks whether the last scattering event was from a perfect specular surface.
  • anyNonSpecularBounces tracks whether any scattering event along the ray’s path has been non-perfect specular. This value is used for path regularization if it is enabled.
  • etaScale is the accumulated product of scaling factors that have been applied to beta due to rays being transmitted between media of different indices of refraction—a detail that is discussed in Section 9.5.2. This value will be used in the Russian roulette computation.
  • Finally, prevIntrCtx stores geometric information about the intersection point from which the sampled ray is leaving. This value is also used in the MIS computation for direct lighting.

<<Update path state variables after surface scattering>>= 
beta *= bs->f * AbsDot(bs->wi, isect.shading.n) / bs->pdf; p_b = bs->pdfIsProportional ? bsdf.PDF(wo, bs->wi) : bs->pdf; specularBounce = bs->IsSpecular(); anyNonSpecularBounces |= !bs->IsSpecular(); if (bs->IsTransmission()) etaScale *= Sqr(bs->eta); prevIntrCtx = si->intr;

<<Declare local variables for PathIntegrator::Li()>>+= 
Float p_b, etaScale = 1; bool specularBounce = false, anyNonSpecularBounces = false; LightSampleContext prevIntrCtx;

The new ray will account for indirect illumination at the intersection point in the following execution of the while loop.

Returning now to the <<Add emitted light at intersection point or from the environment>> fragment at the start of the loop, we can see how the ray from the previous iteration of the while loop can take care of the BSDF sample in Equation (13.10). The ray’s direction was chosen by sampling the BSDF, and so if it happens to hit a light source, then we have everything we need to evaluate the second term of the estimate other than the MIS weight w Subscript normal b Baseline left-parenthesis omega Subscript Baseline right-parenthesis . If the ray does not hit a light source, then that term is zero for the BSDF sample and there is no further work to do.

There are two cases to handle: infinite lights for rays that do not intersect any geometry, and surface emission for rays that do. In the first case, the ray path can terminate once lights have been considered.

<<Add emitted light at intersection point or from the environment>>= 
if (!si) { <<Incorporate emission from infinite lights for escaped ray>> 
for (const auto &light : infiniteLights) { SampledSpectrum Le = light.Le(ray, lambda); if (depth == 0 || specularBounce) L += beta * Le; else { <<Compute MIS weight for infinite light>> 
Float p_l = lightSampler.PMF(prevIntrCtx, light) * light.PDF_Li(prevIntrCtx, ray.d, true); Float w_b = PowerHeuristic(1, p_b, 1, p_l);
L += beta * w_b * Le; } }
break; } <<Incorporate emission from surface hit by ray>> 
SampledSpectrum Le = si->intr.Le(-ray.d, lambda); if (Le) { if (depth == 0 || specularBounce) L += beta * Le; else { <<Compute MIS weight for area light>> 
Light areaLight(si->intr.areaLight); Float lightPDF = lightSampler.PMF(prevIntrCtx, areaLight) * areaLight.PDF_Li(prevIntrCtx, ray.d, true); Float w_l = PowerHeuristic(1, bsdfPDF, 1, lightPDF);
L += beta * w_l * Le; } }

For the initial ray from the camera or after a perfect specular scattering event, emitted radiance should be included in the path without any MIS weighting, since light sampling was not performed at the previous vertex of the path. At this point in execution, beta already includes the BSDF, cosine factor, and PDF value from the previous scattering event, so multiplying beta by the emitted radiance gives the correct contribution.

<<Incorporate emission from infinite lights for escaped ray>>= 
for (const auto &light : infiniteLights) { SampledSpectrum Le = light.Le(ray, lambda); if (depth == 0 || specularBounce) L += beta * Le; else { <<Compute MIS weight for infinite light>> 
Float p_l = lightSampler.PMF(prevIntrCtx, light) * light.PDF_Li(prevIntrCtx, ray.d, true); Float w_b = PowerHeuristic(1, p_b, 1, p_l);
L += beta * w_b * Le; } }

Otherwise, it is necessary to compute the MIS weight w Subscript normal b . p_b gives us the BSDF’s PDF from the previous scattering event, so all we need is the PDF for the ray’s direction from sampling the light. This value is given by the product of the probability of sampling the light under consideration times the probability the light returns for sampling the direction.

Note that the PDF_Li() method is passed a true value for allowIncompletePDF here, again reflecting the fact that because BSDF sampling is capable of sampling all valid directions, it is not required that light sampling do so as well.

<<Compute MIS weight for infinite light>>= 
Float p_l = lightSampler.PMF(prevIntrCtx, light) * light.PDF_Li(prevIntrCtx, ray.d, true); Float w_b = PowerHeuristic(1, p_b, 1, p_l);

The code for the case of a ray hitting an emissive surface is in the fragment <<Incorporate emission from surface hit by ray>>. It is almost the same as the infinite light case, so we will not include it here.

The final issue is Russian roulette–based path termination. As outlined in Section 13.2.1, the task is easy: we compute a termination probability q however we like, make a random choice as to whether to terminate the path, and update beta if the path is not terminated so that all subsequent upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis terms will be scaled appropriately.

However, the details of how q is set can make a big difference. In general, it is a good idea for the termination probability to be based on the path throughput weight; in this way, if the BSDF’s value is small, it is more likely that the path will be terminated. Further, if the path is not terminated, then the scaling factor will generally cause beta to have a value around 1. Thus, all rays that are traced tend to make the same contribution to the image, which improves efficiency.

Another issue is that it is best if the beta value used to compute q does not include radiance scaling due to refraction. Consider a ray that passes through a glass object with a relative index of refraction of 1.5: when it enters the object, beta will pick up a factor of 1 slash 1.5 squared almost-equals 0.44 , but when it exits, that factor will cancel and beta will be back to 1. For ray paths that would exit, to have terminated them after the first refraction would be the wrong decision. Therefore, etaScale tracks those factors in beta so that they can be removed. The image in Figure 13.11 shows the increase in noise if this effect is not corrected for.

Figure 13.11: The Effect of Including Radiance Scaling Due to Transmission in the Russian Roulette Probability q . (a) If etaScale is not included in the probability, then some rays that would have passed through the glass object are terminated unnecessarily, leading to noise in the corresponding parts of the image. (b) Including etaScale in the computation of q fixes this issue. (Transparent Machines scene courtesy of Beeple.)

Finally, note that the termination probability is set according to the maximum component value of rrBeta rather than, for example, its average. Doing so gives better results when surface reflectances are highly saturated and some of the wavelength samples have much lower beta values than others, since it prevents any of the beta components from going above 1 due to Russian roulette.

<<Possibly terminate the path with Russian roulette>>= 
SampledSpectrum rrBeta = beta * etaScale; if (rrBeta.MaxComponentValue() < 1 && depth > 1) { Float q = std::max<Float>(0, 1 - rrBeta.MaxComponentValue()); if (sampler.Get1D() < q) break; beta /= 1 - q; }

Recall that Russian roulette only increases variance. Because it terminates some paths, this must be so, as the final image includes less information when it is applied. However, it can improve efficiency by allowing the renderer to focus its efforts on tracing rays that make the greatest contribution to the final image. Table 13.1 presents measurements of efficiency improvements from Russian roulette for a number of scenes.

Table 13.1: Monte Carlo Efficiency Benefits from Russian Roulette. Measurements of MSE and rendering time when using Russian roulette. All values reported are relative to rendering the same scene without Russian roulette. As expected, MSE increases to varying degrees due to ray termination, but the performance benefit more than makes up for it, leading to an increase in Monte Carlo efficiency.

Scene MSETimeEfficiency
Kroken (Figure 13.4) 1.31 0.261 2.92
Watercolor (Figure 13.6) 1.19 0.187 4.51
San Miguel (Figure 13.7) 1.00 0.239 4.17
BMW M6 (Figure 13.12) 1.00 0.801 1.25

13.4.1 Path Regularization

Scenes with concentrated indirect lighting can pose a challenge to the path-tracing algorithm: the problem is that if the incident indirect radiance at a point has substantial variation but BSDF sampling is being used to generate the direction of indirect rays, then the sampling distribution may be a poor match for the integrand. Variance spikes then occur when the ratio f left-parenthesis x right-parenthesis slash p left-parenthesis x right-parenthesis in the Monte Carlo estimator is large.

Figure 13.12: Image with High Variance Due to Difficult-to-Sample Indirect Lighting. The environment map illuminating the scene includes the sun, which is not only bright but also subtends a small solid angle. When an indirect lighting sample hits a specular surface and reflects to the sun’s direction, variance spikes in the image result because its contribution is not sampled well. (Car model courtesy of tyrant monkey, via Blend Swap.)

Figure 13.12 shows an example of this issue. The car is illuminated by a sky environment map where a bright sun occupies a small number of pixels. Consider sampling indirect lighting at a point on the ground near one of the wheels: the ground material is fairly diffuse, so any direction will be sampled with equal (cosine-weighted) probability. Rarely, a direction will be sampled that both hits the highly specular wheel and then also reflects to a direction where the sun is visible. This is the cause of the bright pixels on the ground. (The lighting in the car interior is similarly difficult to sample, since the glass prevents light source sampling; the variance spikes there follow.)

Figure 13.13: Scene from Figure 13.12 with Roughened BSDFs. (a) Increasing the roughness of all the BSDFs eliminates the variance spikes by allowing the use of MIS at all indirect ray intersection points, though this substantially changes the appearance of the scene. (Note that the car paint is duller and the window glass and headlight covers have the appearance of frosted glass.) (b) Roughening BSDFs only after the first non-specular scattering event along the path preserves visual detail while reducing the error from difficult light paths. (Car model courtesy of tyrant monkey, via Blend Swap.)

Informally, the idea behind path regularization is to blur the function being integrated in the case that it cannot be sampled effectively (or cannot be sampled in the first place). See Figure 13.13, which shows the same scene, but with all the BSDFs made more rough: perfect specular surfaces are glossy specular, and glossy specular surfaces are more diffuse. Although the overall characteristics of the image are quite different, the high variance on the ground has been eliminated: when an indirect lighting ray hits one of the wheels, it is now possible to use a lower variance MIS-based direct lighting calculation in place of following whichever direction is dictated by the law of specular reflection.

Blurring all the BSDFs in this way is an undesirable solution, but there is no need to do so for the camera rays or for rays that have only undergone perfect specular scattering: in those cases, we would like to leave the scene as it was specified. We can consider non-specular scattering itself to be a sort of blurring of the incident light, such that blurring the scene that is encountered after it occurs is less likely to be objectionable—thus the motivation to track this case via the anyNonSpecularBounces variable.

<<Possibly regularize the BSDF>>= 
if (regularize && anyNonSpecularBounces) bsdf.Regularize();

The BSDF class provides a Regularize() method that forwards the request on to its BxDF.

<<BSDF Public Methods>>+= 
void Regularize() { bxdf.Regularize(); }

The BxDF interface in turn requires the implementation of a Regularize() method. For BxDFs that are already fairly broad (e.g., the DiffuseBxDF), the corresponding method implementation is empty.

<<BxDF Interface>>+= 
void Regularize();

However, both the DielectricBxDF and ConductorBxDF can be nearly specular or perfect specular, depending on how smooth their microfacet distribution is. Therefore, their Regularize() method implementations do adjust their scattering properties, through a call to yet one more method named Regularize(), this one implemented by the TrowbridgeReitzDistribution.

<<DielectricBxDF Public Methods>>+= 

Unless the surface is already fairly rough, the TrowbridgeReitzDistribution’s Regularize() method doubles the alpha parameters and then clamps them—to ensure both that perfect specular surfaces with a roughness of zero become non-perfect specular and that surfaces are not excessively roughened.

<<TrowbridgeReitzDistribution Public Methods>>+= 
void Regularize() { if (alpha_x < 0.3f) alpha_x = Clamp(2 * alpha_x, 0.1f, 0.3f); if (alpha_y < 0.3f) alpha_y = Clamp(2 * alpha_y, 0.1f, 0.3f); }