Eurographics Symposium on Rendering 2011
Ravi Ramamoorthi and Erik Reinhard(Guest Editors)

Volume 30 (2011), Number 4

Progressive Expectation–Maximization for
Hierarchical Volumetric Photon Mapping
Wenzel Jakob1,2

Christian Regg1,3

1 Disney Research Zürich

2 Cornell University

Beam Radiance Estimation [Jarosz et al. 2008]

917K Photons (13 + 268 = 281 s)

Per-Pixel Render-time

Wojciech Jarosz1
3 ETH Zürich

Our Method

4K Gaussian fit (54 + 71 = 125 s)

Per-Pixel Render-time

Figure 1: The beam radiance estimate (left) finds all photons along camera rays, which is a performance bottleneck for this C ARS scene due
to high volumetric depth complexity and many photons (917K). Our method (right) fits a hierarchical, anisotropic Gaussian mixture model tothe photons and can render this scene faster, and with higher quality using only (4K) Gaussian components. The listed times denote the costsof the preprocessing stage (including photon tracing and hierarchy construction), as well as the final rendering stage, respectively.

Abstract
State-of-the-art density estimation methods for rendering participating media rely on a dense photon representation of the radiance distribution within a scene. A critical bottleneck of such kernel-based approaches is theexcessive number of photons that are required in practice to resolve fine illumination details, while controllingthe amount of noise. In this paper, we propose a parametric density estimation technique that represents radianceusing a hierarchical Gaussian mixture. We efficiently obtain the coefficients of this mixture using a progressive
and accelerated form of the Expectation–Maximization algorithm. After this step, we are able to create noise-freerenderings of high-frequency illumination using only a few thousand Gaussian terms, where millions of photons
are traditionally required. Temporal coherence is trivially supported within this framework, and the compact
footprint is also useful in the context of real-time visualization. We demonstrate a hierarchical ray tracing-based
implementation, as well as a fast splatting approach that can interactively render animated volume caustics.
Categories and Subject Descriptors (according to ACM CCS): I.3.7 [Computer Graphics]: Three-Dimensional
Graphics and Realism—Ray Tracing; I.6.8 [Simulation and Modeling]: Simulation—Monte Carlo

1. Introduction
Light interactions in participating media are responsible formany subtle but rich lighting effects in fire, water, smoke,
clouds and crepuscular rays. Unfortunately, accurately rendering these effects is very costly. Cerezo et al. [CPP∗ 05]
provide a comprehensive overview of recent techniques.
Though methods such as (bidirectional) path tracing [LW93, VG94, LW96] or Metropolis light transport [PKK00] are general and provide unbiased solutions, it
is costly to obtain high-quality, noise-free solutions except
in the simplest scenes. The most successful approaches
typically combine photon tracing with a biased Monte Carloframework [JC98, WABG06, JZJ08].
Volumetric photon mapping [JC98] traces paths fromlight sources and stores “photons” at the vertices of thesec 2011 The Author(s)This is the author’s version of the work. It is posted here by permission of The Eurographics Association for your personal use. Not for redistribution. The definitive version
is available at http://diglib.eg.org/.

paths. A subsequent rendering pass approximates lighting
using density estimation of the photons. The recent beam
radiance estimate (BRE) [JZJ08] for photon mapping further accelerates this approach and improves the rendering
quality by accounting for all photons that are located along
a ray query.
However, even with such state-of-the-art optimizations,photon mapping suffers from a few remaining sources of
inefficiency: 1) complex high-frequency lighting requires anexcessive number of photons, which increases memory demands and leads to slow rendering times; 2) many photons
may be processed which ultimately contribute very little tothe final pixel color (due to subsequent attenuation or small
projected pixel size); and 3) since photons are generated
stochastically, the reconstructed density is prone to flickering unless many photons are used.

W. Jakob, C. Regg, and W. Jarosz / Progressive Expectation–Maximization for Hierarchical Volumetric Photon Mapping

In this paper, we propose an algorithm that operates ona new representation of volumetric radiance to efficiently
render images. This representation is compact and expressive, solving the aforementioned problems. Instead of storing individual photons and performing non-parametric density estimation [Sil86] during rendering, we progressively fit
a small set of anisotropic Gaussians to the lighting distribution during photon tracing. Since each Gaussian is fit to
the density of a large number of photons, we can accurately
reconstruct complex lighting features during rendering usinga small number (thousands) of Gaussians while providing
better reconstruction than a large number (millions) of photons. Our framework trivially provides temporal coherence
and, to handle the varying contributions of Gaussians to therendered image (due to attenuation and projected pixel size),
we propose a hierarchical approach which provides highquality level-of-detail (LOD). Our compact hierarchical representation provides direct performance benefits which wedemonstrate both in offline and interactive rendering.
2. Related Work
Photon Mapping. Schjøth and colleagues demonstrated
the benefits of anisotropic density estimation [SFES07,
SSO08] to obtain higher quality surface caustics with photon mapping. Hachisuka et al. [HOJ08] progressively refineradiance estimates at fixed measurement points to achieve
convergent bias and variance in a view-dependent manner.
We instead fit to the photon density by progressively refining both the positions and other parameters of an anisotropic
Gaussian mixture, providing reduced variance and faster
rendering from arbitrary viewpoints. Our work has a numberof similarities to both hierarchical photon mapping representations [CB04, SJ09a] which provides LOD for final gather,and photon relaxation [SJ09b] which iteratively moves photon positions to reduce variance. By using Gaussian mixtures
we obtain both of these benefits while drawing upon a set ofprincipled and well-established fitting techniques.
Hierarchical and Point-based Rendering. Level-ofdetail for light transport has been used to accelerate theRadiosity algorithm [HSA91, Sil95] and volume visualization [LH91]. Point-based [DYN04] and mesh-less [LZT∗ 08]
variants allow for more flexible geometric representations
while Lightcuts [WFA∗ 05, WABG06] provides hierarchicalevaluation of lighting from many pointlights.
Our hierarchical rendering approach is related to hierarchical splatting of complex point-sampled geometry [RL00]
and large stellar datasets [HE03, HLE04]; however, we
use an anisotropic Gaussian representation. EWA splatting [ZPvBG02] uses Gaussians to provide improved, aliasfree reconstruction for surface and volume visualization.
While we share similar goals to these approaches, we operate within a parametric density estimation context andprogressively fit the positions and other parameters ofanisotropic Gaussians to the photon density.

