16.1 The Path-Space Measurement Equation

In light of the path integral form of the LTE from Equation (14.16), it’s useful to go back and formally describe the quantity that is being estimated when we compute pixel values for an image. Not only does this let us see how to apply the LTE to a wider set of problems than just computing 2D images (e.g., to precomputing scattered radiance distributions at the vertices of a polygonal model), but this process also leads us to a key theoretical mechanism for understanding the bidirectional path tracing and photon mapping algorithms in this chapter. For simplicity, we’ll use the basic path integral for surfaces rather than the generalized variant from Section 15.1.1, though the conclusions are applicable to both versions of the LTE.

The measurement equation describes the value of an abstract measurement that is found by integrating over some set of rays carrying radiance. For example, when computing the value of a pixel j in the image, we want to integrate over rays starting in the neighborhood of the pixel, with contributions weighted by the image reconstruction filter. Ignoring depth of field for now (so that each point on the film plane corresponds to a single outgoing direction from the camera), we can write the pixel’s value as an integral over points on the film plane of a weighting function times the incident radiance along the corresponding camera rays:

StartLayout 1st Row 1st Column upper I Subscript j 2nd Column equals integral Underscript upper A Subscript normal f normal i normal l normal m Baseline Endscripts integral Underscript script upper S squared Endscripts upper W Subscript normal e Superscript left-parenthesis j right-parenthesis Baseline left-parenthesis normal p Subscript normal f normal i normal l normal m Baseline comma omega Subscript Baseline right-parenthesis upper L Subscript normal i Superscript Baseline left-parenthesis normal p Subscript normal f normal i normal l normal m Baseline comma omega Subscript Baseline right-parenthesis StartAbsoluteValue cosine theta EndAbsoluteValue normal d omega Subscript Baseline normal d upper A Subscript Baseline left-parenthesis normal p Subscript normal f normal i normal l normal m Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column equals integral Underscript upper A Subscript normal f normal i normal l normal m Baseline Endscripts integral Underscript upper A Endscripts upper W Subscript normal e Superscript left-parenthesis j right-parenthesis Baseline left-parenthesis normal p 0 right-arrow normal p 1 right-parenthesis upper L left-parenthesis normal p 1 right-arrow normal p 0 right-parenthesis upper G left-parenthesis normal p 0 left-right-arrow normal p 1 right-parenthesis normal d upper A Subscript Baseline left-parenthesis normal p 1 right-parenthesis normal d upper A Subscript Baseline left-parenthesis normal p 0 right-parenthesis comma EndLayout

where upper I Subscript j is the measurement for the j th pixel and normal p 0 is a point on the film. In this setting, the upper W Subscript normal e Superscript left-parenthesis j right-parenthesis Baseline left-parenthesis normal p 0 right-arrow normal p 1 right-parenthesis term is the product of the filter function around the pixel f Subscript j and a delta function that selects the appropriate camera ray direction of the sample from normal p 0 , omega Subscript normal c normal a normal m normal e normal r normal a Baseline left-parenthesis normal p 1 right-parenthesis :

upper W Subscript normal e Superscript left-parenthesis j right-parenthesis Baseline left-parenthesis normal p 0 right-arrow normal p 1 right-parenthesis equals f Subscript j Baseline left-parenthesis normal p 0 right-parenthesis delta left-parenthesis t left-parenthesis normal p 0 comma omega Subscript normal c normal a normal m normal e normal r normal a Baseline left-parenthesis normal p 1 right-parenthesis right-parenthesis minus normal p 1 right-parenthesis period

This formulation may initially seem gratuitously complex, but it leads us to an important insight. If we expand the upper P left-parenthesis normal p Subscript Baseline overbar Subscript n Baseline right-parenthesis terms of the LTE sum, we have

StartLayout 1st Row 1st Column upper I Subscript j 2nd Column equals integral Underscript upper A Subscript normal f normal i normal l normal m Baseline Endscripts integral Underscript upper A Endscripts upper W Subscript normal e Superscript left-parenthesis j right-parenthesis Baseline left-parenthesis normal p 0 right-arrow normal p 1 right-parenthesis upper L left-parenthesis normal p 1 right-arrow normal p 0 right-parenthesis upper G left-parenthesis normal p 0 left-right-arrow normal p 1 right-parenthesis normal d upper A Subscript Baseline left-parenthesis normal p 1 right-parenthesis normal d upper A Subscript Baseline left-parenthesis normal p 0 right-parenthesis 2nd Row 1st Column Blank 2nd Column equals sigma-summation Underscript i Endscripts integral Underscript upper A Endscripts integral Underscript upper A Endscripts upper W Subscript normal e Superscript left-parenthesis j right-parenthesis Baseline left-parenthesis normal p 0 right-arrow normal p 1 right-parenthesis upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis upper G left-parenthesis normal p 0 left-right-arrow normal p 1 right-parenthesis normal d upper A Subscript Baseline left-parenthesis normal p 1 right-parenthesis normal d upper A Subscript Baseline left-parenthesis normal p 0 right-parenthesis 3rd Row 1st Column Blank 2nd Column equals sigma-summation Underscript i Endscripts ModifyingBelow integral Underscript upper A Endscripts midline-horizontal-ellipsis integral Underscript upper A Endscripts With bottom-brace Underscript i plus 1 times Endscripts upper W Subscript normal e Superscript left-parenthesis j right-parenthesis Baseline left-parenthesis normal p 0 right-arrow normal p 1 right-parenthesis upper T left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis upper L Subscript normal e Baseline left-parenthesis normal p Subscript i plus 1 Baseline right-arrow normal p Subscript i Baseline right-parenthesis upper G left-parenthesis normal p 0 left-right-arrow normal p 1 right-parenthesis 4th Row 1st Column Blank 2nd Column normal d upper A Subscript Baseline left-parenthesis normal p Subscript i plus 1 Baseline right-parenthesis midline-horizontal-ellipsis normal d upper A Subscript Baseline left-parenthesis normal p 0 right-parenthesis comma EndLayout

where upper T left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis is the path throughput function introduced in Equation (14.18).

