16.3 Bidirectional Path Tracing

The path-tracing algorithm described in Section 14.5 was the first fully general light transport algorithm in computer graphics, handling both a wide variety of geometric representations, lights, and BSDF models. Although it works well for many scenes, path tracing can exhibit high variance in the presence of particular tricky lighting conditions. For example, consider the setting shown in Figure 16.12: a light source is illuminating a small area on the ceiling such that the rest of the room is only illuminated by indirect lighting bouncing from that area. If we only trace paths starting from the camera, we will almost never happen to sample a path vertex in the illuminated region on the ceiling before we trace a shadow ray to the light. Most of the paths will have no contribution, while a few of them—the ones that happen to hit the small region on the ceiling—will have a large contribution. The resulting image will have high variance.

Figure 16.12: A Difficult Case for Path Tracing Starting from the Camera. A light source is illuminating a small area on the ceiling (thick line) such that only paths with a second-to-last vertex in the area indicated will be able to find illumination from the light. Bidirectional methods, where a path is started from the light and is connected with a path from the camera, can handle situations like these more robustly.

Difficult lighting settings like this can be handled more effectively by constructing paths that start from the camera on one end and from the light on the other end and are connected in the middle with a visibility ray. The resulting bidirectional path-tracing algorithm (henceforth referred to as BDPT) is a generalization of the standard path-tracing algorithm that can be much more efficient. In contrast to stochastic progressive photon mapping, BDPT is unbiased and does not blur the scene illumination.

BDPT first incrementally constructs a camera subpath starting with a point normal p 0 on the camera. The next vertex, normal p 1 , is found by computing the first intersection along the camera ray. Another vertex is found by sampling the BSDF at normal p 1 and tracing a ray to find a point normal p 2 , and so forth. The resulting path of t vertices is normal p 0 comma normal p 1 comma ellipsis comma normal p Subscript normal t minus 1 Baseline . Following the same process starting from a point on a light source normal q 0 (and then using adjoint BSDFs at each vertex) creates a light subpath of s vertices, normal q 0 comma normal q 1 comma ellipsis comma normal q Subscript normal s minus 1 Baseline .

Given the two subpaths, a complete light-carrying path can be found by connecting a pair of vertices from each path.

normal p Subscript Baseline overbar Subscript Baseline equals normal q 0 comma ellipsis comma normal q Subscript s prime minus 1 Baseline comma normal p Subscript t prime minus 1 Baseline comma ellipsis comma normal p 0 comma

where s prime less-than-or-equal-to s and t prime less-than-or-equal-to t . (Our notation orders the vertices in normal p Subscript Baseline overbar Subscript according to the propagation of light). If a visibility ray between normal q Subscript s prime minus 1 and normal p Subscript t prime minus 1 is unoccluded, then the path contribution can be found by evaluating the BSDFs at the connecting vertices (see Figure 16.13). More generally, these subpaths can be combined using the theory of path-space integration from Section 14.4.4.

Figure 16.13: Bidirectional path tracing is based on generating two subpaths, one starting from a light and the other starting from the camera. Light-carrying paths can be found by attempting to connect pairs of vertices, one from each path. If a visibility ray between the two (dashed line) is unoccluded, the path’s contribution can be added to the radiance estimate.

Superficially, this approach bears some semblance to the two phases of the photon mapper; a key difference is that BDPT computes an unbiased estimate without density estimation. There are also significant differences in how the fact that a given light path could have been constructed in multiple different ways is handled.

There are three refinements to the algorithm that we’ve described so far that improve its performance in practice. The first two are analogous to improvements made to path tracing, and the third is a powerful variance reduction technique.

  • First, subpaths can be reused: given a path normal q 0 comma ellipsis comma normal q Subscript s minus 1 Baseline comma normal p Subscript t minus 1 Baseline comma ellipsis comma normal p 0 , transport can be evaluated over all of the paths given by connecting all the various combinations of prefixes of the two paths together. If two paths have s and t vertices, respectively, then a variety of unique paths can be constructed from them, ranging in length from 2 to s plus t vertices long. Figure 16.14 illustrates these strategies for the case of direct illumination.
  • The second optimization is not to try to connect paths that use only a single vertex from one of the subpaths. It is preferable to generate those paths using optimized sampling routines provided by the camera and light sources; for light sources, these are the direct lighting techniques that were introduced in Section 14.3.
  • The third optimization weights the various strategies for generating paths of a given length more carefully than just averaging all of the strategies that construct paths of the same length. BDPT’s approach of connecting subpaths means that a path containing n scattering events can be generated in n plus 3 different ways. We can expect that some strategies will be a good choice for producing certain types of paths while being quite poor for others. Multiple importance sampling can be applied to combine the set of connection strategies into a single estimator that uses each strategy where it is best. This application of MIS is crucial to BDPT’s efficiency.

Figure 16.14: The four different ways in which bidirectional path tracing can create a direct illumination path. (a) Standard path tracing without direct illumination sampling, where a ray generated using BSDF sampling from a point on a surface happens to intersect a light source. (b) Path tracing with direct illumination sampling, where the explicit shadow ray is indicated with a dashed line. (c) Particle tracing from a light source with an explicit visibility test between a point on a surface and the camera. (d) Particle tracing where a particle happens to intersect the camera’s lens.

One of BDPT’s connection strategies is to directly connect light subpath vertices to the camera: these paths almost always create a contribution whose raster coordinates differ from the current pixel being rendered, which violates the expectations of the SamplerIntegrator interface. Thus, BDPTIntegrator derives from the more general Integrator interface so that it can have more flexibility in how it updates the image. Its implementation is in the files integrators/bdpt.h and integrators/bdpt.cpp.

<<BDPT Declarations>>= 
class BDPTIntegrator : public Integrator { public: <<BDPTIntegrator Public Methods>> 
BDPTIntegrator(std::shared_ptr<Sampler> sampler, std::shared_ptr<const Camera> camera, int maxDepth, bool visualizeStrategies, bool visualizeWeights) : sampler(sampler), camera(camera), maxDepth(maxDepth), visualizeStrategies(visualizeStrategies), visualizeWeights(visualizeWeights) { } void Render(const Scene &scene);
private: <<BDPTIntegrator Private Data>> 
std::shared_ptr<Sampler> sampler; std::shared_ptr<const Camera> camera; const int maxDepth;
};

The BDPTIntegrator constructor, which is straightforward and not included here, initializes member variables with the provided camera, sampler, and the maximum path depth.

<<BDPTIntegrator Private Data>>= 
std::shared_ptr<Sampler> sampler; std::shared_ptr<const Camera> camera; const int maxDepth;

All subpath creation and connection steps are performed in a nested parallel loop over pixels in BDPTIntegrator::Render(). The overall structure of this method is very similar to SamplerIntegrator::Render():

We won’t include this code here, as the details should be familiar now. Instead, we’ll move forward to the fragment responsible for generating and connecting the subpaths for a pixel sample.

<<BDPTIntegrator Public Methods>>= 
void Render(const Scene &scene);

Generating a single BDPT sample involves sampling a pixel position in the image, generating camera and light subpaths, and then connecting them using a variety of specialized connection strategies.

<<Generate a single sample using BDPT>>= 
Point2f pFilm = (Point2f)pPixel + tileSampler->Get2D(); <<Trace the camera and light subpaths>> 
Vertex *cameraVertices = arena.Alloc<Vertex>(maxDepth + 2); Vertex *lightVertices = arena.Alloc<Vertex>(maxDepth + 1); int nCamera = GenerateCameraSubpath(scene, *tileSampler, arena, maxDepth + 2, *camera, pFilm, cameraVertices); int nLight = GenerateLightSubpath(scene, *tileSampler, arena, maxDepth + 1, cameraVertices[0].time(), *lightDistr, lightVertices);
<<Execute all BDPT connection strategies>> 
Spectrum L(0.f); for (int t = 1; t <= nCamera; ++t) { for (int s = 0; s <= nLight; ++s) { int depth = t + s - 2; if ((s == 1 && t == 1) || depth < 0 || depth > maxDepth) continue; <<Execute the left-parenthesis s comma t right-parenthesis connection strategy and update L>> 
Point2f pFilmNew = pFilm; Float misWeight = 0.f; Spectrum Lpath = ConnectBDPT(scene, lightVertices, cameraVertices, s, t, *lightDistr, *camera, *tileSampler, &pFilmNew, &misWeight); if (t != 1) L += Lpath; else film->AddSplat(pFilmNew, Lpath);
} } filmTile->AddSample(pFilm, L);

The Vertex class, which will be introduced in Section 16.3.1, represents a vertex along a subpath. We start by allocating two arrays for vertices for the two subpaths. In addition to a vertex on a surface, a Vertex can represent a vertex for a scattering event in participating media, a vertex on a light source, or a vertex on the camera lens.

For each subpath, one more vertex than the maximum path length must be allocated to store the starting vertex on the light or camera. Camera subpaths get yet again one more vertex, which allows camera paths to randomly intersect light sources—this strategy is important for rendering area lights seen by reflections only from specular surfaces, for example. (The corresponding strategy of allowing a light subpath to randomly intersect the camera lens is less useful in practice.)

The GenerateCameraSubpath() and GenerateLightSubpath() functions, which generate these two subpaths, will be defined in Section 16.3.2, after some preliminaries related to the Vertex representation.

<<Trace the camera and light subpaths>>= 
Vertex *cameraVertices = arena.Alloc<Vertex>(maxDepth + 2); Vertex *lightVertices = arena.Alloc<Vertex>(maxDepth + 1); int nCamera = GenerateCameraSubpath(scene, *tileSampler, arena, maxDepth + 2, *camera, pFilm, cameraVertices); int nLight = GenerateLightSubpath(scene, *tileSampler, arena, maxDepth + 1, cameraVertices[0].time(), *lightDistr, lightVertices);

After the subpaths have been generated, a nested for loop iterates over all pairs of vertices from the two subpaths and attempts to connect them. In these loops, s and t correspond to the number of vertices to use from the corresponding subpath; an index of 0 means that no scattering events are used from the corresponding subpath. In our implementation, this strategy is only supported for the s equals 0 case, which furthermore requires cameraVertices[t] to be a surface intersection involving a light source. Because the dual case—intersecting a camera with t equals 0 —is not supported, the loop over camera subpaths starts at t equals 1 .

A path length of 1 corresponds to connecting a point on the camera lens or a light source to the other subpath. For light endpoints, this is identical to the standard light sampling approach provided by Light::Sample_Li() and first used in Section 14.3.1; our implementation uses this existing functionality. For camera endpoints, we will rely on the symmetric analog Camera::Sample_Wi(). Since Camera::Sample_Wi() and Light::Sample_Li() cannot both be used at the same time, we skip the s equals t equals 1 case.

<<Execute all BDPT connection strategies>>= 
Spectrum L(0.f); for (int t = 1; t <= nCamera; ++t) { for (int s = 0; s <= nLight; ++s) { int depth = t + s - 2; if ((s == 1 && t == 1) || depth < 0 || depth > maxDepth) continue; <<Execute the left-parenthesis s comma t right-parenthesis connection strategy and update L>> 
Point2f pFilmNew = pFilm; Float misWeight = 0.f; Spectrum Lpath = ConnectBDPT(scene, lightVertices, cameraVertices, s, t, *lightDistr, *camera, *tileSampler, &pFilmNew, &misWeight); if (t != 1) L += Lpath; else film->AddSplat(pFilmNew, Lpath);
} } filmTile->AddSample(pFilm, L);

The ConnectBDPT() function attempts to connect the two subpaths with the given number of vertices; it returns the weighted contribution of the radiance carried along the resulting path. (It will be defined shortly, in Section 16.3.3.) In most cases, this contribution is accumulated into the variable L that will be provided to the FilmTile after all of the subpath connections have been attempted. However, the t equals 1 connection connects a vertex of the light subpath directly to the camera and thus will produce different raster positions in every iteration—in this case, the implementation calls Film::AddSplat() to immediately record its sample contribution.

<<Execute the left-parenthesis s comma t right-parenthesis connection strategy and update L>>= 
Point2f pFilmNew = pFilm; Float misWeight = 0.f; Spectrum Lpath = ConnectBDPT(scene, lightVertices, cameraVertices, s, t, *lightDistr, *camera, *tileSampler, &pFilmNew, &misWeight); if (t != 1) L += Lpath; else film->AddSplat(pFilmNew, Lpath);

16.3.1 Vertex Abstraction Layer

A general advantage of path-space rendering techniques is their ability to create paths in a large number of different ways, but this characteristic often leads to cluttered and hard-to-debug implementations. Establishing connections between pairs of vertices on the light and camera subpaths is a simple operation when only surface interactions are involved but quickly becomes unwieldy if one or both of the vertices may represent a scattering event in participating media, for example.