Parameter Fitting. As opposed to previous global illumination techniques, which have relied primarily on nonparametric kernel density estimation such as the nearest
neighbor method [Sil86], our technique leverages a large
body of work on parametric density estimation and parameter fitting. Instead of performing the density estimation
solely during rendering, we fit anisotropic Gaussians to thephoton density during precomputation.
The use of Gaussians for rendering participating media
is not new [ZHG∗ 07, ZRL∗ 08]; however, previous methods
used isotropic Gaussians to approximate the heterogeneous
media density, whereas we fit anisotropic Gaussians to theradiance distribution. We furthermore support hierarchicalevaluation and leverage the well-established Expectation–
Maximization [DLR∗ 77] (EM) method specifically tailoredfor fitting anisotropic Gaussian mixtures.
While EM has seen recent use in image synthesis [TLQ∗ 05, HSRG07, PJJ∗ 11], it has a much richer
history in fields such as machine learning and data mining
where hierarchical techniques [GR05, GNN10] directly
apply to LOD, and incremental techniques [NH98, HG09]
provide parameter fitting within a limited memory footprint
for streaming data. In these fields, however, large Gaussian
mixtures are typically on the order of a few dozen to afew hundred components and their geometric quality isoften only indirectly reproduced (e.g. color segmentation/classification). In contrast, to accurately fit volumetriclighting, we rely on much larger Gaussian mixtures (tensof thousands), and their fitted parameters directly influence
the quality of the rendered image. To obtain the necessary
quality and performance under these conditions we extend
the accelerated EM method [VNV06].
3. Preliminaries
3.1. Radiative Transport in Participating MediaIn the presence of participating media, the light incident at
any point in the scene x (e.g. the camera) from a direction~ω (e.g. the direction through a pixel) as expressed by theradiative transport equation [Cha60] is the sum of two terms:
L(x,~ω) = Tr (x ↔ xb ) Ls (xb ,~ω) + Lm (x,~ω, b),(1)where the first term accounts for the radiance reflected off
surfaces and the second term describes radiance scattered in
the medium. Before reaching x, the radiance from the nearest surface at a distance b, xb = x − b~ω, is reduced by thetransmittance Tr which describes the percentage of light that
propagates between two points (see Dutre et al. [DBB06] fordetails). We summarize our notation in Table 1.
In this paper we are primarily concerned with the media
radiance which integrates the scattered light up to the nearest
surface
Z bZ

Lm (x,~ω, b) =σs (xt ) Tr (x ↔ xt) f (θt ) L(xt ,~ωt ) d~ωt dt.
0 Ω4π

(2)

c 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

W. Jakob, C. Regg, and W. Jarosz / Progressive Expectation–Maximization for Hierarchical Volumetric Photon Mapping
Table 1: The notation used in this paper.

the media using a hierarchical, anisotropic Gaussian mixture.

SymbolDescription{xi }n1inp(x)p(x|Θ)
g(x|Θs )
Θs
Θ
s
k
L(Θ)
F (Θ, q)qi (s)qA (s)
nA
h..iA
Θ̄ j
j

Photon positions used in the EM fit
Index for photons
Total number of photons
True probability density of photon positions
Probability density of GMM approximating p(x)Gaussian density at x of GM component s
Weight, mean, and covariance matrix (ws , µs , Σs ) of component s
Set of all mixture model parameters {Θs }k1
Index for GMM components
Number of GMM components
Log-likelihood function defined in (6)Lower bound on the likelihood function, see (8)Responsibility of point i for GMM component s
Responsibility of cell A for GMM component s
Number of points xi that fall within A
Average of the quantity .. w.r.t. the points xi ∈ A
Parameters (w̄ j , µ̄ j , Σ̄s ) of GMM pyramid component j
Index into GMM pyramid nodes

3.3. Gaussian Mixtures
We briefly review Gaussian mixture models (GMMs), maximum likelihood estimation, and Expectation–Maximizationin terms of our rendering problem.
Given a collection of photons X = {xi }n1 , volumetric photon mapping corresponds to estimating the underlying probability distribution p(x). In general, p(x) could be any distribution obtainable using photon tracing. Photon mapping
estimates this distribution using kernel density estimation,
whereas our approach relies on a GMM. More precisely, we
assume that the underlying probability density can be modeled using a weighted sum of k probability distributions,kHere σs is the scattering coefficient of the medium, f is thephase function, cos θt = ~ωt ·~ω, and xt = x − t~ω. This joint integral along the ray and over the sphere of directions, whichrecursively depends on the incident radiance (1), is whatmakes rendering participating media so costly.
3.2. Volumetric Photon Mapping and the BRE

n

(3)

i=1

where xi , Φi and ~ωi are the position, power and incident direction of photon i respectively, cos θi = ~ωi · ~ω, and Ki is a
2D normalized kernel of radius ri around photon i, evaluatedat xt = x + ((xi − x) ·~ω)~ω. We always use i to index photons.
Note that Jarosz et al. [JZJ08] included an additional σs term.
We account for this during photon instead tracing to remain
consistent with other sources [JC98, JNSJ11].
Unfortunately, there are a few common situations where
the BRE is the bottleneck in rendering. Firstly, complex,
high-frequency lighting may require an excessive number ofphotons, which increases storage costs and also slows down
the beam query during rendering. Secondly, since the BRE
always finds all photons that intersect the query ray, if thescene has a large extent, many photons may be processedwhich ultimately contribute very little to the final pixel color(due to subsequent attenuation or projected pixel size). Figure 1 visualizes the time spent performing the BRE per pixel,
which shows this deficiency. Our approach addresses both ofthese problems by expressing the distribution of photons in
c 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

(4)

s=1

where Θs denotes the parameters for distribution s andΘ = {Θs }k1 defines the model parameters for the entire mixture. Each individual component g is a normalized multivariate Gaussiang(x|Θs ) =Volumetric photon mapping [JC98] accelerates this computation using a two-stage procedure. First a photon map ispopulated by tracing paths from light sources and recording
each scattering event as a “photon.” Then, density estimation is performed over the photons to approximate the media
radiance. To perform this computation efficiently, Jarosz etal. [JZJ08] proposed the beam radiance estimate, which firstassigns a radius ri to each of the n stored photons (using apilot density estimate) and then approximates Equation (2)asLm (x,~ω, s) ≈ ∑ Ki (xt ) Tr (x ↔ xi ) f (θi ) Φi ,p(x) ≈ p(x|Θ) = ∑ ws g(x|Θs ),
exp − 12 (x − µs )| Σ−1(x − µs )ps,(2π)3/2 |Σs |(5)where Θs = (µs , Σs , ws ) with mean µs , covariance matrix Σs ,and weighted by ws . We use |·| to denote the determinant
and always use s to index into the Gaussian components.
Since photons are not actually generated using a Gaussian process, (4) is an approximation, but by increasing ka Gaussian mixture can model any photon distribution with
arbitrary accuracy.
3.4. Maximum Likelihood Estimation
With the assumptions above, our problem is reduced to finding the values of the model parameters Θs for each of the kGaussians, which, when plugged into the model (4), are most
likely to generate the photons. More formally, such maximum
likelihood estimators seek the parameters Θ that maximize ofthe associated data log-likelihoodnL(Θ) = ∑ log (p(xi |Θ)) .

