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 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:

where is the measurement for the th pixel and is a point on the film. In this setting, the term is the product of the filter function around the pixel and a delta function that selects the appropriate camera ray direction of the sample from , :

This formulation may initially seem gratuitously complex, but it leads us to an important insight. If we expand the terms of the LTE sum, we have

where is the path throughput function introduced in Equation (14.18).

Note the nice symmetric way in which the emitted radiance (quantifying the light source emission profile) and the weighting function (quantifying the camera’s sensitivity profile for a pixel ) 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 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 term is known as the importance for the ray between and 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 and direction and evaluates the importance emitted from the point on the camera in a direction . 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 such that attains its maximum value; in practice, the function returns fractional values to specify the raster positions more accurately.

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.

Given the camera-to-world transformation for the provided time, the method checks that the direction points in the same hemisphere as the camera is facing by transforming the camera-space viewing direction to world space and computing the cosine of the angle 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.

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 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 . 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 to get a point along the ray leaving the camera before performing the projection.

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

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 . The following fragment that is part of the PerspectiveCamera constructor uses the RasterToCamera transformation and divides by the coordinate to compute the rectangle’s corner points, which in turn gives the rectangle’s area .

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 (); 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 ; thus, the area-measured PDF for points on the image plane is . 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

where is the angle that makes with the image rectangle normal and is the distance between the point on the lens and the intersection of the ray with the plane. The distance to the point on the image plane is

as is the coordinate of in the local camera coordinate system and . Hence,

By construction, the above density function is normalized when integrated over directions—however, we initially set out to create a normalized importance function that is defined on the space of camera rays , where is the surface region associated with the perspective camera’s lens element. This function must satisfy the ray-space normalization criterion

We can’t directly set equal to due to the extra integration over areas and the additional cosine factor in the above integral.

Note that the lens area of a perspective camera is equal to , where is the lens radius. For a pinhole camera, the lens area is set to 1 and interpreted as a Dirac delta function.

We then define on ray space as

which divides 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 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).

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.

The directional density of the ideal sampling strategy implemented in the method PerspectiveCamera::GenerateRay() was already discussed and is equal to 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.

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.

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

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.

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

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).

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.

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

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.

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 for the spotlight’s illumination distribution is separable with . We thus just need to find a sampling distribution for . The UniformSampleCone() function from Section 13.6.4 provides this functionality.

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.

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.

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

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.

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.

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.

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.

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.

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.

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.

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 from Equation (16.1), which was defined as

The vertices are ordered such that denotes the -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 at vertex , whose only role is to evaluate the original BSDF with swapped arguments:

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 and . 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.

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 ( 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 and transmitted radiance are related by

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

and we can obtain the adjoint BTDF

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).

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 is flipped, hence the expression muI * muO > 0 can be used to check if light is being refracted rather than reflected.

Finally, the transmissive term 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).

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 with an interpolated shading normal 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 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:

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

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

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

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.