14.3 Direct Lighting

Before we introduce the light transport equation in its full generality, we will implement the DirectLightingIntegrator which, unsurprisingly, only accounts for direct lighting—light that has traveled directly from a light source to the point being shaded—and ignores indirect illumination from objects that are not themselves emissive, except for basic specular reflection and transmission effects. Starting out with this integrator allows us to focus on some of the important details of direct lighting without worrying about the full light transport equation. Furthermore, some of the routines developed here will be used again in subsequent integrators that solve the complete light transport equation. Figure 14.12 shows the San Miguel scene rendered with direct lighting only.

<<DirectLightingIntegrator Declarations>>= 
class DirectLightingIntegrator : public SamplerIntegrator { public: <<DirectLightingIntegrator Public Methods>> 
DirectLightingIntegrator(LightStrategy strategy, int maxDepth, std::shared_ptr<const Camera> camera, std::shared_ptr<Sampler> sampler) : SamplerIntegrator(camera, sampler), strategy(strategy), maxDepth(maxDepth) { } Spectrum Li(const RayDifferential &ray, const Scene &scene, Sampler &sampler, MemoryArena &arena, int depth) const; void Preprocess(const Scene &scene, Sampler &sampler);
private: <<DirectLightingIntegrator Private Data>> 
const LightStrategy strategy; const int maxDepth; std::vector<int> nLightSamples;
};

Figure 14.12: Scene Rendered with Direct Lighting Only. Because only direct lighting is considered, some portions of the image are completely black because they are only lit by indirect illumination. (Model courtesy of Guillermo M. Leal Llaguno.)

The implementation provides two different strategies for computing direct lighting. Each computes an unbiased estimate of exitant radiance at a point in a given direction. The LightStrategy enumeration records which approach has been selected. The first strategy, UniformSampleAll, loops over all of the lights and takes a number of samples based on Light::nSamples from each of them, summing the result. (Recall the discussion of splitting in Section 13.7.1—the UniformSampleAll strategy is applying this technique.) The second, UniformSampleOne, takes a single sample from just one of the lights, chosen at random.

<<LightStrategy Declarations>>= 
enum class LightStrategy { UniformSampleAll, UniformSampleOne };

Depending on the scene being rendered, either of these approaches may be more appropriate. For example, if many image samples are being taken for each pixel (e.g., to resolve depth of field without excessive noise), then a single light sample per image sample may be more appropriate: in the aggregate all of the image samples in a pixel will sample the direct lighting well enough to give a high-quality image. Alternatively, if few samples per pixel are being taken, sampling all lights may be preferable to ensure a noise-free result.

The DirectLightingIntegrator constructor, which we won’t include here, just passes a Camera and Sampler to the SamplerIntegrator base class constructor and initializes two member variables. In addition to the direct lighting strategy, the DirectLightingIntegrator stores a maximum recursion depth for rays that are traced to account for specular reflection or specular transmission.

<<DirectLightingIntegrator Private Data>>= 
const LightStrategy strategy; const int maxDepth;

The numbers and types of samples needed by this integrator depend on the sampling strategy used: if a single sample is taken from a single light, then two two-dimensional samples obtained via Sampler::Get2D() suffice—one for selecting a position on a light source and one for sampling a scattered direction from the BSDF.

When multiple samples are taken per light, the integrator requests sample arrays from the sampler before the main rendering process begins. Using sample arrays is preferable to performing many separate calls to Sampler::Get2D() because it gives the sampler an opportunity to use improved sample placement techniques, optimizing the distribution of samples over the entire array and thus, over all of the light source samples for the current point being shaded.

<<DirectLightingIntegrator Method Definitions>>= 
void DirectLightingIntegrator::Preprocess(const Scene &scene, Sampler &sampler) { if (strategy == LightStrategy::UniformSampleAll) { <<Compute number of samples to use for each light>> 
for (const auto &light : scene.lights) nLightSamples.push_back(sampler.RoundCount(light->nSamples));
<<Request samples for sampling all lights>> 
for (int i = 0; i < maxDepth; ++i) { for (size_t j = 0; j < scene.lights.size(); ++j) { sampler.Request2DArray(nLightSamples[j]); sampler.Request2DArray(nLightSamples[j]); } }
} }

