13.2 Path Tracing

Now that we have derived the path integral form of the light transport equation, we will show how it can be used to derive the path-tracing light transport algorithm and will present a path-tracing integrator. Figure 13.4 compares images of a scene rendered with different numbers of pixel samples using the path-tracing integrator. In general, hundreds or thousands of samples per pixel may be necessary for high-quality results.

Figure 13.4: Kroken Scene Rendered with Path Tracing. (a) Rendered with path tracing with 8192 samples per pixel. (b) Rendered with just 8 samples per pixel, giving the characteristic grainy noise that is the hallmark of variance. Although the second image appears darker, the average pixel values of both are actually the same; very large values in some of its pixels cannot be displayed in print. (Scene courtesy of Angelo Ferretti.)

Path tracing was the first general-purpose unbiased Monte Carlo light transport algorithm used in graphics. Kajiya (1986) introduced it in the same paper that first described the light transport equation. Path tracing incrementally generates paths of scattering events starting at the camera and ending at light sources in the scene.

Although it is slightly easier to derive path tracing directly from the basic light transport equation, we will instead approach it from the path integral form, which helps build understanding of the path integral equation and makes the generalization to bidirectional path sampling algorithms easier to understand.

13.2.1 Overview

Given the path integral form of the LTE, we would like to estimate the value of the exitant radiance from the camera ray’s intersection point normal p 1 ,

upper L left-parenthesis normal p 1 right-arrow normal p 0 right-parenthesis equals sigma-summation Underscript i equals 1 Overscript normal infinity Endscripts upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis comma

for a given camera ray from normal p 0 that first intersects the scene at normal p 1 . We have two problems that must be solved in order to compute this estimate:

  1. How do we estimate the value of the sum of the infinite number of upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis terms with a finite amount of computation?
  2. Given a particular upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis term, how do we generate one or more paths normal p Subscript Baseline overbar Subscript in order to compute a Monte Carlo estimate of its multidimensional integral?

For path tracing, we can take advantage of the fact that for physically valid scenes, paths with more vertices scatter less light than paths with fewer vertices overall (this is not necessarily true for any particular pair of paths, just in the aggregate). This is a natural consequence of conservation of energy in BSDFs. Therefore, we will always estimate the first few terms upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis and will then start to apply Russian roulette to stop sampling after a finite number of terms without introducing bias. (Recall that Section 2.2.4 showed how to use Russian roulette to probabilistically stop computing terms in a sum as long as the terms that are not skipped are reweighted appropriately.) For example, if we always computed estimates of upper P left-parenthesis normal p Subscript Baseline overbar Subscript 1 Baseline right-parenthesis , upper P left-parenthesis normal p Subscript Baseline overbar Subscript 2 Baseline right-parenthesis , and upper P left-parenthesis normal p Subscript Baseline overbar Subscript 3 Baseline right-parenthesis but stopped without computing more terms with probability q , then an unbiased estimate of the sum would be

upper P left-parenthesis normal p Subscript Baseline overbar Subscript 1 Baseline right-parenthesis plus upper P left-parenthesis normal p Subscript Baseline overbar Subscript 2 Baseline right-parenthesis plus upper P left-parenthesis normal p Subscript Baseline overbar Subscript 3 Baseline right-parenthesis plus StartFraction 1 Over 1 minus q EndFraction sigma-summation Underscript i equals 4 Overscript normal infinity Endscripts upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis period

Using Russian roulette in this way does not solve the problem of needing to evaluate an infinite sum but has pushed it a bit farther out.

If we take this idea a step further and instead randomly consider terminating evaluation of the sum at each term with probability q Subscript i ,

StartFraction 1 Over 1 minus q 1 EndFraction left-parenthesis upper P left-parenthesis normal p Subscript Baseline overbar Subscript 1 Baseline right-parenthesis plus StartFraction 1 Over 1 minus q 2 EndFraction left-parenthesis upper P left-parenthesis normal p Subscript Baseline overbar Subscript 2 Baseline right-parenthesis plus StartFraction 1 Over 1 minus q 3 EndFraction left-parenthesis upper P left-parenthesis normal p Subscript Baseline overbar Subscript 3 Baseline right-parenthesis plus midline-horizontal-ellipsis comma

we will eventually stop continued evaluation of the sum. Yet, because for any particular value of i there is greater than zero probability of evaluating the term upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis and because it will be weighted appropriately if we do evaluate it, the final result is an unbiased estimate of the sum.

13.2.2 Path Sampling