Note the nice symmetric way in which the emitted radiance upper L Subscript normal e Superscript (quantifying the light source emission profile) and the weighting function upper W Subscript normal e Superscript Baseline Superscript left-parenthesis j right-parenthesis (quantifying the camera’s sensitivity profile for a pixel j ) appear in the above equation: neither term is treated specially, and from this we can infer that the concepts of emission and measurement are mathematically interchangeable.

The implications of this symmetry are important: it says that we can think of the rendering process in two different ways. The first interpretation is that light could be emitted from light sources, bounce around the scene, and arrive at a sensor where upper W Subscript normal e Superscript describes its contribution to the measurement. Alternatively, we can think of the sensor as emitting an imaginary quantity that creates a measurement when it reaches a light source. This idea isn’t just a theoretical construct: it can be applied in practice. A good example is the Dual Photography work of Sen et al. (2005) that showed that it was possible to take photographs from the viewpoint of a video projector by processing input photographs taken with a separate camera—this could be interpreted as turning the projector into a camera, while using the original camera as the “light source” to illuminate the scene.

By simply swapping the role of cameras and light sources in this way, we can create a method known as particle tracing, which traces rays from the light sources to recursively estimate the incident importance arriving on surfaces. This is not a particularly useful rendering technique on its own, but it constitutes an essential ingredient of other methods such as bidirectional path tracing and photon mapping.

The value described by the upper W Subscript normal e Superscript term is known as the importance for the ray between normal p 0 and normal p 1 in the scene. When the measurement equation is used to compute pixel measurements, the importance will often be partially or fully described by delta distributions, as it was in the previous example. Many other types of measurements besides image formation can be described by appropriately constructed importance functions, and thus the formalisms described here can be used to show how the integral over paths described by the measurement equation is also the integral that must be estimated to compute them.

16.1.1 Sampling Cameras

Bidirectional light transport algorithms require the ability to evaluate the value of the importance function for arbitrary points in the scene; this is useful, for example, for computing the importance of a point along a path that started at a light source. The Camera::We() method takes a ray with origin normal p and direction omega and evaluates the importance emitted from the point on the camera normal p Subscript in a direction omega Subscript . If provided, the pRaster2 parameter is used to return the raster position associated with the ray on the film; conceptually this can be understood as the discrete index j such that upper W Subscript normal e Superscript left-parenthesis j right-parenthesis Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis attains its maximum value; in practice, the function returns fractional values to specify the raster positions more accurately.

<<Camera Interface>>+=  
virtual Spectrum We(const Ray &ray, Point2f *pRaster2 = nullptr) const;

The default implementation of this method generates an error message; it’s currently only implemented for the perspective camera model in pbrt. Implementing it for other Camera models is saved for Exercise 16.1 at the end of the chapter.

<<PerspectiveCamera Method Definitions>>+=  
Spectrum PerspectiveCamera::We(const Ray &ray, Point2f *pRaster2) const { <<Interpolate camera matrix and check if omega Subscript is forward-facing>> 
Transform c2w; CameraToWorld.Interpolate(ray.time, &c2w); Float cosTheta = Dot(ray.d, c2w(Vector3f(0, 0, 1))); if (cosTheta <= 0) return 0;
<<Map ray left-parenthesis p comma omega Subscript Baseline right-parenthesis onto the raster grid>> 
Point3f pFocus = ray((lensRadius > 0 ? focalDistance : 1) / cosTheta); Point3f pRaster = Inverse(RasterToCamera)(Inverse(c2w)(pFocus));
<<Return raster position if requested>> 
if (pRaster2) *pRaster2 = Point2f(pRaster.x, pRaster.y);
<<Return zero importance for out of bounds points>> 
Bounds2i sampleBounds = film->GetSampleBounds(); if (pRaster.x < sampleBounds.pMin.x || pRaster.x >= sampleBounds.pMax.x || pRaster.y < sampleBounds.pMin.y || pRaster.y >= sampleBounds.pMax.y) return 0;
<<Compute lens area of perspective camera>> 
Float lensArea = lensRadius != 0 ? (Pi * lensRadius * lensRadius) : 1;
<<Return importance for point on image plane>> 
Float cos2Theta = cosTheta * cosTheta; return Spectrum(1 / (A * lensArea * cos2Theta * cos2Theta));
}

Given the camera-to-world transformation for the provided time, the method checks that the direction omega Subscript points in the same hemisphere as the camera is facing by transforming the camera-space viewing direction left-parenthesis 0 comma 0 comma 1 right-parenthesis to world space and computing the cosine of the angles between them. If these directions are more than 90 degrees apart, then the camera would never return a ray in this direction from its GenerateRay() method, and so an importance value of 0 can be returned immediately.

<<Interpolate camera matrix and check if omega Subscript is forward-facing>>= 
Transform c2w; CameraToWorld.Interpolate(ray.time, &c2w); Float cosTheta = Dot(ray.d, c2w(Vector3f(0, 0, 1))); if (cosTheta <= 0) return 0;

A slightly more involved test next checks if the ray corresponds to one starting from the film area. If its origin is outside of the film’s extent, then the point normal p Subscript is outside of the camera’s viewing volume, and again a zero importance value can be returned.

For a camera with a finite aperture, we have a point on the lens and its direction (Figure 16.1). We don’t yet know the point on the film that this ray corresponds to, but we do know that all rays leaving that point are in focus at the plane z equals normal f normal o normal c normal a normal l normal upper D normal i normal s normal t normal a normal n normal c normal e . Therefore, if we compute the ray’s intersection with the plane of focus, then transforming that point with the perspective projection matrix gives us the corresponding point on the film. For pinhole apertures, we compute the intersection with a plane arbitrarily set at z equals 1 to get a point along the ray leaving the camera before performing the projection.

Figure 16.1: Computing the Focus Point pFocus for a Ray Leaving the Lens. In order to compute the value of the importance function upper W Subscript normal e Superscript for this ray, we need to find the point it corresponds to on the film plane. To do so, we first compute pFocus, the point where the ray intersects the plane of focus. This point can in turn be projected by the camera’s perspective projection matrix to find the corresponding raster-space point on the film.