For the LightStrategy::UniformSampleAll strategy, each light stores a desired number of samples in its Light::nSamples member variable. However, the integrator here only uses that value as a starting point: the Sampler::RoundCount() method is given an opportunity to change that value to a more appropriate one based on its particular sample generation technique. (For example, many samplers only generate collections of samples with power-of-two sizes.) The final number of samples for each light is recorded in the nLightSamples member variable.

<<Compute number of samples to use for each light>>= 
for (const auto &light : scene.lights) nLightSamples.push_back(sampler.RoundCount(light->nSamples));

<<DirectLightingIntegrator Private Data>>+= 
std::vector<int> nLightSamples;

Now the sample arrays can be requested. There are two important details in this fragment: first, although a separate sample request is made for all of the possible ray depths up to the maximum, this doesn’t mean that all intersections will have sample arrays available. Rather, once a sample array is retrieved by a call to Sampler::Get2DArray(), that array won’t be returned again. If both specular reflection and transmission are present, then there may be up to 2 Superscript normal m normal a normal x normal upper D normal e normal p normal t normal h plus 1 Baseline minus 1 intersection points for each camera ray intersection. If the sample arrays are exhausted, the integrator switches to taking a single sample from each light.

Second, note that the arrays are requested in the order they will be consumed by the integrator: at each intersection point, two arrays are used for each light source, in the same order as the lights array.

<<Request samples for sampling all lights>>= 
for (int i = 0; i < maxDepth; ++i) { for (size_t j = 0; j < scene.lights.size(); ++j) { sampler.Request2DArray(nLightSamples[j]); sampler.Request2DArray(nLightSamples[j]); } }

As a SamplerIntegrator, the main method that the DirectLightingIntegrator must implement is Li(). The general form of its implementation here is similar to that of WhittedIntegrator::Li(): the BSDF at the intersection point is computed, emitted radiance is added if the surface is emissive, rays are traced recursively for specular reflection and transmission, and so on. We won’t include the full implementation of DirectLightingIntegrator::Li() here in order to focus on its key fragment, <<Compute direct lighting for DirectLightingIntegrator integrator>>, which estimates the value of the integral that gives the reflected radiance, accumulating it into a value L that will be returned from Li().

Two helper functions, which other integrators will also find useful, take care of the two sampling strategies.

<<Compute direct lighting for DirectLightingIntegrator integrator>>= 
if (strategy == LightStrategy::UniformSampleAll) L += UniformSampleAllLights(isect, scene, arena, sampler, nLightSamples); else L += UniformSampleOneLight(isect, scene, arena, sampler);

To understand the approaches implemented by the two strategies, first recall the scattering equation from Section 5.6, which says that exitant radiance upper L Subscript normal o Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline right-parenthesis from a point normal p Subscript on a surface in direction omega Subscript normal o due to incident radiance at the point is given by an integral of incoming radiance over the sphere times the BSDF for each direction and a cosine term. For the DirectLightingIntegrator, we are only interested in incident radiance directly from light sources, which we will denote by upper L Subscript normal d Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis :

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

This can be broken into a sum over the n lights in the scene

sigma-summation Underscript j equals 1 Overscript n Endscripts integral Underscript script upper S squared Endscripts f Subscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline comma omega Subscript normal i Baseline right-parenthesis upper L Subscript normal d left-parenthesis j right-parenthesis Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal i Baseline right-parenthesis StartAbsoluteValue cosine theta Subscript normal i Baseline EndAbsoluteValue normal d omega Subscript normal i Baseline comma

where upper L Subscript normal d left-parenthesis j right-parenthesis Superscript denotes incident radiance from the j th light and

upper L Subscript normal d Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal i Baseline right-parenthesis equals sigma-summation Underscript j Endscripts upper L Subscript normal d left-parenthesis j right-parenthesis Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal i Baseline right-parenthesis period

One valid approach is to estimate each term of the sum in Equation (14.12) individually, adding the results together. This is the most basic direct lighting strategy and is implemented in UniformSampleAllLights().

In addition to information about the intersection point and additional parameters that are needed to compute the direct lighting, UniformSampleAllLights() also takes a handleMedia parameter that indicates whether the effects of volumetric attenuation should be considered in the direct lighting computation. (Between this parameter and the detail that it takes an Interaction rather than a SurfaceInteraction, this function is actually able to compute reflected radiance at points in participating media; fragments related to this functionality will be defined in the next chapter.)