(6)

i=1

Note that the actual photon density p(x) is not an explicit partof this formulation; instead, it only enters the problem statement via evaluation of the GMM (4) at the photon positions.
3.5. Expectation Maximization
The popular EM algorithm is an iterative technique for finding such a maximum likelihood parameter estimate. Insteadof maximizing L directly (which is usually intractable), EM
starts with an initial guess of the parameters Θ and modifies
them so that a lower bound F ≤ L is strictly increased in every one of its iterations. The bound F is chosen so that its local maxima coincide with those of L; hence, it is guaranteed

W. Jakob, C. Regg, and W. Jarosz / Progressive Expectation–Maximization for Hierarchical Volumetric Photon Mappingthat the algorithm eventually converges to a local maximumof L.

This formulation has several advantages which we will
exploit for photon mapping:

During the optimization, EM makes use of a series of
discrete distributions qi (s) : {1, . . . , k} → [0, 1] (one for eachphoton xi ), which specify the associated photon’s “responsibility” towards each GMM component. The lower boundF (Θ, q) can then be expressed as a function of the parameters
Θ, as well as the auxiliary distributions qi . An intuitive andgeneral interpretation of the EM-type algorithm can then beformulated as [NH98]:• Compact representation: The photons xi can be dis-

E-Step: Given Θ, find q that maximizes F (Θ, q), then
M-Step: Given q, find Θ that maximizes F (Θ, q).
3.6. Accelerated EM
In this paper, we make use of a technique called accelerated
EM [VNV06]. The insight of this approach is to assign thesame distribution qi (s) = qA (s) to all photons xi that fall into
the same cell A of a spatial partition P of R3 . P is simply
a subdivision of space where any point lies in exactly one
cell, e.g. the collection of cells of a Voronoi diagram, or acut through a kd-tree or octree. This allows rephrasing the
EM algorithm purely in terms of operations over cells (asopposed to photons).
For the E-step, Verbeek et al. provide a provably-optimal
assignment of the responsibilities:qA (s) =ws exphlog g(x|Θs )iA∑ka=1 wa exphlog g(x|Θa )iA(7)

where h f (x)iA denotes the average of a function f over all
photons that fall within cell A. They also give an expression
for the lower bound F as a sum over the partition’s cells:


k
p(s)
+ hlog g(x|Θs )iA .
F (Θ, q) = ∑ nA ∑ qA (s) log
qA (s)A∈P s=1(8)

The average log-densities in (7) and (8) can be readily computed in closed-form for Gaussian mixtures [VNV06]. For
this, we only need to store the average photon position hxiAand outer product hxxT iA within each cell.
The accelerated M-step computes the parameters Θs as
weighted averages over cell statistics:
Θs = (ws , µs , Σs ) =1hxiA hxx| iA − µs µ|snA qA (s) 1,
,
, (9)
∑
n A∈P
ws
ws

where n denotes the total number of photons and nA specifies the number of photons that fall within cell A. Note thatthe terms of Θs need to be computed in sequence to satisfy
interdependencies. A naïve accelerated EM implementationis shown in Algorithm 1.
Algorithm 1 ACCELERATED EM-NAÏVE(Θ(0) , P)

4. Overview of our Algorithm
In the remainder of this paper we apply the concepts ofaccelerated EM to the problem of rendering participating
media using photon tracing. At a high-level, our approach
replaces individual photons in the BRE with multivariate
Gaussians which are fit to the photons using EM. Since a
Gaussian is fit to a large collection of photons, each term isinherently less noisy and more expressive than an individual photon. This provides high quality reconstruction using
a relatively low number of terms, which results in faster rendering at the expense of a slightly longer view-independent
precomputation stage.
Our method extends accelerated EM from Section 3.6 inthree simple but important ways. We firstly introduce sparsity into the algorithm to obtain asymptotically faster performance. Secondly, we incorporate progressive refinement
by incrementally shooting photons, updating statistics in thepartitions, and refining the EM fit. This elegantly allows
memory-bounded fitting of extremely large collections ofphotons while simultaneously providing an automatic stopping criterion for photon tracing. We call this process progressive EM. We believe this extension could also find applicability in data mining for incremental fitting of streaming
data. Lastly, we build a full GMM pyramid which allows us
to incorporate hierarchical LOD during rendering.
This algorithm is divided into four key steps (see Figure 2):
1. Photon tracing: Trace photons within the scene (just likein photon mapping), storing the resulting statistics in thepartition P , then discard the photons.
Shoot
Photons
Initial
Guess

E
M

yes

Progressive EM
Shoot More
Photons

RefineCut
noconverged?

1 z=0// iteration counter2 repeat3
Compute the distribution qA (s) for each A ∈ P using (7)4Compute Θ(z+1) using Equation (9)5z = z+1
6 until F (Θ(z) , q(z) )/F (Θ(z−1) , q(z−1) ) − 1 < εcarded, and only a few sufficient statistics of their behavior within each cell A ∈ P must be stored.
• Fast iterations: A small number of statistics succinctly
summarize the photons within a cell. Good results can beobtained with comparatively few cells, where many photons would be needed by non-accelerated EM; this allows
for shorter running times.
• Multiresolution fitting: The partition need not necessarily stay the same over the course of the algorithm. Reminiscent of multigrid methods, we can start with a relatively rough partition and refine it later on.

Build
Hierarchy

seed with converged GMM from previous frame

RenderFigure 2: The program flow of our approach.

c 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

W. Jakob, C. Regg, and W. Jarosz / Progressive Expectation–Maximization for Hierarchical Volumetric Photon Mapping(a) Target density(b) Initial cut(c) Initial GMM(d) Final cut(e) Final GMM

Figure 3: As a simple two-dimensional application of our method, we fit a 1024-term GMM to an image-based target density function (a). Our
method starts with a coarse initial cut (b) and GMM parameters (c) and refines them until convergence (d and e).

2. Progressive GMM fitting: Apply progressive EM (Section 5) to the statistics in P . Should the stored statistics
be too noisy, repeat step 1 by shooting more photons, updating statistics, and continue refining the GMM fit.
3. Level-of-Detail: Construct a full LOD pyramid of
GMMs with decreasingly many components (Section 6).
4. Hierarchical Rendering: For each ray, hierarchicallyprune the LOD pyramid to the set of Gaussians with
sufficient contribution according to a conservative error
heuristic (Section 7).

memory, we use a moderate number of photons (n = 128 · k)
and create a leaf node once the number of photons in a subtree drops below 5. After this step, all memory associatedwith the photons can be released, and only the hierarchy remains.

We focus on practical considerations for implementing thesesteps in the following sections and omit details for the firststage since it is identical to standard photon mapping.