Instead of an inconveniently large number of conditional statements in the core BDPT code, we’ll define the Vertex type, which can represent any kind of path vertex. All of the necessary conditional logic to handle various cases that occur throughout the BDPT implementation will be encapsulated in its methods.

<<BDPT Declarations>>+= 
struct Vertex { <<Vertex Public Data>> 
VertexType type; Spectrum beta; union { EndpointInteraction ei; MediumInteraction mi; SurfaceInteraction si; }; bool delta = false; Float pdfFwd = 0, pdfRev = 0;
<<Vertex Public Methods>> 
Vertex() : ei() {} Vertex(VertexType type, const EndpointInteraction &ei, const Spectrum &beta) : type(type), beta(beta), ei(ei) { } Vertex(const SurfaceInteraction &si, const Spectrum &beta) : type(VertexType::Surface), beta(beta), si(si) { } static inline Vertex CreateCamera(const Camera *camera, const Ray &ray, const Spectrum &beta); static inline Vertex CreateCamera(const Camera *camera, const Interaction &it, const Spectrum &beta); static inline Vertex CreateLight(const Light *light, const Ray &ray, const Normal3f &nLight, const Spectrum &Le, Float pdf); static inline Vertex CreateLight(const EndpointInteraction &ei, const Spectrum &beta, Float pdf); static inline Vertex CreateMedium(const MediumInteraction &mi, const Spectrum &beta, Float pdf, const Vertex &prev); static inline Vertex CreateSurface(const SurfaceInteraction &si, const Spectrum &beta, Float pdf, const Vertex &prev); Vertex(const MediumInteraction &mi, const Spectrum &beta) : type(VertexType::Medium), beta(beta), mi(mi) { } const Interaction &GetInteraction() const { switch (type) { case VertexType::Medium: return mi; case VertexType::Surface: return si; default: return ei; } } const Point3f &p() const { return GetInteraction().p; } Float time() const { return GetInteraction().time; } const Normal3f &ng() const { return GetInteraction().n; } const Normal3f &ns() const { if (type == VertexType::Surface) return si.shading.n; else return GetInteraction().n; } bool IsOnSurface() const { return ng() != Normal3f(); } Spectrum f(const Vertex &next) const { Vector3f wi = Normalize(next.p() - p()); switch (type) { case VertexType::Surface: return si.bsdf->f(si.wo, wi); case VertexType::Medium: return mi.phase->p(mi.wo, wi); } } bool IsConnectible() const { switch (type) { case VertexType::Medium: return true; case VertexType::Light: return (ei.light->flags & (int)LightFlags::DeltaDirection) == 0; case VertexType::Camera: return true; case VertexType::Surface: return si.bsdf->NumComponents( BxDFType(BSDF_DIFFUSE | BSDF_GLOSSY | BSDF_REFLECTION | BSDF_TRANSMISSION)) > 0; } } bool IsLight() const { return type == VertexType::Light || (type == VertexType::Surface && si.primitive->GetAreaLight()); } bool IsDeltaLight() const { return type == VertexType::Light && ei.light && ::IsDeltaLight(ei.light->flags); } bool IsInfiniteLight() const { return type == VertexType::Light && (!ei.light || ei.light->flags & (int)LightFlags::Infinite); } Spectrum Le(const Scene &scene, const Vertex &v) const { if (!IsLight()) return Spectrum(0.f); Vector3f w = Normalize(v.p() - p()); if (IsInfiniteLight()) { <<Return emitted radiance for infinite light sources>> 
Spectrum Le(0.f); for (const auto &light : scene.lights) Le += light->Le(Ray(p(), -w)); return Le;
} else { const AreaLight *light = si.primitive->GetAreaLight(); return light->L(si, w); } } friend std::ostream& operator<<(std::ostream& os, const Vertex &v) { os << "Vertex[" << std::endl << " type = "; switch (v.type) { case VertexType::Camera: os << "camera"; break; case VertexType::Light: os << "light"; break; case VertexType::Surface: os << "surface"; break; case VertexType::Medium : os << "medium"; break; } os << "," << std::endl << " connectible = " << v.IsConnectible() << "," << std::endl << " p = " << v.p() << "," << std::endl << " n = " << v.ng() << "," << std::endl << " pdfFwd = " << v.pdfFwd << "," << std::endl << " pdfRev = " << v.pdfRev << "," << std::endl << " beta = " << v.beta << std::endl << "]" << std::endl; return os; }; Float ConvertDensity(Float pdf, const Vertex &next) const { <<Return solid angle density if next is an infinite area light>> 
if (next.IsInfiniteLight()) return pdf;
Vector3f w = next.p() - p(); Float invDist2 = 1 / w.LengthSquared(); if (next.IsOnSurface()) pdf *= AbsDot(next.ng(), w * std::sqrt(invDist2)); return pdf * invDist2; } Float Pdf(const Scene &scene, const Vertex *prev, const Vertex &next) const { if (type == VertexType::Light) return PdfLight(scene, next); <<Compute directions to preceding and next vertex>> 
Vector3f wp, wn = Normalize(next.p() - p()); if (prev) wp = Normalize(prev->p() - p());
<<Compute directional density depending on the vertex type>> 
Float pdf, unused; if (type == VertexType::Camera) ei.camera->Pdf_We(ei.SpawnRay(wn), &unused, &pdf); else if (type == VertexType::Surface) pdf = si.bsdf->Pdf(wp, wn); else if (type == VertexType::Medium) pdf = mi.phase->p(wp, wn);
<<Return probability per unit area at vertex next>> 
return ConvertDensity(pdf, next);
} Float PdfLight(const Scene &scene, const Vertex &v) const { Vector3f w = v.p() - p(); Float invDist2 = 1 / w.LengthSquared(); w *= std::sqrt(invDist2); Float pdf; if (IsInfiniteLight()) { <<Compute planar sampling density for infinite light sources>> 
Point3f worldCenter; Float worldRadius; scene.WorldBound().BoundingSphere(&worldCenter, &worldRadius); pdf = 1 / (Pi * worldRadius * worldRadius);
} else { <<Get pointer light to the light source at the vertex>>  <<Compute sampling density for non-infinite light sources>> 
Float pdfPos, pdfDir; light->Pdf_Le(Ray(p(), w, time()), ng(), &pdfPos, &pdfDir); pdf = pdfDir * invDist2;
} if (v.IsOnSurface()) pdf *= AbsDot(v.ng(), w); return pdf; } Float PdfLightOrigin(const Scene &scene, const Vertex &v, const Distribution1D &lightDistr) const { Vector3f w = Normalize(v.p() - p()); if (IsInfiniteLight()) { <<Return solid angle density for infinite light sources>> 
return InfiniteLightDensity(scene, lightDistr, w);
} else { <<Return solid angle density for non-infinite light sources>> 
Float pdfPos, pdfDir, pdfChoice = 0; <<Get pointer light to the light source at the vertex>>  <<Compute the discrete probability of sampling light, pdfChoice>> 
for (size_t i = 0; i < scene.lights.size(); ++i) { if (scene.lights[i].get() == light) { pdfChoice = lightDistr.DiscretePDF(i); break; } }
light->Pdf_Le(Ray(p(), w, time()), ng(), &pdfPos, &pdfDir); return pdfPos * pdfChoice;
} }
};

<<Vertex Public Data>>= 
VertexType type;

Altogether, four different types of path vertices are supported in pbrt.

<<BDPT Helper Definitions>>= 
enum class VertexType { Camera, Light, Surface, Medium };

The beta member variable is analogous to the beta Subscript variable in the volumetric path tracer (Section 15.3.1) or the weight carried by particles in the SPPMIntegrator: it contains the product of the BSDF or phase function values, transmittances, and cosine terms for the vertices in the path generated so far, divided by their respective sampling PDFs. For the light subpath, they also include the emitted radiance divided by the density of the emission position and direction. For the camera subpath, radiance is replaced by importance.

<<Vertex Public Data>>+=  
Spectrum beta;

Instances of various types of Interactions represent type-specific data about the vertex. This information is arranged as a space-efficient C++ union since only one of the entries is used at a time.

<<Vertex Public Data>>+=  
union { EndpointInteraction ei; MediumInteraction mi; SurfaceInteraction si; };

EndpointInteraction is a new interaction implementation that is used only by BDPT. It records the position of a path endpoint—i.e., a position on a light source or the lens of the camera—and stores a pointer to the camera or light in question.

<<EndpointInteraction Declarations>>= 
struct EndpointInteraction : Interaction { union { const Camera *camera; const Light *light; }; <<EndpointInteraction Public Methods>> 
EndpointInteraction() : Interaction(), light(nullptr) { } EndpointInteraction(const Interaction &it, const Camera *camera) : Interaction(it), camera(camera) { } EndpointInteraction(const Camera *camera, const Ray &ray) : Interaction(ray.o, ray.time, ray.medium), camera(camera) { } EndpointInteraction(const Light *light, const Ray &r, const Normal3f &nl) : Interaction(r.o, r.time, r.medium), light(light) { n = nl; } EndpointInteraction(const Interaction &it, const Light *light) : Interaction(it), light(light) { } EndpointInteraction(const Ray &ray) : Interaction(ray(1), ray.time, ray.medium), light(nullptr) { n = Normal3f(-ray.d); }
};

There are multiple constructors that initialize the EndpointInteraction contents using a pointer and either an existing Interaction or a sampled ray. For brevity, we only show the constructors for light endpoints.

<<EndpointInteraction Public Methods>>= 
EndpointInteraction(const Light *light, const Ray &r, const Normal3f &nl) : Interaction(r.o, r.time, r.medium), light(light) { n = nl; }

<<EndpointInteraction Public Methods>>+= 
EndpointInteraction(const Interaction &it, const Light *light) : Interaction(it), light(light) { }

A range of static helper functions create Vertex instances for the various types of path vertices. We’ll only include their declarations here, as their implementations are all straightforward. We could instead have provided a range of overloaded constructors that took various Interaction types as parameters, but we think that having the name of the type of vertex being created explicit in a function call makes the following code easier to read.

<<Vertex Public Methods>>= 
static inline Vertex CreateCamera(const Camera *camera, const Ray &ray, const Spectrum &beta); static inline Vertex CreateCamera(const Camera *camera, const Interaction &it, const Spectrum &beta); static inline Vertex CreateLight(const Light *light, const Ray &ray, const Normal3f &nLight, const Spectrum &Le, Float pdf); static inline Vertex CreateLight(const EndpointInteraction &ei, const Spectrum &beta, Float pdf); static inline Vertex CreateMedium(const MediumInteraction &mi, const Spectrum &beta, Float pdf, const Vertex &prev); static inline Vertex CreateSurface(const SurfaceInteraction &si, const Spectrum &beta, Float pdf, const Vertex &prev);

It is often necessary to access the core fields in Interaction that are common to all types of vertices; the Vertex::GetInteraction() method extracts this shared part. Since Vertex::mi, Vertex::si, and Vertex::ei all derive from Interaction and are part of the same union and thus their base Interactions are at the same location in memory, the conditional logic below should be removed by the compiler.

<<Vertex Public Methods>>+=  
const Interaction &GetInteraction() const { switch (type) { case VertexType::Medium: return mi; case VertexType::Surface: return si; default: return ei; } }

The convenience function Vertex::p() returns the vertex position. We omit definitions of Vertex::time(), Vertex::ng(), and Vertex::ns(), which are defined analogously, and return the time, geometric normal, and shading normal, respectively, of the vertex.

<<Vertex Public Methods>>+=  
const Point3f &p() const { return GetInteraction().p; }

The delta attribute is only used by surface interactions and records whether a Dirac delta function was sampled (e.g., when light is scattered by a perfectly specular material).

<<Vertex Public Data>>+=  
bool delta = false;

A simple way to find out whether a vertex (including endpoints) is located on a surface is to check whether Vertex::ng() returns a nonzero result.

<<Vertex Public Methods>>+=  
bool IsOnSurface() const { return ng() != Normal3f(); }

Vertex::f() evaluates the portion of the measurement equation, (16.1), associated with a vertex. This method only needs to handle surface and medium vertices since the BDPT implementation only invokes it in those cases. Note that the next vertex in the path is the only one passed to this method: though the direction to the predecessor vertex is needed to evaluate the BRDF or phase function, this information is already available in Interaction::wo from when the Vertex was first created.

<<Vertex Public Methods>>+=  
Spectrum f(const Vertex &next) const { Vector3f wi = Normalize(next.p() - p()); switch (type) { case VertexType::Surface: return si.bsdf->f(si.wo, wi); case VertexType::Medium: return mi.phase->p(mi.wo, wi); } }

