Chapter 1

Monte Carlo Path Tracing

This Chapter discusses Monte Carlo Path Tracing. Many of these ideas appeared
in James Kajiya’s original paper on the Rendering Equation. Other good original
sources for this material is L. Carter and E. Cashwell’s book Particle-Transport
Simulation with the Monte Carlo Methods and J. Spanier and E. Gelbard’s book
Monte Carlo Principles and Neutron Transport Problems.

1.1 Solving the Rendering Equation

To derive the rendering equation, we begin with the Reflection Equation
Z
Lr (~x, ω
~ r) = ~i → ω
fr (~x, ω ~ r ) Li (~x, ω
~ i ) cos θi dωi .
Ωi

The reflected radiance Lr is computed by integrating the incoming radiance over

a hemisphere centered at a point of the surface and oriented such that its north
pole is aligned with the surface normal. The BRDF fr is a probability distribution
function that describes the probability that an incoming ray of light is scattered in a
random outgoing direction. By convention, the two direction vectors in the BRDF
point outward from the surface.
One of the basic laws of geometric optics is that radiance does not change
as light propagates (assuming there is no scattering or absorption). In free space,
where light travels along straight lines, radiance does not change along a ray. Thus,
assuming that a point ~x on a receiving surface sees a point ~x0 on a source surface,
the incoming radiance is equal to the outgoing radiance:

~ i (~x0 , ~x)) = Lo (~x0 , ω

Li (~x, ω ~ o0 (~x, ~x)).

1
 Li(x,ωi)

Lr(x,ωr)

fr(x,ωi,ωr)

Figure 1.1: Integrating over the upper hemisphere.

Where we use the direction as a function of two points as

~x2 − ~x1
ω ~ (~x1 → ~x2 ) =
~ (~x1 , ~x2 ) = ω .
|~x2 − ~x1 |

The standard convention is to parameterize the incoming radiance Li at ~x with

the direction from the receiver ~x to the source ~x0 . When using this convention,
the incoming radiance is defined on a ray pointing in the direction opposite to the
direction of light propagation.
It is also useful to introduce notation for the two-point radiance

L(~x, ~x0 ) = L(~x → ~x0 ) = Lo (~x, ω

~ (~x, ~x0 )).

and the three-point BRDF

fr (~x, ~x0 , ~x00 ) = fr (~x → ~x0 → ~x00 ) = fr (~x0 , ω

~ (~x0 , ~x), ω
~ (~x00 , ~x0 ))

(Note: the two-point radiance function defined here is different than the two-point
intensity function defined in Kajiya’s original paper.)
Point to point functions are useful since they are intuitive and often clarify the
geometry and physics. For example, if ~x sees ~x0 , then ~x0 sees ~x. This mutual visi-
bility is represented as the two-point visibility function, V (~x, ~x0 ), which is defined

2
 dA(x')
x'
θ'

ω(x',x)

θ ω(x,x')

dA(x)
Figure 1.2: Two-point geometry.

to be 1 if a line segment connecting ~x to ~x0 does not intersect any opaque object,
and 0 otherwise.
The reflection equation involves an integral over the upper hemisphere. This
integral may be converted to an integral over other surfaces by changing of vari-
ables from solid angles to surface areas. This is easily done by relating the solid
angle subtended by the source to its surface area.

cos θo0
dωi = dA(~x0 )
|~x − ~x0 |2

The projected solid angle then becomes

cos θi dωi = G(~x, ~x0 ) dA(~x0 )

where
cos θi cos θo0
G(~x, ~x0 ) = G(~x0 , ~x) = V (~x, ~x0 )
|~x − ~x0 |2
In these equations we are making a distinction between the parameters used to
specify points on the surface (the ~x’s) and the measure that we are using when per-
forming the integral (the differential surface area dA(~x)). Sometimes we will be

3
less rigorous and just use dA or dA0 when we mean dA(~x) and dA(~x0 ). The ge-
ometry factor G is related to the differential form factor by the following equation:
0
F (~x, ~x0 )dA0 = G(~xπ,~x ) dA0 .
Performing this change of variables in the reflection equation leads to the fol-
lowing integral
Z
Lr (~x, ω~) = ~ (~x, ~x0 ), ω
fr (~x, ω ~ ) Lo (~x0 , ω
~ (~x0 , ~x)) G(~x, ~x0 ) V (~x, ~x0 ) dA(~x0 )
A

In this equation, ~x and ω

~ are indendent variables and we are integrating over surface
area which is parameterized by ~x0 ; thus, the incoming direction ω
~ i and the direction
of Lo are functions of these positions and are indicated as such.
The final step in the derivation is to account for energy balance

Lo (~x, ω
~ ) = Le (~x, ω
~ ) + Lr (~x, ω
~ ).

This states that the outgoing radiance is the sum of the emitted and reflected radi-
ances. Placing an emission function on each surface allows us to create area light
sources. Inserting the reflection equation into the energy balance equation results
in the Rendering Equation.
Z
L(~x, ω
~ ) = Le (~x, ω
~) + ~ (~x, ~x0 ), ω
fr (~x, ω ~ ) L(~x0 , ω
~ (~x0 , ~x)) G(~x, ~x0 ) V (~x, ~x0 ) dA0
A

For notational simplicity, we will drop the subscript o on the outgoing radiance.
The rendering equation couples radiance at the receiving surfaces (the left-
hand side) to the radiances of other surfaces (inside the integrand). This equation
applies at all points on all surfaces in the environment. It is important to recognize
the knowns and the unknowns. The emission function Le and the BRDF fr are
knowns since they depends on the scene geometry, the material characteristics, and
the light sources. The unknown is the radiance L on all surfaces. To compute the
radiance we must solve this equation. This equation is an example of an integral
equation, since the unknown L appears inside the integral. Solving this equation in
the main goal of Monte Carlo Path Tracing.
The rendering equation is sometimes written more compactly in operator form.
An operator is a method for mapping a function to another function. In our case,
the function is the radiance.

L = Le + K ◦ L

4
Sometimes it is useful to break the operator K into two operators, T and S. T is
the transfer operator and applied first; it takes outgoing light on one surface and
transfers it to another surface.

~ (~x0 , ~x)) = T ◦ L(~x, ω

Li (~x, ω ~ (~x, ~x0 ))

S is the scattering or reflection operator which takes the incoming light distribution
and computes the outgoing light distribution.

~ 0)
~ ) = S ◦ Li (~x, ω
Lr (~x, ω

Operator equations like the rendering equation may be solved using iteration.

L0 = Le
L1 = Le + K ◦ L0 = Le + K ◦ Le
...
n
X
Ln = Le + K ◦ Ln−1 = K n ◦ Le
i=0

Noting that K 0 = I, where I is the identity operator. This infinite sum is called the
Neumann Series and represents the formal solution (not the computed solution) of
the operator equation.
Another way to interpret the Neumann Series is to draw the analogy between

1
= (1 − x)−1 = 1 + x + x2 ...,
1−x
and
(I − K)−1 = I + K + K 2 ....

The rendering equation

(I − K) ◦ L = Le

then has the following solution

L = (I − K)−1 ◦ Le

Note that (I − K)−1 is just an operator acting on the emission function. This
operator spreads the emitted light over all the surfaces.

5
 fr(x2,x3,x4)
dA(x2)

dA(x0) G(x2,x3) dA(x4)

G(x1,x2)
G(x3,x4)
G(x0,x1) dA(x3)
dA(x1)
fr(x2,x3,x4)
fr(x1,x2,x3)

Figure 1.3: A path from point ~x1 to ~x5 .

It is useful to write out the formal solution

∞
X
L(~x, ω
~) = K i ◦ Le (~x0 , ω
~ 0 ))
i=0

in all its gory detail. Let’s consider one term:

Ln (~x, ω
~ ) = L(~xn , ~xn+1 ) = K n ◦ Le
Z Z
= ... Le (~x0 , ~x1 )G(~x0 , ~x1 )fr (~x0 , ~x1 , ~x2 )G(~x1 , ~x2 )...
A A
~ x)n
G(~xn−1 , ~xn )fr (~xn−1 , ~xn , ~xn+1 ) dA0 dA1 ...dA(~

This integral has a very simple geometric and physical intuition. It represents
a family of light paths. Each path is characterized by the number of bounces or
length n. There are many possible paths of a given length. Paths themselves are
specified by a set vertices. The first vertex is a point on the light source and sub-
sequent vertices are points on reflecting surfaces. The total contribution due to all
paths of a given length is formed by integrating over all possible light and surface
positions. This involves doing n integrals over surface areas. However, we must
weight paths properly when performing the integral. A integrand for a particular
path consists of alternating sequence of geometry and reflection terms. Finally, the
final solution to the equation is the sum of paths of all lengths; or more simply, all
possible light paths.
Note that these are very high dimensional integrals. Specifically, for a path of
length n the integral is over a 2n dimensional space. The integral also involves
very complicated integrands that include visibility terms and complex reflection

6
functions defined over arbitrary shapes. It turns out these complexities is what
makes Monte Carlo the method of choice for solving the rendering equation.
There is one more useful theoretical and practical step, and that is to relate the
solution of the rendering equation to the process of image formation in the camera.
The equation that governs this process is the Measurement Equation
Z Z Z
M= R(~x, ω
~ , t)L(~x, ω
~ , t) dt dω dA.
A Ω T

The response function R(~x, ω

~ , t) depends on the the pixel filter (the ~x dependence),
the aperture (the ω
~ dependence), and the shutter (the t dependence). Other factors
such as transformations of rays by the lens system and spectral sensitivities may
also be included, but we will ignore these factors to simplify the presentation.
As seen above, the pixels value in the image is a function that involves nested
integrals. These integrals are very complicated, but we can easily evaluate the
integrand which corresponds to sampling the function.

• Sampling a pixel over (x, y) prefilters the image and reduces aliasing.

• Sampling the camera aperture (u, v) produces depth of field.

• Sampling in time t (the shutter) produces motion blur.

• Sampling in wavelength λ simulates spectral effects such as dispersion

• Sampling the reflection function produces blurred reflection.

• Sampling the tranmission function produces blurred transmission.

• Sampling the solid angle of the light sources produces penumbras and soft
shadows.

• Sampling paths accounts for interreflection.

Sampling in x, y, u, v and t has been discussed previously. Sampling light sources

and performing hemispherical integration has also been discussed. What remains
is to sample paths.

7
1.2 Monte Carlo Path Tracing
First, let’s introduce some notation for paths. Each path is terminated by the eye
and a light.

E - the eye.

L - the light.

Each bounce involves an interaction with a surface. We characterize the interac-

tion as either reflection or tranmission. There are different types of reflection and
transmission functions. At a high-level, we characterize them as

D - diffuse reflection or transmission

G - glossy reflection or tranmission

S - specular reflection or refraction

Diffuse implies that light is equally likely to be scattered in any direction. Specular
implies that there is a single direction; that is, given an incoming direction there is
a unique outgoing direction. Finally, glossy is somewhere in between.
Particular ray-tracing techniques may be characterized by the paths that they
consider.

Appel Ray casting: E(D|G)L

Whitted Recursive ray tracing: E[S ∗ ](D|G)L

Kajiya Path Tracing: E[(D|G|S)+ (D|G)]L

Goral Radiosity: ED∗ L

The set of traced paths are specified using regular expressions, as was first proposed
by Shirley. Since all paths must involve a light L, the eye E, and at least one
surface, all paths have length at least equal to 3.
A nice thing about this notation is that it is clear when certain types of paths
are not traced, and hence when certain types of light transport is not considered
by the algorithm. For example, Appel’s algorithm only traces paths of length 3,
ignoring longer paths; thus, only direct lighting is considered. Whitted’s algorithm
traces paths of any length, but all paths begin with a sequence of 0 or more mirror

8
reflection and refraction steps. Thus, Whitted’s technique ignores paths such as
the following EDSDSL or E(D|G)∗ L. Distributed ray tracing and path tracing
includes multiple bounces involving non-specular scattering such as E(D|G)∗ L.
However, even these methods ignore paths of the form E(D|G)S ∗ L; that is, multi-
ple specular bounces from the light source as in a caustic. Obviously, any technique
that ignores whole classes of paths will not correctly compute the solution to the
rendering equation.
Let’s now describe the basic Monte Carlo Path Tracing Algorithm:

Step 1. Choose a ray given (x,y,u,v,t)

weight = 1

Step 2. Trace ray to find point of intersection with the nearest surface.

Step 3. Randomly decide whether to compute emitted or reflected light.

Step 3A. If emitted,

return weight * Le
Step 3B. If reflected,
weight *= reflectance
Randomly scatter the ray according to the BRDF pdf
Go to Step 2.

This algorithm will terminate as long as a ray eventually hits a light source. For
simplicity, we assume all light sources are described by emission terms attached to
surfaces. Latter we will discuss how to handle light sources better.
A variation of this algorithm is to trace rays in the opposite direction, from light
sources to the camera. We will assume that reflective surface never absorb light,
and that the camera is a perfect absorber.

Step 1. Choose a light source according to the light source power distribution.
Generate a ray from that light source according to its intensity distribution.
weight = 1

Step 2. Trace ray to find point of intersection.

Step 3. Randomly decide whether to absorb or reflect the ray.

9
 Step 3A. If scattered,
weight *= reflectance
Randomly scatter the ray according to the BRDF.
Go to Step 2.
Step 3B. If the ray is absorbed by the camera film,
Record weight at x, y
Go to Step 1.

The first algorithm is an example of forward ray tracing; in forward ray tracing
rays start from the eye and propagate towards the lights. Forward ray tracing is
also called eye ray tracing. In contrast, in backward ray tracing rays start at the
light and trace towards the eye. As we will discuss in a subsequent section, Both
methods are equivalent because the physics of light transport does not change if
paths are reversed. Both methods have there advantages and disadvantages, and in
fact may be coupled.
The above simple algorithms form the basis of Monte Carlo Path Tracing.
However, we must be more precise. In particular, there are two theoretical and
practical challenges:

Challenge 1 : Sampling an infinite sum of paths in an unbiased way.

Challenge 2 : Finding good estimators with low variance.

1.3 Random Walks and Markov Chains

To understand more about why path tracing works, let’s consider a simpler prob-
lem: a discrete random walk. Instead of a physical system with continuous vari-
ables, such as position and direction, consider a discrete physical system comprised
of n states. Path tracing as described above is an example of a random walk where
we move from sample to sample, or from point to point, where the samples are
drawn from a continuous probability distribution. In a discrete random walk, we
move from state to state, and the samples are drawn from a discrete probability
distribution.
A random walk is characterized by three distributions:

1. Let p0i be the probability of starting in state i.

2. Let pi,j is the probability of moving from state i to state j.

10
 Transition
pi,j

Creation p0i i p*i Termination

Figure 1.4: State transition diagram for a discrete random walk.

3. Let p∗i is the probability of being terminated in state i.

Because the probability of transition and termination must sum to one, p∗i = 1 −
P
j=0 pi,j ; that is, the probability of terminating in state i is equal to the probability
of not moving from state i to j.
A discrete random walk consists of the following steps.

Step 1. Create a random particle in state i with probability p0i .

Step 2. With probability p∗i , terminate in state i.

Score particle in state i by incrementing the counter for state i
Go to Step 1.

Step 3. Randomly select new state according to the transition probability distribu-
tion.
Set i to new j.
Go to Step 2.

Random walks are also called Markov Chains. A Markov Chain is a sequence
of states generated by a random process. We will return to Markov Chains when
we discuss the Metropolis Algorithm. Markov Chains also come up Bayesian rea-
soning and in learning theory (e.g. Hidden Markov Models). Keep your eyes open
for Markov Chains; you will see these techniques used more and more in computer
graphics.
Given a set of particles following random walks, the problem is to compute the
final probability of a particle being terminated in state i. To solve this problem,
we introduce another random variable Pin , which is the probability of being in

11
state i after n transitions. Since each state transition is independent of the previous
transitions, this probability may be computed using a simple recurrence

Pj0 = p0j
X
Pj1 = p0j + pi,j Pi0
i
...
X
Pjn = p0j + pi,j Pin−1 .
i

Defining a matrix M whose entries are Mi,j = pi,j , the above process can be
viewed as the following iterative product of a matrix times a vector

P 0 = p0
P 1 = p0 + M P 0
...
P n = p0 + M P n−1 .

And this procedure may be recognized as the iterative solution of the following
matrix equation
(I − M )P = p0

since then
X
P = (I − M )−1 p0 = p0 + M (p0 + M (p0 ... = M i p0 .
i=0

This process will always converge assuming the matrices are probability distribu-
tions. The basic reason for this is that probabilities are always less than one, and so
a product of probabilities quickly tends towards zero.. Thus, the random walk pro-
vides a means for solving linear systems of equations, assuming that the matrices
are probability transition matrices. Note the similiarity of this discrete iteration of
matrices to the iterative application of the continuous operator when we solve the
rendering equation using Neumann Series.
This method for solving matrix equations using discrete random walks may be
directly applied to the radiosity problem. In the radiosity formulation,
X
Bi = Ei + ρi Fi,j Bj
j

12
where the form factor equals
Z Z
1
Fi,j = G(~x, ~x0 )V (~x, ~x0 ) dA(~x) dA(~x0 )
πAi Ai Aj

Recall that the form factors may be interpreted as a probabilities. The form factor
Fi,j is the percentage of outgoing light leaving surface element Ai that falls on
surface element Aj . In fact, the form factor is the probability that a random ray
leaving surface element Ai makes it to Aj (although one has to be careful about
how one defines a random ray). Thus, form factor matrices may be interprested
as transition matrices. In this equation, ρ is the diffuse reflectance and is equal to
ρ = B/E; The reflectance must be positive and less than 1. The absorption or
termination probability is thus equal to 1 − ρ.
These observations lead to a method for solving the matrix radiosity equation
using a discrete random walk. The major issue, which is always a case with ra-
diosity solution techniques, is computing the form factor matrix. This process is
expensive and error prone because of the complexity of the environment and the
difficulty in doing exact visible surface determination. The form-factor matrix is
also very large. For example, a scene consisting of a million surface elements
would require a million squared matrix. Therefore, in practice, the form factor ma-
trix is often calculated on-the-fly. Assuming a particle is on some surface element
i, an outgoing particle may be sent off in a random direction where the random
direction is chosen from a cosine-weighted distribution (here the cosine is with re-
spect to the surface element normal). The particle is then ray-traced and the closent
point of intersection on surface element j is found. This random process is roughly
equivalent to one generated from a known form-factor matrix.
It is interesting to prove that random walks provide an unbiased estimate of
the solution of the linear system of equations. Although this proof is a bit formal,
it is worthwhile working it through to get a flavor of the mathematical techniques
involved.
The first step is to define a random variable on the space of all paths. Let’s
signify a path of length k as αk = (i1 , i2 , ..., ik ); this path involves a sequence of
transitions from state i1 to state i2 and ending up finally after k transitions in state
ik . The random variable α without the subscript is the set of all paths of length one
to infinity.
The next step is to pick an estimator W (α) for each path. Then we can compute
the expected value of the estimator by weighting W by the probability that a given

13
path is sampled, p(α).
X
E[W ] = p(α)W (α)
α
∞ X
X
= p(αk )W (αk )
k=1 αk
∞ X
X X
= ... p(i1 , ..., ik )W (i1 , ..., ik )
k=1 i1 ik

In the last line we group all paths of the same length together. The sums on each
index i go from 1 to n - the number of states in the system. Thus, there are nk
paths of length k, and of course paths can have infinite length. There are a lot of
paths to consider!
What is the probability p(αk ) of path αk ending in state ik ? Assuming the dis-
crete random walk process described above, this is the probability that the particle
is created in state i1 , times the probability of making the necessary transitions to
arrive at state ik , times the probability of being terminated in state ik

p(αk ) = p0i1 pi1 ,i2 ...pik−1 ,ik p∗k

With these careful definitions, the expected value may be computed

∞ X
X
E[Wj ] = p(αk )Wj (αk )
k=1 αk
X∞ X X
= ... p0i1 pi1 ,i2 ...pik−1 ,ik p∗ik Wj (αk )
k=1 i1 ik

Recall, that our estimator counts the number of counts the number of particles that
terminate in state j. Mathematically, we can describe this counting process with
a delta function, Wj (αk ) = δik ,j /p∗j . This delta function only scores particles
terminating in ik = j. The expected value is then
∞ X
X XX
E[Wj ] = ... p0i1 pi1 ,i2 ...pik−1 ,ik p∗ik δik ,j /p∗j
k=1 i1 ik−1 ik
X∞ X X
= ... p0i1 pi1 ,i2 ...pik−1 ,j .
k=1 i1 ik−1

14
This sum may be recognized as the j component of the matrix equation

E[W ] = p0 + M p0 + M 2 p0 + ...

which is the desired solution of the linear system of equations.

Note that we had to be careful about the estimator. If we hadn’t divided the
count of the particles by the probability of a termination event, the expected value
would not have equaled the right answer. Picking the wrong estimator - that is, an
estimator that results in the wrong expected value - for a complex sampling process
is one of the most common errors when uding Monte Carlo Techniques. Until you
have a lot of experience, it is worthwhile convincing yourself that your estimator
is unbiased.
This technique was originally developed by von Neumann and Ulam, the orig-
inators of the Monte Carlo Method. The estimator they used is often called the
absorption estimator, since only particles that are absorbed are counted. An inter-
esting variation, developed by Wasow, is to count all the number of particles that
pass through state j (including those terminate as well as those that make a tran-
sition). This is called the collision estimator, since it counts all particles colliding
with a surface. It is an interesting exercise to show that the collision estimator also
provides an unbiased estimate of the solution to the linear equation. It is more chal-
lenging, but more interesting, to also derive the conditions when and if the collision
estimator works better than the absorption estimator.
This basic proof technique is easy to generalize to continuous distributions,
but the notation is messy. The details are described in Chapter 3 of the Spanier
and Gelbard book on neutron transport [?], the most authoritative source if you
wish to understand the theory behind Monte Carlo Techniques for solving integral
equations and tranport problems.

1.4 Adjoint Equations and Importance Sampling

Recall, that the pixel response is equal to the sum over paths of length n
Z Z
Mn = ... S(~x0 , ~x1 )G(~x0 , ~x1 )fr (~x0 , ~x1 , ~x2 )G(~x1 , ~x2 )
0 n
~ x)n .
...fr (~xn−2 , ~xn−1 , ~xn )G(xn−1 , xn )R(~xn−1 , ~xn )dA0 dA1 ...dA(~

where we have switched notation and written the source term as S(~x, ~x0 ) = Le (~x, ~x0 ).

15
 As noted above this equation is symmetric under the interchange of lights and
sensors. Switching Le with R, and noting that
Z Z
Mn = ... R(~x0 , ~x1 )G(~x0 , ~x1 )fr (~x0 , ~x1 , ~x2 )G(~x1 , ~x2 )
A A
~ x)n
...fr (~xn−2 , ~xn−1 , ~xn )G(xn−1 , xn )S(~xn , ~xn−1 ) dA0 dA1 ...dA(~
Z Z
= ... S(~xn , ~xn−1 )G(xn , xn−1 )fr (~xn , ~xn−1 , ~xn−2 )G(~xn−1 , ~xn−2 )
A A
~ x)n
...fr (~x2 , ~x1 , ~x0 )G(~x1 , ~x0 )R(~x1 , ~x0 ) dA0 dA1 ...dA(~

In the second step, we noted from the symmetry of the geometry that

G(~xi , ~xj ) = G(~xj , ~xi )

and because of the reciprocity principle the BRDF is also symmetric

fr (~xi , ~xj , ~xk ) = fr (~xk , ~xj , ~xi )

These symmetries implie that we may ray trace from either the light or the eye;
both methods will lead to the same integral.
Suppose now we break the path at some point k. The amount of light that
makes to k is
Z Z
LS (~xk , ~xk+1 ) = ... S(~x0 , ~x1 )G(~x0 , ~x1 )fr (~x0 , ~x1 , ~x2 )G(~x1 , ~x2 )
A A
~ x)k−1
...G(~xk−1 , ~xk )fr (~xk−1 , ~xk , ~xk+1 ) dA0 ...dA(~

In a similar way, treating the sensor as a virtual light source, we can compute the
amount of light coming from the sensor makes it to k.
Z Z
LR (~xk , ~xk+1 ) = ... fr (~xk , ~xk+1 , ~xk+2 )G(~xk+1 , ~xk+2 )
k+2 n
~ x)n
...fr (~xn−2 , ~xn−1 , ~xn )G(xn−1 , xn )R(~xn−1 , ~xn ) dAk+2 ...dA(~

The measured response is then

Z Z
M= LS (~xk , ~xk+1 )G(~xk , ~xk+1 )LR (~xk , ~xk+1 ) dAk dAk+1
A A

Note the use of the notation LS and LR to indicate radiance “cast” from the source
vs. the receiver.

16
 dA(x2)

LS(x2,x3)
R
S G(x2,x3) LR(x2,x3)

dA(x3)

Figure 1.5: A path with both the forward and backward (adjoint) solution of the
transport equation. The forward solution is generated from the source term S and
the backward solution is generated from the received term R. For physical situa-
tions where the transport equation is invariant under path reversal, the forward and
backward equations are the same.

We make two observations about this equation. First, this equation can be
considered the inner product of two radiance functions. If we consider radiance to
be a function on rays r = (~x, ω
~ ), then if we have functions f (r) and g(r), the inner
product of f and g is
Z
< f, g >= f (r)g(r)dµ(r)

where dµ(r) is the appropriate measure on rays. The natural way to measure
the rays between two surface elements A and A0 is dµ(r) = G(x, x0 ) dA dA0 .
Equivalently, considering r to be parameterized by position ~x and direction ω
~ , the
~
dµ(r) = d~ω ◦ dA(~x)(~x).
Second, this integral naturally leads to a method for importance sampling a
path. Suppose we are tracing light and arrive at surface k. To compute the sensor
response, we need to integrate L against R. In this sense, R may be considered
an importance function for sampling the next directions, since we want a sampling
technique that is proportional to R to achieve low variance. But R is the solution
of the reversed transport equation that would be computed if we were to trace rays
from the sensor. R tells us how much light from the sensor would make it to this
point. Thus, the backward solution provides an importance function for the forward
solution, and vice versa. This is the key idea between bidirectional ray tracing.
Manipulating about adjoint equations is easy using the operator notation. Using

17
the operator notation, an integral equation is just
Z
K ◦ f = K(x, y)f (y) dy.

We want to estimate the integral given by the measurement equation, which is just
the inner product of two functions
Z Z 
M =< f, K ◦ g >= f (x) K(x, y)g(y) dy dx.

This of f as the response of the sensor and K ◦ g as the solution of the rendering
equation. This equation may be rearranged
Z Z 
< f, K ◦ g > = f (x) K(x, y)g(y) dy dx
Z 
= f (x)K(x, y) dx g(y) dy

= < K + f, g > .
Note the difference between
Z
K ◦f = K(x, y)f (y) dy

and Z
K+ ◦ f = K(x, y)f (x) dx.
One integral is over the first variable, the other is over the second variable. Of
course, if K(x, y) = K(y, x) these two integrals are the same, in which case
K + = K and the operator is said to be self-adjoint.
This notation provides a succinct way of proving that the forward estimate is
equal to the backward estimate of the rendering equation. Recall
K ◦ LS = S
We can also write a symmetric equation in the other direction
K ◦ LR = R
Then,
< R, LS > = < K ◦ LR , LS >
= < LR , K + ◦ LS >
= < LR , K ◦ LS >
= < LR , S >

18
This result holds even if the operator is not self-adjoint. We will leave the demon-
stration of that fact as an exercise.
This is a beautiful result, but what does it mean in practice. Adjoint equations
have lots of applications in all areas of mathematical physics. What they allow you
to do is create output sensitive algorithms. Normally, when you solve an equation
you solve for the answer everywhere. Think of radiosity; when using the finite
element method you solve for the radiosity on all the surfaces. The same applies to
light ray tracing or the classic discrete random walk; you solve for the probability
of a particle landing in any state. However, in many problems you only want to
find the solution at a few points. In the case of image synthesis, we only need to
compute the radiance that we see, or that falls on the film. Computing the radiance
at other locations only needs to be done if its effects are observable.
We can model the selection of a subset of the solution as the inner product
of the response function times the radiance. If we only want to observe a small
subset of the solution, we make the response function zero in the locations we
don’t care about. Now consider the case when all the surfaces act as sources and
only the film plane contributed a non-zero response. Running a particle tracing
algorithm forward from the sources would be very inefficient, since only rarely
is a particle terminated on the film plane. However, running the algorithm in the
reverse direction is very efficient, since all particles will terminate on sources. Thus
each particle provides useful information. Reversing the problem has led to a much
more efficient algorithm.
The ability to solve for only a subset of the solution is a big advantage of
the Monte Carlo Technique. In fact, in the early days of the development of the
algorithm, Monte Carlo Techniques were used to solve linear systems of equations.
It turns out they are very efficient if you want to solve for only one variable. But
be wary: more conventional techniques like Gaussian elimination are much more
effective if you want to solve for the complete solution.

19