<<Map ray left-parenthesis p comma omega Subscript Baseline right-parenthesis onto the raster grid>>= 
Point3f pFocus = ray((lensRadius > 0 ? focalDistance : 1) / cosTheta); Point3f pRaster = Inverse(RasterToCamera)(Inverse(c2w)(pFocus));

<<Return raster position if requested>>= 
if (pRaster2) *pRaster2 = Point2f(pRaster.x, pRaster.y);

Given the raster-space point, it’s easy to check if it’s inside the image extent.

<<Return zero importance for out of bounds points>>= 
Bounds2i sampleBounds = film->GetSampleBounds(); if (pRaster.x < sampleBounds.pMin.x || pRaster.x >= sampleBounds.pMax.x || pRaster.y < sampleBounds.pMin.y || pRaster.y >= sampleBounds.pMax.y) return 0;

The perspective camera in pbrt is an ideal sensor in the sense that it generates samples with a uniform distribution over the film area. We will now use this fact to derive the corresponding directional sampling distribution. We’ll start by defining a camera space image rectangle that all camera rays pass through and (arbitrarily) choose the one on the plane z equals 1 . The following fragment that is part of the PerspectiveCamera constructor uses the RasterToCamera transformation and divides by the z coordinate to compute the rectangle’s corner points, which in turn gives the rectangle’s area upper A .

<<Compute image plane bounds at z equals 1 for PerspectiveCamera>>= 
Point2i res = film->fullResolution; Point3f pMin = RasterToCamera(Point3f(0, 0, 0)); Point3f pMax = RasterToCamera(Point3f(res.x, res.y, 0)); pMin /= pMin.z; pMax /= pMax.z; A = std::abs((pMax.x - pMin.x) * (pMax.y - pMin.y));

<<PerspectiveCamera Private Data>>+= 
Float A;

Figure 16.2: Deriving the Normalized Importance Function for the PerspectiveCamera. Given a point on the visible region of the image plane at z equals 1 with PDF p left-parenthesis normal p Subscript Baseline right-parenthesis equals 1 slash upper A , we can compute the directional PDF at a point on the lens (filled circle) by applying Equation (5.6) to account for the distance to the point on the lens d and theta , the angle between the vector from normal p Subscript to the point on the lens, and the surface normal of the image plane.

The importance function doesn’t have to obey any normalization constraints (just like emitted radiance from an area light). However, we’ll define the PerspectiveCamera’s importance function as a normalized PDF on ray space, both for convenience in the following and also to be consistent with the weight values of 1 returned from PerspectiveCamera::GenerateRay().

The importance function of PerspectiveCamera varies smoothly over its support ( upper W Subscript normal e Superscript Baseline left-parenthesis omega Subscript Baseline right-parenthesis greater-than 0 ); this variation is defined so that it cancels the vignetting that a real pinhole camera would have and ensures that pixels record values in units of radiance (this is another reason why PerspectiveCamera::GenerateRay() returns weight values of 1).

The camera uniformly generates samples over the image plane area upper A ; thus, the area-measured PDF for points on the image plane is p left-parenthesis normal p Subscript Baseline right-parenthesis equals 1 slash upper A . Consider now the directional PDF at a point on the lens (or the camera’s pinhole) corresponding to a differential area on the image plane (Figure 16.2); applying Equation (5.6) to transform to a directional density gives

p left-parenthesis omega right-parenthesis equals StartLayout Enlarged left-brace 1st Row 1st Column StartFraction d squared Over upper A cosine theta EndFraction comma 2nd Column if omega is within the frustum 2nd Row 1st Column 0 comma 2nd Column otherwise comma EndLayout

where theta is the angle that omega makes with the image rectangle normal and d is the distance between the point on the lens and the intersection of the ray with the z equals 1 plane. The distance to the point on the image plane is

d equals double-vertical-bar StartFraction omega Over cosine theta EndFraction double-vertical-bar equals StartFraction 1 Over cosine theta EndFraction comma

as cosine theta is the z coordinate of omega in the local camera coordinate system and double-vertical-bar omega double-vertical-bar equals 1 . Hence,

p left-parenthesis omega right-parenthesis equals StartLayout Enlarged left-brace 1st Row 1st Column StartFraction 1 Over upper A cosine cubed theta EndFraction comma 2nd Column if omega is within the frustum 2nd Row 1st Column 0 comma 2nd Column otherwise period EndLayout

By construction, the above density function p left-parenthesis omega right-parenthesis is normalized when integrated over directions—however, we initially set out to create a normalized importance function upper W Subscript normal e Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis that is defined on the space of camera rays upper A Subscript normal l normal e normal n normal s Baseline times script upper S squared , where upper A Subscript normal l normal e normal n normal s is the surface region associated with the perspective camera’s lens element. This function must satisfy the ray-space normalization criterion

integral Underscript upper A Subscript normal l normal e normal n normal s Baseline Endscripts integral Underscript script upper S squared Endscripts upper W Subscript normal e Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis StartAbsoluteValue cosine theta EndAbsoluteValue normal d omega Subscript Baseline normal d upper A Subscript Baseline left-parenthesis normal p Subscript Baseline right-parenthesis equals 1 period

We can’t directly set upper W Subscript normal e Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis equal to p left-parenthesis omega right-parenthesis due to the extra integration over areas and the additional cosine factor in the above integral.

Note that the lens area upper A Subscript normal l normal e normal n normal s of a perspective camera is equal to pi r squared , where r is the lens radius. For point camera, the lens area is set to 1 and interpreted as a Dirac delta function.

<<Compute lens area of perspective camera>>= 
Float lensArea = lensRadius != 0 ? (Pi * lensRadius * lensRadius) : 1;

We then define upper W Subscript normal e Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis on ray space as

upper W Subscript normal e Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript Baseline right-parenthesis equals StartFraction p left-parenthesis omega Subscript Baseline right-parenthesis Over pi r squared cosine theta EndFraction equals StartLayout Enlarged left-brace 1st Row 1st Column StartFraction 1 Over upper A pi r squared cosine Superscript 4 Baseline theta EndFraction comma 2nd Column if omega is within the frustum 2nd Row 1st Column 0 comma 2nd Column otherwise comma EndLayout