The Vertex::IsConnectible() method returns a Boolean value that indicates whether a connection strategy involving the current vertex can in principle succeed. If, for example, the vertex is a surface interaction whose BSDF only consists of Dirac delta components, then we can never successfully connect it to a subpath vertex in the other path: there’s zero probability of choosing a direction where the delta distribution is nonzero. The implementation assumes that medium and camera vertices are always connectible (the latter assumption would have to be modified if support for orthographic cameras is added).

<<Vertex Public Methods>>+=  
bool IsConnectible() const { switch (type) { case VertexType::Medium: return true; case VertexType::Light: return (ei.light->flags & (int)LightFlags::DeltaDirection) == 0; case VertexType::Camera: return true; case VertexType::Surface: return si.bsdf->NumComponents( BxDFType(BSDF_DIFFUSE | BSDF_GLOSSY | BSDF_REFLECTION | BSDF_TRANSMISSION)) > 0; } }

A few helper methods are useful for working with lights—these are necessary to deal with the considerable variety of light sources supported by pbrt.

For instance, when the Primitive underlying a surface interaction vertex is itself an area light, the vertex can take on different roles depending on the BDPT connection strategy: it can be re-interpreted as a light source and used as a path endpoint, or it can serve as a normal scattering event to generate paths of greater length. The Vertex::IsLight() method therefore provides a comprehensive test for whether a vertex can be interpreted as a light source.

<<Vertex Public Methods>>+=  
bool IsLight() const { return type == VertexType::Light || (type == VertexType::Surface && si.primitive->GetAreaLight()); }

Light sources that have an emission profile that contains a Dirac delta distribution must be treated specially in the computation of multiple importance sampling weights; Vertex::IsDeltaLight() checks for this case.

<<Vertex Public Methods>>+=  
bool IsDeltaLight() const { return type == VertexType::Light && ei.light && ::IsDeltaLight(ei.light->flags); }

The Vertex::IsInfiniteLight() method indicates whether a vertex is associated with an infinite area light. Such vertices can be created by sampling an emitted ray from an InfiniteAreaLight or by tracing a ray from the camera that escapes into the environment. In the latter case, the vertex is marked with the type VertexType::Light, but ei.light stores nullptr since no specific light source was intersected.

<<Vertex Public Methods>>+=  
bool IsInfiniteLight() const { return type == VertexType::Light && (!ei.light || ei.light->flags & (int)LightFlags::Infinite); }

Finally, Le() can be used to find emitted radiance from an intersected light source toward another vertex.

<<Vertex Public Methods>>+=  
Spectrum Le(const Scene &scene, const Vertex &v) const { if (!IsLight()) return Spectrum(0.f); Vector3f w = Normalize(v.p() - p()); if (IsInfiniteLight()) { <<Return emitted radiance for infinite light sources>> 
Spectrum Le(0.f); for (const auto &light : scene.lights) Le += light->Le(Ray(p(), -w)); return Le;
} else { const AreaLight *light = si.primitive->GetAreaLight(); return light->L(si, w); } }

<<Return emitted radiance for infinite light sources>>= 
Spectrum Le(0.f); for (const auto &light : scene.lights) Le += light->Le(Ray(p(), -w)); return Le;

Probability Densities

BDPT’s multiple importance sampling code requires detailed information about the probability density of light-carrying paths with respect to a range of different path sampling strategies. It is crucial that these densities are expressed in the same probability measure so that ratios of their probabilities are meaningful. The implementation here uses the area product measure for path probabilities. It expresses the density of a path as the product of the densities of its individual vertices, which are in turn given in a simple common (and consistent) measure: probability per unit area. This is the same measure as was initially used to derive the surface form of the LTE in Section 14.4.3.

Recall from Section 5.5 that the Jacobian of the mapping from solid angles to surface area involves the inverse squared distance and the cosine of angle between the geometric normal at next and wn (assuming next is a surface vertex—if it is a point in a participating medium, there is no cosine term (Section 15.1.1)). The ConvertDensity() method returns the product of this Jacobian (computed from the vertex attributes) and the pdf parameter, which should express a solid angle density at the vertex. (Infinite area light sources need special handling here; this case is discussed later, in Section 16.3.5.)

<<Vertex Public Methods>>+=  
Float ConvertDensity(Float pdf, const Vertex &next) const { <<Return solid angle density if next is an infinite area light>> 
if (next.IsInfiniteLight()) return pdf;
Vector3f w = next.p() - p(); Float invDist2 = 1 / w.LengthSquared(); if (next.IsOnSurface()) pdf *= AbsDot(next.ng(), w * std::sqrt(invDist2)); return pdf * invDist2; }

Each vertex has two densities: the first, pdfFwd, stores forward density of the current vertex, which is the probability per unit area of the current vertex as generated by the path sampling algorithm. The second density, pdfRev, is the hypothetical probability density of the vertex if the direction of light transport was reversed—that is, if radiance transport was used in place of importance transport for the camera path and vice versa for the light path. This reverse density will be crucial for computing MIS weights in Section 16.3.4.

<<Vertex Public Data>>+= 
Float pdfFwd = 0, pdfRev = 0;

The Vertex::Pdf() method returns the probability per unit area of the sampling technique associated with a given vertex. Given a preceding vertex prev, it evaluates the density for sampling the vertex next for rays leaving the vertex *this. The prev argument may be equal to nullptr for path endpoints (i.e., cameras or light sources), which have no predecessor. Light sources require some extra care and are handled separately via the PdfLight() method that will be discussed shortly.

<<Vertex Public Methods>>+=  
Float Pdf(const Scene &scene, const Vertex *prev, const Vertex &next) const { if (type == VertexType::Light) return PdfLight(scene, next); <<Compute directions to preceding and next vertex>> 
Vector3f wp, wn = Normalize(next.p() - p()); if (prev) wp = Normalize(prev->p() - p());
<<Compute directional density depending on the vertex type>> 
Float pdf, unused; if (type == VertexType::Camera) ei.camera->Pdf_We(ei.SpawnRay(wn), &unused, &pdf); else if (type == VertexType::Surface) pdf = si.bsdf->Pdf(wp, wn); else if (type == VertexType::Medium) pdf = mi.phase->p(wp, wn);
<<Return probability per unit area at vertex next>> 
return ConvertDensity(pdf, next);
}

For all other vertex types, the function first computes normalized directions to the preceding and next vertex (if present).

<<Compute directions to preceding and next vertex>>= 
Vector3f wp, wn = Normalize(next.p() - p()); if (prev) wp = Normalize(prev->p() - p());

Depending on the vertex type, Pdf() invokes the appropriate PDF method and stores the probability per unit solid angle for sampling the direction to next in the variable pdf.

<<Compute directional density depending on the vertex type>>= 
Float pdf, unused; if (type == VertexType::Camera) ei.camera->Pdf_We(ei.SpawnRay(wn), &unused, &pdf); else if (type == VertexType::Surface) pdf = si.bsdf->Pdf(wp, wn); else if (type == VertexType::Medium) pdf = mi.phase->p(wp, wn);

Finally, the solid angle density is converted to a probability per unit area at next.

<<Return probability per unit area at vertex next>>= 
return ConvertDensity(pdf, next);

Light-emitting vertices can be created in two different ways: by using a sampling routine like Light::Sample_Le(), or by intersecting an emissive surface via ray tracing. To be able to compare these different strategies as part of a multiple importance sampling scheme, it is necessary to know the corresponding probability per unit area for a light vertex. This task is handled by the PdfLight() method.

Its definition resembles that of Vertex::Pdf(): it computes the direction from the current vertex to the provided vertex and invokes Light::Pdf_Le() to retrieve the solid angle density of the underlying sampling strategy, which is subsequently converted into a density per unit area at v. In contrast to Vertex::Pdf(), this method also treats surface vertices located on area lights as if they were light source vertices. Once more, there is a special case for infinite area lights, which we postpone until Section 16.3.5.

<<Vertex Public Methods>>+=  
Float PdfLight(const Scene &scene, const Vertex &v) const { Vector3f w = v.p() - p(); Float invDist2 = 1 / w.LengthSquared(); w *= std::sqrt(invDist2); Float pdf; if (IsInfiniteLight()) { <<Compute planar sampling density for infinite light sources>> 
Point3f worldCenter; Float worldRadius; scene.WorldBound().BoundingSphere(&worldCenter, &worldRadius); pdf = 1 / (Pi * worldRadius * worldRadius);
} else { <<Get pointer light to the light source at the vertex>>  <<Compute sampling density for non-infinite light sources>> 
Float pdfPos, pdfDir; light->Pdf_Le(Ray(p(), w, time()), ng(), &pdfPos, &pdfDir); pdf = pdfDir * invDist2;
} if (v.IsOnSurface()) pdf *= AbsDot(v.ng(), w); return pdf; }

Depending on the vertex type, the pointer to the light source implementation must be obtained from one of two different locations.

<<Get pointer light to the light source at the vertex>>= 

<<Compute sampling density for non-infinite light sources>>= 
Float pdfPos, pdfDir; light->Pdf_Le(Ray(p(), w, time()), ng(), &pdfPos, &pdfDir); pdf = pdfDir * invDist2;

By symmetry, we would now expect a dual routine Vertex::PdfCamera() that applies to camera endpoints. However, cameras in pbrt are never represented using explicit geometry: thus, they cannot be reached by ray intersections, which eliminates the need for a dedicated query function. If desired, a perfectly symmetric implementation could be achieved by instantiating scene geometry that is tagged with an “area camera” analogous to area lights. This increases the set of possible BDPT connection strategies, though their benefit is negligible in most scenes due to the low probability of intersecting the camera.

Note that the Pdf() and PdfLight() methods use the directional probability density of the importance strategy implemented at the current vertex as measured at the location of another given vertex. However, this is not enough to fully characterize the behavior of path endpoints, whose sampling routines generate rays from a 4D distribution. An additional PdfLightOrigin() method fills the gap by providing information about the spatial distribution of samples on the light sources themselves. For the same reason as before, a dedicated PdfCameraOrigin() method for camera endpoints is not needed.

<<Vertex Public Methods>>+= 
Float PdfLightOrigin(const Scene &scene, const Vertex &v, const Distribution1D &lightDistr) const { Vector3f w = Normalize(v.p() - p()); if (IsInfiniteLight()) { <<Return solid angle density for infinite light sources>> 
return InfiniteLightDensity(scene, lightDistr, w);
} else { <<Return solid angle density for non-infinite light sources>> 
Float pdfPos, pdfDir, pdfChoice = 0; <<Get pointer light to the light source at the vertex>>  <<Compute the discrete probability of sampling light, pdfChoice>> 
for (size_t i = 0; i < scene.lights.size(); ++i) { if (scene.lights[i].get() == light) { pdfChoice = lightDistr.DiscretePDF(i); break; } }
light->Pdf_Le(Ray(p(), w, time()), ng(), &pdfPos, &pdfDir); return pdfPos * pdfChoice;
} }

<<Return solid angle density for non-infinite light sources>>= 
Float pdfPos, pdfDir, pdfChoice = 0; <<Get pointer light to the light source at the vertex>>  <<Compute the discrete probability of sampling light, pdfChoice>> 
for (size_t i = 0; i < scene.lights.size(); ++i) { if (scene.lights[i].get() == light) { pdfChoice = lightDistr.DiscretePDF(i); break; } }
light->Pdf_Le(Ray(p(), w, time()), ng(), &pdfPos, &pdfDir); return pdfPos * pdfChoice;

To determine the discrete probability of choosing light among the available light sources, we must find the pointer to the light source and look up the corresponding entry in lightDistr. If there are very many light sources, the linear search here will be inefficient. In that case, this computation could be implemented more efficiently by storing this probability directly in the light source class.

<<Compute the discrete probability of sampling light, pdfChoice>>= 
for (size_t i = 0; i < scene.lights.size(); ++i) { if (scene.lights[i].get() == light) { pdfChoice = lightDistr.DiscretePDF(i); break; } }

16.3.2 Generating the Camera and Light Subpaths

A symmetric pair of functions, GenerateCameraSubpath() and GenerateLightSubpath(), generates the two corresponding types of subpaths. Both do some initial work to get the path started and then call out to a second function, RandomWalk(), which takes care of sampling the following vertices and initializing the path array. Both of these functions return the number of vertices in the subpath.

