Manifold Exploration:
A Markov Chain Monte Carlo technique forrendering scenes with difficult specular transport
Expanded Technical Report
Wenzel Jakob

Steve Marschner

Cornell University
May 7, 2012
Abstract: It is a long-standing problem in unbiased Monte Carlo methods for
rendering that certain difficult types of light transport paths, particularly those
involving viewing and illumination along paths containing specular or glossy
surfaces, cause unusably slow convergence. In this paper we introduce Manifold
Exploration, a new way of handling specular paths in rendering. It is based on
the idea that sets of paths contributing to the image naturally form manifolds inpath space, which can be explored locally by a simple equation-solving iteration.
This paper shows how to formulate and solve the required equations using
only geometric information that is already generally available in ray tracing
systems, and how to use this method in in two different Markov Chain Monte
Carlo frameworks to accurately compute illumination from general familiesof paths. The resulting rendering algorithms handle specular, near-specular,
glossy, and diffuse surface interactions as well as isotropic or highly anisotropic
volume scattering interactions, all using the same fundamental algorithm. An
implementation is demonstrated on a range of challenging scenes and evaluatedagainst previous methods.
Remark: This technical report follows the structure of the main paper. The
main difference is that it contains additional discussion in Sections 3.3 and 6, asecond appendix with C++ code, as well as MEMLT results in Section 8.

1

Figure 1: Two views of an interior scene with complex specular and near-specular transport,
rendered using manifold exploration path tracing: (a) Refractive, reflective, and glossy
tableware, (b) A rough brass luminaire with 24 glass-enclosed bulbs used to light theprevious closeup.

1 Introduction
Certain classes of light paths have traditionally been a source of difficulty in conducting Monte
Carlo simulations of light transport. A well-known example is specular-diffuse-specular paths,
such a tabletop seen through a drinking glass sitting on it, a bottle containing shampoo or othertranslucent liquid, or a shop window viewed and illuminated from outside. Even in scenes wherethese paths do not cause dramatic lighting effects, their presence can lead to unusably slow
convergence in renderers that attempt to account for all transport paths.
Furthermore, wherever ideally specular paths are troublesome, nearly specular paths involving
glossy materials are also troublesome. They can be more problematic, in fact, because they elude
special mechanisms designed to handle specular interactions. These glossy paths have become
more important as material models have evolved, and we would prefer to handle them using
natural generalizations of strategies for specular surfaces, rather than generalizations of strategies
for diffuse surfaces.
Finding these paths efficiently is a key problem of light transport simulations. In this paper,
we show how a Markov Chain Monte Carlo method can be used to efficiently render paths
with illumination and/or viewing through arbitrary chains of specular or glossy reflections or
refractions. The idea is that sets of paths contributing to the image naturally form manifolds
in path space, and using a simple equation-solving iteration it is easy to move around on thesemanifolds. We show how to formulate and solve the equations required to find specular paths,
using only geometric information that is already generally available in ray tracing systems. Thissolution method is then used to define a Markov Chain, which has the right properties to be
applied in a Metropolis-Hastings algorithm to compute lighting through very general families ofpaths that can involve specular, near-specular, glossy, and diffuse surface interactions as well as
isotropic or highly anisotropic volume scattering interactions.
We begin by discussing prior work on which our method is built, then in Section 3 we develop
the theory of the specular manifold, used to handle interactions with ideal specular (polished)surfaces, and offset specular manifolds, which provide a graceful generalization to near-specularmaterials. In Section 5, we derive an algorithm to move from one path to another on a specular
or offset specular manifold, use it to build an algorithm that generates Markov sequences in path2

space, and show how this can be used in the Metropolis Light Transport or Energy Redistribution
Path Tracing frameworks to provide methods for rendering scenes with any kind of light transport.
We go on in the following sections to extend the theory and algorithm to the case of participating
media, and we finish by showing results and comparisons in Section 8.

2 Prior work
Simulating light transport has been a major effort in computer graphics for over 25 years, beginningwith the complementary approaches of finite-element simulation, or radiosity [Goral et al. 1984],
and ray-tracing [Whitted 1980]. The introduction of Monte Carlo methods for ray tracing [Cooket al. 1984], followed by Kajiya’s formulation of global illumination in terms of the Rendering
Equation [Kajiya 1986], established the field of Monte Carlo global illumination. Unbiased
sampling methods, in which each pixel in the image is a random variable with an expected value
exactly equal to the solution of the Rendering Equation, started with Kajiya’s original path tracing
method and continued with bidirectional path tracing [Lafortune and Willems 1993; Veach andGuibas 1994], in which light transport paths can be constructed partly from the light and partly
from the eye, and the seminal Metropolis Light Transport [Veach and Guibas 1997] algorithm,
which uses bidirectional path tracing methods in a Markov Chain Monte Carlo framework.
Various two-pass methods use a particle-tracing pass that sends energy out from light sourcesin the form of “photons” that are traced through the scene [Shirley et al. 1995; Jensen 1996;
Walter et al. 1997]. and stored in a spatial data structure. The second pass then renders theimage using ray tracing, making use of the stored particles to estimate illumination by density
estimation. [Jensen and Christensen 1998; Jarosz et al. 2008; Hachisuka and Jensen 2009].
Photon mapping and other two-pass methods are characterized by storing an approximate
representation of some part of the illumination in the scene, which requires assumptions about
the smoothness of illumination distributions. On one hand, this enables rendering of some modesof transport that are difficult for unbiased methods, since the exact paths by which light travels
do not need to be found; separate paths from the eye and light that end at nearby points sufficeunder assumptions of smoothness. However, this smoothness assumption inherently leads tosmoothing errors in images: the results are biased, in the Monte Carlo sense. Glossy-glossy
transport, without a sufficiently diffuse surface on which to store photons, is challenging tohandle with photon maps, since large numbers of photons must be collected to adequately sample
position-direction space. Some photon mapping variants avoid this by treating glossy materials as
specular, but this means that the the resulting method increasingly resembles path tracing as thenumber of rough surfaces in the input scene grows.

2.1 Specular reflection geometry
Separate from work on global illumination algorithms, various research has examined the propertiesof specular reflection paths. Mitchell and Hanrahan [1992] devised a method to compute irradiancefrom implicitly defined reflectors, using Fermat’s principle with interval Newton’s method to locate
all reflection paths from a source to a point. Walter et al. [2009] proposed a related method thatcomputes the singly scattered radiance within a refractive object with triangle mesh boundaries.
Like these works, our method searches for specular paths. But because it does so within theneighborhood of a given path, it avoids the complexities and constraints entailed by a full global

3

search. Another difference is that our manifold formalism can be used to build a fully general
rendering system that is not limited to the specular paths that prompted its design.
The widely used method of ray differentials for texture antialiasing [Igehy 1999] also reasons
about the local structure of a set of reflected paths—in this case, paths from the eye. Igehy’s
approach requires elementary local differential information only, in the form of derivatives of
surface normals, and does not require global surface descriptions as Mitchell and Hanrahan’smethod does. Manifold exploration requires the same geometric information, and can thus be
implemented in most modern ray tracing sytems.
The analysis of reflection geometry presented by Chen and Arvo [2000a; 2000b] is closest to
the mathematics underlying our proposed methods. Their work relies on a characterization of
specular paths via Fermat’s principle, which asserts that light travels along paths whose opticallength (i.e. the propagation time) constitutes a local extremum amongst neighboring paths.
Using Lagrange multipliers, the authors derive a path Jacobian and path Hessian with respect
to perturbations of the endpoint of a path and use it to accelerate the interactive display ofreflections on curved surfaces.
Our local characterization of specular paths is equivalent to Fermat’s principle, and the path
Jacobian is related to the derivatives that we propose to use to define tangent spaces to the
specular manifold while solving for path transitions. However, the use of this derivative and thegoals of the research are entirely different: in their case, estimating changes to viewing paths,
and in our case, tracking the evolution of specular paths in a very general context, as part of anunbiased rendering system.
A less related but relevant idea is integrating over continuous paths in volume rendering
applications, known as the path integral formulation of radiative transfer [Premoze et al. 2003].
This work, with its implications for the concentration of transport in path space, suggests thepossibility of using MCMC to integrate over multiple-scattering paths in volumes, as discussed in
Section 7.