which divides p by the lens area and a term to cancel out the cosine factor from Equation (16.3). At this point, the implementation has already ensured that omega is within the frustum, hence the only thing left to do is to return the importance value according to the first case of Equation (16.4).

<<Return importance for point on image plane>>= 
Float cos2Theta = cosTheta * cosTheta; return Spectrum(1 / (A * lensArea * cos2Theta * cos2Theta));

In the light of a definition for the perspective camera’s importance function, we can now reinterpret the ray generation function Camera::GenerateRay() as an importance sampling technique for We(). As such, it’s appropriate to define a Camera method that separately returns the spatial and directional PDFs for sampling a particular ray leaving the camera, analogous to the Light::Pdf_Le() method for light sources that will be introduced in Section 16.1.2. As before, this method is currently only implemented for the PerspectiveCamera and the default implementation generates an error message.

<<Camera Interface>>+=  
virtual void Pdf_We(const Ray &ray, Float *pdfPos, Float *pdfDir) const;

The directional density of the ideal sampling strategy implemented in the method PerspectiveCamera::GenerateRay() was already discussed and is equal to p left-parenthesis omega Subscript Baseline right-parenthesis defined in Equation (16.2). The spatial density is the reciprocal of the lens area. Due to this overlap, the first four fragments of PerspectiveCamera::Pdf_We() are therefore either identical or almost identical to the similarly named fragments from PerspectiveCamera::We() with the exception that they return zero probabilities via *pdfPos and *pdfDir upon failure.

<<PerspectiveCamera Method Definitions>>+=  
void PerspectiveCamera::Pdf_We(const Ray &ray, Float *pdfPos, Float *pdfDir) const { <<Interpolate camera matrix and fail if omega Subscript is not forward-facing>> 
Transform c2w; CameraToWorld.Interpolate(ray.time, &c2w); Float cosTheta = Dot(ray.d, c2w(Vector3f(0, 0, 1))); if (cosTheta <= 0) { *pdfPos = *pdfDir = 0; return; }
<<Map ray left-parenthesis p comma omega Subscript Baseline right-parenthesis onto the raster grid>> 
Point3f pFocus = ray((lensRadius > 0 ? focalDistance : 1) / cosTheta); Point3f pRaster = Inverse(RasterToCamera)(Inverse(c2w)(pFocus));
<<Return zero probability for out of bounds points>> 
Bounds2i sampleBounds = film->GetSampleBounds(); if (pRaster.x < sampleBounds.pMin.x || pRaster.x >= sampleBounds.pMax.x || pRaster.y < sampleBounds.pMin.y || pRaster.y >= sampleBounds.pMax.y) { *pdfPos = *pdfDir = 0; return; }
<<Compute lens area of perspective camera>> 
Float lensArea = lensRadius != 0 ? (Pi * lensRadius * lensRadius) : 1;
*pdfPos = 1 / lensArea; *pdfDir = 1 / (A * cosTheta * cosTheta * cosTheta); }

One last additional Camera method completes the symmetry between light sources and cameras: it samples a point on the camera lens and computes an incident direction along which importance is arriving at the given reference position in the scene; it is thus the camera equivalent of Light::Sample_Li().

Like Sample_Li(), the PDF value this method returns is defined with respect to solid angle at the reference point.

<<Camera Interface>>+= 
virtual Spectrum Sample_Wi(const Interaction &ref, const Point2f &u, Vector3f *wi, Float *pdf, Point2f *pRaster, VisibilityTester *vis) const;

The PerspectiveCamera implementation of this method samples a point on the lens to compute the incident importance at the reference point.

<<PerspectiveCamera Method Definitions>>+= 
Spectrum PerspectiveCamera::Sample_Wi(const Interaction &ref, const Point2f &u, Vector3f *wi, Float *pdf, Point2f *pRaster, VisibilityTester *vis) const { <<Uniformly sample a lens interaction lensIntr>> 
Point2f pLens = lensRadius * ConcentricSampleDisk(u); Point3f pLensWorld = CameraToWorld(ref.time, Point3f(pLens.x, pLens.y, 0)); Interaction lensIntr(pLensWorld, ref.time, medium); lensIntr.n = Normal3f(CameraToWorld(ref.time, Vector3f(0, 0, 1)));
<<Populate arguments and compute the importance value>> 
*vis = VisibilityTester(ref, lensIntr); *wi = lensIntr.p - ref.p; Float dist = wi->Length(); *wi /= dist; <<Compute PDF for importance arriving at ref>> 
<<Compute lens area of perspective camera>> 
Float lensArea = lensRadius != 0 ? (Pi * lensRadius * lensRadius) : 1;
*pdf = (dist * dist) / (AbsDot(lensIntr.n, *wi) * lensArea);
return We(lensIntr.SpawnRay(-*wi), pRaster);
}

We can compute an Interaction for a point on the lens by using u to sample the lens and then transforming this point to world space. For pinhole cameras, lensRadius is 0 and pLens always ends up being at the origin.

<<Uniformly sample a lens interaction lensIntr>>= 
Point2f pLens = lensRadius * ConcentricSampleDisk(u); Point3f pLensWorld = CameraToWorld(ref.time, Point3f(pLens.x, pLens.y, 0)); Interaction lensIntr(pLensWorld, ref.time, medium); lensIntr.n = Normal3f(CameraToWorld(ref.time, Vector3f(0, 0, 1)));

Given the point on the lens, most of the output parameters of Sample_We() can be initialized in a straightforward manner.

<<Populate arguments and compute the importance value>>= 
*vis = VisibilityTester(ref, lensIntr); *wi = lensIntr.p - ref.p; Float dist = wi->Length(); *wi /= dist; <<Compute PDF for importance arriving at ref>> 
<<Compute lens area of perspective camera>> 
Float lensArea = lensRadius != 0 ? (Pi * lensRadius * lensRadius) : 1;
*pdf = (dist * dist) / (AbsDot(lensIntr.n, *wi) * lensArea);
return We(lensIntr.SpawnRay(-*wi), pRaster);