<<BDPT Utility Functions>>+=  
int GenerateCameraSubpath(const Scene &scene, Sampler &sampler, MemoryArena &arena, int maxDepth, const Camera &camera, const Point2f &pFilm, Vertex *path) { if (maxDepth == 0) return 0; <<Sample initial ray for camera subpath>> 
CameraSample cameraSample; cameraSample.pFilm = pFilm; cameraSample.time = sampler.Get1D(); cameraSample.pLens = sampler.Get2D(); RayDifferential ray; Spectrum beta = camera.GenerateRayDifferential(cameraSample, &ray); ray.ScaleDifferentials(1 / std::sqrt(sampler.samplesPerPixel));
<<Generate first vertex on camera subpath and start random walk>> 
Float pdfPos, pdfDir; path[0] = Vertex::CreateCamera(&camera, ray, beta); camera.Pdf_We(ray, &pdfPos, &pdfDir); return RandomWalk(scene, ray, sampler, arena, beta, pdfDir, maxDepth - 1, TransportMode::Radiance, path + 1) + 1;
}

A camera path starts with a camera ray from Camera::GenerateRayDifferential(). As in the SamplerIntegrator, the ray’s differentials are scaled so that they reflect the actual pixel sampling density.

<<Sample initial ray for camera subpath>>= 
CameraSample cameraSample; cameraSample.pFilm = pFilm; cameraSample.time = sampler.Get1D(); cameraSample.pLens = sampler.Get2D(); RayDifferential ray; Spectrum beta = camera.GenerateRayDifferential(cameraSample, &ray); ray.ScaleDifferentials(1 / std::sqrt(sampler.samplesPerPixel));

The vertex at position path[0] is initialized with a special endpoint vertex on the camera lens (for cameras with finite apertures) or pinhole. The RandomWalk() function then takes care of generating the rest of the vertices. TransportMode reflects the quantity that is carried back to the origin of the path—hence TransportMode::Radiance is used here. Since the first element of path was already used for the endpoint vertex, RandomWalk() is invoked such that it writes sampled vertices starting at position path[1] with a maximum depth of maxDepth - 1. The function returns the total number of sampled vertices.

<<Generate first vertex on camera subpath and start random walk>>= 
Float pdfPos, pdfDir; path[0] = Vertex::CreateCamera(&camera, ray, beta); camera.Pdf_We(ray, &pdfPos, &pdfDir); return RandomWalk(scene, ray, sampler, arena, beta, pdfDir, maxDepth - 1, TransportMode::Radiance, path + 1) + 1;

The function GenerateLightSubpath() works in a similar fashion, with some minor differences corresponding to the fact that the path starts from a light source.

<<BDPT Utility Functions>>+=  
int GenerateLightSubpath(const Scene &scene, Sampler &sampler, MemoryArena &arena, int maxDepth, Float time, const Distribution1D &lightDistr, Vertex *path) { if (maxDepth == 0) return 0; <<Sample initial ray for light subpath>> 
Float lightPdf; int lightNum = lightDistr.SampleDiscrete(sampler.Get1D(), &lightPdf); const std::shared_ptr<Light> &light = scene.lights[lightNum]; RayDifferential ray; Normal3f nLight; Float pdfPos, pdfDir; Spectrum Le = light->Sample_Le(sampler.Get2D(), sampler.Get2D(), time, &ray, &nLight, &pdfPos, &pdfDir); if (pdfPos == 0 || pdfDir == 0 || Le.IsBlack()) return 0;
<<Generate first vertex on light subpath and start random walk>> 
path[0] = Vertex::CreateLight(light.get(), ray, nLight, Le, pdfPos * lightPdf); Spectrum beta = Le * AbsDot(nLight, ray.d) / (lightPdf * pdfPos * pdfDir); int nVertices = RandomWalk(scene, ray, sampler, arena, beta, pdfDir, maxDepth - 1, TransportMode::Importance, path + 1); <<Correct subpath sampling densities for infinite area lights>> 
if (path[0].IsInfiniteLight()) { <<Set spatial density of path[1] for infinite area light>> 
if (nVertices > 0) { path[1].pdfFwd = pdfPos; if (path[1].IsOnSurface()) path[1].pdfFwd *= AbsDot(ray.d, path[1].ng()); }
<<Set spatial density of path[0] for infinite area light>> 
path[0].pdfFwd = InfiniteLightDensity(scene, lightDistr, ray.d);
}
return nVertices + 1;
}

As usual in this integrator, a specific light is chosen by sampling from the provided Distribution1D. Next, an emitted ray is sampled via the light’s implementation of Light::Sample_Le().

<<Sample initial ray for light subpath>>= 
Float lightPdf; int lightNum = lightDistr.SampleDiscrete(sampler.Get1D(), &lightPdf); const std::shared_ptr<Light> &light = scene.lights[lightNum]; RayDifferential ray; Normal3f nLight; Float pdfPos, pdfDir; Spectrum Le = light->Sample_Le(sampler.Get2D(), sampler.Get2D(), time, &ray, &nLight, &pdfPos, &pdfDir); if (pdfPos == 0 || pdfDir == 0 || Le.IsBlack()) return 0;

The beta variable is initialized with the associated sampling weight, which is given by the emitted radiance multiplied by a cosine factor from the light transport equation and divided by the probability of the sample in ray-space. This step is analogous to Equation (16.15) and the approach implemented in the fragment <<Generate photonRay from light source and initialize beta>> from the particle tracing step of the SPPM integrator.

<<Generate first vertex on light subpath and start random walk>>= 
path[0] = Vertex::CreateLight(light.get(), ray, nLight, Le, pdfPos * lightPdf); Spectrum beta = Le * AbsDot(nLight, ray.d) / (lightPdf * pdfPos * pdfDir); int nVertices = RandomWalk(scene, ray, sampler, arena, beta, pdfDir, maxDepth - 1, TransportMode::Importance, path + 1); <<Correct subpath sampling densities for infinite area lights>> 
if (path[0].IsInfiniteLight()) { <<Set spatial density of path[1] for infinite area light>> 
if (nVertices > 0) { path[1].pdfFwd = pdfPos; if (path[1].IsOnSurface()) path[1].pdfFwd *= AbsDot(ray.d, path[1].ng()); }
<<Set spatial density of path[0] for infinite area light>> 
path[0].pdfFwd = InfiniteLightDensity(scene, lightDistr, ray.d);
}
return nVertices + 1;

RandomWalk() traces paths starting at an initial vertex. It assumes that a position and an outgoing direction at the corresponding path endpoint were previously sampled and that this information is provided via the input arguments ray, a path throughput weight beta, and a parameter pdf that gives the probability of sampling the ray per unit solid angle of ray.d. The parameter mode selects between importance and radiance transport (Section 16.1). The path vertices are stored in the provided path array up to a maximum number of maxDepth vertices, and the actual number of generated vertices is returned at the end.

<<BDPT Utility Functions>>+=  
int RandomWalk(const Scene &scene, RayDifferential ray, Sampler &sampler, MemoryArena &arena, Spectrum beta, Float pdf, int maxDepth, TransportMode mode, Vertex *path) { if (maxDepth == 0) return 0; int bounces = 0; <<Declare variables for forward and reverse probability densities>> 
Float pdfFwd = pdf, pdfRev = 0;
while (true) { <<Attempt to create the next subpath vertex in path>> 
MediumInteraction mi; <<Trace a ray and sample the medium, if any>> 
SurfaceInteraction isect; bool foundIntersection = scene.Intersect(ray, &isect); if (ray.medium) beta *= ray.medium->Sample(ray, sampler, arena, &mi); if (beta.IsBlack()) break; Vertex &vertex = path[bounces], &prev = path[bounces - 1];
if (mi.IsValid()) { <<Record medium interaction in path and compute forward density>> 
vertex = Vertex::CreateMedium(mi, beta, pdfFwd, prev); if (++bounces >= maxDepth) break;
<<Sample direction and compute reverse density at preceding vertex>> 
Vector3f wi; pdfFwd = pdfRev = mi.phase->Sample_p(-ray.d, &wi, sampler.Get2D()); ray = mi.SpawnRay(wi);
} else { <<Handle surface interaction for path generation>> 
if (!foundIntersection) { <<Capture escaped rays when tracing from the camera>> 
if (mode == TransportMode::Radiance) { vertex = Vertex::CreateLight(EndpointInteraction(ray), beta, pdfFwd); ++bounces; }
break; } <<Compute scattering functions for mode and skip over medium boundaries>> 
isect.ComputeScatteringFunctions(ray, arena, true, mode); if (!isect.bsdf) { ray = isect.SpawnRay(ray.d); continue; }
<<Initialize vertex with surface intersection information>> 
vertex = Vertex::CreateSurface(isect, beta, pdfFwd, prev);
if (++bounces >= maxDepth) break; <<Sample BSDF at current vertex and compute reverse probability>> 
Vector3f wi, wo = isect.wo; BxDFType type; Spectrum f = isect.bsdf->Sample_f(wo, &wi, sampler.Get2D(), &pdfFwd, BSDF_ALL, &type); if (f.IsBlack() || pdfFwd == 0.f) break; beta *= f * AbsDot(wi, isect.shading.n) / pdfFwd; pdfRev = isect.bsdf->Pdf(wi, wo, BSDF_ALL); if (type & BSDF_SPECULAR) { vertex.delta = true; pdfRev = pdfFwd = 0; } beta *= CorrectShadingNormal(isect, wo, wi, mode);
ray = isect.SpawnRay(wi);
} <<Compute reverse area density at preceding vertex>> 
prev.pdfRev = vertex.ConvertDensity(pdfRev, prev);
} return bounces; }

The two variables pdfFwd and pdfRev are updated during every loop iteration and satisfy the following invariants: at the beginning of each iteration, pdfFwd records the probability per unit solid angle of the sampled ray direction ray.d. On the other hand, pdfRev denotes the reverse probability at the end of each iteration—that is, the density of the opposite light transport mode per unit solid angle along the same ray segment.

<<Declare variables for forward and reverse probability densities>>= 
Float pdfFwd = pdf, pdfRev = 0;

<<Attempt to create the next subpath vertex in path>>= 
MediumInteraction mi; <<Trace a ray and sample the medium, if any>> 
SurfaceInteraction isect; bool foundIntersection = scene.Intersect(ray, &isect); if (ray.medium) beta *= ray.medium->Sample(ray, sampler, arena, &mi); if (beta.IsBlack()) break; Vertex &vertex = path[bounces], &prev = path[bounces - 1];
if (mi.IsValid()) { <<Record medium interaction in path and compute forward density>> 
vertex = Vertex::CreateMedium(mi, beta, pdfFwd, prev); if (++bounces >= maxDepth) break;
<<Sample direction and compute reverse density at preceding vertex>> 
Vector3f wi; pdfFwd = pdfRev = mi.phase->Sample_p(-ray.d, &wi, sampler.Get2D()); ray = mi.SpawnRay(wi);
} else { <<Handle surface interaction for path generation>> 
if (!foundIntersection) { <<Capture escaped rays when tracing from the camera>> 
if (mode == TransportMode::Radiance) { vertex = Vertex::CreateLight(EndpointInteraction(ray), beta, pdfFwd); ++bounces; }
break; } <<Compute scattering functions for mode and skip over medium boundaries>> 
isect.ComputeScatteringFunctions(ray, arena, true, mode); if (!isect.bsdf) { ray = isect.SpawnRay(ray.d); continue; }
<<Initialize vertex with surface intersection information>> 
vertex = Vertex::CreateSurface(isect, beta, pdfFwd, prev);
if (++bounces >= maxDepth) break; <<Sample BSDF at current vertex and compute reverse probability>> 
Vector3f wi, wo = isect.wo; BxDFType type; Spectrum f = isect.bsdf->Sample_f(wo, &wi, sampler.Get2D(), &pdfFwd, BSDF_ALL, &type); if (f.IsBlack() || pdfFwd == 0.f) break; beta *= f * AbsDot(wi, isect.shading.n) / pdfFwd; pdfRev = isect.bsdf->Pdf(wi, wo, BSDF_ALL); if (type & BSDF_SPECULAR) { vertex.delta = true; pdfRev = pdfFwd = 0; } beta *= CorrectShadingNormal(isect, wo, wi, mode);
ray = isect.SpawnRay(wi);
} <<Compute reverse area density at preceding vertex>> 
prev.pdfRev = vertex.ConvertDensity(pdfRev, prev);

The loop body begins by intersecting the current ray against the scene geometry. If the ray is passing through a participating medium, the call to Medium::Sample() possibly samples a scattering event between the ray and the surface. It returns the medium sampling weight, which is incorporated into the path contribution weight beta.

<<Trace a ray and sample the medium, if any>>= 
SurfaceInteraction isect; bool foundIntersection = scene.Intersect(ray, &isect); if (ray.medium) beta *= ray.medium->Sample(ray, sampler, arena, &mi); if (beta.IsBlack()) break; Vertex &vertex = path[bounces], &prev = path[bounces - 1];