<<Integrator Utility Functions>>= 
Spectrum UniformSampleAllLights(const Interaction &it, const Scene &scene, MemoryArena &arena, Sampler &sampler, const std::vector<int> &nLightSamples, bool handleMedia) { Spectrum L(0.f); for (size_t j = 0; j < scene.lights.size(); ++j) { <<Accumulate contribution of jth light to L>> 
const std::shared_ptr<Light> &light = scene.lights[j]; int nSamples = nLightSamples[j]; const Point2f *uLightArray = sampler.Get2DArray(nSamples); const Point2f *uScatteringArray = sampler.Get2DArray(nSamples); if (!uLightArray || !uScatteringArray) { <<Use a single sample for illumination from light>> 
Point2f uLight = sampler.Get2D(); Point2f uScattering = sampler.Get2D(); L += EstimateDirect(it, uScattering, *light, uLight, scene, sampler, arena, handleMedia);
} else { <<Estimate direct lighting using sample arrays>> 
Spectrum Ld(0.f); for (int k = 0; k < nSamples; ++k) Ld += EstimateDirect(it, uScatteringArray[k], *light, uLightArray[k], scene, sampler, arena, handleMedia); L += Ld / nSamples;
}
} return L; }

This function attempts to retrieve the sample arrays from the Sampler that were previously requested in the Preprocess() method.

<<Accumulate contribution of jth light to L>>= 
const std::shared_ptr<Light> &light = scene.lights[j]; int nSamples = nLightSamples[j]; const Point2f *uLightArray = sampler.Get2DArray(nSamples); const Point2f *uScatteringArray = sampler.Get2DArray(nSamples); if (!uLightArray || !uScatteringArray) { <<Use a single sample for illumination from light>> 
Point2f uLight = sampler.Get2D(); Point2f uScattering = sampler.Get2D(); L += EstimateDirect(it, uScattering, *light, uLight, scene, sampler, arena, handleMedia);
} else { <<Estimate direct lighting using sample arrays>> 
Spectrum Ld(0.f); for (int k = 0; k < nSamples; ++k) Ld += EstimateDirect(it, uScatteringArray[k], *light, uLightArray[k], scene, sampler, arena, handleMedia); L += Ld / nSamples;
}

If all of the requested arrays have been consumed, the code falls back to a single sample estimate via calls to Sampler::Get2D().

<<Use a single sample for illumination from light>>= 
Point2f uLight = sampler.Get2D(); Point2f uScattering = sampler.Get2D(); L += EstimateDirect(it, uScattering, *light, uLight, scene, sampler, arena, handleMedia);

For each light sample, the EstimateDirect() function (which will be defined shortly) computes the value of the Monte Carlo estimator for its contribution. When sample arrays were successfully obtained, all that remains to be done is to average the estimates from each of their sample values.

<<Estimate direct lighting using sample arrays>>= 
Spectrum Ld(0.f); for (int k = 0; k < nSamples; ++k) Ld += EstimateDirect(it, uScatteringArray[k], *light, uLightArray[k], scene, sampler, arena, handleMedia); L += Ld / nSamples;

In a scene with a large number of lights, it may not be desirable to always compute direct lighting from all of the lights at every point that is shaded. The Monte Carlo approach gives a way to do this that still computes the correct result on average. Consider as an example computing the expected value of the sum of two functions upper E left-bracket f left-parenthesis x right-parenthesis plus g left-parenthesis x right-parenthesis right-bracket . If we randomly evaluate just one of f left-parenthesis x right-parenthesis or g left-parenthesis x right-parenthesis and multiply the result by 2, then the expected value of the result will still be f left-parenthesis x right-parenthesis plus g left-parenthesis x right-parenthesis . This idea also generalizes to sums with an arbitrary number of terms; see Ross (2002, p. 102) for a proof. Here we estimate direct lighting for only one randomly chosen light and multiply the result by the number of lights to compensate.

<<Integrator Utility Functions>>+=  
Spectrum UniformSampleOneLight(const Interaction &it, const Scene &scene, MemoryArena &arena, Sampler &sampler, bool handleMedia) { <<Randomly choose a single light to sample, light>> 
int nLights = int(scene.lights.size()); if (nLights == 0) return Spectrum(0.f); int lightNum = std::min((int)(sampler.Get1D() * nLights), nLights - 1); const std::shared_ptr<Light> &light = scene.lights[lightNum];
Point2f uLight = sampler.Get2D(); Point2f uScattering = sampler.Get2D(); return (Float)nLights * EstimateDirect(it, uScattering, *light, uLight, scene, sampler, arena, handleMedia); }