2.2 Markov Chain Monte Carlo in rendering
The Metropolis Light Transport algorithm mentioned above introduced the tools of Markov Chain
Monte Carlo (MCMC) to rendering. A Markov chain is a sequence of points in a state space inwhich the probability of a state appearing at a given position in the sequence depends only on theprevious state. The basic idea of MCMC, first proposed by Metropolis et al. [1953], is to define aMarkov chain that has the function to be integrated as its stationary distribution, meaning that
if the chain is run for a long time the distribution of states it visits will be proportional to thedesired distribution. The following brief introduction follows [Liu 2001].
Defining a Markov chain amounts to defining a transition rule: a process for selecting a newstate x+ randomly, in a way that depends on the current state x. Metropolis et al. provided a way
to take a transition rule that may not produce the desired stationary distribution π(x) and turn it
into one that does. Given a method for sampling a proposal distribution T (x, x0 ), the Metropolistransition rule operates in two steps:
1. Choose x0 according to the probability distribution T (x, x0 ).
(x0 with probability min(1, π(x0 )/π(x))
2. x+ =
x otherwise4

In step 1 we say x0 is proposed as the next state, and in step 2 it is either accepted and becomesthe next state, or it is rejected and the next state repeats the previous one. The probabilitymin(1, π(x0 )/π(x)) is known as the acceptance probability.
If this chain is able to pass from any state in the domain to another using a finite expected
number of steps, and if the length of this journey is in a sense “irregular” (i.e. aperiodic), it isreferred to as ergodic, and the chain’s limiting distribution is guaranteed to converge to π. This
relatively mild criterion is usually satisified by choosing a transition rule with global support.
The original Metropolis algorithm only works when T (x, x0 ) = T (x0 , x). Hastings [1970]
proposed a new acceptance probability:π(x0 )T (x0 , x)0(1)r(x, x ) = min 1,
π(x)T (x, x0 )which relaxes the symmetry restriction to one of symmetric support: T (x, x0 ) must be nonzero
exactly when T (x0 , x) is nonzero. The Metropolis–Hastings algorithm is the starting point for
MCMC rendering methods.
In the rendering context, the state space is the space of all paths through the scene, points inthe space are paths, and the desired probability distribution over paths is proportional to theircontribution to the rendered image (i.e. the amount of illumination they carry to the camera).
The final image is the projection of the path distribution into the image plane.
At the core of an MCMC rendering algorithm is an implementation of a transition rule, andany rule with symmetric support is admissible. But to avoid very low acceptance probabilities,
which lead to poor performance, it is desirable for the transition probability to approximate thecontribution: that is, paths with more light flowing along them should be chosen more often.
Veach’s [1997] transition rule is based on a set of mutations that change the structure of the pathand perturbations that move the vertices by small distances while preserving the structure, bothusing the building blocks of bidirectional path tracing to sample paths. Using his transition rule
to run long chains leads to the Metropolis Light Tranport algorithm (MLT). Kelemen et al. [2002]
later proposed a transition rule based on changing the random numbers used to sample paths bybidirectional path tracing.
Considerable research activity has extended Metropolis light transport in various ways. Paulyet al. [2000] proposed a perturbation rule for rendering participating media with single scattering.
Other projects include Metropolis Instant Radiosity [Segovia et al. 2007], Population Monte Carlorendering [Lai et al. 2007], and Replica Exchange light transport [Kitaoka et al. 2009]. Recently,
two groups [Chen et al. 2011; Hachisuka and Jensen 2011] have combined the transition rule ofKelemen et al. with photon mapping to obtain robust methods based on density estimation.
However, to generate proposals, all of these algorithms ultimately rely on local path samplingstrategies (i.e. path tracing). Specifically, they choose the next interaction vertex along a lightpath by sampling from a directional distribution associated with the current vertex, followed by anintersection search. Our method introduces a new kind of transition rule with different properties.
The original MLT algorithm and subsequent variants all render an image by running a Markov
chain for a long (e.g. > 106 ) sequence of steps, and they guarantee ergodicity using a transition
rule that can generate any path in the domain with some probability. The Energy Redistribution
path tracing (ERPT) technique [Cline et al. 2005], which is readily adapted to work with our
method, is an interesting departure from that approach. It draws on the property that the5Metropolis-Hastings algorithm (1) preserves the stationary distribution of samples even if theunderlying transition rule is not ergodic (e.g. when it cannot reach parts of path space). To obtaincoverage, ERPT first samples a large set of paths via path tracing and runs Markov chains for
short bursts (≈ 103 steps) starting at each sample. The relaxation of the ergodicity requirement
makes it possible to explore paths in a very local fashion.

3 Path space manifolds
Manifold exploration is a technique for integrating the contributions of sets of specular or nearspecular illumination paths to the rendered image of a scene. The general approach applies tosurfaces and volumes and to ideal and non-ideal (glossy) specular surfaces. In this section we
begin by examining the manifold defined by ideal specular reflection or refraction, in the settingof surfaces without participating media. In the following section we develop this theory into a
rendering method for scenes combining ideal specular surfaces with fairly diffuse surfaces. We willthen go on to generalize the method to glossy surfaces, then to generalize the theory to encompass
participating media and to extend the method to handle media with both isotropic and highly
directional scattering.

3.1 Path space
The resulting techniques are all based on the path integral formulation of light transport, describedby Veach [1997] and others, in which the value of each pixel in the image is an integral of acontribution function over path space. Denoting the union of all surfaces in the scene by M, eachsequence x1 . . . xn of at least two points in M is a path along which light may travel. Thus pathspace is
P=

∞
[

Pk

(2)

n=2

Pn = {x1 . . . xn | x1 , . . . , xn ∈ M} .

(3)

As detailed by Veach, the value of pixel j isZIj =fj (x̄)dµ(x̄)
P

where x̄ = x1 . . . xn denotes an element of P and µ is the product measure derived from the area
measure on M. The contribution function fj is a product of terms, one for each vertex and eachedge of the path:fj (x1 · · · xn ) = Le (x1 → x2 )"n−1#YG(xk−1 ↔ xk ) fs (xk−1 → xk → xk+1 ) G(xn−1 ↔ xn ) We(j) (xn−1 → xn ). (4)
k=2

6

S

D

S

Figure 2: The geometry factor (left) and the generalized geometry factor (right) are bothderivatives of projected solid angle at one end with respect to area at the far end.
(j)Here Le is emitted radiance, We is importance for the jth pixel, fs is the BSDF at xk for thegeometry defined by xk−1 , xk , and xk+1 , and G is the geometry factor:G(x ↔ y) =→ |N (y) · −→|N (x) · −xy|yx|
V (x ↔ y)kx − yk2(5)where N (a) is the surface normal at a and V (x ↔ y) is the visibility function.

3.2 Motivating examples
This formulation is the basis for several rendering methods including bidirectional path tracing andMetropolis light transport. However, in the presence of ideal specular reflection, some difficulties
arise, which are normally sidestepped in the transition from theory to algorithm, but which we
prefer to confront directly. When some surface interactions are specular, the entire contribution to
the path space integral is from paths that obey specular reflection or refraction geometry, and theset of such paths is lower in dimension than the full path space. For instance, consider a familyof paths of the form LDSDE (in Heckbert’s [1990] notation) with one specular reflection vertex.
These paths belong to the P5 component of P, but the paths that contribute all have the property−→ −−→
x−
3 x2 + x3 x4 k N (x3 ),that is, the half-vector at x3 is in the direction of the normal. This places two constraints on thepath, meaning that all contributing paths lie on a manifold of dimension 8 embedded in P5 , which
is of dimension 10. The integral is more naturally expressed as an integral over the manifold,rather than as an integral over the whole path space.
To compute illumination due to specular paths, we use a local parameterization of the manifold
in terms of the positions of all nonspecular vertices on the path:
ZZZZ
f (x1 . . . x5 )dx1 dx2 dx4 dx5
M4

Note the missing integral over x3 , the specular vertex. The contribution function f still has

7the same form, a product of terms corresponding to vertices and edges of the path, but the
BSDF values at the specular vertex are replaced by (unitless) specular reflectance values, andthe geometry factors for the two edges involving the specular vertex are replaced by a single
generalized geometry factor that we will denote G(x2 ↔ x3 ↔ x4 ).
The standard geometry factor for a non-specular edge is the derivative of projected solid angle
at one vertex with respect to area at the other vertex, and the generalized geometry factor isdefined analogously: the derivative of solid angle at one end of the specular chain with respect
to area at the other end of the chain, considering the path as a function of the positions of theendpoints. Figure 2 illustrates this for a more complex path involving a chain of three specularvertices. We will explain below how G can be easily computed from the differential geometry of
the specular manifold.

3.3 Specular manifold geometry
In the general case, each path of length k belongs to a class in {D, S}k based on the classificationof each of its vertices. (In this scheme point or orthographic cameras, and point or parallel lights,
are denoted S, while finite-aperture cameras and area lights are D.) Each S surface vertex has an
associated constraint that involves its position and the position of the preceding and following
vertices:
ci (xi−1 , xi , xi+1 ) = 0
The constraint function computes a half-vector at vertex i and projects it into the shading tangent
space; the resulting 2-vector is zero when the half-vector is parallel to the normal. By making useof the generalized half-vector of Walter et al. [2007], both reflection and refraction can be handled
by a single constraint function:−−→ −−−−→
ci (xi−1 , xi , xi+1 ) = T (xi )T h(xi , −x−
i xi−1 , xi xi+1 )η(x, v)v + η(x, w)wh(x, v, w) =kη(x, v)v + η(x, w)wk

(6)
(7)

where T (x) is a matrix whose columns form a basis for the shading tangent plane at x, and η(x, v)
denotes the refractive index associated with the ray (x, v). This generalized geometry factor isrelated to the “extended form factor” discussed by Sillion and Puech [1989].
Specular endpoint vertices also introduce constraints that depend on their type: for instance,
the position of a vertex on a point emitter must remain fixed (e.g. x1 = const). When the emitteris directional, it is the outgoing direction that is constrained, and so on.
From now on, we implicitly identify each vertex xi with an associated point in IR2 using
local parameterizations of M. These parameterizations may be defined on arbitrarily smallneighborhoods, since only their derivatives are relevant to what follows.
With this, the constraints for a length-n path with p specular vertices can be stacked together
into a function C : IR2n → IR2p , and the specular manifold is simply the set
S = {x̄ | C(x̄) = 0}(8)

Expressing S using a constraint in this way makes it convenient to work with neighborhoods ofa particular path. The Implicit Function Theorem [Spivak 1965] guarantees the existence of a

8

D

where

S
D

S
(b) Associated constraints(a) An example path(c) Constraint Jacobian(d) Tangent space

Figure 3: The linear system used to compute the tangent space to the specular manifold, also
known as the derivative of a specular chain with respect to its endpoints.
parameterization of the manifold, in the neighborhood of any path x̄ that is nonsingular (in thesense explained below). This parameterization is a function q : IR2(n−p) → IR2p that determines the
positions of all the specular vertices from the positions of all the nonspecular vertices. Furthermore,
the derivative of q, which gives us the tangent space to the manifold at x̄, is simple to computefrom the derivative of C.
For the specifics we restrict ourselves to the case of a single chain of specular vertices with
non-specular vertices (surfaces, cameras, or light sources) at the ends, which suffices to cover mostcases by applying it separately to mutiple specular chains along a path. Paths with specular
endpoints are handled with simple variations of this scheme. Number the vertices in the chain
x1 , . . . , xk , with x1 and xk being the (non-specular) endpoints of the segment and the remainingk − 2 vertices being specular. In this case C : IR2k → IR2(k−2) , and the derivative ∇C is a matrixof k − 2 by k 2-by-2 blocks, with a tridiagonal structure (Figure 3).
The Implicit Function Theorem gives us a parameterization of the manifold in terms of any
2 vertices, and if we pick x1 and xk this simply says that the path, in a neighborhood of the9current path1 , is a function of the two endpoints. Furthermore, it also tells us the derivative ofthat parameterization, which is to say, the derivative of all the specular vertices’ positions withrespect to the positions of the endpoints. If we block the derivative ∇C, as shown in the figure,
into 2-column matrices B1 and Bk for the first and last vertices and a square matrix A for the
specular chain, then this tangent space to the manifold isTS (x̄) = −A−1 B1 Bk
This matrix is k − 2 by 2 blocks in size, and each block gives the derivative of one vertex (in termsof its own tangent frame) with respect to one endpoint.
The matrix A fails to be invertible when the specular vertices are not even locally unique for
given endpoints, which means that one endpoint is on a caustic for light emitting from the otherendpoint.
We use TS (x̄) for two things: to navigate on the manifold and to compute the generalized
geometry factor. The right two or left two columns of TS (x̄) are useful for updating the specular
chain with x1 or xk held fixed, respectively, as discussed in the next section.
The top-right or bottom-left block of TS (x̄) can be used to compute the generalized geometry
factor as follows. Assuming orthonormal parameterizations2 , the determinant of the top-rightblock gives the ratio of an infinitesimal area at xk to its reflection/refraction, as observed fromx1 , measured on the surface at x2 . To convert this to a ratio of area at xk to solid angle at x1 ,we multiply this determinant by the ordinary geometry factor G(x1 ↔ x2 ); this product is thegeneralized geometry factor G(x1 ↔ · · · ↔ xk ). More succinctly,
G(x1 ↔ · · · ↔ xk ) = P2 A−1 Bk G(x1 ↔ x2 )= Pk−1 A−1(9)

B1 G(xk−1 ↔ xk ),where Pi is a 2 by 2(k − 2) matrix that projects onto the two dimensions associated with vertex i.
A useful property of this framework is its reliance on local information that is easily provided
in ray tracing-based rendering systems. To compute the blocks of the A and B matrices, wemust have access to the partial derivatives of position and shading normal with respect to any
convenient parameterization of the surfaces, along with the refractive indices of all objects. Theseare exactly the same quantities also needed to trace ray differentials through refractive boundaries,
which is part of many mature ray tracing-based rendering systems. A consequence of the simple
form of the constraint (8) is that our technique works with any object that can provide such local
information, including implicitly defined shapes or triangle meshes with shading normals.
These theoretical results about the structure of the specular manifold can be used in an algorithm
to solve for specular paths, which we discuss in Section 4.
A simple example: We now demonstrate how to compute the tangent space and generalized
geometric term associated with a simple example path shown in Figure 4. This is a path withthree vertices; x1 and x3 are planar endpoints used to parameterize the manifold, and x2 lies on a
1

Because it is possible to have several separated specular paths joining two points, the parameterization cannot be
global.
2
If the parameterizations are not orthonormal, two additional determinants are required to account for the change
in area.

10

2

0
-1

S
0

1

Figure 4: An example path with three vertices
specularly reflecting cylinder that stretches out along the z-axis.
Assuming that this path is encountered in a ray tracer similar to PBRT [Pharr and Humphreys
2004] or Mitsuba [Jakob 2010], the rendering system will associate an intersection data record
with each vertex containing (amongst other things) the position and surface normal, as well asderivatives thereof along some parameterization. Suppose that the information provided is asfollows:
x1 = −1, 2,0 , ∂u x1 = −1, 0,0 , ∂v x1 = 0, 0, 1, n1 = 0, −1,0 , ∂u n1 = 0, ∂ v n1 = 0
x2 = 0, 1, 0, ∂u x2 = 1, 0, 0, ∂v x2 = 0, 0, 1, n2 = 0, 1, 0 ,  ∂u n2 = 1, 0, 0 , ∂v n2 = 0
x3 = 1, 2, 0 , ∂u x3 = 1, 0, 0 , ∂v x3 = 0, 0, 1 , n3 = 0, −1, 0 , ∂u n3 = 0,∂v n3 = 0Note that the parameterizations above are locally orthonormal—this is generally preferable, since
it simplifies many computations involving ∇C. In practice, we can simply apply a suitable linear
transformation to the tangents and normal derivatives to to make them correspond to a locallyorthonormal chart. Let us now compute the entries of ∇C.
Recall that the constraint associated with vertex x2 is
x3 − x2
x 1 − x2
x3 − x2
x1 − x2
T
C(x1 , x2 , x2 ) = T (x2 )++kx1 − x2 k kx3 − x2 k
kx1 − x2 k kx3 − x2 k
To compute derivatives of this expression, we can parameterize all needed ingredients to first orderusing the above data records, i.e.
xi (ui , vi ) = xi + ui (∂u xi ) + vi (∂v xi ), ni (ui , vi ) = ni + ui (∂u ni ) + vi (∂v ni )∂u x2 − h∂u x2 , n2 (u2 , v2 )in2 (u2 , v2 )T (u2 , v2 ) =∂v x2 − h∂u x2 , n2 (u2 , v2 )in2 (u2 , v2 )The rest is just a big nested application of the product rule. When differentiating C, we get two
terms: one which accounts for changes of the tangent frame at x2 , with h held constant, and one
which accounts for changes of h, while the tangent frame is held constant. For the former, we
require derivatives of T , e.g:




∂T
∂T
−h∂u x2 , ∂u n2 in2 − h∂u x2 , n2 i∂u n2−h∂u x2 , ∂v n2 in2 − h∂u x2 , n2 i∂v n2==,−h∂v x2 , ∂u n2 in2 − h∂v x2 , n2 i∂u n2−h∂v x2 , ∂v n2 in2 − h∂v x2 , n2 i∂v n2
∂u2
∂v2

11

and for the second term we can repeatedly apply the following vector calculus identity:


∂z(t)
∂ z(t)1 ∂z(t)
z(t)
z(t),
=
−
∂t kz(t)kkz(t)k ∂tkz(t)k3∂twhere z is a vector-valued function that depends on t. The resulting expression is quite messy
and thus we only provide numerical values here. In particular, ∇C is given by 1

− 4 0 − 32 0 41 0∇C =0 12 0 −1 0 12
and therefore−1

TS (x̄) = −A



B1 B3 =



− 61
0

0
1
2

1
6

0

0
1
2


.

Since the parameterizations are orthonormal, the geometric term between x1 and x3 is simply
G(x1 ↔ x3 ) = P2 A−1 B3 G(x1 ↔ x2 ) =

1
6

0

0
1
2

G(x1 ↔ x2 ) =1.
48

For a C++ implementation that computes ∇C, please refer to the appendix A2.

4 Walking on the specular manifold
Our rendering algorithms use MCMC to explore the manifold of specular paths, and for thisthey require some form of local parameterization. With the differential geometry of the specular
manifold in hand, we are now able to develop this extremely useful building block. We propose
an algorithm that moves one of the endpoints of a specular chain and takes all the intermediate
vertices to a valid new configuration. Later, in Section 5, we show how to apply this algorithm torender images.
To simplify the discussion, we will focus on the case where the position of a vertex xn of a
specular chain x1 , . . . xn is adjusted to a given new position x0n , while x1 is held fixed. We shallalso briefly introduce the assumption that xn is located on planar surface of infinite extent.
Our manifold walking algorithm is based on two key insights:
1. The A and B matrices (Section 3.3) may be used to map an infinitesimal in-plane movementof xn to displacements of the vertices x2 , . . . , xn−1 . We can use these displacements toapproximate a finite change to the path simply by adding an offset to each vertex, but thiswill move the path off the specular manifold.
2. Ray tracing provides a deterministic means of projecting an off-specular path back onto thespace of valid configurations based on its first two vertices. Given x1 and x2 , we can trace asequence of rays xi → xi+1 , at each step performing a specular reflection or refraction, and+this leads to corrected positions x+
2 . . . xn .
By combining 1. and 2., we obtain a predictor-corrector type algorithm (Figure 5) that performs
a step according to a local linear model, followed by a projection that restores the specular+configuration, resulting in a new path x1 , x+
2 , . . . , xn , and repeats until convergence. As long as

12the prediction step solves the linear model and moves in the tangent space to the manifold, thisiteration behaves like Newton’s method, exhibiting quadratic convergence near the solution. As
with all Newton-like iterations, it is not guaranteed to converge when started far from the solution,
since the linear model may not be accurate enough to make progress. But since the model isfirst-order accurate, the algorithm is guaranteed to make forward progress when the constraintfunction is differentiable and the partial steps are small enough. Our algorithm uses a simple
heuristic to decrease the step size when progress is not made, then increase back to full steps toget quadratic convergence as it approaches the target configuration. This iteration is illustrated
in Figure 5 and laid out in the following algorithm:
WalkManifold x1 , . . . , xn

x0n )1 Set i = 0 and β = 12 while kxn − x0n k > εL3
p = x2 − β T (x2 )P2 A−1 Bk T (xn )T (x0n − xn )4
Propagate the ray x1 → p through all specular+interactions, producing x+
2 , . . . , xn .
005
if step 4 succeeded and kx+
n − xn k < kxn − xn k+
6x2 , . . . , xn = x+
2 , . . . , xn
7
β = min {1, 2 β}8else
9β = 21 β
10
Set i = i + 1, and fail if i > N .
11 return x2 , . . . , xn−1
Here, i records the number of iterations until a specified maximum N is reached, and ε is arelative error threshold to a scene-scale length L (we use L = maxi kxi k, N = 20 and ε = 10−7 ).
The variable β denotes a step size that is dynamically adjusted to ensure that steps reduce
the target distance. The pseudocode assumes that the plane vertex xn has an orthonormal
parameterization.
To make this algorithm usable for general specular chains, we must remove the previous
assumption that the endpoint xn is located on a plane. However, the actual locations of x+n along
the way from xn to x0n are not needed; all that matters is that x0n is on the surface at convergence.
Therefore, we construct a plane containing xn and x0n , and in step 4 of WalkManifold the lastpropagation step computes an intersection against this hypothetical plane, ignoring the actual
scene geometry. Once the algorithm converges, we must ensure that xn−1 and xn are mutuallyvisible before reporting the manifold walk as successful.
In our implementation, we assign the plane normal using the following symmetric orthogonalization procedure−−−→−−−→ 
nplane := γ γ(n + n0 ) − xn x0n hγ(n + n0 ), xn x0n iwhere n and n0 are the surface normals at xn and x0n , respectively, and γ(v) := v/kvk. This
ensures that the same plane is used for walks xn x0n and x0n xn .

13

(a)

path space(fixed)starting pointnext iterategoal point(b)corrected

predicted
starting
point

path space

Figure 5: One iteration of updating a path using the specular manifold. (a) The path vertices
are modified according to a local linear model, which (b) corresponds to a step along
the tangent plane to the manifold, then (a) a nearby valid specular path is found,which (b) corresponds to projecting back onto the manifold.
To compute the matrix A, we derived symbolic expressions for the manifold constraints in C
(Equation 8). A C++ implementation is given in the appendix. When solving the resulting linear
system in step 3, it is beneficial to exploit the special structure of this matrix, which becomes
important when processing specular chains with more than about ten vertices. We solve for pusing a block tridiagonal LU factorization, and this reduces the time complexity from O(n3 ) toO(n), n being the number of vertices in the chain.
There are several situations under which this algorithm may fail to converge: first, a specularpath between x1 and x0n need not exist at all. Secondly, WalkManifold usually cannot findpaths that lie on a different connected component of the manifold. Thirdly, when the local
structure of the manifold is complex (e.g. due to high-frequency geometric detail of the reflectors
and refractors) and kxn − x0n k is large, the iteration may not converge to a solution. Finally,
the linear system may not be invertible, which happens when the last path lies at the fold of a
wavefront, e.g. a caustic receiving an infinite power per unit area. During the ∼ 1012 manifold
walks performed to produce the results of the paper, this last case only occurred ∼ 1.7 · 105 times,
hence it does not appear to be an issue in practice. However, in the MCMC context it is not an
problem for the iteration to fail occasionally, as explained in the next section.
The manifold walking algorithm works reliably for large chains with over 10 vertices, especiallywhen it is used by the transition rule discussed in the next section, which only moves the endpointsof chains by a small amount. In our scenes, we observe between 92 and 98% successful walks, eachtaking 2-3 iterations on average to converge to the tolerance ε = 10−7 . The failing 2-8% mainlycontain cases where WalkManifold failed for good reasons, because it was asked to walk to apoint for which there is no valid configuration on the manifold.

5 Manifold exploration for surfaces
In this section, we present a new transition rule that proposes steps in path space using manifold
walks. In the context of MCMC, a transition rule, or perturbation, is a random process thatgenerates a proposal state conditioned on the current state of a Markov chain. It provides the basicmeans of navigating through the state space, but to do this correctly the rule must satisfy two basic
criteria: transitions must be reversible (i.e. return to the previous state with nonzero probability
density), and the rule must also supply a function that, up to constant factors, computes the14

probability of proposals conditioned on the current state.
To create an efficient sampling procedure, a transition rule should furthermore propose modifications of an appropriate scale. A rule that takes large steps will tend to leave local maxima
of the target distribution π, and such steps are rejected with high probability. A rule that takes
tiny steps will find most of them accepted, but it will not explore the state space well. Our
perturbation is designed so that its scale naturally adapts to the scene, including the geometry
and material properties.
The new perturbation supports general scenes and can be used both with the ERPT algorithmby Cline et al. and the path space MLT framework proposed by Veach and Guibas, where it
replaces and generalizes the lens, caustic, and multi-chain perturbations. Depending on which
combination is used, we call the resulting algorithm either Manifold Exploration Path Tracing(MEPT) or Manifold Exploration Metropolis Light Transport (MEMLT).

5.1 Manifold perturbation
Given an input path, the manifold perturbation finds a nearby path using a sequence of steps
that can be grouped into sampling and connection phases (Figure 6). The sampling phasechooses a subpath to be modified that consists of three non-specular vertices which are potentiallyseparated by specular chains. We shall denote these non-specular vertices as xa , xb , and xc . Afterestablishing the type of perturbation to be performed, the sampling step generates a perturbed
outgoing direction at the vertex xa and propagates it through the specular chain between xa
and xb (if any) until arriving at a new non-specular vertex x0b in the neighborhood of xb . If
the path configuration (i.e. the arrangement of specular and nonspecular vertices) changed, theperturbation is rejected immediately.
Up to this point, the proposed scheme is similar to the set of perturbations proposed by Veach
and Guibas. However, recall that in their work, perturbations must propagate through the pathuntil arriving at a pair of adjacent non-specular vertices (“DD” in Heckbert’s [1990] notation)that can be used to establish a connection edge. Any attempt to connect two sampled subpaths
that involves a specular vertex must fail, since the probability of creating a valid path in thismanner is zero.
In comparison, our perturbation can stop at vertex xb and use the WalkManifold algorithm
to update the configuration of the specular chain between xb and xc (Figure 7). This seemingly
subtle difference has major repercussions on the types of scenes that can be rendered efficiently. In
particular, the resulting method can systematically explore large classes of specular paths insteadof having to rely on finding them by chance.

Connect

SampleFigure 6: The manifold perturbation samples a perturbed outgoing direction from a vertex xa
and propagates it through a specular chain (if any) using ray tracing until arriving at a
non-specular vertex xb . To connect the vertices xb and xc , the perturbation performsa manifold walk to determine the positions of intermediate specular vertices (if any).

15

S
ed

trac

ted

a
upd

S

half-vector equal
to surface normal

D

Figure 7: Manifold perturbation example: a slightly perturbed outgoing direction at xa is
propagated until encountering the non-specular vertex x0b . Previously, it was not clear
how to “connect” x0b to xc through multiple specular interactions. Our method canfind this connection given knowledge about the previous path.
In the following, we will discuss the perturbation in more detail; first for the ideally specularcase and then in a more general form that extends to rough surfaces.
Strategy sampling: Motivated by the desire to attempt a large range of different types of pathmodifications with high probability, the sampling step first chooses among the possible perturbation
strategies for a given path by selecting three vertices as follows. Given a path x1 , . . . , xk , uniformly
select a non-specular initial vertex xa , as well as a perturbation direction (i.e. towards the lightsource or towards the camera). Walk along the path in this direction until the first non-specularvertex is encountered, and continue until a second non-specular vertex is found. This path traversal
may fail by walking past the end of the path, in which case the strategy sampling phase is simply
restarted from scratch. This determines xa , xb and xc . For notational convenience, assume thata < b < c.
Perturbation sampling: With the overall strategy established, the sampling phase now perturbs
the path segment xa+1 , . . . , xb . The goal here is to produce a new subpath x0a+1 , . . . , x0b thatis “nearby”. When the vertex xa denotes a surface scattering event with incident and exitant
−−→
−−−−→
0
directions ωi = −x−a xa−1 and ωo = xa xa+1 , the perturbation determines xa+1 by tracing a ray in a0
direction ωo that is sampled from a suitable spherical distribution D(ωo0 ) concentrated aroundωo . It is absolutely critical that this distribution generates direction changes of the appropriate
scale: for instance, when xa is a diffuse material, relative large perturbations are in order. On theother hand, when xa is a glossy material that only reflects into a small cone of directions, large
perturbations will almost always be rejected, reducing performance.
Observe that a useful hint about the right scale can be obtained directly from the scatteringmodel at xa , in particular from the associated importance sampling density p(ωi → ωo ). When
this sampling density is high, ωo is likely located on a sharp peak of the scattering function, andsmall steps are appropriate.
Our strategy is to sample from a distribution D(ωo0 ) centered at ωo whose concentration is set
so that D(ωo ) equals λ2 p(ωi → ωo ). For D(ωo0 ), we use the spherical von Mises-Fisher distribution.

16

(a) λ = 1(b) λ = 3(c) λ = 10(d) λ = 30(e) λ = 90(f ) λ = 270

Figure 8: The effects of the λ parameter. This scene was rendered with relatively short Markov
Chains (100 steps), hence the range of “good” parameter values is lower than in ourother examples scenes. Modeled after a scene by Eric Veach.
17

The parameter λ (generally set between 50 and 500) specifies how large the perturbations are
relative to standard BSDF sampling. This is the main parameter of our technique, and it affects
how far perturbations will move in path space. When λ is set to an inappropriately low or highvalue, the amount of noise present in the output renderings increases. In the first case, too few
mutations are accepted, causing the chain to become “stuck” in certain paths for many iterations.
In the latter case, the steps taken by the chain are too small to effectively explore path space, andthis results in the typical coherent noise patterns that are known from other MLT-type algorithms.
A comparison involving different settings is shown in Figure 8. We currently set this parameter
manually to achieve a desired acceptance ratio, but this could in theory be automated using
adaptive MCMC.
When xa is a camera or light source, we choose a new outgoing direction in much the same way,
but query the underlying model for the directional density of the associated sampling method (e.g.
uniformly choosing pixels in screen-space). To further enlarge the space of possible perturbations,following Veach, we split camera and light source endpoints into two separate vertices corresponding
to the position and direction components. When xa is such a position vertex, we perturb its
location on the aperture or light source by sampling a tangential displacement from a 2D normal
distribution with variance ρ/(2πλ2 ), where ρ is the surface area.
After x0a+1 has been determined in this manner, the perturbation is propagated through the
specular chain until reaching x0b . This process is deterministic.
Connection: When there is no specular chain between x0b and xc , the connection step only entails
checking that the vertices are mutually visible, and that their scattering models carry illumination
along the connection edge. When there is a manifold, we first set
x0c−1 , . . . x0b+1 = WalkManifold(xc , . . . xb → x0b )and then perform the same verification.
Recall that a key requirement of the Metropolis-Hastings algorithm discussed in Section 2.2
was that a nonzero transition probability T (x, x0 ) > 0 also implies that T (x0 , x) > 0. This createsa potential issue when walking on the manifold, because WalkManifold can be non-reversible.
It might succeed in moving from x to x0 but fail to move from x0 to x. Even when the reverse
iteration converges, the manifold can contain bifurcations so that it may converge to a different
solution. Therefore, we always perform another manifold walk in the reverse direction and reject
the perturbation if the path did not return to its original configuration. In the example scenes, weobserved between 0% and 0.4% non-reversible walks.
Transition probability computation Finally, the change in the contribution function is computed
and, together with the transition probabilities, used to randomly accept or reject the proposal
with probability
fj (x̄0 )T (x̄0 , x̄)0r(x̄, x̄ ) = min 1,(10)
fj (x̄)T (x̄, x̄0 )where fj is the contribution function (Equation 4) and x̄ and x̄0 , denote the original and proposal
path respectively. Note that many factors cancel in the above ratio, particularly all of thosein fj that are associated with the unchanged path segment, or common terms in the transition18

probability. For instance, the probability of choosing a particular sampling strategy cancels, since
it only depends on the (unchanged) path configuration.
We require that T (x̄, x̄0 ) and T (x̄0 , x̄) express the density of forward and reverse proposals in
a common measure. Observe that sampling an outgoing direction ωo0 from D(ωo0 ) at xa , andpropagating it through the first specular chain, produces area density D⊥ (ωo0 ) G(xa ↔ · · · ↔ x0b )on the surface at x0b (where D⊥ (ωo0 ) = D(ωo0 )/| cos(na , ωo )| denotes probability with respect tothe projected solid angle measure at xa ). This is the needed transition probability T (x̄, x̄0 ).
In a MCMC-based rendering system, we will generally want to sample paths based on theircontribution to the entire image rather than to a single pixel. This is accomplished by replacing(j)
We in (4) with an importance function that measures the overall luminance received on theimage plane.

5.2 Extension to glossy materials
The method just presented can be used for scenes with both specular and non-specular transport,
and the two classes are treated separately. This is unfortunate, since a near-specular chain throughan almost-smooth dielectric object cannot be explored as effectively as a perfectly specular one.
However, it turns out that a simple generalization suffices to encompass these materials as well.
The path space integrand corresponding to a chain of glossy interactions has its energy concentrated in a thin “band” near the specular manifold, and the key idea of how we handle glossy
materials is to take steps along a family of offset manifolds that are parallel to the specular
manifold, so that path space near the specular manifold can be explored without stepping outof this thin band of near-specular transport. In this section, we add a simple extension thatendows the perturbation with the ability to walk on offset manifolds and to recognize when this isappropriate.
For this, we first replace Equation (8) with the offset manifold
So = {x̄ | C(x̄) = o} ,(11)where o captures the offset from ideal specular transport. Inspecting C (Equation 6) yields an
intuitive explanation for the contents of the vector o in terms of microfacet theory. The two
entries associated with each vertex xi record the microsurface normal mi responsible for thereflection or refraction (which is now different from the shading normal ni ), projected into theshading tangent space T (xi ). Our extended perturbation then preserves the projection of thismicrogeometry normal mi as an invariant during the manifold traversal.
Changes to the manifold perturbation: Since the differential geometry of the offset manifold isidentical to that of the ordinary manifold, the only required change in WalkManifold affects
the ray tracing step on line four, where the algorithm now reflects and refracts using mi instead ofni . Similarly, the deterministic phase of the manifold perturbation responsible for propagating thesampled direction at xa to a position x0b uses these normals instead. Note that it is straightforward
to handle both cases, near-specular and specular perturbations, using the same implementation.
Recognizing near-specular transport: An important aspect about the treatment of generalscattering is the decision of whether the surface associated with a scattering event is “smooth

19

S
ed

trac

d
ate
d
p

u

“S”

half-vector preservedin surface frame

D

Figure 9: Perturbation of path with near-specular surface interactions: instead of requiring thatthe half-vectors agree with the surface normals, their direction is preserved in thesurface frame.
enough” to be classified as part of a specular chain. We make this decision randomly by assigninga specular probability ψ(xi ) to each vertex that takes on values 0 and 1 when xi is diffuse andspecular, respectively, and values in (0, 1) when xi is at a rough interface. This avoids the issuesof “hard” classifications that are commonly used in rendering algorithms. For specifics on our
specular probability function ψ, please refer to the appendix.
Transition probability: In the purely specular case discussed earlier, the proposal and target
distributions were both supported on the same space S. This property permitted computing transition probabilities under an arbitrary projection (e.g. onto the vertex xb ), since the determinantof the associated change of variables canceled when considering ratios of the form fj /T .
In the glossy case, this does not hold anymore: T S
is a distribution on an offset manifold So ,whereas fj is defined on the higher-dimensional space o So . To compute transition probabilities,
we must perform a change of variables to separate out these “perpendicular” dimensions in fj ,(a) always non-specular(b) always specular(c) probabilistic

Figure 10: Consistently classifying glossy materials as non-specular or specular produces unsatisfactory results. Instead, our method makes this decision randomly wheneverencountering a rough object (modeled after a scene by Cline et al.)20

A perturbation with 3 glossy
and 1 specular vertices“S”

invert“S”(full matrix)

“S”
S

“S”

discardentries ofspecularvertices

Figure 11: An illustration of the matrices involved in the glossy transition probability computation.
and this causes a determinant to appear in the final transition probability:
T (x̄, x̄0 ) = Tspec (x̄, x̄0 )c−1
Y∂ [xb , xi1 , . . . , xik ]ψ(xi ).
(1 − ψ(xb ))(1 − ψ(xc )))∂ [xb , oi1 , . . . , oik ](12)i=a+1,i6=bHere, Tspec is the transition probability of the purely specular case and the indices i1 , . . . , ik ∈{a + 1, . . . , c − 1} \ {b} refer to glossy vertices that were classified as specular. The terms involving
ψ represent the discrete probability of performing this classification; roughness values may change
during a perturbation, and hence we must account for them to maintain detailed balance.
The determinant in (12) is not hard to compute. Recall that the A-matrix (Section 3.3) maps
perturbations of the vertex positions to changes of the half-vectors. Here, we seek the opposite:
how the non-specular vertices move as a function of o and xb . We therefore compute the matrix
A over the vertex range xa+1 , . . . , xc−1 with one small adjustment: the two rows associated withvertex xb are set so that there is a 2 × 2 identity matrix on the diagonal and zeroes everywhere
else. (Figure 11) We invert this matrix and discard all rows and columns in this inverse thatare associated with specular vertices, and the determinant of the result is the sought change ofvariables factor. Note that when the two chains only consists of glossy vertices, it is equivalent(and considerably faster) to compute the inverse determinant of the block tridiagonal matrix.