When Medium::Sample() generates a medium scattering event, the corresponding Interaction is stored in a Vertex and appended at the end of the path array. The Vertex::CreateMedium() method converts the solid angle density in pdfFwd to a probability per unit area and stores the result in Vertex::pdfFwd.

<<Record medium interaction in path and compute forward density>>= 
vertex = Vertex::CreateMedium(mi, beta, pdfFwd, prev); if (++bounces >= maxDepth) break;

If the maximum path depth has not yet been exceeded, a scattered direction is sampled from the phase function and used to spawn a new ray that will be processed by the next loop iteration.

At this point, we could evaluate the phase function with swapped arguments to obtain the sampling density at the preceding vertex for a hypothetical random walk that would have produced the same scattering interactions in reverse order. Since phase functions are generally symmetric with respect to their arguments, instead we simply reuse the value computed for pdfFwd.

<<Sample direction and compute reverse density at preceding vertex>>= 
Vector3f wi; pdfFwd = pdfRev = mi.phase->Sample_p(-ray.d, &wi, sampler.Get2D()); ray = mi.SpawnRay(wi);

For surfaces, the overall structure is similar, though some extra care is required to deal with non-symmetric scattering and surfaces that mark transitions between media.

<<Handle surface interaction for path generation>>= 
if (!foundIntersection) { <<Capture escaped rays when tracing from the camera>> 
if (mode == TransportMode::Radiance) { vertex = Vertex::CreateLight(EndpointInteraction(ray), beta, pdfFwd); ++bounces; }
break; } <<Compute scattering functions for mode and skip over medium boundaries>> 
isect.ComputeScatteringFunctions(ray, arena, true, mode); if (!isect.bsdf) { ray = isect.SpawnRay(ray.d); continue; }
<<Initialize vertex with surface intersection information>> 
vertex = Vertex::CreateSurface(isect, beta, pdfFwd, prev);
if (++bounces >= maxDepth) break; <<Sample BSDF at current vertex and compute reverse probability>> 
Vector3f wi, wo = isect.wo; BxDFType type; Spectrum f = isect.bsdf->Sample_f(wo, &wi, sampler.Get2D(), &pdfFwd, BSDF_ALL, &type); if (f.IsBlack() || pdfFwd == 0.f) break; beta *= f * AbsDot(wi, isect.shading.n) / pdfFwd; pdfRev = isect.bsdf->Pdf(wi, wo, BSDF_ALL); if (type & BSDF_SPECULAR) { vertex.delta = true; pdfRev = pdfFwd = 0; } beta *= CorrectShadingNormal(isect, wo, wi, mode);
ray = isect.SpawnRay(wi);

The fragment <<Capture escaped rays when tracing from the camera>> is necessary to support infinite area lights. It will be discussed in Section 16.3.5. The following fragment, <<Compute scattering functions for mode and skip over medium boundaries>>, is analogous to <<Compute scattering functions and skip over medium boundaries>> from the basic path tracer except that scattering functions are requested for the current light transport mode (radiance or importance transport) using the mode parameter.

<<Compute scattering functions for mode and skip over medium boundaries>>= 
isect.ComputeScatteringFunctions(ray, arena, true, mode); if (!isect.bsdf) { ray = isect.SpawnRay(ray.d); continue; }

Given a valid intersection, the current path vertex is initialized with the corresponding surface intersection vertex, where, again, the solid angle density pdfFwd is converted to an area density before being stored in Vertex::pdfFwd.

<<Initialize vertex with surface intersection information>>= 
vertex = Vertex::CreateSurface(isect, beta, pdfFwd, prev);

If the maximum path depth has not yet been exceeded, a scattered direction is sampled from the BSDF and the path contribution in beta is updated. For the surface case, we can’t generally assume that BSDF::Pdf() is symmetric; hence we must re-evaluate the sampling density with swapped arguments to obtain pdfRev. In case of a specular sampling event, we mark the vertex using the flag Vertex::delta and set pdfFwd and pdfRev to 0 to indicate that the underlying interaction has no continuous density function. Finally, we correct for non-symmetry related to the use of shading normals (see Section 16.1.3 for details).

<<Sample BSDF at current vertex and compute reverse probability>>= 
Vector3f wi, wo = isect.wo; BxDFType type; Spectrum f = isect.bsdf->Sample_f(wo, &wi, sampler.Get2D(), &pdfFwd, BSDF_ALL, &type); if (f.IsBlack() || pdfFwd == 0.f) break; beta *= f * AbsDot(wi, isect.shading.n) / pdfFwd; pdfRev = isect.bsdf->Pdf(wi, wo, BSDF_ALL); if (type & BSDF_SPECULAR) { vertex.delta = true; pdfRev = pdfFwd = 0; } beta *= CorrectShadingNormal(isect, wo, wi, mode);

The loop wraps up by converting the reverse density pdfRev to a probability per unit area and storing it in the Vertex data structure of the preceding vertex.

<<Compute reverse area density at preceding vertex>>= 
prev.pdfRev = vertex.ConvertDensity(pdfRev, prev);

16.3.3 Subpath Connections

The ConnectBDPT() function takes the light and camera subpaths and the number of vertices s and t to use from each one, respectively. It returns the corresponding strategy’s contribution.

The connection strategy with t equals 1 uses only a single camera vertex, the camera’s position; the raster position of the path’s contribution is then based on which pixel the last vertex of the light subpath is visible in (if any). In this case, the resulting position is returned via the pRaster argument.

<<BDPT Method Definitions>>= 
Spectrum ConnectBDPT(const Scene &scene, Vertex *lightVertices, Vertex *cameraVertices, int s, int t, const Distribution1D &lightDistr, const Camera &camera, Sampler &sampler, Point2f *pRaster, Float *misWeightPtr) { Spectrum L(0.f); <<Ignore invalid connections related to infinite area lights>> 
if (t > 1 && s != 0 && cameraVertices[t - 1].type == VertexType::Light) return Spectrum(0.f);
<<Perform connection and write contribution to L>> 
Vertex sampled; if (s == 0) { <<Interpret the camera subpath as a complete path>> 
const Vertex &pt = cameraVertices[t - 1]; if (pt.IsLight()) L = pt.Le(scene, cameraVertices[t - 2]) * pt.beta;
} else if (t == 1) { <<Sample a point on the camera and connect it to the light subpath>> 
const Vertex &qs = lightVertices[s - 1]; if (qs.IsConnectible()) { VisibilityTester vis; Vector3f wi; Float pdf; Spectrum Wi = camera.Sample_Wi(qs.GetInteraction(), sampler.Get2D(), &wi, &pdf, pRaster, &vis); if (pdf > 0 && !Wi.IsBlack()) { <<Initialize dynamically sampled vertex and L for t equals 1 case>> 
sampled = Vertex::CreateCamera(&camera, vis.P1(), Wi / pdf); L = qs.beta * qs.f(sampled) * vis.Tr(scene, sampler) * sampled.beta; if (qs.IsOnSurface()) L *= AbsDot(wi, qs.ns());
} }
} else if (s == 1) { <<Sample a point on a light and connect it to the camera subpath>> 
const Vertex &pt = cameraVertices[t-1]; if (pt.IsConnectible()) { Float lightPdf; VisibilityTester vis; Vector3f wi; Float pdf; int lightNum = lightDistr.SampleDiscrete(sampler.Get1D(), &lightPdf); const std::shared_ptr<Light> &light = scene.lights[lightNum]; Spectrum lightWeight = light->Sample_Li(pt.GetInteraction(), sampler.Get2D(), &wi, &pdf, &vis); if (pdf > 0 && !lightWeight.IsBlack()) { EndpointInteraction ei(vis.P1(), light.get()); sampled = Vertex::CreateLight(ei, lightWeight / (pdf * lightPdf), 0); sampled.pdfFwd = sampled.PdfLightOrigin(scene, pt, lightDistr); L = pt.beta * pt.f(sampled) * vis.Tr(scene, sampler) * sampled.beta; if (pt.IsOnSurface()) L *= AbsDot(wi, pt.ns()); } }
} else { <<Handle all other bidirectional connection cases>> 
const Vertex &qs = lightVertices[s - 1], &pt = cameraVertices[t - 1]; if (qs.IsConnectible() && pt.IsConnectible()) { L = qs.beta * qs.f(pt) * pt.f(qs) * pt.beta; if (!L.IsBlack()) L *= G(scene, sampler, qs, pt); }
}
<<Compute MIS weight for connection strategy>> 
Float misWeight = L.IsBlack() ? 0.f : MISWeight(scene, lightVertices, cameraVertices, sampled, s, t, lightDistr); L *= misWeight; if (misWeightPtr) *misWeightPtr = misWeight;
return L; }

A number of cases must be considered when handling connections; special handling is needed for those involving short subpaths with only zero or one vertex. Some strategies dynamically sample an additional vertex, which is stored in the temporary variable sampled.

<<Perform connection and write contribution to L>>= 
Vertex sampled; if (s == 0) { <<Interpret the camera subpath as a complete path>> 
const Vertex &pt = cameraVertices[t - 1]; if (pt.IsLight()) L = pt.Le(scene, cameraVertices[t - 2]) * pt.beta;
} else if (t == 1) { <<Sample a point on the camera and connect it to the light subpath>> 
const Vertex &qs = lightVertices[s - 1]; if (qs.IsConnectible()) { VisibilityTester vis; Vector3f wi; Float pdf; Spectrum Wi = camera.Sample_Wi(qs.GetInteraction(), sampler.Get2D(), &wi, &pdf, pRaster, &vis); if (pdf > 0 && !Wi.IsBlack()) { <<Initialize dynamically sampled vertex and L for t equals 1 case>> 
sampled = Vertex::CreateCamera(&camera, vis.P1(), Wi / pdf); L = qs.beta * qs.f(sampled) * vis.Tr(scene, sampler) * sampled.beta; if (qs.IsOnSurface()) L *= AbsDot(wi, qs.ns());
} }
} else if (s == 1) { <<Sample a point on a light and connect it to the camera subpath>> 
const Vertex &pt = cameraVertices[t-1]; if (pt.IsConnectible()) { Float lightPdf; VisibilityTester vis; Vector3f wi; Float pdf; int lightNum = lightDistr.SampleDiscrete(sampler.Get1D(), &lightPdf); const std::shared_ptr<Light> &light = scene.lights[lightNum]; Spectrum lightWeight = light->Sample_Li(pt.GetInteraction(), sampler.Get2D(), &wi, &pdf, &vis); if (pdf > 0 && !lightWeight.IsBlack()) { EndpointInteraction ei(vis.P1(), light.get()); sampled = Vertex::CreateLight(ei, lightWeight / (pdf * lightPdf), 0); sampled.pdfFwd = sampled.PdfLightOrigin(scene, pt, lightDistr); L = pt.beta * pt.f(sampled) * vis.Tr(scene, sampler) * sampled.beta; if (pt.IsOnSurface()) L *= AbsDot(wi, pt.ns()); } }
} else { <<Handle all other bidirectional connection cases>> 
const Vertex &qs = lightVertices[s - 1], &pt = cameraVertices[t - 1]; if (qs.IsConnectible() && pt.IsConnectible()) { L = qs.beta * qs.f(pt) * pt.f(qs) * pt.beta; if (!L.IsBlack()) L *= G(scene, sampler, qs, pt); }
}

The first case ( s equals 0 ) applies when no vertices on the light subpath are used and can only succeed when the camera subpath normal p 0 comma normal p 1 comma ellipsis comma normal p Subscript t minus 1 Baseline is already a complete path—that is, when vertex normal p Subscript t minus 1 can be interpreted as a light source. In this case, L is set to the product of the path throughput weight and the emission at normal p Subscript t minus 1 .

<<Interpret the camera subpath as a complete path>>= 
const Vertex &pt = cameraVertices[t - 1]; if (pt.IsLight()) L = pt.Le(scene, cameraVertices[t - 2]) * pt.beta;

The second case applies when t equals 1 —that is, when a prefix of the light subpath is directly connected to the camera (Figure 16.15). To permit optimized importance sampling strategies analogous to direct illumination routines for light sources, we will ignore the actual camera vertex normal p 0 and sample a new one using Camera::Sample_Wi()—this optimization corresponds to the second bullet listed at the beginning of Section 16.3. This type of connection can only succeed if the light subpath vertex normal q Subscript s minus 1 supports sampled connections; otherwise the BSDF at normal q Subscript s minus 1 will certainly return 0 and there’s no reason to attempt a connection.

