Manifold Exploration: A Markov Chain Monte Carlo techniquefor rendering scenes with difficult specular transport
Wenzel Jakob

Steve Marschner

Cornell University

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

Abstract

1

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 ideathat sets of paths contributing to the image naturally form manifolds
in path space, which can be explored locally by a simple equationsolving 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 thismethod in in two different Markov Chain Monte Carlo frameworks
to accurately compute illumination from general families of 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 rangeof challenging scenes and evaluated against previous methods.

Certain classes of light paths have traditionally been a source ofdifficulty in conducting Monte Carlo simulations of light transport.
A well-known example is specular-diffuse-specular paths, such atabletop seen through a drinking glass sitting on it, a bottle containing shampoo or other translucent liquid, or a shop window viewed
and illuminated from outside. Even in scenes where these paths do
not cause dramatic lighting effects, their presence can lead to unusably slow convergence in renderers that attempt to account for alltransport paths.

CR Categories: I.3.7 [Computer Graphics]: Three-Dimensional
Graphics and Realism
Keywords: Rendering, Specular manifold, MCMC
Links:

DL

PDF

W EB

V IDEOC ODEIntroduction

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 naturalgeneralizations 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 asimple 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 thatis already generally available in ray tracing systems. This solution method is then used to define a Markov Chain, which has theright properties to be applied in a Metropolis-Hastings algorithm
to compute lighting through very general families of paths that caninvolve specular, near-specular, glossy, and diffuse surface interactions as well as isotropic or highly anisotropic volume scatteringinteractions.
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-specular materials. 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 path 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 wefinish by showing results and comparisons in Section 7.

2

Prior work