Given this method for evaluating only a finite number of terms of the infinite sum, we also need a way to estimate the contribution of a particular term upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis . We need i plus 1 vertices to specify the path, where the last vertex normal p Subscript i is on a light source and the first vertex normal p 0 is a point on the camera film or lens (Figure 13.5). Looking at the form of upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis , a multiple integral over surface area of objects in the scene, the most natural thing to do is to sample vertices normal p Subscript i according to the surface area of objects in the scene, such that all points on surfaces in the scene are sampled with equal probability. (We do not actually use this approach in the integrator implementations in this chapter for reasons that will be described later, but this sampling technique could possibly be used to improve the efficiency of our basic implementation and helps to clarify the meaning of the path integral LTE.)

Figure 13.5: A path normal p Subscript Baseline overbar Subscript i of i plus 1 vertices from the camera at normal p Subscript , intersecting a series of positions on surfaces in the scene, to a point on the light normal p Subscript i . Scattering according to the BSDF occurs at each path vertex from normal p 1 to normal p Subscript i minus 1 such that the radiance estimate at the camera due to this path is given by the product of the path throughput upper T left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis and the emitted radiance from the light divided by the path sampling weights.

With this sampling approach, we might define a discrete probability over the n objects in the scene. If each has surface area upper A Subscript i , then the probability of sampling a path vertex on the surface of the i th object should be

p Subscript i Baseline equals StartFraction upper A Subscript i Baseline Over sigma-summation Underscript j Endscripts upper A Subscript j Baseline EndFraction period

Then, given a method to sample a point on the i th object with uniform probability, the probability density function (PDF) for sampling any particular point on object i is 1 slash upper A Subscript i . Thus, the overall probability density for sampling the point is

StartFraction upper A Subscript i Baseline Over sigma-summation Underscript j Endscripts upper A Subscript j Baseline EndFraction StartFraction 1 Over upper A Subscript i Baseline EndFraction comma

and all samples normal p Subscript i have the same PDF value:

p Subscript upper A Baseline left-parenthesis normal p Subscript i Baseline right-parenthesis equals StartFraction 1 Over sigma-summation Underscript j Endscripts upper A Subscript j Baseline EndFraction period

It is reassuring that they all have the same weight, since our intent was to choose among all points on surfaces in the scene with equal probability.

Given the set of vertices normal p 0 comma normal p 1 comma ellipsis comma normal p Subscript i minus 1 Baseline sampled in this manner, we can then sample the last vertex normal p Subscript i on a light source in the scene, defining its PDF in the same way. Although we could use the same technique used for sampling path vertices to sample points on lights, this would usually lead to high variance, since for all the paths where normal p Subscript i was not on the surface of an emitter, the path would have zero value. The expected value would still be the correct value of the integral, but convergence would be extremely slow. A better approach is to sample over the areas of only the emitting objects with probabilities updated accordingly. Given a complete path, we have all the information we need to compute the estimate of upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis ; it is just a matter of evaluating each of the terms.

It is easy to be more creative about how we set the sampling probabilities with this general approach. For example, if we knew that indirect illumination from a few objects contributed to most of the lighting in the scene, we could assign a higher probability to generating path vertices normal p Subscript i on those objects, updating the sample weights appropriately.

However, there are two interrelated problems with sampling paths in this manner. The first can lead to high variance, while the second can lead to incorrect results. The first problem is that many of the paths will have no contribution if they have pairs of adjacent vertices that are not mutually visible. Consider applying this area sampling method in a complex building model: adjacent vertices in the path will almost always have a wall or two between them, giving no contribution for the path and high variance in the estimate.

The second problem is that if the integrand has delta functions in it (e.g., a point light source or a perfect specular BSDF), this sampling technique will never be able to choose path vertices such that the delta distributions are nonzero. Even if there are no delta distributions, as the BSDFs become increasingly glossy almost all the paths will have low contributions since the points in f Subscript Baseline left-parenthesis normal p Subscript i plus 1 Baseline right-arrow normal p Subscript i Baseline right-arrow normal p Subscript i minus 1 Baseline right-parenthesis will cause the BSDF to have a small or zero value, and again we will suffer from high variance.

13.2.3 Incremental Path Construction

A solution that solves both of these problems is to construct the path incrementally, starting from the vertex at the camera normal p 0 . At each vertex, the BSDF is sampled to generate a new direction; the next vertex normal p Subscript i plus 1 is found by tracing a ray from normal p Subscript i in the sampled direction and finding the closest intersection. We are effectively trying to find a path with a large overall contribution by making a series of choices that find directions with important local contributions. While one can imagine situations where this approach could be ineffective, it is generally a good strategy.