6 A unified path space framework for surfaces and volumes
We have presented Manifold Exploration in the context of surfaces, but the extension to volumes
poses no fundamental difficulties. Before we can proceed, however, we require a description of
how volumetric interactions are modeled in the path-space setting.
The purpose of this section is to find an operator-based description of surface and volume lighttransport similar in style to the frameworks presented by Arvo [1995] and Veach [1997], which islater used to derive a path-space integral. Such a unified surface and volume path space has been
proposed in prior work by Pauly et al. [2000], but the context of that work was computational,
and hence it only definitions were supplied.
The goal here is to provide a complete derivation that provides solid theoretical underpinnings

21

for the development of new path-space rendering algorithms. We roughly follow the approachof Chapter 4 in Eric Veach’s Ph.D. thesis [Veach 1997] and conclude with an extended set ofoperators that can describe interactions with surfaces and volumes.
6.0.1 Definitions
It will be convenient to establish some notation: as before, M denotes the set of surfaces. We
define the distance to the next surface as seen by the ray (x, ω) ⊆ R3 × S 2 asdM (x, ω) := inf {d > 0 | x + dω ∈ M}
where inf ∅ = ∞. Based on this distance, we can define a ray-casting function xM :xM (x, ω) = x + dM (x, ω)ω.
We shall also make use of the following notation to express integration along rays
Z xM (x,ω)
Z dM (x,ω)
f (z) dz :=f (x + tω) dt.
x0Note that that while xM (x, ω) is generally undefined for “points at infinity” (dM (x, ω) = ∞), it
has a well-defined meaning when encountered in integrals of the above form. We will also writeZ yf (z) dz :=Z kx−yk→ dt,
f (x + t −
xy)