5.2. Progressive Update and Cut Refinement

5. Progressive GMM Fitting
Our progressive EM algorithm is outlined in Algorithm 2.
The core of the approach is based on accelerated EM, but
we introduce sparsity to lines 4–5, progressively update thestatistics (line 8), and use a more efficient cut refinement
procedure (line 10).
Algorithm 2 P ROGRESSIVE EM(Θ(0) , P)
1 z=0// iteration counter2 repeat3
repeat4
Compute the distr. qA (s) for all A ∈ P using (7)5Compute Θ(z+1) using Equation (9)6z = z+1
7
until F (Θ(z) , q(z) )/F (Θ(z−1) , q(z−1) ) − 1 < ε8
if Equation (10) indicates an unsatisfactory fit
9Shoot more photons and update statistics
10
Expand cut along 50% “best” nodes according to F (8)11 until the relative change in the expanded F drops < ε5.1. Partition Construction
Our algorithm assumes the availability of a partition P . To
obtain a spatial partition of R3 , we follow an approach similar to the one taken by Verbeek et al.: we create an entire
family of possible partitions by constructing a balanced kdtree over the photon positions, where each leaf and innernode is augmented with the aforementioned sufficient statistics; this can be done in O(n log n) time. Any cut through thetree then forms a valid partition. To ensure that the tree fits in
c 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

We initialize the cut to the level of the kd-tree where thenumber of cells equals eight times the number of GMM
components (we found this to be a good compromise interms of the total running time).

After running E- and M-steps until convergence, we checkthe effective number of photons that contribute to GMM
component s:
neff (s) = ∑ nA qA (s).
(10)
A∈P

If any Gaussian is fit to fewer than 64 photons, we deem
these statistics to be unreliable, so we shoot more photons,
and progressively update the statistics. Shooting more photons reduces the variance of the statistics, converging to theproper average of these quantities within each cell in thelimit.
After this step, we use a greedy criterion to push down
the partition cut. This entails removing nodes and replacing
them by their children in the kd-tree. Verbeek et al. proposed
expanding the cut by one node in each outer iteration, but
we found that this results in slow convergence. Instead, we
determine the most suitable set of nodes by computing theincrease in F (8) for each possible replacement and executethe top 50% of them. We repeat this entire procedure until
fitting parameters converge. Figure 3 illustrates this procedure on an image-based density function. The progressive
EM iterations and cut refinement, as well as the computation of Equation 10, can all be done efficiently using ideas
discussed in the next section.
5.3. Exploiting Sparsity
Unfortunately, the complexity of this algorithm renders it
impractical for large-scale problem instances, which are themain focus of this paper. In particular, although the computation of (7) and (9) no longer depends on the numberof photons, the cost is linear in the product of |P| and k.
To obtain a reasonable fit, the final cut will usually satisfy

W. Jakob, C. Regg, and W. Jarosz / Progressive Expectation–Maximization for Hierarchical Volumetric Photon Mapping|P| = c · k, for some constant c > 1, hence the running time isin O(k2 ).
To work around the complexity issue, our fitting pipeline
exploits the inherent sparsity of this problem: due to the exponential falloff in g, it is reasonable to expect that a cellA ∈ P will excert little influence on a distant Gaussian component s. In the context of accelerated EM, this means thatthe value of qA (s) will be negligible. We may thus want toskip such cell-GMM component pairs in the calculation of(9) or (10).
There is one caveat to the previous idea: suppose thata cell A ∈ P is located relatively far away from all of the
GMM components. Although the numerator in qA (s) (7) will
be very small in this case, qA (s) will nonetheless exert a significant influence on distant GMM components due to thenormalization term in its denominator. This means that a
pruning criterion cannot be based on spatial distance alone.
Our method therefore proceeds in two stages: first, we determine the normalization term in the denominator of qA (s)(7) for each cell A ∈ P . Following this, we recompute the
GMM parameters Θ and effective photon counts neff (s) using
Equation (9) and (10). All of these steps exploit sparsity andobey a user-specified relative error tolerance τ, such that thealgorithm becomes faster and more approximate with higher
τ and degenerates to the O(n2 ) approach when τ = 0. Specifically, we loosely bound qA as a function of spatial distanceand drop Gaussian-cell pairs with the help of a bounding volume hierarchy when this does not cause the result to exceed
the specified error bound. All our results use τ = 10−4 .
5.4. Initial Guess and Temporal Coherence
EM requires an initial guess which is iteratively refined.
Since EM will, in general, only converge to a local maximum of L, the starting parameters are an important factor in
the final quality of the EM fit.
For still images we obtain an initial guess by extracting
a suitably-sized horizontal cut through the partition kd-tree.
Each cell in this cut already stores the sufficient statistics
that are needed to analytically fit a Gaussian to its contents,
hence this step takes negligible time.
For animations we could apply the same procedure foreach frame independently. However, like any photon mapping method, the stochastic positions of photons can introduce temporal flicking. We can achieve significantly better
quality by seeding EM with the converged GMM parameters from the previous frame (a similar idea was used by
Zhou et al. [ZRL∗ 08] for temporally coherent non-linear optimization). This optional step does introduce dependence
between frames, but we found it to be worthwhile since itcompletely eliminated temporal flicker in our results. Furthermore, since the previous frame is often a remarkably
good match to the current lighting, this significantly reduces
the computation time needed for fitting (for instance, fitting

the B UMPY S PHERE animation is 4.35× faster compared toframe-independent initialization).
6. Level-of-Detail Construction
To render the GMM at various levels of detail, we wish toconstruct a hierarchical representation of the Gaussian mixture. In our approach, the leaves are the components of ourGMM, interior levels approximate the GMM using decreasing numbers of Gaussian components, up to the root which
expresses the entire distribution as a single Gaussian. Notethat we cannot simply compute a multi-resolution representation (where each level of the tree is an independent GMM
with 2l components). Since we will dynamically select subsets of this tree for rendering, there needs to be a proper
nesting between the levels of the hierarchy such that a node
is a representative for its entire subtree.
To accomplish this efficiently, we apply a modified version of Walter et al.’s heap-based agglomerative clustering
algorithm [WBKP08] as shown in Algorithm 3. We denotenodes of the GMM pyramid as Θ̄ and index a specific nodeas Θ̄ j . We initialize the leaves of the pyramid (line 1) to theGaussian mixture, {Θs }k1 , constructed in Section 5. To form
the next higher level, we start with a heap of all minimum
Gaussian distances in the current level (lines 6–8) and thenmerge (and remove) the two “most similar” Gaussians until
no Gaussians are left in the heap (lines 9–15). Repeating thisentire process log2 k times creates a full binary tree.
Algorithm 3 C ONSTRUCT G AUSSIAN P YRAMID(Θ = {Θs }k1 )1 Θ̄ = Θ// construct tree starting with leaf nodes
2 level = log2 k
3 while level > 0 // loop until we reach the root4Θ̄0 = ∅// initialize next level to empty set
5H =∅// initialize min-heap to empty set
6for each Gaussian j1 ∈ Θ̄7Find closest Gaussian j2 to j1 according to (11).
8
Insert edge ( j1 , j2 ) into H.
9while H 6= ∅10
Pop an edge e from the heap11
if neither endpoint of e has been previously merged
12Merge the two Gaussians of e according toEquations (12–14) and insert into Θ̄013
for all vertices j1 connected to either end of e
14
Find the closest Gaussian j2 to j1 .
15
Insert edge ( j1 , j2 ) into H.
16Θ̄ = Θ̄017level = level − 1
18 return Θ̄We calculate the distance between Gaussians using a symmetric dissimilarity functiond( j1 , j2 ) = w j1 · w j2 · min(DKL (Θ̄ j1 k Θ̄ j2 ), DKL (Θ̄ j2 kΘ̄ j1 )), (11)where DKL is the Kullback-Leibler divergence between component j1 and j2 , and takes on a simple analytical form in thecase of multivariate Gaussians [GNN10]. Kullback-Leibler
divergence is a fundamental statistical measure of relativec 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