Because this approach constructs the path by sampling BSDFs according to solid angle, and because the path integral LTE is an integral over surface area in the scene, we need to apply the correction to convert from the probability density according to solid angle p Subscript omega Sub Subscript to a density according to area p Subscript upper A (Section 4.2). If omega Subscript i minus 1 is the normalized direction sampled at normal p Subscript i minus 1 , it is:

p Subscript upper A Baseline left-parenthesis normal p Subscript i Baseline right-parenthesis equals p Subscript omega Sub Subscript Subscript Baseline left-parenthesis omega Subscript i minus 1 Baseline right-parenthesis StartFraction StartAbsoluteValue cosine theta Subscript i Baseline EndAbsoluteValue Over double-vertical-bar normal p Subscript i minus 1 Baseline minus normal p Subscript i Baseline double-vertical-bar squared EndFraction period

This correction causes all the factors of the corresponding geometric function upper G left-parenthesis normal p Subscript i plus 1 Baseline left-right-arrow normal p Subscript i Baseline right-parenthesis to cancel out of upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis except for the cosine theta Subscript i plus 1 term. Furthermore, we already know that normal p Subscript i minus 1 and normal p Subscript i must be mutually visible since we traced a ray to find normal p Subscript i , so the visibility term is trivially equal to 1. An alternative way to think about this is that ray tracing provides an operation to importance sample the visibility component of upper G .

With path tracing, the last vertex of the path, which is on the surface of a light source, gets special treatment. Rather than being sampled incrementally, it is sampled from a distribution that is just over the surfaces of the lights. (Sampling the last vertex in this way is often referred to as next event estimation (NEE), after a Monte Carlo technique with that name.) For now we will assume there is such a sampling distribution p Subscript normal e over the emitters, though in Section 13.4 we will see that a more effective estimator can be constructed using multiple importance sampling.

With this approach, the value of the Monte Carlo estimate for a path is

StartLayout 1st Row 1st Column upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis almost-equals 2nd Column StartFraction upper L Subscript normal e Baseline left-parenthesis normal p Subscript i Baseline right-arrow normal p Subscript i minus 1 Baseline right-parenthesis f Subscript Baseline left-parenthesis normal p Subscript i Baseline right-arrow normal p Subscript i minus 1 Baseline right-arrow normal p Subscript i minus 2 Baseline right-parenthesis upper G left-parenthesis normal p Subscript i Baseline left-right-arrow normal p Subscript i minus 1 Baseline right-parenthesis Over p Subscript normal e Baseline left-parenthesis normal p Subscript i Baseline right-parenthesis EndFraction 2nd Row 1st Column Blank 2nd Column times left-parenthesis product Underscript j equals 1 Overscript i minus 2 Endscripts StartFraction f Subscript Baseline left-parenthesis normal p Subscript j plus 1 Baseline right-arrow normal p Subscript j Baseline right-arrow normal p Subscript j minus 1 Baseline right-parenthesis StartAbsoluteValue cosine theta Subscript j Baseline EndAbsoluteValue Over p Subscript omega Sub Subscript Subscript Baseline left-parenthesis omega Subscript j Baseline right-parenthesis EndFraction right-parenthesis period EndLayout
(13.7)

Because this sampling scheme reuses vertices of the path of length i minus 1 (except the vertex on the emitter) when constructing the path of length  i , it does introduce correlation among the upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis terms. This does not affect the unbiasedness of the Monte Carlo estimator, however. In practice this correlation is more than made up for by the improved efficiency from tracing fewer rays than would be necessary to make the upper P left-parenthesis normal p Subscript Baseline overbar Subscript i Baseline right-parenthesis terms independent.

Relationship to the RandomWalkIntegrator

With this derivation of the foundations of path tracing complete, the implementation of the RandomWalkIntegrator from Chapter 1 can now be understood more formally: at each path vertex, uniform spherical sampling is used for the distribution p Subscript omega Sub Subscript —and hence the division by 1 slash left-parenthesis 4 pi right-parenthesis , corresponding to the uniform spherical PDF. The factor in parentheses in Equation (13.7) is effectively computed via the product of beta values through recursive calls to RandomWalkIntegrator::LiRandomWalk(). Emissive surfaces contribute to the radiance estimate whenever a randomly sampled path hits a surface with a nonzero upper L Subscript normal e ; because directions are sampled with respect to solid angle, the p Subscript normal e Baseline left-parenthesis normal p Subscript i Baseline right-parenthesis factor in Equation (13.7) is not over emissive geometry but is the uniform directional probability p Subscript omega Sub Subscript . Most of the remaining upper G factor then cancels out due to the change of variables from integrating over area to integrating over solid angle.