x

where

0

−
→ := y − x .
xyky − xk

6.1 Light transport model
The radiative transfer equation (RTE) is commonly used to describe the steady-state light transport
on a three-dimensional domain V ⊆ R3 . It is given by the following energy balance equation:(ω · ∇) L(x, ω) + σt (x) L(x, ω) =Z
σs (x)
fp (x, ω 0 → ω) L(x, ω 0 ) dσx (ω 0 ) + Le (x, ω),x ∈ V ◦.

(13)

S2

where σs and σt are the scattering and extinction coefficients, fp is the phase function, L denotes
the volume radiance, and Le represents emitted radiance. dσx denotes differential solid angles at
x. Please note that due to the volumetric setting, the emission term Le has units of radiance perunit length at interior points x ∈ V ◦ .
We will later specify boundary conditions on the set of surfaces M ⊆ V, hence the above
equation only holds on the interior of the domain, i.e. V ◦ := V \ M.
Restricted to a line of direction ω parameterized by s ∈ R, Equation (13) can be rewritten as a
one-dimensional ordinary differential equation of the form
L0 (s) + σt (s)L(s) = Z(s)where Z(s) represents the emitted and in-scattered radiance at s. Given an initial condition at22s = 0, this ODE has the solution
 Z s

Z t


Z s
L(s) = exp −σt (r) dr L(0) +exp
σt (r) dr Z(t) dt
 Z0 s Z0 sZ s0= exp −σt (r) dr L(0) +exp −σt (r) dr Z(t) dt