Which of the nLights to sample illumination from is determined using a 1D sample from Sampler::Get1D().

<<Randomly choose a single light to sample, light>>= 
int nLights = int(scene.lights.size()); if (nLights == 0) return Spectrum(0.f); int lightNum = std::min((int)(sampler.Get1D() * nLights), nLights - 1); const std::shared_ptr<Light> &light = scene.lights[lightNum];

It’s possible to be even more creative in choosing the individual light sampling probabilities than the uniform method used in UniformSampleOneLight(). In fact, we’re free to set the probabilities any way we like, as long as we weight the result appropriately and there is a nonzero probability of sampling any light that contributes to the reflection at the point. The better a job we do at setting these probabilities to reflect the relative contributions of lights to reflected radiance at the reference point, the more efficient the Monte Carlo estimator will be, and the fewer rays will be needed to lower variance to an acceptable level. (This is just the discrete instance of importance sampling.)

One widely used approach to this task is to base the sample distribution on the total power of each light. In a similar manner, we could take more than one light sample with this approach; indeed, any number of samples can be taken in the end, as long as they are weighted appropriately.

14.3.1 Estimating the Direct Lighting Integral

Having chosen a particular light to estimate direct lighting from, we need to estimate the value of the integral

integral Underscript script upper S squared Endscripts f Subscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline comma omega Subscript normal i Baseline right-parenthesis upper L Subscript normal d Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal i Baseline right-parenthesis StartAbsoluteValue cosine theta Subscript normal i Baseline EndAbsoluteValue normal d omega Subscript normal i

for that light. To compute this estimate, we need to choose one or more directions omega Subscript normal r and apply the Monte Carlo estimator:

StartFraction 1 Over upper N EndFraction sigma-summation Underscript j equals 1 Overscript upper N Endscripts StartFraction f Subscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline comma omega Subscript normal r Baseline right-parenthesis upper L Subscript normal d Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal r Baseline right-parenthesis StartAbsoluteValue cosine theta Subscript j Baseline EndAbsoluteValue Over p left-parenthesis omega Subscript normal r Baseline right-parenthesis EndFraction period

To reduce variance, we will use importance sampling to choose the directions omega Subscript normal r . Because both the BSDF and the direct radiance terms are individually complex, it can be difficult to find sampling distributions that match their product well. (However, see the “Further Reading” section as well as Exercise 14.8 at the end of this chapter for approaches that sample their product directly.) Here we will use the BSDF’s sampling distribution for some of the samples and the light’s for the rest. Depending on the characteristics of each of them, one of these two sampling methods may be far more effective than the other. Therefore, we will use multiple importance sampling to reduce variance for the cases where one or the other is more effective.

Figure 14.13 shows cases where one of the 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 increasing size. Figures 14.13(1) and (2) show the BSDF and light sampling strategies on their own, while Figure 14.13(3) shows their combination computed using multiple importance sampling. As the example illustrates, sampling the BSDF can be 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 14.13: 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.

By applying multiple importance sampling, not only can we use both of the two sampling methods, but we can also do so in a way that eliminates the extreme variance from the cases where a sampling method unexpectedly finds a high-contribution direction, since the weighting terms from MIS reduce these contributions substantially.

EstimateDirect() implements this approach and computes a direct lighting estimate for a single light source sample. Its handleMedia parameter indicates whether the effect of attenuation from participating media should be accounted for, and its specular parameter indicates whether or not perfectly specular lobes should be considered in the direct illumination estimate. The default value for both specular and handleMedia arguments is set to false in the function declaration, which is not shown here.