The PDF of the sample is the probability of sampling a point on the lens (1 / lensArea), converted into a probability per unit solid angle at the reference point. For pinhole cameras, there is an implied delta distribution in both the PDF and the importance function value that cancel out later. (This is following the same convention as was used for BSDFs and light sources in Sections 14.1.3 and 14.2.1).

<<Compute PDF for importance arriving at ref>>= 
<<Compute lens area of perspective camera>> 
Float lensArea = lensRadius != 0 ? (Pi * lensRadius * lensRadius) : 1;
*pdf = (dist * dist) / (AbsDot(lensIntr.n, *wi) * lensArea);

16.1.2 Sampling Light Rays

For bidirectional light transport algorithms, it’s also necessary to add a light sampling method, Sample_Le(), that samples a ray from a distribution of rays leaving the light, returning the ray in *ray and the surface normal at the point on the light source in *nLight (effectively, the analog of Camera::GenerateRay()). A total of four sample values are passed to this method in the u1 and u2 parameters so that two are available to sample the ray’s origin and two are available for its direction. Not all light implementations need all of these values—for example, the origin of all rays leaving a point light is the same.

This method returns two PDF values: the ray origin’s probability density with respect to surface area on the light and its direction’s probability density with respect to solid angle. The joint probability of sampling the ray is the product of these two probabilities.

<<Light Interface>>+=  
virtual Spectrum Sample_Le(const Point2f &u1, const Point2f &u2, Float time, Ray *ray, Normal3f *nLight, Float *pdfPos, Float *pdfDir) const = 0;

So that multiple importance sampling can be applied, there is also a method to return the position and direction PDFs for a given ray.

<<Light Interface>>+= 
virtual void Pdf_Le(const Ray &ray, const Normal3f &nLight, Float *pdfPos, Float *pdfDir) const = 0;

Point Lights

The sampling method for generating rays leaving point lights is straightforward. The origin of the ray must be the light’s position; this part of the density is described by a delta distribution. Directions are uniformly sampled over the sphere, and the overall sampling density is the product of these two densities. As usual, we’ll ignore the delta distribution that is in the actual PDF because it is canceled out by a (missing) corresponding delta term in the radiance value in the Spectrum returned by the sampling routine.

<<PointLight Method Definitions>>+=  
Spectrum PointLight::Sample_Le(const Point2f &u1, const Point2f &u2, Float time, Ray *ray, Normal3f *nLight, Float *pdfPos, Float *pdfDir) const { *ray = Ray(pLight, UniformSampleSphere(u1), Infinity, time, mediumInterface.inside); *nLight = (Normal3f)ray->d; *pdfPos = 1; *pdfDir = UniformSpherePdf(); return I; }

<<PointLight Method Definitions>>+= 
void PointLight::Pdf_Le(const Ray &, const Normal3f &, Float *pdfPos, Float *pdfDir) const { *pdfPos = 0; *pdfDir = UniformSpherePdf(); }

Spotlights

The method for sampling an outgoing ray with a reasonable distribution for the spotlight is more interesting. While it could just sample directions uniformly on the sphere as was done for the point light, this distribution is likely to be a bad match for the spotlight’s actual distribution. For example, if the light has a very narrow beam angle, many samples will be taken in directions where the light doesn’t cast any illumination. Instead, we will sample from a uniform distribution over the cone of directions in which the light casts illumination. Although the sampling distribution does not try to account for the falloff toward the edges of the beam, this is only a minor shortcoming in practice.

The PDF p left-parenthesis theta comma phi right-parenthesis for the spotlight’s illumination distribution is separable with p left-parenthesis phi right-parenthesis equals 1 slash left-parenthesis 2 pi right-parenthesis . We thus just need to find a sampling distribution for theta . The UniformSampleCone() function from Section 13.6.4 provides this functionality.

<<SpotLight Method Definitions>>+=  
Spectrum SpotLight::Sample_Le(const Point2f &u1, const Point2f &u2, Float time, Ray *ray, Normal3f *nLight, Float *pdfPos, Float *pdfDir) const { Vector3f w = UniformSampleCone(u1, cosTotalWidth); *ray = Ray(pLight, LightToWorld(w), Infinity, time, mediumInterface.inside); *nLight = (Normal3f)ray->d; *pdfPos = 1; *pdfDir = UniformConePdf(cosTotalWidth); return I * Falloff(ray->d); }

The SpotLight’s Pdf_Le() method for sampled rays must check to see if the direction is inside the cone of illuminated directions before returning the cone sampling PDF.

<<SpotLight Method Definitions>>+= 
void SpotLight::Pdf_Le(const Ray &ray, const Normal3f &, Float *pdfPos, Float *pdfDir) const { *pdfPos = 0; *pdfDir = (CosTheta(WorldToLight(ray.d)) >= cosTotalWidth) ? UniformConePdf(cosTotalWidth) : 0; }

The sampling routines for ProjectionLights and GonioPhotometricLights are essentially the same as the ones for SpotLights and PointLights, respectively. For sampling outgoing rays, ProjectionLights sample uniformly from the cone that encompasses their projected image map (hence the need to compute ProjectionLight::cosTotalWidth in the constructor), and those for GonioPhotometricLights sample uniformly over the unit sphere. Exercise 16.2 at the end of this chapter discusses improvements to these sampling methods that better account for the directional variation of these lights.

Area Lights

The method for sampling a ray leaving an area light is also easily implemented in terms of the shape sampling methods from Section 14.2.2.