Figure 16.15: t equals 1 Sampling Strategy for BDPT. We’d like to connect a subset of the light subpath to the camera. Given the last vertex of the light subpath, Camera::Sample_Wi() samples a vertex on the lens normal p 0 corresponding to a ray leaving the camera to the light path vertex (if there is such a ray that intersects the film).

<<Sample a point on the camera and connect it to the light subpath>>= 
const Vertex &qs = lightVertices[s - 1]; if (qs.IsConnectible()) { VisibilityTester vis; Vector3f wi; Float pdf; Spectrum Wi = camera.Sample_Wi(qs.GetInteraction(), sampler.Get2D(), &wi, &pdf, pRaster, &vis); if (pdf > 0 && !Wi.IsBlack()) { <<Initialize dynamically sampled vertex and L for t equals 1 case>> 
sampled = Vertex::CreateCamera(&camera, vis.P1(), Wi / pdf); L = qs.beta * qs.f(sampled) * vis.Tr(scene, sampler) * sampled.beta; if (qs.IsOnSurface()) L *= AbsDot(wi, qs.ns());
} }

If the camera vertex was generated successfully, pRaster is initialized and vis holds the connection segment. Following Equation (16.1), we can compute the final contribution as the product of the subpath weights, the transmittance over the connecting segment, the BRDF or phase function, and a cosine factor when normal q Subscript s minus 1 is a surface vertex.

<<Initialize dynamically sampled vertex and L for t equals 1 case>>= 
sampled = Vertex::CreateCamera(&camera, vis.P1(), Wi / pdf); L = qs.beta * qs.f(sampled) * vis.Tr(scene, sampler) * sampled.beta; if (qs.IsOnSurface()) L *= AbsDot(wi, qs.ns());

We omit the next case, s equals 1 , here. It corresponds to performing a direct lighting calculation at the last vertex of the camera subpath. Its implementation is similar to the t equals 1 case—the main differences are that roles of lights and cameras are exchanged and that a light source must be chosen using lightDistr before a light sample can be generated.

The last case, <<Handle all other bidirectional connection cases>>, is responsible for most types of connections: it applies whenever the camera and light subpath prefixes are long enough so that no special cases are triggered (i.e., when s comma t greater-than 1 ). If we consider the generalized path contribution equation from Section 15.1.1, we have constructed camera and light subpaths with the incremental path construction approach used in Section 14.5.3 for regular path tracing. Given the throughput of these paths up to the current vertices, ModifyingAbove upper T With caret left-parenthesis normal q overbar Subscript s Baseline right-parenthesis and ModifyingAbove upper T With caret left-parenthesis normal p Subscript Baseline overbar Subscript t Baseline right-parenthesis , respectively, where

normal p overbar Subscript t Baseline equals normal p 0 comma normal p 1 comma ellipsis comma normal p Subscript t minus 1 Baseline comma

and similarly for normal q overbar Subscript s , we can find that the contribution of a path of t light vertices and s camera vertices is given by

StartLayout 1st Row 1st Column ModifyingAbove upper P With caret left-parenthesis normal q overbar Subscript s Baseline normal p overbar Subscript t Baseline right-parenthesis equals upper L Subscript normal e Superscript Baseline ModifyingAbove upper T With caret left-parenthesis normal q overbar Subscript s Baseline right-parenthesis 2nd Column left-bracket ModifyingAbove f With caret left-parenthesis normal q Subscript Baseline Subscript s minus 2 Baseline right-arrow normal q Subscript Baseline Subscript s minus 1 Baseline right-arrow normal p Subscript Baseline Subscript t minus 1 Baseline right-parenthesis ModifyingAbove upper G With caret left-parenthesis normal q Subscript s minus italic 1 Baseline left-right-arrow normal p Subscript t minus italic 1 Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column ModifyingAbove f With caret left-parenthesis normal q Subscript Baseline Subscript s minus 1 Baseline right-arrow normal p Subscript Baseline Subscript t minus 1 Baseline right-arrow normal p Subscript Baseline Subscript t minus 2 Baseline right-parenthesis right-bracket ModifyingAbove upper T With caret left-parenthesis normal p Subscript Baseline overbar Subscript t Baseline right-parenthesis upper W Subscript normal e Superscript Baseline period EndLayout

The first and last products involving the emission, importance and generalized throughput terms, upper L Subscript normal e Superscript Baseline ModifyingAbove upper T With caret left-parenthesis normal q overbar Subscript s Baseline right-parenthesis for the light path and ModifyingAbove upper T With caret left-parenthesis normal p Subscript Baseline overbar Subscript t Baseline right-parenthesis upper W Subscript normal e Superscript for the camera path, are already available in the Vertex::beta fields of the connection vertices, so we only need to compute the term in brackets to find the path’s overall contribution. The symmetric nature of BDPT is readily apparent: the final contribution is equal to the product of the subpath weights, the BRDF or phase functions and a (symmetric) generalized geometry term. Note that this strategy cannot succeed when one of the connection vertices is marked as not connectible—in this case, no connection attempt is made.

The product of subpath weights and the two BSDFs is often 0; this case happens, for example, if the connecting segment requires that light be transmitted through one of the two surfaces but the corresponding surface isn’t transmissive. In this case, it’s worth avoiding the unnecessary call to the G() function, which traces a shadow ray to test visibility.

<<Handle all other bidirectional connection cases>>= 
const Vertex &qs = lightVertices[s - 1], &pt = cameraVertices[t - 1]; if (qs.IsConnectible() && pt.IsConnectible()) { L = qs.beta * qs.f(pt) * pt.f(qs) * pt.beta; if (!L.IsBlack()) L *= G(scene, sampler, qs, pt); }

The generalized geometry term, Equation (15.5), is computed in a separate function G().

<<BDPT Utility Functions>>+=  
Spectrum G(const Scene &scene, Sampler &sampler, const Vertex &v0, const Vertex &v1) { Vector3f d = v0.p() - v1.p(); Float g = 1 / d.LengthSquared(); d *= std::sqrt(g); if (v0.IsOnSurface()) g *= AbsDot(v0.ns(), d); if (v1.IsOnSurface()) g *= AbsDot(v1.ns(), d); VisibilityTester vis(v0.GetInteraction(), v1.GetInteraction()); return g * vis.Tr(scene, sampler); }

The computation of the multiple importance sampling weight for the connected path is implemented as a separate function MISWeight(), which we discuss next.

<<Compute MIS weight for connection strategy>>= 
Float misWeight = L.IsBlack() ? 0.f : MISWeight(scene, lightVertices, cameraVertices, sampled, s, t, lightDistr); L *= misWeight; if (misWeightPtr) *misWeightPtr = misWeight;

16.3.4 Multiple Importance Sampling

Recall the example of a light pointed up at the ceiling, indirectly illuminating a room. Even without multiple importance sampling, bidirectional path tracing will do much better than path tracing by reducing the number of paths with no contribution, since the paths from the light provide camera path vertices more light-carrying targets to hit with connection segments (see Figure 16.17, which shows the effectiveness of various types of bidirectional connections). However, the image will still suffer from variance caused by paths with unexpectedly large contributions due to vertices on the camera subpaths that happen to find the bright spot on the ceiling. MIS can be applied to address this issue; it automatically recognizes that connection strategies that involve at least one scattering event on the light subpath lead to superior sampling strategies in this case. This ability comes at the cost of having to know the probabilities for constructing a path according to all available strategies and is the reason for caching the Vertex::pdfFwd and Vertex::pdfRev values earlier.

In this section, we will explain the MISWeight() function that computes the multiple importance sampling weight associated with a particular BDPT sampling strategy. It takes a light and camera subpath and an integer pair left-parenthesis s comma t right-parenthesis identifying the prefixes used by a successful BDPT connection attempt, producing a complete path normal q 0 comma ellipsis comma normal q Subscript s minus 1 Baseline comma normal p Subscript t minus 1 Baseline comma ellipsis comma normal p 0 . It iterates over all alternative strategies that could hypothetically have generated the same input path but with an earlier or later crossover point between the light and camera subpaths. Figure 16.16 shows the basic idea and Figure 16.17 shows images from the various strategies individually.

Figure 16.16: Multiple Importance Sampling in the Context of BDPT. Given a specific connection strategy left-parenthesis s equals 2 comma t equals 2 right-parenthesis shown at the top, MISWeight() considers other strategies that could have produced the same path (bottom). The t equals 0 case is omitted for simplicity. (It only makes sense in systems where the sensor of the camera can be intersected by rays.)

Figure 16.17: The Individual BDPT Strategies. Each row corresponds to light paths of a certain length. Note how almost every sampling strategy has deficiencies of some kind, evident in the form of high variance in these images. (Regular path tracing only samples the s equals 1 paths.) Applying multiple importance sampling to path contributions is an effective way to reduce this variance.

The function then reweights the path contribution using the balance heuristic from Section 13.10.1, taking all of these possible sampling strategies into account. It is straightforward to switch to other MIS variants (e.g., based on the power heuristic) if desired.

Weighted versions of the images in Figure 16.17 are shown in Figure 16.18.

Figure 16.18: Variance Reduction Due to Multiple Importance Sampling. The same sampling strategies as in Figure 16.17, but now weighted using multiple importance sampling—effectively “turning off” each strategy where it does not perform well. The final result is computed by summing all of these images.

Let left-parenthesis s comma t right-parenthesis denote the currently considered connection strategy, which concatenates a prefix normal q 0 comma ellipsis comma normal q Subscript s minus 1 Baseline from the light subpath and a (reversed) prefix normal p Subscript t minus 1 Baseline comma ellipsis comma normal p 0 from the camera subpath, producing a path normal x Subscript Baseline overbar of length n equals s plus t with vertices that we will refer to as normal x Subscript i left-parenthesis 0 less-than-or-equal-to i less-than n right-parenthesis :

normal x Subscript Baseline overbar equals left-parenthesis normal x 0 comma ellipsis comma normal x Subscript n minus italic 1 Baseline right-parenthesis equals left-parenthesis normal q 0 comma ellipsis comma normal q Subscript s minus 1 Baseline comma normal p Subscript t minus 1 Baseline comma ellipsis comma normal p 0 right-parenthesis period

Suppose that the probability per unit area of vertex normal x Subscript i is given by p Superscript right-arrow Baseline left-parenthesis normal x Subscript i Baseline right-parenthesis and p Superscript left-arrow Baseline left-parenthesis normal x Subscript i Baseline right-parenthesis for sampling strategies based on importance and radiance transport, respectively. Then the area product density of the current path is simply the product of the importance transport densities up to vertex normal x Subscript s minus 1 and the radiance transport densities for the remainder:

p Subscript s Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis equals p Superscript right-arrow Baseline left-parenthesis normal x 0 right-parenthesis midline-horizontal-ellipsis p Superscript right-arrow Baseline left-parenthesis normal x Subscript s minus 1 Baseline right-parenthesis dot p Superscript left-arrow Baseline left-parenthesis normal x Subscript s Baseline right-parenthesis midline-horizontal-ellipsis p Superscript left-arrow Baseline left-parenthesis normal x Subscript n minus 1 Baseline right-parenthesis period

Implementation-wise, the above expression is straightforward to evaluate: the importance transport densities p Superscript right-arrow Baseline left-parenthesis normal x Subscript i Baseline right-parenthesis are already cached in the Vertex::pdfFwd fields of the light subpath, and the same holds true for the radiance transport densities on the camera subpath.

More generally, we are also interested in the path density according to other connection strategies left-parenthesis i comma j right-parenthesis that could in theory have created this path. This requires that they generate paths of a compatible length—i.e., that i plus j equals s plus t . Their corresponding path density is given by

p Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis equals p Superscript right-arrow Baseline left-parenthesis normal x 0 right-parenthesis midline-horizontal-ellipsis p Superscript right-arrow Baseline left-parenthesis normal x Subscript i minus 1 Baseline right-parenthesis dot p Superscript left-arrow Baseline left-parenthesis normal x Subscript i Baseline right-parenthesis midline-horizontal-ellipsis p Superscript left-arrow Baseline left-parenthesis normal x Subscript n minus 1 Baseline right-parenthesis comma

where 0 less-than-or-equal-to i less-than-or-equal-to n . Evaluating these will also involve the reverse probabilities in Vertex::pdfRev.

Recall from Section 13.10.1 that the balance heuristic weight for strategy s out of a set of n sampling strategies with uniform sample allocation was given by

w Subscript s Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis equals StartFraction p Subscript s Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis Over sigma-summation Underscript i Endscripts p Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis EndFraction period

This is the expression we would like to evaluate in MISWeight(), though there are two practical issues that must first be addressed.