W. Jakob, C. Regg, and W. Jarosz / Progressive Expectation–Maximization for Hierarchical Volumetric Photon Mapping

1
2

3

1

3

1

1

(b) Intersect bounding box

1

Tr < ɛTr

Ω(r)(a) Hierarchy usedin (b) - (e)

2

Ω1(r)

(c) Criterion 1: small
projected size(d) Criterion 2: lowtransmittance(e) Otherwise: recurse

Figure 4: We (a) traverse our hierarchy of Gaussians from top to bottom, (b) intersecting the current node’s AABB. We then (c) test the solid
angle of the representative and also check the extinction (d) before traversing down the tree.

entropy between two distributions and has been used extensively for clustering Gaussian distributions [GR05, GNN10].
To merge Gaussians Θ̄ j1 and Θ̄ j2 into Θ̄ j , we analytically
fit parameters as follows [GR05]:w̄ j = w̄ j1 + w̄ j2 ,(12)µ̄ j = w̄−1j (w̄ j1 µ̄ j1 + w̄ j2 µ̄ j2 ) ,Σ̄ j = w̄−1w̄kjk∈{ j1 , j2 }

∑

(13)
|

Σ̄k + (µ̄k − µ̄ j )(µ̄k − µ̄ j )



.

(14)

This procedure eventually forms a full binary tree that we
use for hierarchical pruning during rendering.
7. Rendering
During rendering, for a given query ray passing through themedium, we are interested in approximating the medium radiance (2) using our GMM. If we ignored the hierarchy constructed in the last section, we could directly adapt the beam
radiance estimate (3) by replacing photons xi and kernels Ki
with the components g(x|Θs ) of our Gaussian mixture.
All the Gaussian terms, in theory, contribute to the media
radiance along a query ray, but in reality some components
may contribute very little (due to distance to the ray, significant attenuation, or small projected solid angle). The hierarchy constructed in the last section allows us to efficiently
incorporate LOD. The main idea behind our approach is
to traverse the GMM pyramid top-down and hierarchicallyprune the traversal when the interior nodes are deemed areasonable approximation for their subtree. To facilitate thishierarchical pruning, we augment the GMM pyramid withan axis-aligned bounding-box (AABB) hierarchy (BBH).
7.1. Bounding-Box Hierarchy Construction
We spatially bound each Gaussian by computing the distanceat which its remaining energy is a small percentage of thetotal energy of the Gaussian mixture. This truncates eachanisotropic Gaussian to an ellipsoid. We initialize the bounding box of each node to tightly fit its truncated Gaussian.
We then complete the BBH by propagating the boundingbox information up the tree: each node’s bounding box isunioned with the bounding box of its two child nodes. This
construction can be performed efficiently in O(k log k) using
two sweeps of the Gaussian pyramid.
c 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

7.2. Hierarchical Traversal
For fast hierarchical rendering, we prune entire subtrees of
the GMM pyramid during ray traversal by accounting for the
expected contribution of each Gaussian to the rendered image. For a given query ray r, we start at the root. At each node
we decide between three possible outcomes: 1) the entire
subtree can be skipped, 2) the subtree can be well approximated by adding the contribution of the current node to theBRE, or 3) we need to descend further.
We first intersect the ray with the bounding box of thenode. If there is no hit, we skip the entire subtree. Otherwise,
we compute the projected solid angle of the node’s bounding
sphere Ω j (r) and compare it to the solid angle of a pixel ofthe ray Ω(r). If Ω j (r) < εΩ Ω(r), for a user-specified threshold εΩ , we evaluate the contribution of the node’s Gaussian.
If this test fails, we additionally compute the transmittance
between the origin of the ray and the first hit point with thebounding box, Tr (r, Θ j ). If this is less than a user-specified
constant, εTr , we evaluate the node’s Gaussian. Only if allthese tests fail do we descend further down the BBH. These
cases are illustrated in Figure 4 and outlined in Algorithm 4.
Algorithm 4 H IERARCHICAL BRE(r, Θ j )1 if r does not intersect AABB(Θ j )2
return 0
3 if (Ω j (r) < εΩ Ω(r)) or (Tr (r, Θ j ) < εTr )4return Gaussian contribution using Equation (15)
5 return H IERARCHICAL BRE(r, Θ2 j ) +H IERARCHICAL BRE(r, Θ2 j+1 )To evaluate the contribution of a Gaussian, we need tointegrate it along r(t) = x + t~ω while accounting for transmittance.
Lmj (x,~ω, b) =

1
w̄ j
4π

Z b

g(r(t)|Θ̄ j ) Tr (x ↔ r(t)) dt.

(15)

0

This corresponds to a single summand in (3) for isotropic
scattering (we discuss anisotropic scattering in Section 9).
As in Equation (3), neither σs nor the albedo are used here
because we account for these terms during photon tracing.
For homogeneous media, the integral in Equation (15) ac-

W. Jakob, C. Regg, and W. Jarosz / Progressive Expectation–Maximization for Hierarchical Volumetric Photon Mapping
Beam Radiance Estimation [Jarosz et al. 2008]

Our Method

3.6M Photons (133 + 344 = 477 s)
Per-Pixel Render-time
16K Gaussian fit (226 + 100 = 326 s)
Per-Pixel Render-time
Figure 5: In the M USEUM scene, the high volumetric depth complexity and large number of photons (3.1M) makes the BRE (left) inefficient.
Our method addresses these problems by using a hierarchical, anisotropic Gaussian mixture with only 16K Gaussian components.

cepts the closed-form solution: 
C3 + 2C2 b√g(r(t)|Θ̄ j ) e
dt = C0 erf
−
2 C2
a