<<DiffuseAreaLight Method Definitions>>+= 
Spectrum DiffuseAreaLight::Sample_Le(const Point2f &u1, const Point2f &u2, Float time, Ray *ray, Normal3f *nLight, Float *pdfPos, Float *pdfDir) const { <<Sample a point on the area light’s Shape, pShape>> 
Interaction pShape = shape->Sample(u1); pShape.mediumInterface = mediumInterface; *pdfPos = shape->Pdf(pShape); *nLight = pShape.n;
<<Sample a cosine-weighted outgoing direction w for area light>> 
Vector3f w = CosineSampleHemisphere(u2); *pdfDir = CosineHemispherePdf(w.z); Vector3f v1, v2, n(pShape.n); CoordinateSystem(n, &v1, &v2); w = w.x * v1 + w.y * v2 + w.z * n;
*ray = pShape.SpawnRay(w); return L(pShape, w); }

The surface area–based variant of Shape::Sample() is used to find the ray origin, sampled from some density over the surface.

<<Sample a point on the area light’s Shape, pShape>>= 
Interaction pShape = shape->Sample(u1); pShape.mediumInterface = mediumInterface; *pdfPos = shape->Pdf(pShape); *nLight = pShape.n;

The ray’s direction is sampled from a cosine-weighted distribution about the surface normal at the sampled point. Incorporating this cosine weighting means that rays leaving the light carry uniform differential power, which is preferable for bidirectional light transport algorithms. Because the direction returned by CosineSampleHemisphere() is in the canonical coordinate system, it must be transformed to the coordinate system about the surface normal at the sampled point here.

<<Sample a cosine-weighted outgoing direction w for area light>>= 
Vector3f w = CosineSampleHemisphere(u2); *pdfDir = CosineHemispherePdf(w.z); Vector3f v1, v2, n(pShape.n); CoordinateSystem(n, &v1, &v2); w = w.x * v1 + w.y * v2 + w.z * n;

Distant Lights

Sampling a ray from the DistantLight’s distribution of outgoing rays is a more interesting problem. The ray’s direction is determined in advance by a delta distribution; it must be the same as the light’s negated direction. For its origin, there are an infinite number of 3D points where it could start. How should we choose an appropriate one, and how do we compute its density?

The desired property is that rays intersect points in the scene that are illuminated by the distant light with uniform probability. One way to do this is to construct a disk that has the same radius as the scene’s bounding sphere and has a normal that is oriented with the light’s direction and then choose a random point on this disk, using the ConcentricSampleDisk() function (Figure 16.3). Once this point has been chosen, if the point is displaced along the light’s direction by the scene’s bounding sphere radius and used as the origin of the light ray, the ray origin will be outside the bounding sphere of the scene.

Figure 16.3: To sample an outgoing ray direction for a distant light source, the DistantLight::Sample_Le() method finds the disk oriented in the light’s direction that is large enough so that the entire scene can be intersected by rays leaving the disk in the light’s direction. Ray origins are sampled uniformly by area on this disk, and ray directions are given directly by the light’s direction.

This is a valid sampling approach, since by construction it has nonzero probability of sampling all incident rays into the sphere due to the directional light. The area component of the sampling density is uniform and therefore equal to the reciprocal of the area of the disk that was sampled. The directional density is given by a delta distribution based on the light’s direction.

<<DistantLight Method Definitions>>+= 
Spectrum DistantLight::Sample_Le(const Point2f &u1, const Point2f &u2, Float time, Ray *ray, Normal3f *nLight, Float *pdfPos, Float *pdfDir) const { <<Choose point on disk oriented toward infinite light direction>> 
Vector3f v1, v2; CoordinateSystem(wLight, &v1, &v2); Point2f cd = ConcentricSampleDisk(u1); Point3f pDisk = worldCenter + worldRadius * (cd.x * v1 + cd.y * v2);
<<Set ray origin and direction for infinite light ray>> 
*ray = Ray(pDisk + worldRadius * wLight, -wLight, Infinity, time);
*nLight = (Normal3f)ray->d; *pdfPos = 1 / (Pi * worldRadius * worldRadius); *pdfDir = 1; return L; }

Choosing the point on the oriented disk is a simple application of vector algebra. We construct a coordinate system with two vectors perpendicular to the disk’s normal (the light’s direction); see Figure 16.4. Given a random point on the canonical unit disk, computing the offsets from the disk’s center with respect to its coordinate vectors gives the corresponding point.

Figure 16.4: Given sample points left-parenthesis d 1 comma d 2 right-parenthesis on the canonical unit disk, points on an arbitrarily oriented and sized disk with normal bold n Subscript can be found by computing an arbitrary coordinate system left-parenthesis bold v 1 comma bold v 2 comma bold n Subscript Baseline right-parenthesis and then computing points on the disk with the offset d 1 bold v 1 plus d 2 bold v 2 from the disk’s center.

<<Choose point on disk oriented toward infinite light direction>>= 
Vector3f v1, v2; CoordinateSystem(wLight, &v1, &v2); Point2f cd = ConcentricSampleDisk(u1); Point3f pDisk = worldCenter + worldRadius * (cd.x * v1 + cd.y * v2);

Finally, the point is offset along the light direction and the ray can be initialized. Recall from Section 12.4 that DistantLights can’t be embedded in any medium other than a vacuum. Therefore, no medium needs to be specified for the ray.

<<Set ray origin and direction for infinite light ray>>= 
*ray = Ray(pDisk + worldRadius * wLight, -wLight, Infinity, time);

Infinite Area Lights

Generating a random ray leaving an infinite light source can be done by sampling a direction with the same approach as the earlier InfiniteAreaLight::Sample_Li() method. The sampled ray’s origin is then set using the same approach as was used for DistantLights, where a disk that covers the scene’s bounding sphere is oriented along the ray’s direction (recall Figure 16.3). We therefore won’t include the <<Compute direction for infinite light sample ray>> or <<Compute origin for infinite light sample ray>> fragments here.