<<Integrator Utility Functions>>+=  
Spectrum EstimateDirect(const Interaction &it, const Point2f &uScattering, const Light &light, const Point2f &uLight, const Scene &scene, Sampler &sampler, MemoryArena &arena, bool handleMedia, bool specular) { BxDFType bsdfFlags = specular ? BSDF_ALL : BxDFType(BSDF_ALL & ~BSDF_SPECULAR); Spectrum Ld(0.f); <<Sample light source with multiple importance sampling>> 
Vector3f wi; Float lightPdf = 0, scatteringPdf = 0; VisibilityTester visibility; Spectrum Li = light.Sample_Li(it, uLight, &wi, &lightPdf, &visibility); if (lightPdf > 0 && !Li.IsBlack()) { <<Compute BSDF or phase function’s value for light sample>> 
Spectrum f; if (it.IsSurfaceInteraction()) { <<Evaluate BSDF for light sampling strategy>> 
const SurfaceInteraction &isect = (const SurfaceInteraction &)it; f = isect.bsdf->f(isect.wo, wi, bsdfFlags) * AbsDot(wi, isect.shading.n); scatteringPdf = isect.bsdf->Pdf(isect.wo, wi, bsdfFlags);
} else { <<Evaluate phase function for light sampling strategy>> 
const MediumInteraction &mi = (const MediumInteraction &)it; Float p = mi.phase->p(mi.wo, wi); f = Spectrum(p); scatteringPdf = p;
}
if (!f.IsBlack()) { <<Compute effect of visibility for light source sample>> 
if (handleMedia) Li *= visibility.Tr(scene, sampler); else if (!visibility.Unoccluded(scene)) Li = Spectrum(0.f);
<<Add light’s contribution to reflected radiance>> 
if (!Li.IsBlack()) { if (IsDeltaLight(light.flags)) Ld += f * Li / lightPdf; else { Float weight = PowerHeuristic(1, lightPdf, 1, scatteringPdf); Ld += f * Li * weight / lightPdf; } }
} }
<<Sample BSDF with multiple importance sampling>> 
if (!IsDeltaLight(light.flags)) { Spectrum f; bool sampledSpecular = false; if (it.IsSurfaceInteraction()) { <<Sample scattered direction for surface interactions>> 
BxDFType sampledType; const SurfaceInteraction &isect = (const SurfaceInteraction &)it; f = isect.bsdf->Sample_f(isect.wo, &wi, uScattering, &scatteringPdf, bsdfFlags, &sampledType); f *= AbsDot(wi, isect.shading.n); sampledSpecular = sampledType & BSDF_SPECULAR;
} else { <<Sample scattered direction for medium interactions>> 
const MediumInteraction &mi = (const MediumInteraction &)it; Float p = mi.phase->Sample_p(mi.wo, &wi, uScattering); f = Spectrum(p); scatteringPdf = p;
} if (!f.IsBlack() && scatteringPdf > 0) { <<Account for light contributions along sampled direction wi>> 
Float weight = 1; if (!sampledSpecular) { lightPdf = light.Pdf_Li(it, wi); if (lightPdf == 0) return Ld; weight = PowerHeuristic(1, scatteringPdf, 1, lightPdf); } <<Find intersection and compute transmittance>> 
SurfaceInteraction lightIsect; Ray ray = it.SpawnRay(wi); Spectrum Tr(1.f); bool foundSurfaceInteraction = handleMedia ? scene.IntersectTr(ray, sampler, &lightIsect, &Tr) : scene.Intersect(ray, &lightIsect);
<<Add light contribution from material sampling>> 
Spectrum Li(0.f); if (foundSurfaceInteraction) { if (lightIsect.primitive->GetAreaLight() == &light) Li = lightIsect.Le(-wi); } else Li = light.Le(ray); if (!Li.IsBlack()) Ld += f * Li * Tr * weight / scatteringPdf;
} }
return Ld; }

First, one sample is taken from the light’s sampling distribution using Sample_Li(), which also returns the light’s emitted radiance and the value of the PDF for the sampled direction.

<<Sample light source with multiple importance sampling>>= 
Vector3f wi; Float lightPdf = 0, scatteringPdf = 0; VisibilityTester visibility; Spectrum Li = light.Sample_Li(it, uLight, &wi, &lightPdf, &visibility); if (lightPdf > 0 && !Li.IsBlack()) { <<Compute BSDF or phase function’s value for light sample>> 
Spectrum f; if (it.IsSurfaceInteraction()) { <<Evaluate BSDF for light sampling strategy>> 
const SurfaceInteraction &isect = (const SurfaceInteraction &)it; f = isect.bsdf->f(isect.wo, wi, bsdfFlags) * AbsDot(wi, isect.shading.n); scatteringPdf = isect.bsdf->Pdf(isect.wo, wi, bsdfFlags);
} else { <<Evaluate phase function for light sampling strategy>> 
const MediumInteraction &mi = (const MediumInteraction &)it; Float p = mi.phase->p(mi.wo, wi); f = Spectrum(p); scatteringPdf = p;
}
if (!f.IsBlack()) { <<Compute effect of visibility for light source sample>> 
if (handleMedia) Li *= visibility.Tr(scene, sampler); else if (!visibility.Unoccluded(scene)) Li = Spectrum(0.f);
<<Add light’s contribution to reflected radiance>> 
if (!Li.IsBlack()) { if (IsDeltaLight(light.flags)) Ld += f * Li / lightPdf; else { Float weight = PowerHeuristic(1, lightPdf, 1, scatteringPdf); Ld += f * Li * weight / lightPdf; } }
} }