C3 + 2C2 a
√
erf
,
(16)
2 C2
where the constants C0 through C3 are:C2
exp 4C3 −C1
2
1
q, C1 = (x − µ̄ j )| Σ̄−1
C0 =j (x − µ̄ j ) − σt b,
3
2
2
Σ̄ j
(2π)
Z b

C2 =

−σt t

1 | −1~ω Σ̄ j ~ω,2

C3 = ~ω| Σ̄−1j (x − µ̄ j ) + σt .

(17)

Though we don’t currently support this in our implementation, for heterogeneous media we could take an approachsimilar to the BRE and compute this integral numericallyusing ray marching.
8. Implementation and Results
We implemented our hierarchical GMM rendering approachin a Monte Carlo ray tracer and also implemented the BRE
in the same system to facilitate comparison. All timings arecaptured on a dual-core Intel Core i7 2.80GHz CPU with 6
GB memory and an NVIDIA GeForce GT240.
8.1. GPU Splatting
In addition to the general hierarchical ray tracing solution
described in Section 7.2, we have also implemented a GPU
splatting approach which allows us to interactively preview
arbitrary levels of the GMM mixture.
Each Gaussian component is rendered as an elliptical splat
corresponding to the screen-space projection of the truncatedGaussian ellipsoid. Each splat is rendered as a triangle strip
and Equation (15) is evaluated in a fragment shader with additive blending. This gives correct results for homogeneous
media. While our GPU implementation ignores surface interactions, it nonetheless provides a useful tool to interactively preview the media radiance (refer to the video).

In Figure 1 and 5 we highlight the problems that ourapproach addresses. Both the C ARS and M USEUM scenes
contains several light sources and have a relatively large
depth-extent. We show a rendering with the BRE, and with
our approach, as well as visualizations of the rendering
time spent in each pixel. With the BRE, regions of the image that contain large concentrations of photons take longer
to render, even if these photons contribute little to the final pixel color. Our hierarchical rendering technique solves
this problem by detecting and pruning such low contribution subtrees in the GMM pyramid. The resulting images are higher quality and render faster. Furthermore, since
our pre-computation is view-independent, fly-through animations render very quickly once the Gaussian mixture isconstructed.
Figure 6 shows a progression of renderings of the
B UMPY S PHERE scene with increasing numbers of Gaussiancomponents. Though our compact GMM representation
uses only a few thousand Gaussian components (ranging
from 1K to 64K) we are able to faithfully capture highfrequency caustic structures because we progressively fit
to a large number of photons (up to 18 million). We also
show BRE results which use all the available photons
during rendering. Each BRE result takes longer to renderthan the corresponding GMM rendering (between 7.8 and9.2× slower for the rendering pass, and 1.2 to 6.0× slower
including the precomputation stage). Since the BRE uses
isotropic density estimation, even with all the photons it
produces both noisier and blurrier results. Our anisotropic
GMM fitting procedure acts like an intelligent noisereduction and sharpening filter that automatically finds
structures in the underlying density in a way not possible
with on-the-fly density estimation.

8.2. Results

Figure 7 shows the O CEAN and S PHERE C AUSTIC scenes
depicting detailed volume caustics. These high-frequency
structures are difficult to reconstruct even when using
over 4M isotropic photon points with the BRE, while ouranisotropic fitting process produces slightly better results
1.2× and 1.6× faster using only 16K Gaussian components.

Parameters and rendering statistics for all scenes are listed inTable 2. Note that the preprocessing costs listed for the BRE
and progressive EM methods include the time spent tracing
photons.

The accompanying video shows animations for the
B UMPY S PHERE and O CEAN scenes using both the offline
renderer and the GPU splatting approach. The videos also
highlight the benefit of our temporal coherence as comparedc 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

W. Jakob, C. Regg, and W. Jarosz / Progressive Expectation–Maximization for Hierarchical Volumetric Photon Mapping
BRE All 1.0M photons

All 294K photons(6 + 126 = 132 s)
1K Gaussians(23 + 192 = 215 s)

Our Method 4K Gaussians(7 + 15 = 22 s)

BRE All 4.7M photons(109 + 337 = 446 s)

Our Method 16K Gaussians(35 + 24 = 59 s)

BRE All 18M photons(507 + 609 = 1116 s)

Our Method 64K Gaussians(184 + 39 = 223 s)

BRE

Our Method(868 + 66 = 934 s)

Figure 6: B UMPY S PHERE renderings using our method (bottom) and corresponding BRE renderings (top) using all photons in the Gaussian
fit. As before, the listed times correspond to the preprocessing and rendering stages, respectively.

to the BRE. When animated, we also obtain increased
performance since the previous frame is a closer match tothe converged solution and less EM iterations are required.
The total fitting time for the B UMPY S PHERE animation is
4.35× lower using the temporally coherent method.
9. Discussion and Future Work
Anisotropic Media. A significant limitation of our current
implementation is that we only handle isotropic media. We
believe our general methodology, however, could be naturally extended to handle this case by additionally storing
directional information during fitting. One option is to perform the fitting procedure entirely in a 5D spatio-directional
domain where each component is the tensor product of a spatial Gaussian and a von Mises-Fischer distribution. Anotherpossibility would be to treat the two domains separately andfit von Mises-Fischer distributions to the directional statistics under the support of each 3D Gaussian in a similar
fashion to previous multi-resolution methods for surface reflectance [TLQ∗ 05, HSRG07].
Performance, Quality, and Scalability. Our render passis considerably faster than the BRE due to the reduced
number of terms and hierarchical pruning during traversal.
Yet, we are able to obtain higher quality results with suchfew terms since our Gaussians are anisotropic and inherently less noisy than individual photons. Even including precomputation, our total time to render is faster than the BRE
for the scenes we tried. Furthermore, since our fitting isview-independent, we provide an even greater performance
benefit for walk-throughs of static or pre-scripted lighting.
GMMs are also an attractive low-memory representation
c 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

BRE (73 + 195 = 268 s)

Our Method (178 + 40 = 218 s)

BRE (89 + 638 = 727 s)

Our Method (330 + 127 = 457 s)

Figure 7: The O CEAN and S PHERE C AUSTIC scenes comparing the
BRE with 2M photons to our method with 16K Gaussians.

for volumetric lighting: For instance, a 4096-term GMM requires only 240 kilobytes of storage.
Though our pre-computation stage is slightly slower, thestatistics in Table 2 indicate that progressive EM is actuallyquite competitive (in speed and scalability) with the precomputation needed for the BRE. Nonetheless, our fitting
does incur a performance penalty which may not be warranted in situations where other aspects of our algorithm are
not beneficial (e.g. in very low-frequency illumination, or
scenes with a small depth extent).
Fitting to Other Representations. Recently, Jarosz etal. [JNSJ11] proposed the photon beams method which performs density estimation over photon path segments instead