Simulating light transport has been a major effort in computer
graphics for over 25 years, beginning with 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 [Cook et 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 inthe 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 pathtracing [Lafortune and Willems 1993; Veach and Guibas 1994], inwhich light transport paths can be constructed partly from the lightand 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 sources in 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 thenrenders the image using ray tracing, making use of the stored particles to estimate illumination by density estimation. [Jensen andChristensen 1998; Jarosz et al. 2008; Hachisuka and Jensen 2009].
Photon mapping and other two-pass methods are characterized bystoring an approximate representation of some part of the illumination in the scene, which requires assumptions about the smoothnessof illumination distributions. On one hand, this enables renderingof some modes of 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 suffice under assumptions of smoothness. However, this smoothness
assumption inherently leads to smoothing 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, ischallenging to handle with photon maps, since large numbers ofphotons 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 the number 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 properties of specular reflection paths.
Mitchell and Hanrahan [1992] devised a method to compute irradiance from implicitly defined reflectors, using Fermat’s principlewith interval Newton’s method to locate all reflection paths from asource to a point. Walter et al. [2009] proposed a related methodthat computes 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 the neighborhood of a given path, it avoids the complexities and constraints
entailed by a full global search. Another difference is that our manifold formalism can be used to build a fully general rendering systemthat 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 ofreflected paths—in this case, paths from the eye. Igehy’s approachrequires elementary local differential information only, in the formof derivatives of surface normals, and does not require global surface descriptions as Mitchell and Hanrahan’s method does. Manifold exploration requires the same geometric information, and canthus be implemented in most modern ray tracing sytems.
The analysis of reflection geometry presented by Chen andArvo [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 optical length (i.e. the propagation time) constitutes a local extremum amongst neighboring paths. Using Lagrange
multipliers, the authors derive a path Jacobian and path Hessianwith respect to perturbations of the endpoint of a path and use it toaccelerate the interactive display of reflections 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 thisderivative and the goals of the research are entirely different: intheir case, estimating changes to viewing paths, and in our case,tracking the evolution of specular paths in a very general context,
as part of an unbiased 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 the possibility of using MCMC to integrate over multiplescattering paths in volumes, as discussed in Section 6.

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 spacein which the probability of a state appearing at a given position in
the sequence depends only on the previous 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 timethe distribution of states it visits will be proportional to the desired
distribution. The following brief introduction follows [Liu 2001].
Defining a Markov chain amounts to defining a transition rule: aprocess for selecting a new state x+ randomly, in a way that depends on the current state x. Metropolis et al. provided a way totake 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 Metropolis transition 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 otherwise
In step 1 we say x0 is proposed as the next state, and in step

2 it is either accepted and becomes the 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 anotherusing a finite expected number of steps, and if the length of thisjourney is in a sense “irregular” (i.e. aperiodic), it is referred toas ergodic, and the chain’s limiting distribution is guaranteed toconverge to π. This relatively mild criterion is usually satisified bychoosing 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)(1)r(x, x0 ) = 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 in the space are paths, and the desired
probability distribution over paths is proportional to their contribution 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 implementationof a transition rule, and any rule with symmetric support is admissible. But to avoid very low acceptance probabilities, which leadto poor performance, it is desirable for the transition probability toapproximate the contribution: 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 thepath and perturbations that move the vertices by small distances
while preserving the structure, both using the building blocks ofbidirectional 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 basedon changing the random numbers used to sample paths by bidirectional path tracing.
Considerable research activity has extended Metropolis light transport in various ways. Pauly et 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 Carlo rendering [Lai et al. 2007], andReplica Exchange light transport [Kitaoka et al. 2009]. Recently,
two groups [Chen et al. 2011; Hachisuka and Jensen 2011] havecombined the transition rule of Kelemen et al. with photon mapping to obtain robust methods based on density estimation.

distribution of samples even if the underlying 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 near-specular illumination paths to therendered image of a scene. The general approach applies to surfaces
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 setting of surfaces without
participating media. In the following section we develop this theory into a rendering method for scenes combining ideal specularsurfaces with fairly diffuse surfaces. We will then go on to generalize the method to glossy surfaces, then to generalize the theory toencompass 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, described by Veach [1997] and others, inwhich the value of each pixel in the image is an integral of a contribution function over path space. Denoting the union of all surfaces
in the scene by M, each sequence x1 . . . xn of at least two points
in M is a path along which light may travel. Thus path space 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 productmeasure derived from the area measure on M. The contributionfunction fj is a product of terms, one for each vertex and each edgeof the path:fj (x1 · · · xn ) = Le (x1 → x2 )"n−1#YG(xk−1 ↔ xk ) fs (xk−1 → xk → xk+1 )
k=2

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

(4)

(j)

However, to generate proposals, all of these algorithms ultimately
rely on local path sampling strategies (i.e. path tracing). Specifically, they choose the next interaction vertex along a light path bysampling from a directional distribution associated with the current
vertex, followed by an intersection search. Our method introduces
a new kind of transition rule with different properties.

Here Le is emitted radiance, We is importance for the jth pixel,
fs is the BSDF at xk for the geometry defined by xk−1 , xk , andxk+1 , and G is the geometry factor:→ |N (y) · −→|N (x) · −xy|yx|G(x ↔ y) =V (x ↔ y)(5)kx − yk2

The original MLT algorithm and subsequent variants all render animage by running a Markov chain for a long (e.g. > 106 ) sequenceof steps, and they guarantee ergodicity using a transition rule thatcan 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 aninteresting departure from that approach. It draws on the property
that the Metropolis-Hastings algorithm (1) preserves the stationary

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 and Metropolis 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 somesurface interactions are specular, the entire contribution to the pathspace integral is from paths that obey specular reflection or refraction geometry, and the set of such paths is lower in dimension than
the full path space. For instance, consider a family of paths of theform 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−→+−→ k N (x ),
x−3−
x
x−3−
x
2
4
3
that is, the half-vector at x3 is in the direction of the normal. This
places two constraints on the path, meaning that all contributing
paths lie on a manifold of dimension 8 embedded in P5 , which isof dimension 10. The integral is more naturally expressed as anintegral 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 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 reflectancevalues, and the 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 theother vertex, and the generalized geometry factor is defined 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 the endpoints. Figure 3
illustrates this for a more complex path involving a chain of three
specular vertices. We will explain below how G can be easily computed from the differential geometry of the specular manifold.
This generalized geometry factor is related to the “extended form
factor” discussed by Sillion and Puech [1989].

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 classification of each of its vertices. (In thisscheme 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 andprojects it into the shading tangent space; the resulting 2-vector iszero 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:

S

D

S

Figure 3: The geometry factor (left) and the generalized geometry
factor (right) are both derivatives of projected solid angle at one
end with respect to area at the far end.
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).
Specular endpoint vertices also introduce constraints that depend on
their type: for instance, the position of a vertex on a point emittermust remain fixed (e.g. x1 = const). When the emitter is 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 small neighborhoods,
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 , andthe specular manifold is simply the set
S = {x̄ | C(x̄) = 0}(8)

Expressing S using a constraint in this way makes it convenient towork with neighborhoods of a particular path. The Implicit Function Theorem [Spivak 1965] guarantees the existence of a parameterization of the manifold, in the neighborhood of any path x̄ thatis nonsingular (in the sense explained below). This parameterization is a function q : IR2(n−p) → IR2p that determines the positionsof all the specular vertices from the positions of all the nonspecular vertices. Furthermore, the derivative of q, which gives us thetangent space to the manifold at x̄, is simple to compute from thederivative of C.
For the specifics we restrict ourselves to the case of a single chainof specular vertices with non-specular vertices (surfaces, cameras,
or light sources) at the ends, which suffices to cover most cases
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 theremaining k − 2 vertices being specular. In this case C : IR2k →IR2(k−2) , and the derivative ∇C is a matrix of k − 2 by k 2-by-2
blocks, with a tridiagonal structure (Figure 2).

−−→ −−−−→
ci (xi−1 , xi , xi+1 ) = T (xi )T h(xi , −x−
i xi−1 , xi xi+1 )(6)

The Implicit Function Theorem gives us a parameterization of the
manifold in terms of any 2 vertices, and if we pick x1 and xk thissimply says that the path, in a neighborhood of the current path1 ,is a function of the two endpoints. Furthermore, it also tells us thederivative of that parameterization, which is to say, the derivativeof all the specular vertices’ positions with respect to the positionsof the endpoints. If we block the derivative ∇C, as shown in theη(x, v)v + η(x, w)wh(x, v, w) =kη(x, v)v + η(x, w)wk

(7)

1 Because it is possible to have several separated specular paths joining

two points, the parameterization cannot be global.

D

where

S
D

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

Figure 2: 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.
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 isTS (x̄) = −A−1 B1

Bk

These theoretical results about the structure of the specular manifold can be used in an algorithm to solve for specular paths, which
we discuss next.

4

This matrix is k − 2 by 2 blocks in size, and each block gives thederivative of one vertex (in terms of its own tangent frame) withrespect 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 other endpoint.
We use TS (x̄) for two things: to navigate on the manifold and tocompute the generalized geometry factor. The right two or left twocolumns of TS (x̄) are useful for updating the specular chain withx1 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-right block gives theratio of an infinitesimal area at xk to its reflection/refraction, as observed from x1 , measured on the surface at x2 . To convert this to aratio of area at xk to solid angle at x1 , we multiply this determinantby the ordinary geometry factor G(x1 ↔ x2 ); this product is thegeneralized geometry factor G(x1 ↔ · · · ↔ xk ). More succinctly,

Walking on the specular manifold

Our rendering algorithms use MCMC to explore the manifold ofspecular paths, and for this they 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 endpointsof a specular chain and takes all the intermediate vertices to a valid
new configuration. Later, in Section 5, we show how to apply thisalgorithm to render images.
To simplify the discussion, we will focus on the case where theposition of a vertex xn of a specular chain x1 , . . . xn is adjusted
to a given new position x0n , while x1 is held fixed. We shall also
briefly introduce the assumption that xn is located on planar surface
of infinite extent.
Our manifold walking algorithm is based on two key insights:(9)

1. The A and B matrices (Section 3.3) may be used to map an
infinitesimal in-plane movement of xn to displacements ofthe vertices x2 , . . . , xn−1 . We can use these displacements
to approximate a finite change to the path simply by adding
an offset to each vertex, but this will move the path off the
specular manifold.

where Pi is a 2 by 2(k − 2) matrix that projects onto the two dimensions associated with vertex i. Please see the supplementary
technical report on how to compute TS (x̄) and the generalized geometry factor for a simple example path.

2. Ray tracing provides a deterministic means of projecting an
off-specular path back onto the space of valid configurations
based on its first two vertices. Given x1 and x2 , we can trace
a sequence of rays xi → xi+1 , at each step performing a specular reflection or refraction, and this leads to corrected posi+tions x+2 . . . xn .

G(x1 ↔ · · · ↔ xk ) = P2 A−1 Bk G(x1 ↔ x2 )−1

= Pk−1 A

B1 G(xk−1 ↔ xk ),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, we must have access to the partial derivatives of position and shading normal withrespect to any convenient parameterization of the surfaces, alongwith the refractive indices of all objects. These are 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.
2 If the parameterizations are not orthonormal, two additional determi-

nants are required to account for the change in area.

By combining 1. and 2., we obtain a predictor-corrector type algorithm (Figure 4) that performs a step according to a local linear
model, followed by a projection that restores the specular config+
uration, resulting in a new path x1 , x+
2 , . . . , xn , and repeats until convergence. As long as the prediction step solves the linear
model and moves in the tangent space to the manifold, this iteration
behaves like Newton’s method, exhibiting quadratic convergencenear the solution. As with all Newton-like iterations, it is not guaranteed to converge when started far from the solution, since thelinear model may not be accurate enough to make progress. But
since the model is first-order accurate, the algorithm is guaranteed
to make forward progress when the constraint function is differentiable and the partial steps are small enough. Our algorithm uses
a simple heuristic to decrease the step size when progress is notmade, then increase back to full steps to get quadratic convergence(a)path space(fixed)starting pointnext iterategoal point(b)corrected

predicted
starting
point

path space

Figure 4: 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 tothe manifold, then (a) a nearby valid specular path is found, which(b) corresponds to projecting back onto the manifold.
as it approaches the target configuration. This iteration is illustrated
in Figure 4 and laid out in the following algorithm:
WALK M ANIFOLD x1 , . . . , xn

x0n )Set i = 0 and β = 1
while kxn − x0n k > εLp = x2 − β T (x2 )P2 A−1 Bk T (xn )T (x0n − xn )Propagate the ray x1 → p through all specular+interactions, producing x+
2 , . . . , xn .
+05
if step 4 succeeded and kxn − xn k < kxn − x0n k+
6x2 , . . . , xn = x+
2 , . . . , xn
7
β = min {1, 2 β}8else
9β = 12 β
10
Set i = i + 1, and fail if i > N .
11 return x2 , . . . , xn−11234Here, i records the number of iterations until a specified maximum
N is reached, and ε is a relative 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 ensurethat steps reduce the target distance. The pseudocode assumes thatthe 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 ona plane. However, the actual locations of x+n along the way fromxn 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 WALK M ANIFOLD the last propagation step
computes an intersection against this hypothetical plane, ignoring
the actual scene geometry. Once the algorithm converges, we mustensure that xn−1 and xn are mutually visible 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 forwalks xn x0n and x0n xn .
To compute the matrix A, we derived symbolic expressions for the
manifold constraints in C (Equation 8). A C++ implementation isprovided in the supplementary technical report. When solving theresulting 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 timecomplexity from O(n3 ) to O(n), n being the number of vertices in
the chain.
There are several situations under which this algorithm may fail toconverge: first, a specular path between x1 and x0n need not exist
at all. Secondly, WALK M ANIFOLD usually cannot find paths thatlie on a different connected component of the manifold. Thirdly,
when the local structure of the manifold is complex (e.g. due tohigh-frequency geometric detail of the reflectors and refractors) andkxn −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 forthe iteration to fail occasionally, as explained in the next section.
The manifold walking algorithm works reliably for large chainswith over 10 vertices, especially when 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, each taking 2-3 iterations on average
to converge to the tolerance ε = 10−7 . The failing 2-8% mainlycontain cases where WALK M ANIFOLD failed for good reasons, because it was asked to walk to a point 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 that generates a proposal state conditioned on the current state of a Markov
chain. It provides the basic means 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 statewith nonzero probability density), and the rule must also supply afunction that, up to constant factors, computes the 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 rulethat 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 sothat its scale naturally adapts to the scene, including the geometry
and material properties.
The new perturbation supports general scenes and can be used bothwith the ERPT algorithm by Cline et al. and the path space MLT
framework proposed by Veach and Guibas, where it replaces andgeneralizes 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 pathusing a sequence of steps that can be grouped into sampling andconnection phases (Figure 5). The sampling phase chooses a subpath to be modified that consists of three non-specular vertices
which are potentially separated by specular chains. We shall denoteIn comparison, our perturbation can stop at vertex xb and use
the WALK M ANIFOLD algorithm to update the configuration of the
specular chain between xb and xc (Figure 6). This seemingly subtle
difference has major repercussions on the types of scenes that canbe rendered efficiently. In particular, the resulting method can systematically explore large classes of specular paths instead of having
to rely on finding them by chance.
In the following, we will discuss the perturbation in more detail;first for the ideally specular case 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 path modifications 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 initialvertex 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-specular vertex is encountered, and continue until
a second non-specular vertex is found. This path traversal may failby 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 that a < b < c.

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 that is“nearby”. When the vertex xa denotes a surface scattering event−−→with incident and exitant directions ωi = −x−a xa−1 and ωo =
−
−
−
−
→
0
xa xa+1 , the perturbation determines xa+1 by tracing a ray in a
direction ωo0 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 the other hand, when xa is a glossy material
that only reflects into a small cone of directions, large perturbations
Perturbation sampling:

SampleConnect

Figure 5: 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-specularvertex xb . To connect the vertices xb and xc , the perturbation performs a manifold walk to determine the positions of intermediate
specular vertices (if any).

ed

Up to this point, the proposed scheme is similar to the set of perturbations proposed by Veach and Guibas. However, recall that intheir work, perturbations must propagate through the path until 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 this manner is zero.

S
trac

these non-specular vertices as xa , xb , and xc . After establishing thetype 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 thepath configuration (i.e. the arrangement of specular and nonspecular vertices) changed, the perturbation is rejected immediately.

d

ate

upd

S

half-vector equal
to surface normal

D

Figure 6: Manifold perturbation example: a slightly perturbed
outgoing direction at xa is propagated until encountering the nonspecular 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.
will almost always be rejected, reducing performance.
Observe that a useful hint about the right scale can be obtained directly from the scattering model 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 scatteringfunction, and small steps are appropriate.
Our strategy is to sample from a distribution D(ωo0 ) centered at ωowhose concentration is set so that D(ωo ) equals λ2 p(ωi → ωo ). ForD(ωo0 ), we use the spherical von Mises-Fisher distribution. 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 high value, 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 toeffectively explore path space, and this results in the typical coherent noise patterns that are known from other MLT-type algorithms.
We currently set this parameter manually to achieve a desired acceptance ratio, but this could in theory be automated using adaptive
MCMC. Please refer to the supplemental material for comparison
involving different values of λ.
When xa is a camera or light source, we choose a new outgoingdirection in much the same way, but query the underlying modelfor the directional density of the associated sampling method (e.g.
uniformly choosing pixels in screen-space). To further enlarge thespace 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 bysampling 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 andxc , the connection step only entails checking that the vertices aremutually visible, and that their scattering models carry illumination
along the connection edge. When there is a manifold, we first set