<<InfiniteAreaLight Method Definitions>>+= 
Spectrum InfiniteAreaLight::Sample_Le(const Point2f &u1, const Point2f &u2, Float time, Ray *ray, Normal3f *nLight, Float *pdfPos, Float *pdfDir) const { <<Compute direction for infinite light sample ray>> 
Point2f u = u1; <<Find left-parenthesis u comma v right-parenthesis sample coordinates in infinite light texture>> 
Float mapPdf; Point2f uv = distribution->SampleContinuous(u, &mapPdf); if (mapPdf == 0) return Spectrum(0.f);
Float theta = uv[1] * Pi, phi = uv[0] * 2.f * Pi; Float cosTheta = std::cos(theta), sinTheta = std::sin(theta); Float sinPhi = std::sin(phi), cosPhi = std::cos(phi); Vector3f d = -LightToWorld(Vector3f(sinTheta * cosPhi, sinTheta * sinPhi, cosTheta)); *nLight = (Normal3f)d;
<<Compute origin for infinite light sample ray>> 
Vector3f v1, v2; CoordinateSystem(-d, &v1, &v2); Point2f cd = ConcentricSampleDisk(u2); Point3f pDisk = worldCenter + worldRadius * (cd.x * v1 + cd.y * v2); *ray = Ray(pDisk + worldRadius * -d, d, Infinity, time);
<<Compute InfiniteAreaLight ray PDFs>> 
*pdfDir = sinTheta == 0 ? 0 : mapPdf / (2 * Pi * Pi * sinTheta); *pdfPos = 1 / (Pi * worldRadius * worldRadius);
return Spectrum(Lmap->Lookup(uv), SpectrumType::Illuminant); }

The PDFs for these rays are the PDF for sampling the direction (as derived in Section 14.2.4) and the PDF for sampling a point on the disk.

<<Compute InfiniteAreaLight ray PDFs>>= 
*pdfDir = sinTheta == 0 ? 0 : mapPdf / (2 * Pi * Pi * sinTheta); *pdfPos = 1 / (Pi * worldRadius * worldRadius);

The Pdf_Le() method applies the same formulas, so we won’t include its implementation here.

16.1.3 Non-symmetric Scattering

Certain aspects in the input scene specification of materials and geometry can lead to non-symmetric behavior in light transport simulations, where incident radiance and importance are scattered in different ways at a point. If these differences aren’t accounted for, rendering algorithms based on radiance and importance transport will produce different and inconsistent results when rendering the same input scene. Bidirectional techniques that combine radiance and importance transport are particularly affected, since their design is fundamentally based on the principle of symmetry.

In this section, we will briefly enumerate cases that give rise to non-symmetry and explain how they can be addressed to arrive at a consistent set of bidirectional estimators.

Recall the path throughput term upper T left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis from Equation (16.1), which was defined as

upper T left-parenthesis normal p Subscript Baseline overbar Subscript n Baseline right-parenthesis equals product Underscript i equals 1 Overscript n minus 1 Endscripts f Subscript Baseline left-parenthesis normal p Subscript Baseline Subscript i plus 1 Baseline right-arrow normal p Subscript Baseline Subscript i Baseline right-arrow normal p Subscript Baseline Subscript i minus 1 Baseline right-parenthesis upper G left-parenthesis normal p Subscript Baseline Subscript i plus 1 Baseline left-right-arrow normal p Subscript Baseline Subscript i Baseline right-parenthesis period

The vertices are ordered such that normal p Subscript Baseline Subscript i denotes the i -th scattering event as seen from the camera.

Sampling techniques based on finding importance-carrying paths trace rays starting at the light sources to estimate the incident importance at the light, which means that the vertices will be generated in reverse compared to the above ordering. As such, the incident and outgoing direction arguments of the BSDFs will be (incorrectly) reversed unless special precautions are taken. We thus define the adjoint BSDF f Superscript asterisk at vertex normal p Subscript Baseline Subscript i , whose only role is to evaluate the original BSDF with swapped arguments:

f Superscript asterisk Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline comma omega Subscript normal i Baseline right-parenthesis equals f left-parenthesis normal p Subscript Baseline comma omega Subscript normal i Baseline comma omega Subscript normal o Baseline right-parenthesis period

All sampling steps based on importance transport will then use the adjoint form of the BSDF rather than its original version. Most BSDFs in pbrt are symmetric so that there is no actual difference between f and f Superscript asterisk . However, certain cases related to shading normals and light refracting into media with a different index of refraction require additional attention.

The TransportMode enumeration is used to inform such non-symmetric BSDFs about the transported quantity so that they can correctly switch between the adjoint and non-adjoint forms.

<<TransportMode Declarations>>= 
enum class TransportMode { Radiance, Importance };

Non-symmetry Due to Refraction

When light refracts into a material with a higher index of refraction than the incident medium’s index of refraction, the energy is compressed into a smaller set of angles. This is easy to see yourself, for instance, by looking at the sky from underwater in a quiet outdoor swimming pool. Because no light can be refracted above the critical angle ( tilde 48.6 Superscript ring for water), the incident hemisphere of light is squeezed into a considerably smaller subset of the hemisphere, which covers the remaining set of angles. Radiance along rays that do refract must thus increase so that energy is preserved when light passes through the interface. More precisely, the incident left-parenthesis upper L Subscript normal i Baseline right-parenthesis and transmitted left-parenthesis upper L Subscript normal t Baseline right-parenthesis radiance are related by

upper L Subscript normal i Baseline equals StartFraction eta Subscript normal i Superscript 2 Baseline Over eta Subscript normal t Superscript 2 Baseline EndFraction upper L Subscript normal t Baseline comma

where eta Subscript normal i and eta Subscript normal t are the refractive indices on the incident and transmitted sides, respectively. The symmetry relationship satisfied by a BTDF is

eta Subscript normal t Superscript 2 Baseline f left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline comma omega Subscript normal i Baseline right-parenthesis equals eta Subscript normal i Superscript 2 Baseline f left-parenthesis normal p Subscript Baseline comma omega Subscript normal i Baseline comma omega Subscript normal o Baseline right-parenthesis comma

and we can obtain the adjoint BTDF

f Superscript asterisk Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline comma omega Subscript normal i Baseline right-parenthesis equals f left-parenthesis normal p Subscript Baseline comma omega Subscript normal i Baseline comma omega Subscript normal o Baseline right-parenthesis equals StartFraction eta Subscript normal t Superscript 2 Baseline Over eta Subscript normal i Superscript 2 Baseline EndFraction f left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline comma omega Subscript normal i Baseline right-parenthesis comma