First, path densities can easily under- or overflow the range of representable values in single or even double precision. Area densities of individual vertices can be seen to be inversely proportional to the square of the scene dimensions: for instance, uniformly scaling the scene to half its size will quadruple the vertex densities. When computing the area product density of a path with 10 vertices, the same scaling operation will cause an increase in path density of approximately one million times. When working with very small or large scenes (compared to a box of unit volume), the floating point exponent of p Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis can quickly exceed the valid range.

Second, a naive MIS implementation has time complexity of upper O left-parenthesis n Superscript 4 Baseline right-parenthesis , where n is the maximum path length. Evaluation of p Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis based on Equation (16.17) involves a linear sweep over n vertices, and the MIS weight in Equation (16.18) requires another sweep over n strategies. Given that this must be done once per connection strategy for a number of strategies that is proportional to the square of the subpath length, we are left with an algorithm of quartic complexity.

We will avoid both of these issues by using a more efficient incremental computation that works with ratios of probability densities to achieve better numerical and run-time behavior.

Dividing both the numerator and denominator of Equation (16.18) by p Subscript s Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis yields

w Subscript s Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis equals StartStartFraction 1 OverOver sigma-summation Underscript i Endscripts StartFraction p Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis Over p Subscript s Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis EndFraction EndEndFraction equals left-parenthesis sigma-summation Underscript i equals 0 Overscript s minus 1 Endscripts StartFraction p Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis Over p Subscript s Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis EndFraction plus 1 plus sigma-summation Underscript i equals s plus 1 Overscript n Endscripts StartFraction p Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis Over p Subscript s Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis EndFraction right-parenthesis Superscript negative 1 Baseline period

The two sums above consider alternative strategies that would have taken additional steps on the camera or light subpath, respectively. Let us define a more concise notation for the individual summand terms:

r Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis equals StartFraction p Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis Over p Subscript s Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis EndFraction period

These satisfy the following recurrence relations:

StartLayout 1st Row 1st Column r Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis 2nd Column equals StartFraction p Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis Over p Subscript i plus 1 Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis EndFraction StartFraction p Subscript i plus 1 Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis Over p Subscript s Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis EndFraction equals StartFraction p Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis Over p Subscript i plus 1 Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis EndFraction r Subscript i plus 1 Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis left-parenthesis i less-than s right-parenthesis comma 2nd Row 1st Column r Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis 2nd Column equals StartFraction p Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis Over p Subscript i minus 1 Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis EndFraction StartFraction p Subscript i minus 1 Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis Over p Subscript s Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis EndFraction equals StartFraction p Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis Over p Subscript i minus 1 Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis EndFraction r Subscript i minus 1 Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis left-parenthesis i greater-than s right-parenthesis period EndLayout

The recurrence weights in the above equations are ratios of path densities of two adjacent sampling strategies, which differ only in how a single vertex is generated. Thus, they can be reduced to probability ratios of the affected vertex:

StartLayout 1st Row 1st Column StartFraction p Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis Over p Subscript i plus 1 Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis EndFraction 2nd Column equals StartFraction p Superscript right-arrow Baseline left-parenthesis normal x 0 right-parenthesis midline-horizontal-ellipsis p Superscript right-arrow Baseline left-parenthesis normal x Subscript i minus 1 Baseline right-parenthesis dot p Superscript left-arrow Baseline left-parenthesis normal x Subscript i Baseline right-parenthesis dot p Superscript left-arrow Baseline left-parenthesis normal x Subscript i plus 1 Baseline right-parenthesis midline-horizontal-ellipsis p Superscript left-arrow Baseline left-parenthesis normal x Subscript n minus 1 Baseline right-parenthesis Over p Superscript right-arrow Baseline left-parenthesis normal x 0 right-parenthesis midline-horizontal-ellipsis p Superscript right-arrow Baseline left-parenthesis normal x Subscript i minus 1 Baseline right-parenthesis dot p Superscript right-arrow Baseline left-parenthesis normal x Subscript i Baseline right-parenthesis dot p Superscript left-arrow Baseline left-parenthesis normal x Subscript i plus 1 Baseline right-parenthesis midline-horizontal-ellipsis p Superscript left-arrow Baseline left-parenthesis normal x Subscript n minus 1 Baseline right-parenthesis EndFraction equals StartFraction p Superscript left-arrow Baseline left-parenthesis normal x Subscript i Baseline right-parenthesis Over p Superscript right-arrow Baseline left-parenthesis normal x Subscript i Baseline right-parenthesis EndFraction comma 2nd Row 1st Column StartFraction p Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis Over p Subscript i minus 1 Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis EndFraction 2nd Column equals StartFraction p Superscript right-arrow Baseline left-parenthesis normal x 0 right-parenthesis midline-horizontal-ellipsis p Superscript right-arrow Baseline left-parenthesis normal x Subscript i minus 2 Baseline right-parenthesis dot p Superscript right-arrow Baseline left-parenthesis normal x Subscript i minus 1 Baseline right-parenthesis dot p Superscript left-arrow Baseline left-parenthesis normal x Subscript i Baseline right-parenthesis midline-horizontal-ellipsis p Superscript left-arrow Baseline left-parenthesis normal x Subscript n minus 1 Baseline right-parenthesis Over p Superscript right-arrow Baseline left-parenthesis normal x 0 right-parenthesis midline-horizontal-ellipsis p Superscript right-arrow Baseline left-parenthesis normal x Subscript i minus 2 Baseline right-parenthesis dot p Superscript left-arrow Baseline left-parenthesis normal x Subscript i minus 1 Baseline right-parenthesis dot p Superscript left-arrow Baseline left-parenthesis normal x Subscript i Baseline right-parenthesis midline-horizontal-ellipsis p Superscript left-arrow Baseline left-parenthesis normal x Subscript n minus 1 Baseline right-parenthesis EndFraction equals StartFraction p Superscript right-arrow Baseline left-parenthesis normal x Subscript i minus 1 Baseline right-parenthesis Over p Superscript left-arrow Baseline left-parenthesis normal x Subscript i minus 1 Baseline right-parenthesis EndFraction period EndLayout

Combining this result with Equation (16.20), we obtain the following recursive expression for r Subscript i :

r Subscript i Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis equals StartLayout Enlarged left-brace 1st Row 1st Column 1 comma 2nd Column if i equals s 2nd Row 1st Column StartFraction p Superscript left-arrow Baseline left-parenthesis normal x Subscript i Baseline right-parenthesis Over p Superscript right-arrow Baseline left-parenthesis normal x Subscript i Baseline right-parenthesis EndFraction r Subscript i plus 1 Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis comma 2nd Column if i less-than s period 3rd Row 1st Column StartFraction p Superscript right-arrow Baseline left-parenthesis normal x Subscript i minus 1 Baseline right-parenthesis Over p Superscript left-arrow Baseline left-parenthesis normal x Subscript i minus 1 Baseline right-parenthesis EndFraction r Subscript i minus 1 Baseline left-parenthesis normal x Subscript Baseline overbar right-parenthesis comma 2nd Column if i greater-than s period EndLayout

The main portion of the MISWeight() function accumulates these probability ratios in a temporary variable sumRi using an incremental evaluation scheme based on Equation (16.21). The last line returns the reciprocal of the r Subscript i terms according to Equation (16.19). There is also a special case at the beginning, which directly returns a weight of 1 for paths with two vertices, which can only be generated by a single strategy.

<<BDPT Utility Functions>>+= 
Float MISWeight(const Scene &scene, Vertex *lightVertices, Vertex *cameraVertices, Vertex &sampled, int s, int t, const Distribution1D &lightPdf) { if (s + t == 2) return 1; Float sumRi = 0; <<Define helper function remap0 that deals with Dirac delta functions>> 
auto remap0 = [](float f) -> float { return f != 0 ? f : 1; };
<<Temporarily update vertex properties for current strategy>> 
<<Look up connection vertices and their predecessors>> 
Vertex *qs = s > 0 ? &lightVertices[s - 1] : nullptr, *pt = t > 0 ? &cameraVertices[t - 1] : nullptr, *qsMinus = s > 1 ? &lightVertices[s - 2] : nullptr, *ptMinus = t > 1 ? &cameraVertices[t - 2] : nullptr;
<<Update sampled vertex for s equals 1 or t equals 1 strategy>> 
ScopedAssignment<Vertex> a1; if (s == 1) a1 = { qs, sampled }; else if (t == 1) a1 = { pt, sampled };
<<Mark connection vertices as non-degenerate>> 
ScopedAssignment<bool> a2, a3; if (pt) a2 = { &pt->delta, false }; if (qs) a3 = { &qs->delta, false };
<<Update reverse density of vertex normal p Subscript Baseline Subscript t minus 1 >> 
ScopedAssignment<Float> a4; if (pt) a4 = { &pt->pdfRev, s > 0 ? qs->Pdf(scene, qsMinus, *pt) : pt->PdfLightOrigin(scene, *ptMinus, lightPdf) };
<<Update reverse density of vertex normal p Subscript Baseline Subscript t minus 2 >> 
ScopedAssignment<Float> a5; if (ptMinus) a5 = { &ptMinus->pdfRev, s > 0 ? pt->Pdf(scene, qs, *ptMinus) : pt->PdfLight(scene, *ptMinus) };
<<Update reverse density of vertices normal q Subscript Baseline Subscript s minus 1 and normal q Subscript Baseline Subscript s minus 2 >> 
<<Consider hypothetical connection strategies along the camera subpath>> 
Float ri = 1; for (int i = t - 1; i > 0; --i) { ri *= remap0(cameraVertices[i].pdfRev) / remap0(cameraVertices[i].pdfFwd); if (!cameraVertices[i].delta && !cameraVertices[i - 1].delta) sumRi += ri; }
<<Consider hypothetical connection strategies along the light subpath>> 
ri = 1; for (int i = s - 1; i >= 0; --i) { ri *= remap0(lightVertices[i].pdfRev) / remap0(lightVertices[i].pdfFwd); bool deltaLightvertex = i > 0 ? lightVertices[i - 1].delta : lightVertices[0].IsDeltaLight(); if (!lightVertices[i].delta && !deltaLightvertex) sumRi += ri; }
return 1 / (1 + sumRi); }

A helper function remap0() returns its argument while mapping 0-valued arguments to 1. It is used to handle the special case of Dirac delta functions in the path, which have a continuous density of 0. Such degenerate vertices cannot be joined using any deterministic connection strategy, and their discrete probability cancels when iterating over the remaining set of strategies because it occurs both in the numerator and denominator of the summands in Equation (16.19). The purpose of the helper function is to temporarily map these densities to a nonzero value to make sure that this cancellation occurs without causing a division by 0.

<<Define helper function remap0 that deals with Dirac delta functions>>= 
auto remap0 = [](float f) -> float { return f != 0 ? f : 1; };

To avoid an excessively large number of calls to the various Vertex PDF functions, the weight computation uses the cached probabilities in Vertex::pdfFwd and Vertex::pdfRev. Since these values only capture information about the original camera and light subpaths, they must still be updated to match the full path configuration near the crossover point—specifically normal q Subscript s minus 1 and normal p Subscript t minus 1 and their predecessors. This is implemented in the somewhat technical fragment <<Temporarily update vertex properties for current strategy>>, which we discuss last.

We iterate over hypothetical strategies that would have taken additional steps from the light direction, using a temporary variable ri to store the current iterate r Subscript i . The fragment name makes reference to the camera subpath, since these extra steps involve vertices that were in reality sampled from the camera side. All vertex densities are passed through the function remap0(), and the ratio is only added to a running sum when endpoints of the current hypothetical connection strategy are marked as non-degenerate. The loop terminates before reaching the left-parenthesis n comma 0 right-parenthesis strategy, which shouldn’t be considered since the camera cannot be intersected.

<<Consider hypothetical connection strategies along the camera subpath>>= 
Float ri = 1; for (int i = t - 1; i > 0; --i) { ri *= remap0(cameraVertices[i].pdfRev) / remap0(cameraVertices[i].pdfFwd); if (!cameraVertices[i].delta && !cameraVertices[i - 1].delta) sumRi += ri; }

The next step considers additional steps along the light subpath and largely resembles the previous case. A special case arises when the current strategy would involve intersecting a light source (i.e., when s equals 0 ): this will fail when the endpoint involves a Dirac delta distribution, hence the additional test below.

<<Consider hypothetical connection strategies along the light subpath>>= 
ri = 1; for (int i = s - 1; i >= 0; --i) { ri *= remap0(lightVertices[i].pdfRev) / remap0(lightVertices[i].pdfFwd); bool deltaLightvertex = i > 0 ? lightVertices[i - 1].delta : lightVertices[0].IsDeltaLight(); if (!lightVertices[i].delta && !deltaLightvertex) sumRi += ri; }