x0c−1 , . . . x0b+1 = WALK M ANIFOLD(xc , . . . xb → x0b )and then perform the same verification.


fj (x̄0 )T (x̄0 , x̄)r(x̄, x̄0 ) = min 1,
fj (x̄)T (x̄, x̄0 )(10)where fj is the contribution function (Equation 4) and x̄ and x̄0 ,denote the original and proposal path respectively. Note that manyfactors cancel in the above ratio, particularly all of those in fj thatare associated with the unchanged path segment, or common terms
in the transition probability. For instance, the probability of choosing a particular sampling strategy cancels, since it only depends onthe (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 thatsampling an outgoing direction ωo0 from D(ωo0 ) at xa , and propagating it through the first specular chain, produces area density D⊥ (ωo0 ) G(xa ↔ · · · ↔ x0b ) on the surface at x0b (whereD⊥ (ωo0 ) = D(ωo0 )/| cos(na , ωo )| denotes probability with respect
to the 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 tosample paths based on their contribution to the entire image rather(j)than to a single pixel. This is accomplished by replacing We
in (4) with an importance function that measures the overall luminance received on the image 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 through analmost-smooth dielectric object cannot be explored as effectively asa 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 totake steps along a family of offset manifolds that are parallel to the
specular manifold, so that path space near the specular manifold canbe explored without stepping out of this thin band of near-speculartransport. In this section, we add a simple extension that endows
the perturbation with the ability to walk on offset manifolds and torecognize when this is appropriate.
For this, we first replace Equation (8) with the offset manifold
So = {x̄ | C(x̄) = o} ,

(11)

ed

Finally, the change in thecontribution function is computed and, together with the transition
probabilities, used to randomly accept or reject the proposal withprobability
Transition probability computation

S
trac

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
creates a potential issue when walking on the manifold, because
WALK M ANIFOLD can be non-reversible. It might succeed in moving from x to x0 but fail to move from x0 to x. Even when thereverse 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 andreject the perturbation if the path did not return to its original configuration. In the example scenes, we observed between 0% and0.4% non-reversible walks.

ted

a
pd

u

“S”

half-vector preservedin surface frame

D

Figure 7: Perturbation of path with near-specular surface interactions: instead of requiring that the half-vectors agree with thesurface normals, their direction is preserved in the surface frame.
where o captures the offset from ideal specular transport. Inspecting C (Equation 6) yields an intuitive explanation for the contentsof the vector o in terms of microfacet theory. The two entries associated with each vertex xi record the microsurface normal mi
responsible for the reflection or refraction (which is now different
from the shading normal ni ), projected into the shading tangent
space T (xi ). Our extended perturbation then preserves the projection of this microgeometry normal mi as an invariant during the
manifold traversal.
Changes to the manifold perturbation: Since the differential
geometry of the offset manifold is identical to that of the ordinary
manifold, the only required change in WALK M ANIFOLD affects theray tracing step on line four, where the algorithm now reflects andrefracts using mi instead of ni . Similarly, the deterministic phaseof the manifold perturbation responsible for propagating the sampled direction at xa to a position x0b uses these normals instead.
Note that it is straightforward to handle both cases, near-specularand specular perturbations, using the same implementation.
Recognizing near-specular transport: An important aspect
about the treatment of general scattering is the decision of whether
the surface associated with a scattering event is “smooth enough” tobe classified as part of a specular chain. We make this decision randomly by assigning a specular probability ψ(xi ) to each vertex thattakes on values 0 and 1 when xi is diffuse and specular, respectively, and values in (0, 1) when xi is at a rough interface. Thisavoids the issues of “hard” classifications that are commonly usedin rendering algorithms. For specifics on our specular probability
function ψ, please refer to the appendix.
Transition probability: In the glossy case, T is a distribution
on an offset manifold So , whereas fj is defined on the higher-(a) always non-specular(b) always specular(c) probabilistic

Figure 8:
Consistently classifying glossy materials as nonspecular or specular produces unsatisfactory results. Instead, ourmethod makes this decision randomly whenever encountering arough object (modeled after a scene by Cline et al.)The path space over which the contribution function is integrated
now consists of all possible arrangements of surface and medium
vertices in paths of length n. Each component has a measure thatis the Cartesian product of surface area and volume measures at thevertices. Using this formulation, our implementation can cleanly
abstract away the differences between surface and volume vertices
when operating on paths.

Figure 9: It is possible to extend the path space specular integration framework to volumes, for instance to render caustics from refractive objects, such as this dodecahedron-shaped luminaire withtinted glass inlays.
Sdimensional space o So . To compute transition probabilities, wemust perform a change of variables to separate out these “perpendicular” dimensions in fj , and this causes a determinant to appear
in the final transition probability:
T (x̄, x̄0 ) = Tspec (x̄, x̄0 )∂ [xb , xi1 , . . . , xik ]∂ [xb , oi1 , . . . , oik ](1 − ψ(xb ))(1 − ψ(xc )))c−1Yψ(xi ).

(12)

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 from arough mirror. The fast-varying part ofthe mirror BRDF is a function of thehalf-vector, and hence our method preserves it during manifold walks. In themedium case, we are interested in being able to handle highly peaked phasefunctions that are a function of the scat- Figure 10: Medium constraint
tering 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 )(15)

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