0

0

(14)

t

where the first term corresponds to attenuated radiance from surfaces and the second term accounts
for in-scattered radiance.
In this solution, the parameterization and boundary condition at s = 0 were arbitrarily chosen.
In general, the boundary condition will be specified by the nearest surface visible along the ray(x, −ω). To formalize this constraint, we will first define the reduced surface radiance at (x, ω) as(
Lo (xM (x, −ω), ω) τ (xM (x, −ω) ↔ x), dM (x, −ω) < ∞Lred (x, ω) :=0,otherwise
where

 Z y

τ (x ↔ y) := exp −σt (z) dzxand Lo represents the surface radiance scattered into direction ω (we will consider this termin more detail later). Now, we can substitute a suitable parameterization into Equation (14)and replace the first term by the reduced surface radiance, which yields the integral form of theradiative transfer equation:
Z dM (x,−ω)L(x, ω) = Lred (x, ω) +
τ (x ↔ x − tω)
0


Z
0
0
0
σs (x − tω) fp (x − tω, ω → ω)L(x − tω, ω ) dσx−tω (ω ) + Le (x − tω, ω) dt
S2

which can be written more compactly using the introduced notation:

Z xM (x,−ω) Z000L(x, ω) = Lred (x, ω) + τ (x ↔ y) σs (y) fp (y, ω → ω)L(y, ω ) dσy (ω ) +Le (y, ω) dy. (15)xS2

6.2 Boundary conditions
By specifying boundary conditions on the surfaces M, we will invariably also introduce discontinuities in the volumetric radiance function L. These discontinuities are problematic when referringto the surface radiance field, since it will generally not be identical on both sides.
Analogous to a one-sided limit in the 1D case, we shall therefore distinguish between the frontand back-facing components L+ (x, ω) and L− (x, ω) as determined by the surface normal N (x)for x ∈ M. Based on this separation, the more intuitive incident and exitant radiance functions

23

can then be defined as(
L+ (x, −ω),Li (x, ω) :=L− (x, −ω),(L+ (x, ω), ω · N (x) > 0
and Lo (x, ω) :=.
L− (x, ω), ω · N (x) < 0ω · N (x) > 0
ω · N (x) < 0

With the help of these definitions, we can introduce the usual boundary condition known as therendering equation [Kajiya 1986]:
Z
Lo (x, ω) =fs (x, ω 0 → ω)Li (x, ω 0 ) dσx⊥ (ω 0 ) + Le (x, ω), x ∈ M(16)
S2

where fs is the BSDF of the surface.
Away from surfaces, the distinction between Li and Lo is technically not necessary (involving
only a change of direction) but we introduce it here for consistency. Specifically, we will set
Li (x, ω) := L(x, −ω) and Lo (x, ω) := L(x, ω) for x ∈ V ◦ and to analogously express the integral
form of the RTE (15) in terms of Li and Lo :Lo (x, ω) = Lred (x, ω)+Z xM (x,−ω) Zτ (x ↔ y) σs (y)fp (y, −ω → ω)Li (y, ω ) dσy + Le (y, ω) dy. (17)
0

0

S2

x

Note the negated ω 0 term in the spherical integral, which has a positive correspondence in thesurface boundary condition (16). This apparent discrepancy is due to the the different parameter
conventions of surface and volume scattering models.3

6.3 Operator notation
Volumetric light transport can be seen as the alternation of two steps: scattering on a surface
or within the medium, followed by propagation and attenuation. Analogous to [Veach 1997] and[Arvo 1995], these two steps can be formulated as linear operators defined on the space of radiancefunctions. The goal of this section is to partition Equations (16) and (17) so that they can be
expressed in this manner.
We define the scattering operator K as one of the in-scattering integrals in Equations (16) or(17) dependent on whether or not x ∈ M:Zfs (x, ω 0 → ω) h(x, ω 0 ) dσx⊥ (ω 0 ),
x∈M

2
S
Z
:=
(Kh)(x, ω)(18)σs (x)fp (x, −ω 0 → ω) h(x, ω 0 ) dσx (ω 0 ), otherwise
S2

On surfaces, K turns incident radiance into outgoing radiance. On the interior of the domain, it
turns incident radiance into outgoing radiance per unit length.
In comparison, the transport operator G has a single definition on the whole domain. Its role
is to transform outgoing radiance per unit length from the volume and outgoing radiance from3

As a somewhat unfortunate consequence of notations in different fields, i n volume scattering models, the incident
direction arguments points towards the scattering location, whereas it points away in the surface case.

24

surfaces into incident radiance.
Z xM (x,ω)(Gh)(x, ω) := τ (x ↔ y) h(y, −ω) dy + τ (x ↔ xM (x, ω)) h(xM (x, ω), −ω)

(19)

x

This leads to interesting differences in comparison to previous work: to match equations (16) and(17), we arrive at an equilibrium equation of the form
Li = G (KLi + Le )(20)

whereas Lo = Le + KGLo was used by [Veach 1997]. The differences arise since volume lighttransport forces us to deal with both radiance and radiance per unit length. To express theequilibrium condition purely in terms of radiance, the order of the G and K operators must thus
be reversed. Assuming invertibility, the solution operator S can now be found:
Li = G (KLi + Le )⇔ (I − GK) Li = GLe⇔ Li = (I − GK)−1 G Le .
|
{z
}
=: S

When the following Neumann series converges in operator norm, the solution operator can also be
expressed as∞
XS=(GK)k G.
(21)
k=0

6.4 Unified path integral formulation
Having defined the necessary operators, we can now proceed to find the desired path integral
formulation. Note that in the following, we will often switch between parameterizing functions interms of directions and positions. The arrow direction indicates the flow of light, which leads tothe following conventions:→
Le (x → y) := Le (x, −xy)→
W (x ← y) := W (x, −
xy)
e

e

→ −
→
fs (x → y → z) := fs (y, −
yx,
yz)
→ −
→
f (x → y → z) := f (y, −
xy,
yz)
p

p

Here, We denotes emitted importance. When dealing with rays that do not intersect any surfaces,
we consider the ray-casting operator xM (x, ω) to return special ‘points at infinity’. This ensuresthat the following is well-defined even if dM (x, ω) = ∞:Le (x → xM (x, ω)) = Le (x, ω).

25

6.4.1 Change of variables
When dealing with light transport on surfaces, it is often convenient to switch between integration
over a sphere and integration over all surfaces of the scene, e.g.
ZZ⊥
f (xM (x, ω)) dσx (ω) =f (y) Gsurf (x ↔ y) dA(y).
S2MThis change of variables involves the geometric term
Gsurf (x ↔ y) = V (x ↔ y) ·|cos θ1 cos θ2 |kx − yk2→ makes with the surface normals N (x) and N (y),xywhere θ1 and θ2 denote the angles that −respectively, and V is a visibility function defined as(
1, if {αx + (1 − α)y | α ∈ (0, 1)} ∩ M = ∅V (x ↔ y) :=0, otherwise
In the volume setting, a similar change of variables is possible: we consider a volume integral overrays radially emanating from a point x. Fubini’s theorem and some manipulation leads toZ xM (x,ω)

Z

f (y) dy dω
S2

Z

Zx∞

f (x + rω)V (x ↔ x + rω) dr dω
2
ZS∞ 0 Z1=f (y)V (x ↔ y) dydt
r2 S 2 (x,r)
0Z
V (x ↔ y)=f (y) dy.
2
R3 kx − yk=

To handle all possible combinations of volume and surface endpoints, we redefine the basic
geometry term as as follows:→ D −
→
Dx −
xy
y yx
G(x ↔ y) := V (x ↔ y) ·kx − yk2
where(|N (a) · ω| , a ∈ M
Da (ω) :=1,
otherwise

6.4.2 Path space measurement integral
Based on the generalized operators G and K, we can proceed to find a path integral formulationof light transport similar to what is done in Chapter 8 of Eric Veach’s thesis [Veach 1997].
We begin by considering a measurement Ij , which is defined as the ray-space convolution of the

26

(j)

importance We

emitted by sensor j with the incident radiance Li (see [Veach 1997] for details):
Ij = hWe(j) , Li i
Z Z=
We(j) (x, ω)Li (x, ω) dσx⊥ (ω) dA(x).
M

(22)

S2

For simplicity, we restrict the ray origins to the boundary, hence volumetric sensors are notsupported. The goal of the path space approach is to express the measurement integral (22) overthe set of transport paths PZIj =fj (x̄) dµ(x̄),
P

where µ is a measure on P, and fj is a contribution weighting function specific to the measurement.
To render the definitions precise, let us define the path space P as the union of all fixed-length
paths formed by concatenating vertices from M and V ◦ according to a configuration vectorc ∈ {0, 1}k , i.e.
∞
[
[
P :=
Pkc
k=1 c∈{0,1}kwhere
Pkc := {x0 x1 · · · xk | (xi ∈ M if ci = 0) and (xi ∈ V ◦ if ci = 1) (i = 1, . . . , k)} .
Using the Lebesgue measures for area and volume on M and V ◦ , we can define a combined productmeasure on P:
(
Z Y
∞
k
X
X
dA(xi ), ci = 0cccµ (D) :=µk (D ∩ Pk ) where µk (D) :=D i=1 dV (xi ), ci = 1k
k=1c∈{0,1}

To find the path space formulation, let us insert the operator form of the equilibrium equation(20) into the measurement integral (22):Z Z
Ij =We(j) (x, ω)(G(KLi + Le ))(x, ω) dσx⊥ (ω) dA(x).
M

S2Expanding the G operator leads to
Z

Z

=
M

S2

"Z

We(j) (x, ω)xM (x,ω)τ (x ↔ y) [KLi + Le ](y, −ω) dyx#+ τ (x ↔ xM (x, ω)) [KLi + Le ](xM (x, ω), −ω) dσx⊥ (ω) dA(x).

27Using a change of variables (c.f. Section 6.4.1), we can turn the two summands into a volume anda surface integral, respectively:
Z

"Z

=
V◦

M

We(j) (x ← y) τ (x ↔ y) G(x ↔ y) [KLi + Le ](y → x) dV (y)#Z

We(j) (x ← y)τ (x ↔ y) G(x ↔ y)[KLi + Le ](y → x) dA(y) dA(x).

+M

Using the fact that KLi = KSLe =
=

∞ Z
X
k=0

"Z
V◦

M

P∞

k
k=1 (KG) Le , this can be rewritten as

We(j) (x ← y) Ḡ(x ↔ y) [(KG)k Le ](y → x) dV (y)#Z

We(j) (x ← y) Ḡ(x ↔ y)[(KG)k Le ](y → x) dA(y) dA(x).

+
M

(23)

where we have defined the attenuated geometric term Ḡ(x ↔ y) := G(x ↔ y)τ (x ↔ y).
In the above geometric sum, the operators G and K are reversed in comparison previous
encounters (e.g. Equation (21)). This makes it possible to perform additional simplifications:consider the concatenated operator KG for y ∈ M:
"ZZ
xM (y,ω 0 )0(KGh)(y → x) =fs xM (y, ω ) → y → xτ (y ↔ z) h(z, −ω 0 ) dz
S2y#+ τ (y ↔ xM (y, ω 0 )) h(xM (y, ω 0 ), −ω 0 ) dσy⊥ (ω 0 ).
As before, we can perform a change of variables to volume and surface integrals:Z=VZ◦

fs (z → y → x) Ḡ(y ↔ z) h(z → y) dV (z)fs (z → y → x) Ḡ(y ↔ z) h(z → y) dA(z)+M

For y ∈ V ◦ , we get
Z(KGh)(y → x) = σsZ+fp (z → y → x) Ḡ(y ↔ z) h(z → y) dV (z)
fp (z → y → x) Ḡ(y ↔ z) h(z → y) dA(z)

V◦

M

28

To unify the two cases, we define the generalized scattering function f¯ as(
σs fp (z → y → x), y ∈ V ◦f¯(z → y → x) :=fs (z → y → x),
y∈M

(24)

which leads to a single definition on the whole domain:Z(KGh)(y → x) =f¯ (z → y → x) Ḡ(y ↔ z) h(z → y) dV (z)VZ◦+f¯ (z → y → x) Ḡ(y ↔ z) h(z → y) dA(z).
M

Inserting increasing powers of KG into Equation (23) leads to a nested integral, which has thefollowing compact representation on path space:ZIj =fj (x̄) dµ(x̄)
(25)
P

where
fj (x̄) = fj (x1 · · · xn ) = Le (x1 → x2 )

"n−1
Y

#
¯
f (xk−1 → xk → xk+1 ) Ḡ(xk−1 ↔ xk )k=2

· Ḡ(xn−1 ↔ xn ) We(j) (xn−1 → xn ).

(26)

The seeming direction reversal in the importance emission function is due to the convention that“→” indicates the flow of light.

7 Manifold exploration for volumes
There are no specular reflections per se in the volume, butfrom a purely mathematical standpoint the phase function of astrongly forward-scattering volume is not unlike reflection froma rough mirror. The fast-varying part of the mirror BRDF is afunction of the half-vector, and hence our method preserves it
during manifold walks. In the medium case, we are interested
in being able to handle highly peaked phase functions that are
Figure 12: Medium constraint
a function of the scattering angle. For this purpose we treat
the scattering angle as analogous to the half vector, introducing the specular manifold constraint−−→ T −−−−→c(xi−1 , xi , xi+1 ) = T (−x−
i−1 xi ) xi xi+1 ( = oi )(27)

where T (v) is a basis for the plane orthogonal to the direction v (Figure 12). This still leaves
one DOF per vertex, which we remove by preserving the distance to the scattering event (i.e.
kxi−1 − xi k = const.). While computing the entries of the constraint Jacobian ∇C, we use thesetwo constraints in place of the previous definition (Section 3.3) whenever a vertex describes a
medium interaction.

29

Figure 13: It is possible to extend the path space specular integration framework to volumes, forinstance to render caustics from refractive objects, such as this dodecahedron-shaped
luminaire with tinted glass inlays.

7.1 Medium manifold perturbation
From an algorithmic perspective, manifold exploration for volumes is almost identical to thesurface case. Our implementation handles both cases jointly, which permits constructing specular
chains that contain both surface and medium interactions.
Apart from the new type of constraint (27), the computation of offset manifold (11) tangent
vectors is unchanged. In the ray tracing step 4 of WalkManifold when encountering a vertex
xi−1 that is followed by a medium interaction vertex xi , we set x0i = xi−1 +kxi−1 −xi k d, where d isthe outgoing direction at xi−1 (this enforces the length constraint mentioned earlier). Afterwards,
the manifold offset oi is transformed into an outgoing direction in the new frame at x0i (Figure 12).
As in the glossy surface case, we require a criterion that clarifies when treating a medium vertex
xi as non-specular is in order, and when it is better handled by the manifold. Again, this decision
is made probabilistically, based on a modified specular probability function ψ(xi ) described in the
appendix.
When computing the generalized geometric term through a specular chain with medium
endpoints, we consider them to be parameterized on a surface perpendicular to the direction
facing the chain.

8 Results
We have implemented the proposed technique and prior work as extension modules to the Mitsuba
renderer [Jakob 2010]. All techniques operate on top of a newly added bidirectional abstraction layerthat exposes cameras, light sources, scattering models, and participating media as generalized pathvertex and edge objects with a common basic interface. This greatly simplified the implementation
effort, as bidirectional rendering algorithms can usually be stated much more succinctly in termsof operations on vertices that are oblivious to whether they contain e.g. a camera model or a
medium scattering event.

30

SceneSeed generator

Perturb.

samples/px. chains/px.
Torus (MEPT)
Chandelier (MEPT)
Chandelier (MEMLT)
Table (MEPT)
Table (MEMLT)GlassEgg (MEPT)
GlassEgg (MEMLT)

32
64
—
32
—
128
—

2
3
—
1
—
2
—

λ
100
160
160
300
300
90
90

Mutations

Manifold walks

total accepted

total

1178 M
1216 M
3626 M
1074 M
1921 M
1533 M
2774 M

Manif. size

conv. avg. iter. avg.

78.3% 1002 M 96.7%
73.6%
975 M 97.6%
34.1%
376 M 98.4%
77.5%
868 M 95.5%
35.9% 772.4 M 94.3%
72.4% 1246 M 92.2%
40.3% 1101 M 92.5%

2.3
2.4
2.1
2.8
3.2
3.4
3.3

3.4
4.6
4.5
4.4
4.4
4.1
4.1

max.
7
13
13
14
14
14
14

Table 1: Listing of seed generator and perturbation parameters, as well as performance statistics.
We compare the following algorithms:• Primary sample space MLT by Kelemen et al., implemented on top of bidirectional pathtracing (PSSMLT).
• Path space MLT by Veach and Guibas (MLT).
• An extended form of energy redistribution path tracing by Cline et al. (ERPT), which isseeded by bidirectional rather than unidirectional path tracing. The ERPT implementation
shares the caustic, lens, and multi-chain perturbation with the previous algorithm. Since
they introduce bias, we did not use the post-processing filters proposed in the original paper.
• Manifold exploration path tracing (MEPT), which is structured similarly to ERPT. We
modified the original algorithm by replacing its highly specialized caustic, lens, and multichain perturbations with the manifold perturbation. Due to its general design, the new
perturbation subsumes and extends the capabilities of the original set.
• Manifold Metropolis Light Transport (MEMLT), which corresponds to MLT with ourperturbation (i.e. the bidirectional mutation and manifold perturbation, but none of theoriginal perturbations from MLT).
Due to the aforementioned abstraction layer, all techniques transparently support participating
media even if this was originally not part of their description. We found that MEPT generally
performs better that MEMLT due to certain limitations of the bidirectional mutation that are
discussed later.
The rendering of result images was conducted on Amazon EC2 cc1.4xlarge instances, which
are eight-core Intel Xeon X5570 machines. A single machine was used per image. To exploit
the local parallelism, our implementation runs a separate Markov chain on each core, and theresulting buffers are averaged together when exposing the image.
We have rendered three views of a challenging interior scene containing approximately 2 million
triangles with shading normals and a mixture of glossy, diffuse, and specular surfaces and somescattering volumes. One hour of processing time was allocated to each rendering technique,and a comparison of the resulting images is shown in Figures 14, 15, and 16. The converged
images in Figure 1 were rendered using 48 hours. The one hour renderings are intentionallyunconverged to permit a visual analysis of the convergence behavior. Table 1 lists parametersand statistics collected during these renderings. The path generator columns refer to the seeding
scheme used by ERPT and MEPT, which samples and subsequently resamples a number of paths

31per pixel before launching Markov Chains. The statistics include the total number of mutations
and acceptance ratio, as well as the the convergence behavior of the manifold walks and vertex
count of encountered manifolds.
Chandelier: In this set of results, the poor performance of MLT is most apparent and iscaused by the ineffectiveness of the bidirectional mutation in finding long specular paths. Because
it must decide up front on the configuration of a path before generating it, most of the time themutation fails, resulting in acceptance rates under 1%. Consequently, too few jumps between
disjoint connected components of path space occur, causing parts of the image have an incorrect
relative brightness. This weakness is inherited by MEMLT, which also builds upon the bidirectional
mutation. It is more successful at exploring some of the diffuse-specular-diffuse paths throughthe bulbs but overall does not work well on this scene. Seeding the same perturbations usingbidirectional path tracing, which does not suffer from this disadvantage, performs much better, as
can be seen in the ERPT and MEPT renderings.
Table: This scene is lit by the chandelier, with its glass-enclosed sources, so all illumination isby specular paths. By reasoning about the geometry of the specular and offset specular manifolds
for the paths it encounters, our perturbation rule is more successful at rendering paths—such as
illumination that refracts from the bulbs into the butter dish, then to the camera (6 specularvertices)—that the other methods struggle with. The MLT rendering looks too dark, because it
did not find enough of these paths and mainly captures diffuse illumination from the walls. The
noise in the ERPT result reveals that the underlying bidirectional path tracer encountered someof those paths but the Veach–Guibas perturbations are not able to explore path space aroundthem effectively. The primary sample space MLT variant also has difficulties rendering this scene,
because it has no knowledge about the underlying path geometry. MEMLT produces a clean
result, but the relative brightness of different parts of the image is incorrect due to bidirectional
mutation’s difficulty in performing sufficiently many jumps between them.
GlassEgg: In this scene, our technique’s ability to create a specular chain containing both
medium and surface interactions leads to fast convergence when rendering the forward-scattering
medium inside the glass egg. MLT and ERPT perform poorly here, since they do not havesuitable perturbations for exploring this space. Because the MLT perturbations treat glossy anddiffuse materials identically, they have difficulty rendering the near-specular tabletop, producing
streak-like artifacts in the output rendering.

9 Conclusions
We have presented a new type of Markov Chain Monte Carlo rendering method that models
the space of valid specular and near-specular light paths using high-dimensional manifolds. The
differential geometry of these manifolds provides a powerful tool to efficiently explore these paths,
which can be very hard for previous methods to find. Our technique applies in the frameworks of
MLT or ERPT, producing rendering algorithms with support for specular paths fundamentally
built in at the core. In equal-time comparisons on very challenging scenes, the new Manifold
Exploration Path Tracing algorithm compares favorably to previous work in Monte Carlo and
MCMC rendering.
This new algorithm does still share certain limitations with its predecessors. Most important,it needs well distributed seed paths, because it can only explore connected components of the

32

MLT

ERPT

PSSMLT

Figure 14: Chandelier: This view contains a brass chandlier with 24 light bulbs, each surrounded by a glass enclosure. The chandlier uses a realistic metal material based on
microfacet theory and is attached to the ceiling using specular metal cylinders. This
scene is challenging, as certain important light paths are found with low probability,
particularly those involving interreflection between the bulbs and the body of thechandelier. In this and the following comparisons, one hour of processing time wasallocated to each rendering technique.
33

MEMLT

MEPT

MEPT (48h)Figure 14: Chandelier (continued)

34

MLT

ERPT

PSSMLT

Figure 15: Table: This view of our room scene shows chinaware (using a BRDF with both
diffuse and specular components), a teapot containing an absorbing medium, anda butter dish on a glossy silver tray. Illumination comes from the chandelier inFigure 14.

35

MEMLT

MEPT

MEPT (48h)Figure 15: Table (continued)

36

MLT

ERPT

PSSMLT

Figure 16: GlassEgg: This view of a different part of the room, now lit through windows using
a spherical environment map surrounding the scene, contains a forward-scattering
participating medium (g = 0.8) inside the glass egg.
manifold for which seed paths are provided. Bidirectional Path Tracing is reasonably effective
but still has trouble finding many components, and this problem fundamentally becomes more
and more difficult as the number of path types increases. Ultimately, as the number of different
path types exceeds the number of samples that can be generated, local exploration of path spacebecomes ineffective; future algorithms could be designed to attempt exploration only in sufficiently
large path space components.
Unlike many methods for caustics and other specular phenomena, we have shown how to
generalize Manifold Exploration almost trivially to handle glossy surfaces and volumes. Similarrefinements can let the same method handle perfectly anisotropic reflections, strongly oriented
volume scattering media, and other kinds of problems with exactly or approximately constrained
paths.
While MCMC rendering is a natural match for our methods of dealing with specular paths,
our predictor-corrector iteration can be used in other kinds of algorithms as well, including
deterministic ones to map out specular paths, for instance in design of luminaires or opticalsystems.

37

MEMLT

MEPT (48h)

MEPT

Figure 16: GlassEgg (continued)

10 Acknowledgements
This research was conducted in conjunction with the Intel Science and Technology Center for
Visual Computing. Additional funding was provided by the National Science Foundation under
grant IIS-1011919. The authors are indebted to Olesya Isaenko, who meticulously crafted theexample scenes used in the evaluation.

Appendix A1Specular probability function for surfaces: We found the following following heuristic basedon microfacet theory to work well: when the microsurface normals at xi follow a distribution
Dα (mi ) with roughness parameter α, the BSDF at xi will take on small values when mi moves
into a region where Dα (mi ) has low density. We thus set ψ(xi ) by computing the expected
probability that treating vertex xi as non-specular during a manifold perturbation would move its
microsurface normal mi from a region of high density to one of low density, and we choose the90th-percentile to classify the support of Dα into such regions.

38

perturbed by

90%
percentile

90%
percentile

To obtain the specular probability, our implementations must know the expected angular change∆θ of microsurface normals during a perturbation, which is found by briefly running the Markov
chain before rendering starts. During rendering, ψ(xi ) is computed as the area ratio of the twohighlighted regions on the sphere:ψ(xi ) =1 − cos θq (α(xi )),1 − cos (θq (α(xi )) + ∆θ)(28)

where θq is the aforementioned percentile (with q set to 0.9). For the Beckmann distribution, thisis given byθq (α) := tan−1 (−α2 log(1 − q)).
Implementation-wise, this heuristic requires material models to be able to compute their Beckmann
distribution-equivalent roughness or provide a custom quantile function.
Specular probability function for participating media: In the medium case, we use the sameprobability (28), but now with a percentile that is suitable for volumetric scattering. We useθq (g) = cos−1(1+|g|)2 −2(1+|g|)(1+g 2 )q+2|g|(1+g 2 )q 2
(1+|g|−2|g|q)2

where g is the mean cosine of the phase function, q is set to 0.5, and θq was derived from theHenyey-Greenstein phase function. In Ψ(xi ) (Equation 28) we must also replace ∆θ with theaverage change in scattering angle at medium vertices, again determined in a brief phase before
rendering.

Appendix A2
Let us assume the existence of a Vertex data structure with the following contents:
struct Vertex {Point p;Vector dpdu, dpdv;
Normal n;
Vector dndu, dndv;
float eta;
Matrix2x2 A, B, C;// Vertex position// Tangent vectors// Normal vector// Normal derivatives// Relative index of refraction// Matrix blocks of ∇C in the row// associated with the current vertex};

39

Then the following C++ code computes the derivatives of ∇C associated with a single surface
interaction vertex. It assumes that the function is called with a pointer into a list of vertices.
void computeDerivatives(Vertex *v) {/* Compute relevant directions and a few useful projections */Vector wi = v[-1].p - v[0].p;
Vector wo = v[ 1].p - v[0].p;
float ili = 1/wi.length();float ilo = 1/wo.length();wi *= ili; wo *= ilo;
Vector H = wi + v[0].eta * wo;
float ilh = 1/H.length();H *= ilh;
float dot_H_n= dot(v[0].n, H),
dot_H_dndu = dot(v[0].dndu, H),
dot_H_dndv = dot(v[0].dndv, H),
dot_u_n= dot(v[0].dpdu, v[0].n),
dot_v_n= dot(v[0].dpdv, v[0].n);/* Local shading tangent frame */Vector s = v[0].dpdu - dot_u_n * v[0].n;
Vector t = v[0].dpdv - dot_v_n * v[0].n;
ilo *= v[0].eta * ilh; ili *= ilh;/* Derivatives of C with respect to x_{i-1} */VectordH_du = (v[-1].dpdu - wi * dot(wi, v[-1].dpdu)) * ili,
dH_dv = (v[-1].dpdv - wi * dot(wi, v[-1].dpdv)) * ili;
dH_du -= H * dot(dH_du, H);
dH_dv -= H * dot(dH_dv, H);
v[0].A = Matrix2x2(
dot(dH_du, s), dot(dH_dv, s),dot(dH_du, t), dot(dH_dv, t));/* Derivatives of C with respect to x_i */dH_du = -v[0].dpdu * (ili + ilo) + wi * (dot(wi, v[0].dpdu) * ili)
+ wo * (dot(wo, v[0].dpdu) * ilo);
dH_dv = -v[0].dpdv * (ili + ilo) + wi * (dot(wi, v[0].dpdv) * ili)
+ wo * (dot(wo, v[0].dpdv) * ilo);
dH_du -= H * dot(dH_du, H);
dH_dv -= H * dot(dH_dv, H);
v[0].B = Matrix2x2(
dot(dH_du, s) - dot(v[0].dpdu, v[0].dndu) * dot_H_n - dot_u_n * dot_H_dndu,
dot(dH_dv, s) - dot(v[0].dpdu, v[0].dndv) * dot_H_n - dot_u_n * dot_H_dndv,
dot(dH_du, t) - dot(v[0].dpdv, v[0].dndu) * dot_H_n - dot_v_n * dot_H_dndu,

40dot(dH_dv, t) - dot(v[0].dpdv, v[0].dndv) * dot_H_n - dot_v_n * dot_H_dndv);/* Derivatives of C with respect to x_{i+1} */dH_du = (v[1].dpdu - wo * dot(wo, v[1].dpdu)) * ilo;
dH_dv = (v[1].dpdv - wo * dot(wo, v[1].dpdv)) * ilo;
dH_du -= H * dot(dH_du, H);
dH_dv -= H * dot(dH_dv, H);
v[0].C = Matrix2x2(
dot(dH_du, s), dot(dH_dv, s),dot(dH_du, t), dot(dH_dv, t));}References
Arvo, J. R. 1995. Analytic Methods for Simulated Light Transport. PhD thesis, Yale University.
Chen, M., and Arvo, J. 2000. Perturbation methods for interactive specular reflections. IEEE
Trans. Vis. and Comp. Graph. 6, 3 (July/Sept.), 253–264.
Chen, M., and Arvo, J. 2000. Theory and application of specular path perturbation. ACM
Trans. Graph. 19, 4 (Oct.), 246–278.
Chen, J., Wang, B., and Yong, J.-H. 2011. Improved stochastic progressive photon mappingwith metropolis sampling. Computer Graphics Forum 30, 4, 1205–1213.
Cline, D., Talbot, J., and Egbert, P. 2005. Energy redistribution path tracing. ACM Trans.
Graph. 24, 3 (Aug.), 1186–1195.
Cook, R. L., Porter, T., and Carpenter, L. 1984. Distributed ray tracing. In Computer
Graphics (Proceedings of SIGGRAPH 84), 137–145.
Goral, C. M., Torrance, K. E., Greenberg, D. P., and Battaile, B. 1984. Modeling theinteraction of light between diffuse surfaces. In Computer Graphics (Proceedings of SIGGRAPH84), 213–222.
Hachisuka, T., and Jensen, H. W. 2009. Stochastic progressive photon mapping. ACM Trans.
Graph. 28, 5 (Dec.).
Hachisuka, T., and Jensen, H. W. 2011. Robust adaptive photon tracing using photon pathvisibility. ACM Trans. Graph. 30, 5 (Oct.), 114:1–114:11.
Hastings, W. K. 1970. Monte carlo sampling methods using markov chains and their applications.
Biometrika 57, 1, 97–109.
Heckbert, P. S. 1990. Adaptive radiosity textures for bidirectional ray tracing. In Computer
Graphics (Proceedings of SIGGRAPH 90), 145–154.

41

Igehy, H. 1999. Tracing ray differentials. In Computer Graphics (Proceedings of SIGGRAPH99), 179–186.
Jakob, W., 2010. Mitsuba renderer. http://www.mitsuba-renderer.org.
Jarosz, W., Zwicker, M., and Jensen, H. W. 2008. The beam radiance estimate forvolumetric photon mapping. Computer Graphics Forum 27, 2 (Apr.), 557–566.
Jensen, H. W., and Christensen, P. H. 1998. Efficient simulation of light transport inscenes with participating media using photon maps. In Computer Graphics (Proceedings of
SIGGRAPH 98), 311–320.
Jensen, H. W. 1996. Global illumination using photon maps. In Eurographics Rendering
Workshop 1996, 21–30.
Kajiya, J. T. 1986. The rendering equation. In Computer Graphics (Proceedings of SIGGRAPH86), 143–150.
Kelemen, C., Szirmay-Kalos, L., Antal, G., and Csonka, F. 2002. A simple and robust
mutation strategy for the metropolis light transport algorithm. Computer Graphics Forum 21,
3, 531–540.
Kitaoka, S., Kitamura, Y., and Kishino, F. 2009. Replica exchange light transport. Computer
Graphics Forum 28, 8 (Dec.), 2330–2342.
Lafortune, E. P., and Willems, Y. D. 1993. Bi-directional path tracing. In Proceedings ofCompugraphics 93.
Lai, Y.-C., Fan, S. H., Chenney, S., and Dyer, C. 2007. Photorealistic image rendering withpopulation monte carlo energy redistribution. In Rendering Techniques 2007: 18th Eurographics
Workshop on Rendering, 287–296.
Liu, J. S. 2001. Monte Carlo strategies in scientific computing. Springer.
Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller,
E. 1953. Equation of state calculations by fast computing machines. J. Chem. Phys. 21, 6.
Mitchell, D. P., and Hanrahan, P. 1992. Illumination from curved reflectors. In Computer
Graphics (Proceedings of SIGGRAPH 92), 283–291.
Pauly, M., Kollig, T., and Keller, A. 2000. Metropolis light transport for participating
media. In Rendering Techniques 2000: 11th Eurographics Workshop on Rendering, 11–22.
Pharr, M., and Humphreys, G. 2004. Physically Based Rendering: From Theory to Implementation. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA.
Premoze, S., Ashikhmin, M., and Shirley, P. 2003. Path integration for light transport involumes. In Eurographics Symposium on Rendering: 14th Eurographics Workshop on Rendering.
Segovia, B., Iehl, J., and Péroche, B. 2007. Metropolis instant radiosity. Computer Graphics
Forum 26, 3 (Sept.), 425–434.

42

Shirley, P. S., Wade, B., Hubbard, P., Zareski, D., Walter, B., and Greenberg, D. P.
1995. Global illumination via density estimation. In Eurographics Rendering Workshop 1995,
219–231.
Sillion, F. X., and Puech, C. 1989. A general two-pass method integrating specular anddiffuse reflection. In Computer Graphics (Proceedings of SIGGRAPH 89), 335–344.
Spivak, M. 1965. Calculus on Manifolds. Addison-Wesley.
Veach, E., and Guibas, L. 1994. Bidirectional estimators for light transport. In Fifth
Eurographics Workshop on Rendering.
Veach, E., and Guibas, L. J. 1997. Metropolis light transport. In Computer Graphics(Proceedings of SIGGRAPH 97), 65–76.
Veach, E. 1997. Robust Monte Carlo Methods for Light Transport Simulation. PhD thesis,
Stanford University.
Walter, B., Hubbard, P. M., Shirley, P. S., and Greenberg, D. F. 1997. Global
illumination using local linear density estimation. ACM Trans. Graph. 16, 3 (July), 217–259.
Walter, B., Marschner, S. R., Li, H., and Torrance, K. E. 2007. Microfacet models forrefraction through rough surfaces. In Rendering Techniques 2007: 18th Eurographics Workshop
on Rendering, 195–206.
Walter, B., Zhao, S., Holzschuch, N., and Bala, K. 2009. Single scattering in refractive
media with triangle mesh boundaries. ACM Trans. Graph. 28, 3 (July), 92:1–92:8.
Whitted, T. 1980. An improved illumination model for shaded display. Communications of theACM 23, 6 (June), 343–349.

43