Only if the light successfully samples a direction and returns nonzero emitted radiance does EstimateDirect() go ahead and evaluate the BSDF or phase function at the provided Interaction; otherwise, there’s no reason to go through the computational expense. (Consider, for example, a spotlight, which returns no radiance for points outside its illumination cone.)

<<Compute BSDF or phase function’s value for light sample>>= 
Spectrum f; if (it.IsSurfaceInteraction()) { <<Evaluate BSDF for light sampling strategy>> 
const SurfaceInteraction &isect = (const SurfaceInteraction &)it; f = isect.bsdf->f(isect.wo, wi, bsdfFlags) * AbsDot(wi, isect.shading.n); scatteringPdf = isect.bsdf->Pdf(isect.wo, wi, bsdfFlags);
} else { <<Evaluate phase function for light sampling strategy>> 
const MediumInteraction &mi = (const MediumInteraction &)it; Float p = mi.phase->p(mi.wo, wi); f = Spectrum(p); scatteringPdf = p;
}

<<Evaluate BSDF for light sampling strategy>>= 
const SurfaceInteraction &isect = (const SurfaceInteraction &)it; f = isect.bsdf->f(isect.wo, wi, bsdfFlags) * AbsDot(wi, isect.shading.n); scatteringPdf = isect.bsdf->Pdf(isect.wo, wi, bsdfFlags);

(We postpone the medium-specific fragment that evaluates the phase function at the interaction point, <<Evaluate phase function for light sampling strategy>>, to Chapter 15.)

If participating media are to be accounted for, the radiance from the light to the illuminated point is scaled by the beam transmittance between the two points to account for attenuation due to participating media. Otherwise, a call to the VisibilityTester’s Unoccluded() method traces a shadow ray to determine if the sampled point on the light source is visible. (This step and the next are skipped if the BSDF or phase function returned a black SPD.)

<<Compute effect of visibility for light source sample>>= 
if (handleMedia) Li *= visibility.Tr(scene, sampler); else if (!visibility.Unoccluded(scene)) Li = Spectrum(0.f);

The light sample’s contribution can now be accumulated. Recall from Section 14.2.1 that if the light is described by a delta distribution then there is an implied delta distribution in both the emitted radiance value returned from Sample_Li() as well as the PDF and that they are expected to cancel out when the estimator is evaluated. In this case, we must not try to apply multiple importance sampling and should compute the standard estimator instead. If this isn’t a delta distribution light source, then the BSDF’s PDF value for sampling the direction omega Subscript normal i , which was returned by BSDF::Pdf(), is used with the MIS estimator, where the weight is computed here with the power heuristic.

<<Add light’s contribution to reflected radiance>>= 
if (!Li.IsBlack()) { if (IsDeltaLight(light.flags)) Ld += f * Li / lightPdf; else { Float weight = PowerHeuristic(1, lightPdf, 1, scatteringPdf); Ld += f * Li * weight / lightPdf; } }

Next, a sample is generated using the BSDF’s sampling distribution. This step should be skipped if the light source’s emission profile involves a delta distribution because, in that case, there’s no chance that sampling the BSDF will give a direction that receives light from the source. Otherwise, the BSDF can be sampled.