W. Jakob, C. Regg, and W. Jarosz / Progressive Expectation–Maximization for Hierarchical Volumetric Photon Mapping
Table 2: Parameters and performance statistics for the scenes used in this paper. For BRE the preprocess includes photon tracing, kd-tree andBBH construction, and radius estimation. For our method this includes progressive photon shooting, GMM fitting, and hierarchy construction.
BRE [Jarosz et al. 2008]
Scene

k

n

|P|

EM iter.

C ARS

4K

917,504

87,615M USEUM

16K

3,670,016

B UMPY S PHERE
B UMPY S PHERE
B UMPY S PHERE
B UMPY S PHERE

1K
4K
16K
64K

O CEANS PHERE C AUSTIC

Our Method

Preprocess

Render

Total

Preprocess

RenderTotal

12

13.48 s268.20 s

281.69 s

54.43 s71.21 s (3.8×)

125.64 s (2.2×)

381,60210

133.90 s

343.95 s

477.89 s

225.73 s

100.39 s (3.4×)

326.12 s (1.7×)

294,9121,085,621
4,718,592
18,874,368

30,264
121,910489,602
1,961,891

19
19
20
19

5.73 s23.46 s109.30 s
507.30 s

126.35 s
192.30 s
337.42 s
609.25 s

132.08 s
215.76 s
446.72 s
1116.55 s

6.81 s35.13 s184.06 s
868.48 s

15.19 s (8.3×)
24.52 s (7.8×)
39.35 s (8.6×)
66.03 s (9.2×)

22.00 s (6.0×)
59.65 s (3.6×)
223.41 s (2.0×)
934.51 s (1.2×)

16K

4,718,592

485,33918

73.47 s195.46 s

268.93 s

178.25 s

40.67 s (4.8×)

218.92 s (1.2×)

16K

4,194,304

458,36611

89.55 s638.39 s

727.90 s

329.83 s127.30 s (5.0×)

457.13 s (1.6×)of photon points. They showed that this results in improved
rendering quality compared to the BRE. Their method shares
similar goals to ours (more compact representation, andhigher quality reconstruction); however, their benefits is actually orthogonal to ours. Our Gaussian fitting pipeline could
be applied in the context of photon beams as well as photon
points. A straightforward combination could update statistics in all partition cells intersected by a photon beam. This
could significantly improve the quality of the statistics (andthe EM fit) compared to photon points. We believe that inhomogeneous media the statistical contributions of a beam
could be evaluated analytically and plan to investigate thisin future work. Another, more ambitious, avenue would beto develop a fitting procedure analogous to EM which operates more explicitly on photon beams, by clustering intohigher-dimensional density distributions such as blurred 2Dline segments.
Convergence and Automatic Termination. Though thename of our progressive EM algorithm is inspired by progressive photon mapping (PPM) [HOJ08] there are important differences. Our method progressively fits a Gaussian mixture to photon statistics whereas PPM progressively
refines radiance estimates by updating statistics at fixed
measurement points. Our Gaussian components could beroughly interpreted as “measurement points” which not only
adapt in intensity, but also location, to the underlying photon density. However, whereas PPM provides guarantees forconvergence of both bias and noise, our solution eliminates
noise in the limit but is always biased since we use a fixed
number of Gaussian terms. An interesting possibility for future work would be to progressively split Gaussian terms
when sufficient statistics indicate a single Gaussian cannot
represent the underlying distribution, or to automatically deduce the required number of Gaussian terms.
Automatic termination by specifying a desired error in
the image has been investigated previously in the context
of PPM [HJJ10]. While our progressive EM algorithm also
provides an automatic stopping criterion for shooting photons, this tolerance does not have an intuitive relation to theerror in the image, but rather describes the convergence of
the GMM fit to the underlying density distribution. Error isfairly well understood in the context of kernel-based density
estimation used in photon mapping. A more rigorous investigation of the error and bias introduced by our method is aninteresting avenue of future work.
10. Conclusion
We have presented a progressive, hierarchical method for
rendering participating media. Our approach can be viewed
as a variant of volumetric photon mapping which replaces
the usual on-the-fly non-parametric density estimate witha parametric density model based on a Gaussian mixture
fit to the photons. The fitting is done using a sparse andprogressive form of the accelerated EM algorithm. We construct a hierarchical representation which is used to provide
estimates of scattered radiance at appropriate scales during
rendering. Our approach shows improved quality and temporal coherence as compared to the BRE and also supports
interactive preview using splatting.
Acknowledgments. We thank Steve Marschner, PeterPike Sloan, and Derek Nowrouzezahrai for helpful discussions and comments. The M USEUM scene is courtesy of Alvaro Luna Bautista and Joel Anderson; the B UMPY S PHERE
scene is courtesy of Bruce Walter.
References[CB04] C HRISTENSEN P. H., BATALI D.: An irradiance atlas for
global illumination in complex production scenes. In Rendering
Techniques (June 2004), pp. 133–142. 2[Cha60] C HANDRASEKHAR S.: Radiative Transfer. Dover Publications, New York, 1960. 2[CPP∗ 05] C EREZO E., P ÉREZ F., P UEYO X., S ERON F. J., S IL LION F. X.: A survey on participating media rendering techniques. The Visual Computer 21, 5 (2005). 1[DBB06] D UTRÉ P., BALA K., B EKAERT P.: Advanced global
illumination. AK Peters Ltd, 2006. 2[DLR∗ 77] D EMPSTER A., L AIRD N., RUBIN D., ET AL .: Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal Statistical Society. 39, 1 (1977), 1–38. 2[DYN04] D OBASHI Y., YAMAMOTO T., N ISHITA T.: Radiosity
for point-sampled geometry. In 12th Pacific Conference on Computer Graphics and Applications (Oct. 2004), pp. 152–159. 2[GNN10] G ARCIA V., N IELSEN F., N OCK R.: HierarchicalGaussian mixture model. In Proceedings of the IEEE Interna-

c 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