i=a+1,i6=bHere, Tspec is the transition probability of the purely specular caseand the indices i1 , . . . , ik ∈ {a + 1, . . . , c − 1} \ {b} refer to glossy
vertices that were classified as specular. The second line accounts
for the discrete probability of performing this classification; roughness values may change during a perturbation, and hence we mustaccount for them to maintain detailed balance. The determinant iseasy to compute from the entries of A; for details, please refer tothe supplementary material.

6

Manifold exploration for volumes

We have presented Manifold Exploration first in the context of surfaces, but the extension to volumes poses no fundamental difficulties. We briefly sketch the required changes to the mathematicalframework, following Pauly et al. [2000]. To our knowledge, noformal proof of these modifications has been published, and wetherefore include a derivation using a generalized operator theory
for media and surfaces as supplemental material of the paper.
Separate 3-space into M, the surfaces, and V, the volume between
them, replace the definition in (5) withG(x ↔ y) = τ (x ↔ y)D(x, y)D(y, x)
kx − yk2

and replace fs in (4) by the function(
fs (z → y → x)y∈M¯f (z → y → x) =σs (y)fp (z → y → x) y ∈ V

(13)

6.1

From an algorithmic perspective, manifold exploration for volumes
is almost identical to the surface 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 (15), the computation of
offset manifold (11) tangent vectors is unchanged. In the ray
tracing step 4 of WALK M ANIFOLD 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 is the outgoing direction
at xi−1 (this enforces the length constraint mentioned earlier). Afterwards, the manifold offset oi is transformed into an outgoingdirection in the new frame at x0i (Figure 10).
As in the glossy surface case, we require a criterion that clarifies
when treating a medium vertex xi as non-specular is in order, andwhen 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.

7
(14)

−
→
where D(a, b) = |N (a) · ab| if a ∈ M and 1 otherwise. The
function σs (y) denotes the scattering coefficient of the medium at
y, fp is the medium’s phase function, and τ (x ↔ y) is the transmittance between x and y. For more detail on these quantities, referto the supplemental material.

Medium manifold perturbation

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
layer that exposes cameras, light sources, scattering models, andparticipating media as generalized path vertex 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 terms of operations on vertices

SceneT ORUSC HANDELIERTABLEG LASS E GG

Seed path generator parameters

MEPT parameters

Mutations

Manifold walks

Manifold size

samples / pixel

chains / pixel

λ

total

accepted

total

converged

avg. iter.

avg.

max.

32
64
32
128

2
3
1
2

100
160
300
90

1178 M
1216 M
1074 M
1533 M

78.3%
73.6%
77.5%
72.4%

1002 M
975 M
868 M
1246 M

96.7%
97.6%
95.5%
92.2%

2.3
2.4
2.8
3.4

3.4
4.6
4.4
4.1

7
13
14
14

Table 1: Listing of seed generator and perturbation parameters, as well as captured performance statistics.
that are oblivious to whether they contain e.g. a camera model or a
medium scattering event.
We compare the following algorithms:• Primary sample space MLT by Kelemen et al., implemented
on top of bidirectional path tracing (PSSMLT).
• Path space MLT by Veach and Guibas (MLT).
• An extended form of energy redistribution path tracing by
Cline et al. (ERPT), which is seeded by bidirectional rather
than unidirectional path tracing. The ERPT implementation
shares the caustic, lens, and multi-chain perturbation with theprevious 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.
Due to the aforementioned abstraction layer, all techniques transparently support participating media even if this was originally notpart of their description. We also implemented Manifold Exploration Metropolis Light Transport (MEMLT), which corresponds
to MLT with our perturbation (i.e. the bidirectional mutation andmanifold perturbation, but none of the original perturbations). We
found that MEPT generally performs better that MEMLT and
therefore do not present results in the main paper. They can befound in the supplementary technical report along with convergedreference images.
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 the resulting 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 anda mixture of glossy, diffuse, and specular surfaces and some scattering volumes. One hour of processing time was allocated to eachrendering technique, and a comparison of the resulting images isshown in Figures 11, 12, and 13. The converged images in Figure 1 were rendered using 48 hours. The one hour renderings areintentionally unconverged to permit a visual analysis of the convergence behavior. Table 1 lists parameters and statistics collectedduring these renderings. The path generator columns refer to theseeding scheme used by ERPT and MEPT, which samples and subsequently resamples a number of paths per pixel before launching
Markov Chains. The statistics include the total number of mutations and acceptance ratio, as well as the the convergence behaviorof the manifold walks and vertex count of encountered manifolds.
C HANDELIER: In this set of results, the poor performance of MLT
is most apparent and is caused 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 the mutation 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. Seeding the same perturbations using bidirectional 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 is by 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 thebulbs into the butter dish, then to the camera (6 specular vertices)—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 some of those paths but the Veach–Guibas perturbations are
not able to explore path space around them effectively. The primary
sample space MLT variant also has difficulties rendering this scene,
because it has no knowledge about the underlying path geometry.
G LASS E GG: In this scene, our technique’s ability to create a specular chain containing both medium and surface interactions leads tofast convergence when rendering the forward-scattering medium inside the glass egg. MLT and ERPT perform poorly here, since they
do not have suitable perturbations for exploring this space. Because
the MLT perturbations treat glossy and diffuse materials identically,
they have difficulty rendering the near-specular tabletop, producing
streak-like artifacts in the output rendering.

8

Conclusions

We have presented a new type of Markov Chain Monte Carlo rendering method that models the space of valid specular and nearspecular light paths using high-dimensional manifolds. The differential geometry of these manifolds provides a powerful tool toefficiently 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 workin 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 manifold
for which seed paths are provided. Bidirectional Path Tracing isreasonably 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 numberof different path types exceeds the number of samples that can begenerated, local exploration of path space becomes ineffective; fu-

(a) MLT(b) ERPT(c) PSSMLT(d) MEPT

Figure 11: C HANDELIER: 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 ischallenging, as certain important light paths are found with low probability, particularly those involving interreflection between the bulbs
and the body of the chandelier. In this and the following comparisons, one hour of processing time was allocated to each rendering technique.

(a) MLT(b) ERPT(c) PSSMLT(d) MEPT

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

(a) MLT(b) ERPT(c) PSSMLT(d) MEPT

Figure 13: G LASS E GG: 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.
ture algorithms could be designed to attempt exploration only insufficiently 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. Similar refinements
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 beused in other kinds of algorithms as well, including deterministic
ones to map out specular paths, for instance in design of luminaires
or optical systems.

9

Acknowledgements

move its microsurface normal mi from a region of high density to
one of low density, and we choose the 90th-percentile to classify
the support of Dα into such regions.

90%
percentile

perturbed by

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 thearea ratio of the two highlighted regions on the sphere:ψ(xi ) =1 − cos θq (α(xi )),1 − cos (θq (α(xi )) + ∆θ)(16)This research was conducted in conjunction with the Intel Scienceand Technology Center for Visual Computing. Additional funding
was provided by the National Science Foundation under grant IIS1011919. The authors are indebted to Olesya Isaenko, who meticulously crafted the example scenes used in the evaluation.

where θq is the aforementioned percentile (with q set to 0.9). For
the Beckmann distribution, this is given by

Appendix

Implementation-wise, this heuristic requires material models to beable to compute their Beckmann distribution-equivalent roughness
or provide a custom quantile function.

We found the following following heuristic based on 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θq (α) := tan−1 (−α2 log(1 − q)).

Specular probability function for surfaces:

In themedium case, we use the same probability (16), but now with apercentile that is suitable for volumetric scattering. We use
Specular probability function for participating media:θ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 the Henyey-Greenstein phase function. In
Ψ(xi ) (Equation 16) we must also replace ∆θ with the averagechange in scattering angle at medium vertices, again determined ina brief phase before rendering.

K ITAOKA , S., K ITAMURA , Y., AND K ISHINO , F. 2009. Replica
exchange light transport. Computer Graphics Forum 28, 8(Dec.), 2330–2342.

References

L AI , Y.-C., FAN , S. H., C HENNEY, S., AND DYER , C. 2007. Photorealistic image rendering with population monte carlo energy
redistribution. In Rendering Techniques 2007: 18th Eurographics Workshop on Rendering, 287–296.

C HEN , M., AND A RVO , J. 2000. Perturbation methods for interactive specular reflections. IEEE Trans. Vis. and Comp. Graph.
6, 3 (July/Sept.), 253–264.
C HEN , M., AND A RVO , J. 2000. Theory and application of specular path perturbation. ACM Trans. Graph. 19, 4 (Oct.), 246–278.
C HEN , J., WANG , B., AND YONG , J.-H. 2011. Improved stochastic progressive photon mapping with metropolis sampling. Computer Graphics Forum 30, 4, 1205–1213.
C LINE , D., TALBOT, J., AND E GBERT, P. 2005. Energy redistribution path tracing. ACM Trans. Graph. 24, 3 (Aug.), 1186–1195.
C OOK , R. L., P ORTER , T., AND C ARPENTER , L. 1984. Distributed ray tracing. In Computer Graphics (Proceedings of SIGGRAPH 84), 137–145.
G ORAL , C. M., T ORRANCE , K. E., G REENBERG , D. P., AND
BATTAILE , B. 1984. Modeling the interaction of light between
diffuse surfaces. In Computer Graphics (Proceedings of SIGGRAPH 84), 213–222.
H ACHISUKA , T., AND J ENSEN , H. W. 2009. Stochastic progressive photon mapping. ACM Trans. Graph. 28, 5 (Dec.).
H ACHISUKA , T., AND J ENSEN , H. W. 2011. Robust adaptive
photon tracing using photon path visibility. ACM Trans. Graph.
30, 5 (Oct.), 114:1–114:11.
H ASTINGS , W. K. 1970. Monte carlo sampling methods using
markov chains and their applications. Biometrika 57, 1, 97–109.
H ECKBERT, P. S. 1990. Adaptive radiosity textures for bidirectional ray tracing. In Computer Graphics (Proceedings of SIGGRAPH 90), 145–154.
I GEHY, H. 1999. Tracing ray differentials. In Computer Graphics(Proceedings of SIGGRAPH 99), 179–186.
JAKOB , W., 2010.
renderer.org.

Mitsuba renderer.

L AFORTUNE , E. P., AND W ILLEMS , Y. D. 1993. Bi-directional
path tracing. In Proceedings of Compugraphics 93.

L IU , J. S. 2001. Monte Carlo strategies in scientific computing.
Springer.
M ETROPOLIS , N., ROSENBLUTH , A. W., ROSENBLUTH , M. N.,
T ELLER , A. H., AND T ELLER , E. 1953. Equation of statecalculations by fast computing machines. J. Chem. Phys. 21, 6.
M ITCHELL , D. P., AND H ANRAHAN , P. 1992. Illumination fromcurved reflectors. In Computer Graphics (Proceedings of SIGGRAPH 92), 283–291.
PAULY, M., KOLLIG , T., AND K ELLER , A. 2000. Metropolis
light transport for participating media. In Rendering Techniques
2000: 11th Eurographics Workshop on Rendering, 11–22.
P REMOZE , S., A SHIKHMIN , M., AND S HIRLEY, P. 2003. Path integration for light transport in volumes. In Eurographics Symposium on Rendering: 14th Eurographics Workshop on Rendering.
S EGOVIA , B., I EHL , J., AND P ÉROCHE , B. 2007. Metropolis
instant radiosity. Computer Graphics Forum 26, 3 (Sept.), 425–
434.
S HIRLEY, P. S., WADE , B., H UBBARD , P., Z ARESKI , D., WAL TER , B., AND G REENBERG , D. P. 1995. Global illumination
via density estimation. In Eurographics Rendering Workshop
1995, 219–231.
S ILLION , F. X., AND P UECH , C. 1989. A general two-pass method
integrating specular and diffuse reflection. In Computer Graphics (Proceedings of SIGGRAPH 89), 335–344.
S PIVAK , M. 1965. Calculus on Manifolds. Addison-Wesley.
V EACH , E., AND G UIBAS , L. 1994. Bidirectional estimators forlight transport. In Fifth Eurographics Workshop on Rendering.
V EACH , E., AND G UIBAS , L. J. 1997. Metropolis light transport.
In Computer Graphics (Proceedings of SIGGRAPH 97), 65–76.

http://www.mitsuba-

V EACH , E. 1997. Robust Monte Carlo Methods for Light Transport
Simulation. PhD thesis, Stanford University.

JAROSZ , W., Z WICKER , M., AND J ENSEN , H. W. 2008. The
beam radiance estimate for volumetric photon mapping. Computer Graphics Forum 27, 2 (Apr.), 557–566.

WALTER , B., H UBBARD , P. M., S HIRLEY, P. S., AND G REEN BERG , D. F. 1997. Global illumination using local linear density
estimation. ACM Trans. Graph. 16, 3 (July), 217–259.

J ENSEN , H. W., AND C HRISTENSEN , P. H. 1998. Efficient simulation of light transport in scenes with participating media using photon maps. In Computer Graphics (Proceedings of SIGGRAPH 98), 311–320.

WALTER , B., M ARSCHNER , S. R., L I , H., AND T ORRANCE ,K. E. 2007. Microfacet models for refraction through roughsurfaces. In Rendering Techniques 2007: 18th Eurographics
Workshop on Rendering, 195–206.

J ENSEN , H. W. 1996. Global illumination using photon maps. In
Eurographics Rendering Workshop 1996, 21–30.

WALTER , B., Z HAO , S., H OLZSCHUCH , N., AND BALA , K. 2009.
Single scattering in refractive media with triangle mesh boundaries. ACM Trans. Graph. 28, 3 (July), 92:1–92:8.

K AJIYA , J. T. 1986. The rendering equation. In Computer Graphics (Proceedings of SIGGRAPH 86), 143–150.
K ELEMEN , C., S ZIRMAY-K ALOS , L., A NTAL , G., AND
C SONKA , F. 2002. A simple and robust mutation strategy for themetropolis light transport algorithm. Computer Graphics Forum21, 3, 531–540.

W HITTED , T. 1980. An improved illumination model for shaded
display. Communications of the ACM 23, 6 (June), 343–349.