Finally, we will define the missing fragment <<Temporarily update vertex properties for current strategy>>, which modifies Vertex attributes with new values specific to the current connection strategy left-parenthesis s comma t right-parenthesis . To reduce the amount of code needed for both the update and the subsequent cleanup operations, we will introduce a helper class ScopedAssignment that temporarily modifies a given variable and then reverts its change when program execution leaves the scope it was defined in. It stores a pointer ScopedAssignment::target to a memory location of arbitrary type (specified via the Type template parameter) and a snapshot of the original value in ScopedAssignment::backup.

<<BDPT Helper Definitions>>+=  
template <typename Type> class ScopedAssignment { public: <<ScopedAssignment Public Methods>> 
ScopedAssignment(Type *target = nullptr, Type value = Type()) : target(target) { if (target) { backup = *target; *target = value; } } ~ScopedAssignment() { if (target) *target = backup; } ScopedAssignment &operator=(ScopedAssignment &&other) { target = other.target; backup = other.backup; other.target = nullptr; return *this; }
private: Type *target, backup; };

The ScopedAssignment constructor takes a pointer to a target memory location and overwrites it with the value parameter after making a backup copy. The destructor simply reverts any changes.

<<ScopedAssignment Public Methods>>= 
ScopedAssignment(Type *target = nullptr, Type value = Type()) : target(target) { if (target) { backup = *target; *target = value; } } ~ScopedAssignment() { if (target) *target = backup; }

The main update operation then consists of finding the connection vertices and their predecessors and updating vertex probabilities and other attributes so that the two Vertex arrays reflect the chosen connection strategy.

<<Temporarily update vertex properties for current strategy>>= 
<<Look up connection vertices and their predecessors>> 
Vertex *qs = s > 0 ? &lightVertices[s - 1] : nullptr, *pt = t > 0 ? &cameraVertices[t - 1] : nullptr, *qsMinus = s > 1 ? &lightVertices[s - 2] : nullptr, *ptMinus = t > 1 ? &cameraVertices[t - 2] : nullptr;
<<Update sampled vertex for s equals 1 or t equals 1 strategy>> 
ScopedAssignment<Vertex> a1; if (s == 1) a1 = { qs, sampled }; else if (t == 1) a1 = { pt, sampled };
<<Mark connection vertices as non-degenerate>> 
ScopedAssignment<bool> a2, a3; if (pt) a2 = { &pt->delta, false }; if (qs) a3 = { &qs->delta, false };
<<Update reverse density of vertex normal p Subscript Baseline Subscript t minus 1 >> 
ScopedAssignment<Float> a4; if (pt) a4 = { &pt->pdfRev, s > 0 ? qs->Pdf(scene, qsMinus, *pt) : pt->PdfLightOrigin(scene, *ptMinus, lightPdf) };
<<Update reverse density of vertex normal p Subscript Baseline Subscript t minus 2 >> 
ScopedAssignment<Float> a5; if (ptMinus) a5 = { &ptMinus->pdfRev, s > 0 ? pt->Pdf(scene, qs, *ptMinus) : pt->PdfLight(scene, *ptMinus) };
<<Update reverse density of vertices normal q Subscript Baseline Subscript s minus 1 and normal q Subscript Baseline Subscript s minus 2 >> 

We begin by obtaining pointers to the affected connection vertices normal q Subscript s minus 1 and normal p Subscript t minus 1 and their predecessors.

<<Look up connection vertices and their predecessors>>= 
Vertex *qs = s > 0 ? &lightVertices[s - 1] : nullptr, *pt = t > 0 ? &cameraVertices[t - 1] : nullptr, *qsMinus = s > 1 ? &lightVertices[s - 2] : nullptr, *ptMinus = t > 1 ? &cameraVertices[t - 2] : nullptr;

Recall that strategies with s equals 1 or t equals 1 perform camera and light source sampling and thus generate a new endpoint. The implementation accounts for this by temporarily overriding *qs or *pt with the sampled vertex provided via the sampled argument of MISWeight().

<<Update sampled vertex for s equals 1 or t equals 1 strategy>>= 
ScopedAssignment<Vertex> a1; if (s == 1) a1 = { qs, sampled }; else if (t == 1) a1 = { pt, sampled };

Certain materials in pbrt (e.g., the UberMaterial) instantiate both specular and non-specular BxDF lobes, which requires some additional consideration at this point: suppose the specular lobe of such a material is sampled while generating the camera or light subpath. In this case, the associated Vertex will have its Vertex::delta flag set to true, causing MISWeight() to (correctly) ignore it as a hypothetical connection vertex when comparing the densities of different strategies. On the other hand, it is possible that a BDPT strategy later connects this degenerate vertex to a vertex on the other subpath using its non-specular component. In this case, its “personality” must temporarily change to that of a non-degenerate vertex. We always force the Vertex::delta attribute of the connection vertices to false to account for this possibility.

<<Mark connection vertices as non-degenerate>>= 
ScopedAssignment<bool> a2, a3; if (pt) a2 = { &pt->delta, false }; if (qs) a3 = { &qs->delta, false };

Next, we will update the reverse sampling densities of the connection vertices and their predecessors, starting with normal p Subscript t minus 1 . This vertex was originally sampled on the camera subpath, but it could also have been reached using an extra step from the light side ( normal q Subscript s minus 2 Baseline right-arrow normal q Subscript s minus 1 Baseline right-arrow normal p Subscript t minus 1 in three-point form); the resulting density at normal p Subscript t minus 1 is evaluated using Vertex::Pdf().

The case s equals 0 is special: here, normal p Subscript t is an intersection with a light source found on the camera subpath. The alternative reverse sampling strategy generates a light sample using Light::Sample_Le(), and we evaluate its spatial density with the help of Vertex::PdfLightOrigin().

<<Update reverse density of vertex normal p Subscript Baseline Subscript t minus 1 >>= 
ScopedAssignment<Float> a4; if (pt) a4 = { &pt->pdfRev, s > 0 ? qs->Pdf(scene, qsMinus, *pt) : pt->PdfLightOrigin(scene, *ptMinus, lightPdf) };

The next fragment initializes the pdfRev field of normal p Subscript t minus 2 with the density of the reverse strategy normal q Subscript s minus 1 Baseline right-arrow normal p Subscript t minus 1 Baseline right-arrow normal p Subscript t minus 2 . Once more, there is a special case for s equals 0 , where the alternative reverse strategy samples an emitted ray via Light::Sample_Le() and intersects it against the scene geometry; the corresponding density is evaluated using Vertex::PdfLight().

<<Update reverse density of vertex normal p Subscript Baseline Subscript t minus 2 >>= 
ScopedAssignment<Float> a5; if (ptMinus) a5 = { &ptMinus->pdfRev, s > 0 ? pt->Pdf(scene, qs, *ptMinus) : pt->PdfLight(scene, *ptMinus) };

The last fragment, <<Update reverse density of vertices normal q Subscript Baseline Subscript s minus 1 and normal q Subscript Baseline Subscript s minus 2 >>, is not included here. It is analogous except that it does not require a similar special case for t equals 0 .

16.3.5 Infinite Area Lights and BDPT

The infinite area light, first introduced in Section 12.6, provides a convenient way of illuminating scenes with realistic captured illumination. Unfortunately, its definition as an infinitely distant directional source turns out to be rather difficult to reconcile with BDPT’s path integral formulation, which expresses probabilities in terms of area densities, which in turn requires finite-sized emitting surfaces.

Through some gymnastics, we could represent infinite area lights with finite shapes. For example, an infinite area light’s radiance emission distribution could be described by a large emitting sphere that surrounded the scene, where the directional distribution of emitted radiance at each point on the interior of the sphere was the same as the infinite area light’s emitted radiance for the same direction. This approach would require a significantly more complex implementation of InfiniteAreaLight with no practical benefits apart from BDPT compatibility.

Instead of changing the functionality of InfiniteAreaLight, we will instead make infinite area lights a special case in BDPT. Since illumination from these lights is most naturally integrated over solid angles, our approach will be to add support for solid angle integration to the vertex abstraction layer. In practice, scenes may contain multiple infinite area lights; we will follow the convention of treating them as one combined light when evaluating the emitted radiance or determining sample probabilities.

First, we will create a special endpoint vertex any time a camera path ray escapes the scene. The <<Capture escaped rays when tracing from the camera>> fragment is invoked by RandomWalk() whenever no surface intersection could be found while generating the camera subpath. In the implementation of the Vertex::CreateLight() method called here, the pdfFwd variable recording the probability per unit solid angle is stored directly in Vertex::pdfFwd without conversion by Vertex::ConvertDensity().

<<Capture escaped rays when tracing from the camera>>= 
if (mode == TransportMode::Radiance) { vertex = Vertex::CreateLight(EndpointInteraction(ray), beta, pdfFwd); ++bounces; }

The existence of light vertices on the camera subpath leads to certain nonsensical connection strategies. For instance, we couldn’t possibly connect a light vertex on the camera subpath to another vertex on the light subpath. The following check in ConnectBDPT() detects and ignores such connection attempts.

<<Ignore invalid connections related to infinite area lights>>= 
if (t > 1 && s != 0 && cameraVertices[t - 1].type == VertexType::Light) return Spectrum(0.f);

Some parts of the code may still attempt to invoke ConvertDensity() with a next Vertex that refers to an infinite area light. The following fragment detects this case at the beginning of ConvertDensity() and directly returns the supplied solid angle density without conversion in that case.

<<Return solid angle density if next is an infinite area light>>= 
if (next.IsInfiniteLight()) return pdf;

Next, we need to adapt the light subpath sampling routine to correct the probability values returned by the ray sampling function InfiniteAreaLight::Sample_Le(). This case is detected in an additional fragment at the end of GenerateLightSubpath().

<<Correct subpath sampling densities for infinite area lights>>= 
if (path[0].IsInfiniteLight()) { <<Set spatial density of path[1] for infinite area light>> 
if (nVertices > 0) { path[1].pdfFwd = pdfPos; if (path[1].IsOnSurface()) path[1].pdfFwd *= AbsDot(ray.d, path[1].ng()); }
<<Set spatial density of path[0] for infinite area light>> 
path[0].pdfFwd = InfiniteLightDensity(scene, lightDistr, ray.d);
}

Recall that InfiniteAreaLight::Sample_Le() samples a ray direction (with a corresponding density pdfDir) and a ray origin on a perpendicular disk (with a corresponding density pdfPos) that touches the scene’s bounding sphere. Due to foreshortening, the resulting ray has a corresponding spatial density of monospace p monospace d monospace f monospace upper P monospace o monospace s StartAbsoluteValue cosine theta EndAbsoluteValue at its first intersection with the scene geometry, where theta is the angle between monospace r monospace a monospace y monospace period monospace d and the geometric normal.

<<Set spatial density of path[1] for infinite area light>>= 
if (nVertices > 0) { path[1].pdfFwd = pdfPos; if (path[1].IsOnSurface()) path[1].pdfFwd *= AbsDot(ray.d, path[1].ng()); }

Following our new convention, the spatial density of infinite area light endpoints is now expressed as a probability per unit solid angle. We will create a helper function InfiniteLightDensity() that determines this value while also accounting for the presence of other infinite area lights.

<<Set spatial density of path[0] for infinite area light>>= 
path[0].pdfFwd = InfiniteLightDensity(scene, lightDistr, ray.d);

This function performs a weighted sum of the directional densities of all infinite area lights using the light probabilities in lightDistr.

<<BDPT Helper Definitions>>+= 
inline Float InfiniteLightDensity(const Scene &scene, const Distribution1D &lightDistr, const Vector3f &w) { Float pdf = 0; for (size_t i = 0; i < scene.lights.size(); ++i) if (scene.lights[i]->flags & (int)LightFlags::Infinite) pdf += scene.lights[i]->Pdf_Li(Interaction(), -w) * lightDistr.func[i]; return pdf / (lightDistr.funcInt * lightDistr.Count()); }

The remaining two changes are analogous and address the probability computation in the PdfLightOrigin() and PdfLight() methods. For the former, we similarly return the combined solid angle density when an infinite area light is detected.

<<Return solid angle density for infinite light sources>>= 
return InfiniteLightDensity(scene, lightDistr, w);

In PdfLight(), we compute the probability of sampling a ray origin on a disk whose radius is equal to the scene’s bounding sphere. The remainder of PdfLight() already accounts for the necessary cosine foreshortening factor; hence we do not need to multiply by it here.

<<Compute planar sampling density for infinite light sources>>= 
Point3f worldCenter; Float worldRadius; scene.WorldBound().BoundingSphere(&worldCenter, &worldRadius); pdf = 1 / (Pi * worldRadius * worldRadius);