<<Sample BSDF with multiple importance sampling>>= 
if (!IsDeltaLight(light.flags)) { Spectrum f; bool sampledSpecular = false; if (it.IsSurfaceInteraction()) { <<Sample scattered direction for surface interactions>> 
BxDFType sampledType; const SurfaceInteraction &isect = (const SurfaceInteraction &)it; f = isect.bsdf->Sample_f(isect.wo, &wi, uScattering, &scatteringPdf, bsdfFlags, &sampledType); f *= AbsDot(wi, isect.shading.n); sampledSpecular = sampledType & BSDF_SPECULAR;
} else { <<Sample scattered direction for medium interactions>> 
const MediumInteraction &mi = (const MediumInteraction &)it; Float p = mi.phase->Sample_p(mi.wo, &wi, uScattering); f = Spectrum(p); scatteringPdf = p;
} if (!f.IsBlack() && scatteringPdf > 0) { <<Account for light contributions along sampled direction wi>> 
Float weight = 1; if (!sampledSpecular) { lightPdf = light.Pdf_Li(it, wi); if (lightPdf == 0) return Ld; weight = PowerHeuristic(1, scatteringPdf, 1, lightPdf); } <<Find intersection and compute transmittance>> 
SurfaceInteraction lightIsect; Ray ray = it.SpawnRay(wi); Spectrum Tr(1.f); bool foundSurfaceInteraction = handleMedia ? scene.IntersectTr(ray, sampler, &lightIsect, &Tr) : scene.Intersect(ray, &lightIsect);
<<Add light contribution from material sampling>> 
Spectrum Li(0.f); if (foundSurfaceInteraction) { if (lightIsect.primitive->GetAreaLight() == &light) Li = lightIsect.Le(-wi); } else Li = light.Le(ray); if (!Li.IsBlack()) Ld += f * Li * Tr * weight / scatteringPdf;
} }

Once more, we postpone the medium-related code to Chapter 15. Given a surface interaction, the implementation samples a scattered direction and records whether or not a delta distribution was sampled.

<<Sample scattered direction for surface interactions>>= 
BxDFType sampledType; const SurfaceInteraction &isect = (const SurfaceInteraction &)it; f = isect.bsdf->Sample_f(isect.wo, &wi, uScattering, &scatteringPdf, bsdfFlags, &sampledType); f *= AbsDot(wi, isect.shading.n); sampledSpecular = sampledType & BSDF_SPECULAR;

One important detail is that the light’s PDF and the multiple importance sampling weight are only computed if the BSDF component used for sampling omega Subscript normal i is non-specular; in the specular case, MIS shouldn’t be applied since there is no chance of the light sampling the specular direction.

<<Account for light contributions along sampled direction wi>>= 
Float weight = 1; if (!sampledSpecular) { lightPdf = light.Pdf_Li(it, wi); if (lightPdf == 0) return Ld; weight = PowerHeuristic(1, scatteringPdf, 1, lightPdf); } <<Find intersection and compute transmittance>> 
SurfaceInteraction lightIsect; Ray ray = it.SpawnRay(wi); Spectrum Tr(1.f); bool foundSurfaceInteraction = handleMedia ? scene.IntersectTr(ray, sampler, &lightIsect, &Tr) : scene.Intersect(ray, &lightIsect);
<<Add light contribution from material sampling>> 
Spectrum Li(0.f); if (foundSurfaceInteraction) { if (lightIsect.primitive->GetAreaLight() == &light) Li = lightIsect.Le(-wi); } else Li = light.Le(ray); if (!Li.IsBlack()) Ld += f * Li * Tr * weight / scatteringPdf;

Given a direction sampled by the BSDF or a medium’s phase function, we need to find out if the ray along that direction intersects this particular light source and if so, how much radiance from the light reaches the surface. When participating media are being accounted for, the transmittance up to the intersection point on the light is recorded.

<<Find intersection and compute transmittance>>= 
SurfaceInteraction lightIsect; Ray ray = it.SpawnRay(wi); Spectrum Tr(1.f); bool foundSurfaceInteraction = handleMedia ? scene.IntersectTr(ray, sampler, &lightIsect, &Tr) : scene.Intersect(ray, &lightIsect);

The code must account for both regular area lights, with geometry associated with them, as well as lights like the InfiniteAreaLight that don’t have geometry but need to return their radiance for the sample ray via the Light::Le() method.

<<Add light contribution from material sampling>>= 
Spectrum Li(0.f); if (foundSurfaceInteraction) { if (lightIsect.primitive->GetAreaLight() == &light) Li = lightIsect.Le(-wi); } else Li = light.Le(ray); if (!Li.IsBlack()) Ld += f * Li * Tr * weight / scatteringPdf;