W. Jakob, C. Regg, and W. Jarosz / Progressive Expectation–Maximization for Hierarchical Volumetric Photon Mapping
tional Conference on Acoustics, Speech, and Signal Processing(2010), pp. 4070–4073. 2, 6, 7
[GR05] G OLDBERGER J., ROWEIS S.: Hierarchical clustering ofa mixture model. Advances in Neural Information Processing
Systems 17, 505-512 (2005), 2–4. 2, 7
[HE03] H OPF M., E RTL T.: Hierarchical splatting of scattereddata. In Proceedings of VIS (2003). 2[HG09] H ASAN B. A. S., G AN J.: Sequential EM for unsupervised adaptive Gaussian mixture model based classifier. In
Machine Learning and Data Mining in Pattern Recognition,
vol. 5632. 2009, pp. 96–106. 2[HJJ10] H ACHISUKA T., JAROSZ W., J ENSEN H. W.: A progressive error estimation framework for photon density estimation.
ACM Transactions on Graphics. 29 (December 2010). 10[HLE04] H OPF M., L UTTENBERGER M., E RTL T.: Hierarchical splatting of scattered 4d data. IEEE Computer Graphics andApplications 24 (July 2004), 64–72. 2[HOJ08] H ACHISUKA T., O GAKI S., J ENSEN H. W.: Progressive photon mapping. ACM Transactions on Graphics 27, 5 (Dec.
2008), 130:1–130:8. 2, 10[HSA91] H ANRAHAN P., S ALZMAN D., AUPPERLE L.: A rapid
hierarchical radiosity algorithm. In Computer Graphics (Proceedings of SIGGRAPH 91) (July 1991), pp. 197–206. 2[HSRG07] H AN C., S UN B., R AMAMOORTHI R., G RINSPUN
E.: Frequency domain normal map filtering. ACM Transactions
on Graphics. 26, 3 (2007). 2, 9
[JC98] J ENSEN H. W., C HRISTENSEN P. H.: Efficient simulationof light transport in scenes with participating media using photon
maps. In Proceedings of SIGGRAPH. (July 1998). 1, 3
[JNSJ11] JAROSZ W., N OWROUZEZAHRAI D., S ADEGHI I.,J ENSEN H. W.: A comprehensive theory of volumetric radianceestimation using photon points and beams. ACM Transactions on
Graphics 30, 1 (Jan. 2011), 5:1–5:19. 3, 9
[JZJ08] JAROSZ W., Z WICKER M., J ENSEN H. W.: The beam
radiance estimate for volumetric photon mapping. Computer
Graphics Forum 27, 2 (Apr. 2008). 1, 3
[LH91] L AUR D., H ANRAHAN P.: Hierarchical splatting: A progressive refinement algorithm for volume rendering. In Proceedings of SIGGRAPH. (July 1991), pp. 285–288. 2[LW93] L AFORTUNE E. P., W ILLEMS Y. D.: Bi-directional pathtracing. In Compugraphics (1993). 1[LW96] L AFORTUNE E. P., W ILLEMS Y. D.: Rendering participating media with bidirectional path tracing. In Eurographics
Rendering Workshop (June 1996), pp. 91–100. 1[LZT∗ 08] L EHTINEN J., Z WICKER M., T URQUIN E., KONTKA NEN J., D URAND F., S ILLION F. X., A ILA T.: A meshlesshierarchical representation for light transport. ACM Transactions
on Graphics 27, 3 (Aug. 2008), 37:1–37:9. 2[NH98] N EAL R., H INTON G.: A view of the EM algorithm
that justifies incremental, sparse, and other variants. Learningin graphical models 89 (1998), 355–368. 2, 4
[PJJ∗ 11] PAPAS M., JAROSZ W., JAKOB W., RUSINKIEWICZ S.,M ATUSIK W., W EYRICH T.: Goal-based caustics. Computer
Graphics Forum (Proceedings of Eurographics ’11) 30, 2 (June2011). 2[PKK00] PAULY M., KOLLIG T., K ELLER A.: Metropolis lighttransport for participating media. In Eurographics Workshop on
Rendering (June 2000), pp. 11–22. 1[RL00] RUSINKIEWICZ S., L EVOY M.: Qsplat: A multiresolution point rendering system for large meshes. In Proceedings of
ACM SIGGRAPH. (July 2000), pp. 343–352. 2
c 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

[SFES07] S CHJØTH L., F RISVAD J. R., E RLEBEN K.,
S PORRING J.: Photon differentials. In GRAPHITE. (2007),
pp. 179–186. 2[Sil86] S ILVERMAN B. W.: Density Estimation for Statistics andData Analysis. Monographs on Statistics and Applied Probability. Chapman and Hall, New York, NY, 1986. 2[Sil95] S ILLION F. X.: A unified hierarchical algorithm for
global illumination with scattering volumes and object clusters.
IEEE Transactions on Visualization and Computer Graphics 1, 3
(1995). 2[SJ09a] S PENCER B., J ONES M. W.: Hierarchical photon mapping. IEEE Transactions on Visualization and Computer Graphics 15, 1 (Jan./Feb. 2009), 49–61. 2[SJ09b] S PENCER B., J ONES M. W.: Into the blue: Better caustics through photon relaxation. Computer Graphics Forum 28, 2(Apr. 2009), 319–328. 2[SSO08] S CHJØTH L., S PORRING J., O LSEN O. F.: Diffusion
based photon mapping. Computer Graphics Forum 27, 8 (Dec.
2008), 2114–2127. 2[TLQ∗ 05] TAN P., L IN S., Q UAN L., G UO B., S HUM H.-Y.:
Multiresolution reflectance filtering. In Proceedings Eurographics Symposium on Rendering 2005 (2005), pp. 111–116. 2, 9
[VG94] V EACH E., G UIBAS L.: Bidirectional estimators for lighttransport. In Eurographics Workshop on Rendering (1994). 1[VNV06] V ERBEEK J. J., N UNNINK J. R., V LASSIS N.: Accelerated EM-based clustering of large data sets. Data Mining andKnowledge Discovery 13, 3 (November 2006), 291–307. 2, 4
[WABG06] WALTER B., A RBREE A., BALA K., G REENBERG
D. P.: Multidimensional lightcuts. ACM Transactions on Graphics. 25, 3 (2006). 1, 2
[WBKP08] WALTER B., BALA K., K ULKARNI M., P INGALI
K.: Fast agglomerative clustering for rendering. In IEEE Symposium on Interactive Ray Tracing (August 2008), pp. 81–86. 6[WFA∗ 05] WALTER B., F ERNANDEZ S., A RBREE A., BALA
K., D ONIKIAN M., G REENBERG D. P.: Lightcuts: a scalable
approach to illumination. ACM Transactions on Graphics. 24, 3(Aug. 2005). 2[ZHG∗ 07] Z HOU K., H OU Q., G ONG M., S NYDER J., G UO B.,
S HUM H.-Y.: Fogshop: Real-time design and rendering of inhomogeneous, single-scattering media. In Pacific Graphics. (2007),
pp. 116–125. 2[ZPvBG02] Z WICKER M., P FISTER H., VAN BAAR J., G ROSS
M.: EWA splatting. IEEE Transactions on Visualization andComputer Graphics 8, 3 (July/Sept. 2002), 223–238. 2[ZRL∗ 08] Z HOU K., R EN Z., L IN S., BAO H., G UO B., S HUM
H.-Y.: Real-time smoke rendering using compensated ray
marching. ACM Transactions on Graphics. 27, 3 (2008). 2, 6