which effectively cancels out the scale factor in Equation (16.5). With these equations, we can now define the last missing piece in the implementation of SpecularTransmission::Sample_f(). Whenever radiance is transported over a refractive boundary, we apply the scale factor from Equation (16.5). For importance transport, we use the adjoint BTDF, which lacks the scaling factor due to the combination of Equations (16.5) and  (16.6).

<<Account for non-symmetry with transmission to different medium>>= 
if (mode == TransportMode::Radiance) ft *= (etaI * etaI) / (etaT * etaT);

A similar adjustment is also needed for the FourierBSDF::f() method in the case of refraction. In this case, FourierBSDFTable::eta provides the relative index of refraction. Recall that this model uses a convention where the sign of mu Subscript normal i Baseline equals cosine theta Subscript normal i is flipped, hence the expression muI * muO > 0 can be used to check if light is being refracted rather than reflected.

<<Update scale to account for adjoint light transport>>= 
if (mode == TransportMode::Radiance && muI * muO > 0) { float eta = muI > 0 ? 1 / bsdfTable.eta : bsdfTable.eta; scale *= eta * eta; }

Finally, the transmissive term upper S Subscript omega of the SeparableBSSRDF requires a similar correction when light leaves the medium after a second refraction (the first one being handled by the material’s BSDF).

<<Update BSSRDF transmission term to account for adjoint light transport>>= 
if (bssrdf->mode == TransportMode::Radiance) f *= bssrdf->eta * bssrdf->eta;

Non-symmetry Due to Shading Normals

Shading normals are another cause of non-symmetric scattering. As previously discussed in Section 3.6.3, shading normals are mainly used to make polygonal surfaces appear smoother than their actual discretization. This entails replacing the “true” geometric normal bold n Subscript normal g with an interpolated shading normal bold n Subscript normal s whenever the BSDF or the cosine term in the light transport equation are evaluated. Bump or normal mapping can be interpreted as another kind of shading normal, where bold n Subscript normal s is obtained from a texture map.

This kind of modification to the normal of a surface interaction causes a corresponding change in the underlying reflectance model, producing an effective BSDF that is generally non-symmetric. Without additional precautions, this non-symmetry can lead to visible artifacts in renderings based on adjoint techniques, including discontinuities in shading effects resembling flat-shaded polygons that interpolated normals were originally meant to avoid.

Recall the light transport equation, (14.13), which relates the incident and outgoing radiance on surfaces:

StartLayout 1st Row 1st Column upper L Subscript normal o Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline right-parenthesis 2nd Column equals upper L Subscript normal e Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline right-parenthesis plus integral Underscript script upper S squared Endscripts f Subscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline comma omega Subscript normal i Baseline right-parenthesis upper L Subscript normal i Superscript Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal i Baseline right-parenthesis StartAbsoluteValue bold n Subscript normal g Baseline dot omega Subscript normal i Baseline EndAbsoluteValue normal d omega Subscript normal i Baseline period EndLayout

Here, the cosine factor is expressed as an inner product involving omega Subscript normal i and the true normal of the underlying geometry. Suppose now that we’d like to replace bold n Subscript normal g with the shading normal bold n Subscript normal s . Instead of modifying the scattering equation, another mathematically equivalent way of expressing this change entails switching to a new BSDF f Subscript shade defined as

f Subscript shade Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline comma omega Subscript normal i Baseline right-parenthesis equals StartFraction StartAbsoluteValue bold n Subscript normal s Baseline dot omega Subscript normal i Baseline EndAbsoluteValue Over StartAbsoluteValue bold n Subscript normal g Baseline dot omega Subscript normal i Baseline EndAbsoluteValue EndFraction f left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline comma omega Subscript normal i Baseline right-parenthesis period

The first factor in the above expression makes this BSDF non-symmetric with respect to the arguments omega Subscript normal i and omega Subscript normal o . To avoid artifacts and inconsistencies in bidirectional rendering algorithms, the adjoint BSDF f Subscript shade Superscript asterisk should be used in simulations whenever importance transport is used. It is given by

f Subscript shade Superscript asterisk Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline comma omega Subscript normal i Baseline right-parenthesis equals StartFraction StartAbsoluteValue bold n Subscript normal s Baseline dot omega Subscript normal o Baseline EndAbsoluteValue Over StartAbsoluteValue bold n Subscript normal g Baseline dot omega Subscript normal o Baseline EndAbsoluteValue EndFraction f Superscript asterisk Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline comma omega Subscript normal i Baseline right-parenthesis period

Rather than integrating this special case into all BxDF subclasses, we find it cleaner to detect this case in the integrator and apply a correction factor

upper C Subscript shade Baseline left-parenthesis normal p Subscript Baseline comma omega Subscript normal o Baseline comma omega Subscript normal i Baseline right-parenthesis equals StartLayout Enlarged left-brace 1st Row 1st Column StartFraction StartAbsoluteValue bold n Subscript normal s Baseline dot omega Subscript normal o Baseline EndAbsoluteValue StartAbsoluteValue bold n Subscript normal g Baseline dot omega Subscript normal i Baseline EndAbsoluteValue Over StartAbsoluteValue bold n Subscript normal g Baseline dot omega Subscript normal o Baseline EndAbsoluteValue StartAbsoluteValue bold n Subscript normal s Baseline dot omega Subscript normal i Baseline EndAbsoluteValue EndFraction 2nd Column if importance is being transported 2nd Row 1st Column 1 2nd Column if radiance is being transported comma EndLayout

which corrects the normal dependence of the non-adjoint version into that of the adjoint when importance transport is indicated by the mode parameter. This adjustment is implemented by the helper function CorrectShadingNormal() below.

<<BDPT Utility Functions>>= 
Float CorrectShadingNormal(const SurfaceInteraction &isect, const Vector3f &wo, const Vector3f &wi, TransportMode mode) { if (mode == TransportMode::Importance) return (AbsDot(wo, isect.shading.n) * AbsDot(wi, isect.n)) / (AbsDot(wo, isect.n) * AbsDot(wi, isect.shading.n)); else return 1; }