Pacific Graphics 2011
Bing-Yu Chen, Jan Kautz, Tong-Yee Lee, and Ming C. Lin(Guest Editors)

Volume 30 (2011), Number 7

Toward Optimal Space Partitioning for Unbiased, AdaptiveFree Path Sampling of Inhomogeneous Participating Media
Yonghao Yue1 , Kei Iwasaki2 , Bing-Yu Chen3 , Yoshinori Dobashi4 and Tomoyuki Nishita1
1 The University of Tokyo 2 Wakayama University 3 National Taiwan University 4 Hokkaido University
1 {yonghao,nis}@nis-lab.is.s.u-tokyo.ac.jp 2 iwasaki@sys.wakayama-u.ac.jp 3 robin@ntu.edu.tw 4 doba@ime.ist.hokudai.ac.jp

Abstract
Photo-realistic rendering of inhomogeneous participating media with light scattering in consideration is important in computer graphics, and is typically computed using Monte Carlo based methods. The key technique insuch methods is the free path sampling, which is used for determining the distance (free path) between successive
scattering events. Recently, it has been shown that efficient and unbiased free path sampling methods can be constructed based on Woodcock tracking. The key concept for improving the efficiency is to utilize space partitioning(e.g., kd-tree or uniform grid), and a better space partitioning scheme is important for better sampling efficiency.
Thus, an estimation framework for investigating the gain in sampling efficiency is important for determining how
to partition the space. However, currently, there is no estimation framework that works in 3D space. In this paper,
we propose a new estimation framework to overcome this problem. Using our framework, we can analyticallyestimate the sampling efficiency for any typical partitioned space. Conversely, we can also use this estimation
framework for determining the optimal space partitioning. As an application, we show that new space partitioning schemes can be constructed using our estimation framework. Moreover, we show that the differences in theperformances using different schemes can be predicted fairly well using our estimation framework.
Categories and Subject Descriptors (according to ACM CCS): I.3.7 [Computer Graphics]: Three-Dimensional
Graphics and Realism—I.3.3 [Computer Graphics]: Picture/Image Generation—G.3 [Probability and Statistics]:
Probabilistic Algorithms—

1. Introduction
Steam, water, fire, smoke, explosions, volcanic eruptions,clouds, atmosphere, mist due to waterfalls and splashes due
to ocean waves are common participating media aroundus. These participating media are usually inhomogeneous.
Photo-realistic rendering of these inhomogeneous participating media is important in computer graphics, as the renderedresults are usually impressive. Indeed, the demand for highquality rendering of those media is dramatically increasingin, e.g., film industry.
Light scattering is essential in participating media rendering. Typically, Monte Carlo based methods [LW96, JC98,
PKK00, RSK06] are used for solving the radiative transport
equation [Cha50] to account for the light transport in participating media. A Monte Carlo based method generates anumber of light paths to simulate the light transfer in thescene. The technique at the heart of such a method is thec 2011 The Author(s)Computer Graphics Forum c 2011 The Eurographics Association and Blackwell Publishing Ltd. Published by Blackwell Publishing, 9600 Garsington Road, Oxford OX4 2DQ,
UK and 350 Main Street, Malden, MA 02148, USA.

free path sampling, which is used for determining the distance (free path) between successive scattering events. The
computation efficiency of participating media rendering andthe quality of the resulting images highly depend on the freepath sampling technique used.
In the computer graphics field, ray marching is generally
used as the technique for free path sampling. However, as
figured out by previous works, e.g., Raab et al. [RSK06], raymarching results in a stochastically biased solution which
does not converge to the exact solution. This bias depends
on the sampling interval used in ray marching, and different
sampling intervals result in different pixel colors.
On the other hand, an unbiased free path sampling technique called Woodcock tracking [WMHL65] was proposedin the nuclear science field. However, Woodcock tracking
is known to become less efficient for more inhomogeneousparticipating media [Lep07]. To overcome this problem, re-Y. Yue et al. / Toward Optimal Space Partitioning for Unbiased, Adaptive Free Path Sampling

cently, unbiased acceleration techniques based on space partitioning were introduced. Yue et al. [YIC∗ 10] used a kdtree, and Szirmay-Kalos et al. [SKTM11] used a uniform
grid. The key concept of these methods is to partition thespatial domain according to the spatial variation of the mean
free path of the participating medium and to sample the freepath in an adaptive manner.
In these unbiased, adaptive free path sampling techniques,
the structure of the partitioned space influences the sampling efficiency, and an estimation framework for investigating the gain in sampling efficiency is important in order toobtain a good structure. However, currently, there is no estimation framework that works in the 3D space. For example, in [YIC∗ 10], the space partitioning scheme they used isbased on a heuristic approach that approximately solves the
3D problem by a set of disjoint 1D problems. As a result, theresulting sampling efficiency is not always optimal.
In this paper, we show for the first time an estimation
framework that works in 3D space. Using our new framework, we can analytically estimate the sampling efficiency
given a particular space partitioning. Conversely, we can also
use this estimation framework for determining the optimal
space partitioning. Our framework is derived from the formulation of the 1D problem presented in [YIC∗ 10], and canbe regarded as a generalization of it. Thanks to the generality, our estimation framework is the first one that can beused to develop automatic space partitioning schemes basedon any typical spatial subdivision structures. In this paper,
we show new automatic space partitioning schemes aiming
at optimal sampling efficiency, using a uniform grid, octree
and kd-tree. Moreover, we show that the differences in theperformances using different schemes can be predicted fairly
well using our estimation framework.
2. Related WorkIn this section, we focus on previous work related to freepath sampling. For previous research on participating media rendering, the reader may refer to the surveys [CPP∗ 05,
GJJD09].
The techniques for free path sampling can be classified
into the following two categories: 1) ray marching or its
variants [PH89, JC98, PKK00, BM03], and 2) Woodcock
tracking or its variants [WMHL65, CCW72, RSK06, Lep07,
YIC∗ 10, SKTM11].
Although ray marching is widely used in the computer
graphics field, it produces a biased solution which does notconverge to the exact solution. The bias due to ray marching
could be problematic because 1) different sampling intervals
would result in different results and 2) the error between abiased solution and the exact solution is easily noticeable
but is difficult to be estimated in advance.
An alternative that is free from the aforementioned bias

induced problems is Woodcock tracking [WMHL65], which
was proposed in the nuclear science field and was first introduced to the computer graphics field by Raab et al. [RSK06].
Woodcock tracking is proven to be unbiased by Coleman [Col68], thus will produce exact solutions. A shortcoming of Woodcock tracking is that it becomes less efficient for more inhomogeneous participating media [Lep07].
Badal and Badano [BB09] accelerated Woodcock tracking
by implementing it on the GPU. Another improvement wasdone by Leppänen [Lep07], in which a two-level spacepartitioning was used. Very recently, much more efficienttechniques using more general space partitioning schemes
were proposed. Yue et al. [YIC∗ 10] and Szirmay-Kalos etal. [SKTM11] proposed to use a kd-tree and a uniform grid,
respectively.
In these space partitioning based methods, it is important
to obtain a good structure of the partitioned space as thespace partitioning influences the sampling efficiency. Moreover, we claim that an estimation technique to investigate thegain in sampling efficiency is important for determining how
to partition the space. However, an estimation frameworkthat works in 3D space is lacking. For example, when using
uniform grids, the user needs to determine the resolution of
the grid manually. The kd-tree based method offers an automatic space partitioning, but is inherently based on a heuristic approach that utilizes 1D formulations [YIC∗ 10]. As aresult, the method is imperfect, and the partitioned space is
not guaranteed to be optimal. Thus, in this paper, we focus
on deriving an estimation framework that works in 3D space,aiming at optimal space partitioning.
3. Free Path Sampling
In this section, we first state the free path sampling problem and briefly review Woodcock tracking and the conceptbehind previous methods utilizing space partitioning.
To generate a light path in a Monte Carlo based rendering method, scattering events are generated successively. To
generate the (i + 1)-th scattering event, the free path di andthe scattering direction ~ωi need to be determined throughstochastic sampling. The scattering direction can be obtainedusing a conventional importance sampling technique. For thefree path, we need to sample di according to the following
probability density function [PKK00]:
pd f f p (xi+1 = xi + di~ωi ) = e−τ(xi ,xi+1 ) k(xi+1 ),(1)

where xi and xi+1 are the Rlocations of the scattering eventsi and (i + 1), τ(xi , xi+1 ) = xxii+1 k(x′ )dx′ is the optical depth
between xi and xi+1 , and k is the extinction coefficient.
Woodcock tracking samples the free path by employinga rejection sampling technique. First, a majorant extinction
coefficient kM which is never less than the extinction coefficient of the participating medium is computed. Then, Woodcock tracking samples pseudo scattering events by regarding the medium as a uniform medium with the extinction
c 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

Y. Yue et al. / Toward Optimal Space Partitioning for Unbiased, Adaptive Free Path Sampling

coefficient being kM . For unbiasedness, such pseudo scattering events are only accepted as ‘real’ scattering events withthe probability k(xi + di~ωi )/kM . When generating a newpseudo scattering event, the free path di is incremented by− ln(1 − rand())/kM , whose expectation value is 1/kM . The
reader may also refer to [RSK06] for implementing Woodcock tracking.
Woodcock tracking becomes less efficient in more inhomogeneous participating medium for the following reasons.
In an inhomogeneous participating medium, k(xi + di~ωi ) isoften much smaller than the majorant extinction coefficientkM , and the ratio k(xi + di~ωi )/kM becomes small where themedium is sparse. Therefore, pseudo scattering generation
will be iterated many times until a real scattering event is
detected.
Previous methods utilizing space partitioning improved
Woodcock tracking by subdividing the spatial domain into
multiple domains and used different majorant extinction coefficients for different partitioned domains. Unbiasedness is
ensured as follows: when going across a partition, the sampling location is rewound back to the intersection point between the ray and the partition. Note that this rewinding operation requires an additional iteration.
The way toward finding the optimal space partitioning(i.e., the space partitioning that makes the sampling efficiency optimal) is to minimize the number of iterations
needed to generate a ‘real’ scattering event, taking into account the trade-off between the following two aspects: 1)
if we partition the space more finely, the ratio of k(xi +di~ωi )/kM in each partitioned domain would be closer to its
upper bound (i.e., 1), and the number of iterations insideeach partition would be decreased; 2) the number of the partitions should be kept small as we need additional iterations
to go across these partitions due to the rewind process.
To account for the above trade-off, we need an estimation
framework of the number of iterations, given a space partitioning. Currently, the only analytic estimation framework is
the one formulated in 1D space [YIC∗ 10]. In this paper, we
present a generalized version of their estimation framework.
Our framework is fully formulated in 3D space.
4. Evaluating the Partitioned SpaceIn this section, we first show an evaluation framework for 1D
space, which is a slightly modified version of the formulation
given by Yue et al. [YIC∗ 10]. Then, we generalize it to 3D
space.

M,I4
M,I1
M,I2

M,I5
M,I3

I1

I2
P1

I3
P2

I4
P3

I5
P4

Figure 1: Example distributions of a participating medium
in 1D. The horizontal axis shows the locations x, and thevertical axis shows the value of k(x).

In 1D space (assuming the x axis), the rays travel only
along the x axis, as shown in Figure 1. Suppose that wewant to estimate the expected number of iterations when aray travels through the interval (s,t], and k(x) is given at anarbitrary location x. For simplicity, we assume that a ‘real’
scattering event does not happen in the interval (s,t].
For the interval (s,t], we assume that we are given a spacepartitioning, which is represented by a set of n subintervals,
I j , where j = 1, ..., n, and the majorant extinction coefficientin a subinterval I j is given as kM,I j . Between any two adjacent subintervals I j and I j+1 , there is a partition Pj . Whenthe ray travels through the interval (s,t], the ray will passthrough the subintervals and the partitions in alternate order.
To estimate the expected number of iterations, we accountfor the expected numbers of iterations needed for passing
through each subinterval and partition as follows.
Let |I j | be the length of the subinterval I j . Since the majorant extinction coefficient in this subinterval is kM,I j , theexpected distance that the ray will proceed in a single iteration is 1/kM,I j . Thus, the expected number of iterations totravel through this subinterval is given by |I j |kM,I j . To go
across a partition, we need to perform a single rewind process. Thus, the expected number of iterations to go acrossa partition is 1. Combining these two results, the expected
number of iterations N̄ needed to travel the interval (s,t] canbe formulated asnN̄ = ∑ |I j |kM,I j + (n − 1).

(2)

j=1

Similarly, we may also account for the expected computation
time. Let titer and trewind be the average computation time
for a single iteration inside a subinterval and for a rewinding
process, respectively. Then, the expected computation timeT̄ for the ray to travel the interval (s,t] can be formulated asnT̄ = titer ∑ |I j |kM,I j + trewind (n − 1).

(3)

j=1

4.1. Estimation Framework in 1D SpaceThe efficiency of the free path sampling is tightly coupledwith the average numbers of the iterations needed before therays encounter ‘real’ scattering events. Therefore, it is important to estimate the expected number of iterations.
c 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

4.2. Estimation Framework in 3D SpaceIn 3D space, we need to account for any rays traveling in arbitrary directions. Thus, the formulation in 3D space is much
more difficult than that in 1D space.

Y. Yue et al. / Toward Optimal Space Partitioning for Unbiased, Adaptive Free Path SamplingRΓ(E j ) kM,E j lE j (r)dr. For simplicity, let us omit thesubscription E j in kM,E j and lE j . By calculating the integral
first according to the directions of the rays, we obtain
Z
Z

Z Z
l(~ω, u, v)dudv dω, (6)kMkM l(r)dr =

(b)

(a)
j

j
j

S2

Γ(E j )D j (~ω)2

Figure 2: (a): An example showing the space E j and its projection onto the direction ~ω. (b): Another view of (a) (theviewing ray is parallel to H(~ω)).

To account for the rays in all the directions, we consider
the average of the expected numbers of iterations for all therays passing through a subspace or a partition. Since a typical space partitioning scheme uses axis-aligned planes, we
assume the subspaces and the partitions are rectangular parallelepipeds and planes, respectively.
Let the bounding box of the participating medium be E.
We assume that a space partitioning is given so that E is subdivided into ns subspaces E j ( j = 1, ..., ns ) by n p partitions
Pk (k = 1, ..., n p ). We estimate the average of the expected
numbers of iterations N̄E j for the rays passing through a subspace E j asZ.ZN̄E j =
kM,E j lE j (r)drdr,(4)
Γ(E j )

Γ(E)and the average of the expected numbers of iterations N̄Pk forthe rays passing through a partition Pk as
Z
.Z
dr,
(5)
1 · drN̄Pk =
Γ(Pk )

Γ(E)where Γ(E j ), Γ(E) and Γ(Pk ) indicate the sets of rays passing through E j , E and Pk , respectively. kM,E j is the majorant extinction coefficient of the subspace E j , lE j (r) denotes
the length of the intersection between the ray r and the subspace E j , and dr is the measure of the ray r. The meaningsof Eqs.(4) and (5) are as follows. The numerators in themdescribe the summation of the numbers of iterations for all
the rays, and the denominators describe the size of the set ofrays passing through E. Thus, by the division, we obtain the
averages of the numbers of iterations for a single ray passing
through E.

One of the important contributions of this paper is that the
integrals appearing in Eqs.(4) and (5) can be calculated analytically by assuming 1) all the rays pass through the space
E and do not encounter ‘real’ scattering events and 2) therays are distributed uniformly in the space. We show the analytical solutions below.
we

show

ZZ

D j (~ω)l(~ω, u, v)dudv = |E j |,theanalyticalsolutionof(7)because l(~ω, u, v)dudv represents an infinitesimal pillar
which is sliced from E j in the direction ~ω as shown in Figure 2(b), and the integral of such pillars is identical to thevolume of E j . Therefore, we obtainZΓ(E j )kM l(r)dr = 4πkM |E j |.

(8)R

Next, we show the analytical solution of Γ(E) dr. By calculating the integral first according to the directions of therays, we obtainZ
Z Z ZZ|D(~ω)|dω, (9)dudv dω =
dr =
S2

Γ(E)

D(~ω)

S2

where |D(~ω)| denotes the area of D(~ω), which is the projected region of the space E onto a plane H(~ω) perpendicular to the direction ~ω. |D(~ω)| can be computed as follows.
Let ~ex , ~ey and ~ez be the unit vectors along the x, y and z axes,
respectively, and let Sx , Sy and Sz be the areas of the faces ofE perpendicular to the x, y and z axes, respectively. Then,|D(~ω)| = Sx |~ex ·~ω| + Sy |~ey ·~ω| + Sz |~ez ·~ω|.

(10)

Let us use the polar coordinates (θ, φ) to describe ~ω, i.e.,~ω = (sin θ cos φ, cos θ, sin θ sin φ). Then,
Z

dr =

Z

=8

Z πZ π(Sx sin θ cos φ + Sy cos θ + Sz sin θ sin φ) sin θdφdθ

Γ(E)

2

0

2

0

S2

|D(~ω)|dω= 2π(Sx + Sy + Sz ) = πS(E),4.3. Analytic Solution

First,

where S denotes the set of all the directions, dω is themeasure of the solid angle for the direction ~ω. As shownin Figure 2(a), D j (~ω) is the projected region of the subspace E j onto a plane H(~ω) which is perpendicular to thedirection ~ω. u and v are the orthogonal coordinates in H(~ω).
l(ω, u, v) = l(r) is the length of the intersection between E j
and the ray in the direction ~ω passing through a point (u, v)
on H(~ω). Let |E j | be the volume of E j , then the following
formula holds,(11)

where S(E) denotes the surface area of E.
R

Finally, we show the analytical solution of Γ(Pk ) 1 · dr. By
calculating the integral first according to the directions of therays, we obtainZ
Z Z ZZ
1 · dr =dudv dω =|Pk,⊥ (~ω)|dω,Γ(Pk )

S2

Pk,⊥ (~ω)

S2

(12)

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

Y. Yue et al. / Toward Optimal Space Partitioning for Unbiased, Adaptive Free Path Samplingwhere Pk,⊥ (~ω) represents the projected region of Pk ontoa plane perpendicular to ~ω, and |Pk,⊥ (~ω)| denotes its area.
Then, using the polar coordinates for ~ω to compute the integral, we obtainZ

Γ(Pk )

1 · dr =Z

S2|Pk,⊥ (~ω)|dω = 2π|Pk |,(13)

minimum from all possible m. Typically, we choose m as a
power of two.
To determine whether the partitioning according to thisresolution is beneficial, we also estimate the average of expected numbers of iterations N̄nopart without partitioning as,

where |Pk | is the area of Pk .

N̄nopart =Then, by putting Eqs.(4), (5), (8), (11) and (13) all together, we obtain an estimator for the average of the expected
numbers of iterations N̄ for a given space partitioning as
nsnpN̄ = ∑ N̄E j + ∑ N̄Pk =
j=1

k=1

np

 ns

∑ 4kM,E |E j |+ ∑ 2|Pk |
j

j=1

k=1

.

S(E),

(14)
where ns and n p are the numbers of the partitioned domains(subspaces) and partitions, respectively. Note that N̄ is a dimensionless value because kM,E j has the dimension m−1 ,|E j | has the dimension m3 , and |Pk | and S(E) have the dimension m2 , respectively. The average computation time T̄
can be estimated similarly as
np

.
ns
T̄ = titer · ∑ 4kM,E j |E j | +trewind · ∑ 2|Pk | S(E). (15)
j=1

k=1

An important contribution of the above formulations is thatthey do not contain any heuristic user specified parameters.
They can be computed strictly.
5. On Optimal Space Partitioning Schemes
In this section, we apply our estimation framework to thedevelopment of new schemes aiming at optimal partitioning,
using the uniform grid, octree and kd-tree, respectively.
5.1. Uniform Grid based Scheme
When using a uniform grid, we need to decide its resolution.
That is, from a variety of the possible resolutions, we need tocheck for 1) the best resolution and 2) whether partitioning
the space according to this resolution is beneficial.
First, we show that we can use Eq.(14) to estimate the average of the numbers of iterations. For simplicity, we assume
that the simulation space (the bounding box) is a cube, andthe length of each side is w. We also assume that each side
is subdivided into m slices so that we have m × m × m gridcells. Then, in Eq.(14), the volumes of the grid-cells are thesame, thus |E j | = (w/m)3 , the summation of the areas of thenp
2|Pk | = 3(m − 1) · 2w2 , and S(E) = 6w2 .
partitions is ∑k=1
Thus, we have3N̄ part =2w 1 m
·∑ kM,E j + (m − 1).
3 m3 j=1(16)

To find the best resolution, we find m which makes N̄ partc 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

2w
kM,E .
3(17)

If N̄ part < N̄nopart , then it is beneficial to partition the simulation space according to the best resolution. Note that theabove scheme can be straightforwardly extended to handle asimulation space represented by a rectangular parallelepipedor uneven number of slices for each side.
5.2. Octree based Scheme
We consider a partitioning scheme in a top-down and recursive style. We let the root node represent the simulation
space (the boundary box). For the root node, first, we estimate whether it is beneficial to partition this node. If we
decide to partition this node, then we recursively continue
the partitioning for the eight child nodes.
For simplicity, again, we assume the simulation space isa cube, with the length of each side being w. We let l be thelevel (or depth) of a node, e.g., the root node is in the level
0. We partition a node by the axis-aligned planes containing
the center of the node, so that all the eight child nodes havethe same shape and size. Therefore, a node in the level l is acube, with the length of each side being w/2l .
First, we estimate the average of the expected numbers of
iterations N̄ part for the case with partitioning. Assume thatthe parent node is in the level l, then the volume of eachchild node is (w/2l+1 )3 , the summation of the areas of thepartitioning planes is 3(w/2l )2 , and the surface area of thesimulation space S(E) is 6w2 . Thus, we haveN̄ part =1  2w 1 8
+
1
,
k
·
M,E
∑
j
4l 3 · 2l 8 j=1(18)

where kM,E j denotes the majorant extinction coefficient for
the eight child nodes. Next, we estimate the average for thecase without partitioning as
1 2w
(19)
kM,E p ,4l 3 · 2l
where kM,E p denotes the majorant extinction coefficient forthe parent node. If N̄ part < N̄nopart , then it is beneficial topartition the node.
N̄nopart =Sometimes (e.g., a participating medium with the high extinction coefficient portion coagulated around the center of anode), even if N̄ part < N̄nopart is not satisfied, the average ofthe numbers of iterations can be decreased if we continue
partitioning. To handle such a case, we introduce the failcount, which is used when deciding the spatial subdivision
structure for standard ray tracing.

Y. Yue et al. / Toward Optimal Space Partitioning for Unbiased, Adaptive Free Path Sampling

(a)

(b)

Figure 3: (a) and (b): Before and after removing redundantpartitioning. For simplicity, the number of child nodes is reduced to 2 instead of 8 in this illustration. Blue and red circles indicate the leaf nodes and the nodes with the fail countincremented, respectively. If the child nodes of a parent node(marked by the green arrow) are all leaf nodes and the failcount of the parent node is incremented, it is beneficial tomerge the leaf nodes because N̄nopart ≥ N̄ part . The merge
step can be processed recursively, e.g., for the node marked
by the navy arrow. Conversely, if the child nodes of a parent
node (marked by the purple arrow) are not all leaf nodes, it
may be beneficial not to merge the child nodes, as there isbenefit in the subsequent partitioning.

A fail count is an integer value kept in each node during
the partitioning process. Initially, the fail count of the root
node is set to 0. During partitioning, if N̄ part < N̄nopart is notsatisfied, we increment the fail count by 1 instead of terminating the partitioning process. The fail counts of the childnodes are set equal to that of the parent node. If the fail countof a node exceeds a threshold, we terminate partitioning.
After the octree is constructed, we try to remove redundant space partitioning (see Figure 3). That is, if the leaf
nodes of a parent node are generated with the fail countincreased, we can merge these leaf nodes. We perform thismerge process recursively in a bottom-up fashion to remove
all the redundant partitionings.

a fairly good sampling efficiency as demonstrated in Section 6, the partitioning scheme is currently time consuming
because of the following two reasons. First, the number ofpossible subspaces is large. Assume that we use an auxiliary grid with n grid-cells to search for the subspaces, thenwe have approximately O(n2 ) possible subspaces. Second,
the number of possible partitioning for a single subspace isalso large (e.g. if a subspace is fully contained in the space
E p , there are 426 different ways for partitioning. The number 426 can be obtained as follows. Distinct partitioning canbe enumerated by considering how the space is subdivided
by the partitioning planes. For each axis, there are two suchplanes, which we describe as x+ , x− , y+ , y− , z+ and z− .
Then, a partitioning can be described as a permutation ofthese labels, like y+ x− x+ z− y− z+ . The possible cases in the
permutation, however, include duplicated cases, which canbe omitted. That is, if two planes for the same axis, e.g., x+
and x− , are neighboring, transpose the order of them would
result in the same partitioning.)
Currently, we are limited to using an auxiliary grid which
has equal to or less than 163 grid-cells. This is because thecomputation time for space partitioning becomes expensive
for higher resolutions. The computation times are approximately 2 seconds for 83 grid-cells, 2.5 minutes for 163 gridcells, and 3 hours for 323 grid-cells. We believe we can overcome this limitation in future research.
6. Evaluations
We evaluate the sampling efficiencies for the partitioned
spaces obtained using the uniform grid based scheme (abbreviated as UG), the octree based scheme (OC), the kd-tree
based scheme using the heuristic approach [YIC∗ 10] (hKD)
and the kd-tree based scheme presented in this paper (KD).

Our kd-tree based scheme is also in a top-down and recursive style. Given a space (let us call the space E p ), we checkfor each subspace in E p and estimate the benefit when partitioning the space using the planes containing the faces of thesubspace. The average of the expected numbers of iterations
with partitioning N̄ part can be estimated using Eq.(14). From
all the possible subspaces, we find one that makes N̄ part minimum. Note that for a single subspace, there are many different ways for the partitioning. The summation of the areas ofnp
2|Pk |) and the summation with respect tothe partitions (∑k=1s
4kM,E j |E j |) would
the majorant extinction coefficient (∑nj=1
differ according to the partitioning. So for a single subspace,
we also check for all the possible partitioning. Then, we estimate the average of the expected numbers of iterations without partitioning N̄nopart . If N̄ part < N̄nopart , we partition thespace E p accordingly, and continue the partitioning process
for the partitioned spaces recursively.

As we did not take into account the termination of iteration inside the analytical space in our estimation framework,
the prediction of the numbers of iterations in actual situations is not accurate. Nevertheless, we show that our estimation framework can predict the ratio of the performances
between any two schemes fairly well. Let us shortly explain
the reason. For any straight line in the analytical space, there
are a collection of random rays closely aligned with this line.
The random rays can be assumed to be uniformly distributedin the collection, and we can connect some of these rays togenerate a collection of rays that pass through the analytical space without termination of iteration inside the space.
The average number of iterations for these connected rays
can be assumed to be very close to that of the rays considered in our estimation framework. Thus, the average number of iterations for the random rays in the whole analyticalspace can be assumed to be the product of a constant andthe estimated average number of iterations using our framework. When computing the acceleration ratio between any
two schemes, such constants will be canceled out.

Although using this kd-tree based scheme will result in

The evaluation is conducted by both comparing the es-

5.3. kd-Tree based Scheme

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

timated acceleration ratios over Woodcock tracking (WT)
(i.e., the case without partitioning) and comparing the actual
acceleration ratios measured by shooting random rays. Here,the acceleration ratio is defined as the ratio between the averages of the expected numbers of iterations. We did not use
the computation time, as it is dependent on the optimizationof the code and the computing architecture.

Higher contrast

Y. Yue et al. / Toward Optimal Space Partitioning for Unbiased, Adaptive Free Path Sampling
8
7
6
5
4
3
2
1
1

For the evaluation, we prepared the following 168 types ofparticipating media, shown in Figure 4, with different spatial
distributions of the extinction coefficient. Using 3D Perlin
noise [Per02], we first prepared 3 types of base media (represented by 8 × 8 × 8 grid-cells) by setting the frequency
parameter from 0 to 2. After that, we normalized the extinction coefficient in the range [0, 1]. Then for each of thebase medium, we prepared 8 types of different variations bychanging the contrast to 2 j , ( j = 1, ..., 8). The contrast is defined as the ratio of the maximum difference of the extinction
coefficient to the average extinction coefficient. To obtain thedesired contrast, we applied a power r to the normalized extinction coefficients, where an appropriate r is found using abisection method. Finally, for each of these 24 types of participating media, we multiplied the extinction coefficient by
1/8, 1/4, 1/2, 1, 2, 4, 8 times of a baseline value to create 7
sets of participating media with different maximum extinction coefficients. In these 168 types of participating media,
a medium with higher contrast parameter, higher maximum
extinction coefficient and lower frequency parameter is more
inhomogeneous, because the portion with higher extinction
coefficient is coagulated in a smaller region.
First, we compare the estimated acceleration ratios in Figures 5(a) to (d). We can see that the tendency of the performance is basically KD ≥ hKD ≥ OC ≥ UG ≥ W T . Especially, for highly inhomogeneous participating media, kdtree based schemes (hKD and KD) performed much better.
To investigate this tendency in more detail, we show inTable 1(a) the win-loss standings between any two schemes,
and in Table 2(a) a quantitative evaluation of the estimatedacceleration ratios. From Table 1(a), we can see that although hKD is basically superior to UG and OC, there areabout 30 to 40% cases where hKD is inferior to UG and
OC. Conversely, KD is always equal or superior to UG andOC, and outperforms hKD for almost all the cases except
for only 7 cases. From Table 2(a), we can see that KD performs more than 2 times better than UG or OC for highly
inhomogeneous participating media.
Next, we show the actual acceleration ratios of theschemes over Woodcock tracking in Figures 5(e) to (h).
These ratios were obtained by shooting random rays andmeasuring the actual numbers of iterations. The computation times for the space partitioning are shown in Figure 6
and are not taken into account when computing the acceleration ratios. We can see that the tendency in the results
matches our estimation fairly well. The win-loss standings
and the qualitative evaluation of the acceleration ratio are
c 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

Higher extinction coefficient Higher extinction coefficient Higher extinction coefficient
Higher frequency parameter

Figure 4: Rendered results of the participating media used
for the evaluation of the sampling efficiency. The horizontal
axis corresponds to 21 sets of media with varying frequency
parameter and multiplier to the extinction coefficient. The
multipliers are 1/8 in the sets 1, 8 and 15; 1/4 in the sets 2,
9 and 16; 1/2 in the sets 3, 10 and 17; 1 in the sets 4, 11 and18; 2 in the sets 5, 12 and 19; 4 in the sets 6, 13 and 20; 8 in
the sets 7, 14 and 21. The frequency parameters are 0, 1 and
2 in the sets 1 to 7, 8 to 14 and 15 to 21, respectively. The
vertical axis shows the contrast parameter j.

Table 1: Win-loss standings for (a) the estimated acceleration ratios and (b) the measured acceleration ratios. For
example, in (a), UG is equal or superior to WT in 168 cases
and inferior in 0 case; hKD is equal or superior to OC in
102 cases and inferior in 66 cases.
(a) WTUG
OC
hKD
UG 168-0
OC 168-0 168-0
hKD 168-0 120-48 102-66
KD 168-0 168-0 168-0 161-7(b) WTUG
OC
hKD
UG 168-0
OC 168-0 165-3
hKD 168-0 124-44 102-66
KD 168-0 168-0 167-1 162-6

shown in Tables 1(b) and 2(b), respectively. By comparing toTables 1(a) and 2(a), we can see that our estimation framework works fairly good (e.g., the estimation error is within
10% for the average acceleration ratios). One may be curious
in the fact that hKD is always superior to WT in our result,while in [YIC∗ 10], it is shown that hKD becomes slightly
inferior to WT for nearly homogeneous participating media(media with low contrast values and high frequency parameters). This is due to the difference in the settings of the media. In this paper, the frequency parameters of the media arelower, which is due to the low resolution (83 ) of the participating media. We would like to conduct an evaluation forhigher resolution participating media in the future.
Finally, we show an evaluation using a volume data of asmoke (Figure 7). The smoke is lit by an environment lightsource. We calculated multiple scattering inside the smoke

Y. Yue et al. / Toward Optimal Space Partitioning for Unbiased, Adaptive Free Path Sampling

30
20
10

st
tra

t

t
as

t
as

0
19 21
13 15 177 9 111
er ext.
gh
Hi
1 3 5
t.
ex
t. Higher
q. param.
Higher ex
Higher fre

8

n
Co

ntr

ntr

0
19 21
13 15 177 9 111
er ext.
gh
Hi
1 3 5
t.
ex
t. Higher
q. param.
Higher ex
Higher fre(f) OC

50
40

Acceleration ratio

10

as
ntr
Co

Co

Co

(e) UG

30
20

8

Acceleration ratio

10
0
19 21
13 15 177 9 11
1
ext.
er
gh
1 3 5
Hi
t.
t. Higher ex
q. param.
Higher ex
Higher fre

8

t

30
20

40

Acceleration ratio

10
0
19 21
13 15 177 9 11
1
ext.
er
gh
1 3 5
Hi
t.
t. Higher ex
q. param.
Higher ex
Higher fre

40(d) KD50

Acceleration ratio

3020(c) hKD
50

Acceleration ratio

40

as
ntr
Co

t
as

t
as

(b) OC
50

70
60
50
40
30
20
10
0
8
19 21
17
15
13
7 9 111
ext.
1 3 5
t. Higher
t. Higher ex
q. param.
Higher ex
Higher fre

Acceleration ratio

70
60
50
40
30
20
10
0
8
19 21
17
15
13
7 9 111
ext.
1 3 5
t. Higher
t. Higher ex
q. param.
Higher ex
Higher fre

ntr

ntr

(a) UG

8

Co

Co

70
60
50
40
30
20
10
0
8
19 21
17
15
13
7 9 111
ext.
1 3 5
t. Higher
t. Higher ex
q. param.
Higher ex
Higher fre

Acceleration ratio

Acceleration ratio

t
as
ntr

Co

70
60
50
40
30
20
10
0
8
19 21
17
15
13
7 9 111
ext.
1 3 5
t. Higher
t. Higher ex
q. param.
Higher ex
Higher fre(g) hKD(h) KD

Figure 5: (a) to (d) and (e) to (h): Comparisons of the estimated and measured acceleration ratios, respectively. (a) and (e):Uniform grid based scheme, (b) and (f): octree based scheme, (c) and (g): kd-tree based scheme using the heuristic approach,and (d) and (h): kd-tree based scheme presented in this paper. The vertical axis shows the acceleration ratio over Woodcock
tracking. The left-right and anterior-posterior axes correspond to the horizontal and vertical axes in Figure 4, respectively.

1

Construction time (s)(c) hKD

0.4
8

t
as

13
7 9 11er
1 3 5
Set numb

1.2
0.8

ntr

0.000015 17 19 21

2.0
1.6

Co

t

(b) OC1

as

7 9 11 13er
1 3 5
Set numb

0.0004
8

ntr

t
as

(a) UG

10.000
15 17 19 21

0.00120.0008

Co

8

ntr

13 15
7 9 11er
1 3 5
Set numb

0.003

Co

t
as

ntr

Co

1

0.000017 19 21

0.009
0.006

0.00200.0016

Construction time (s)0.00048

0.012

Construction time (s)0.00120.00080.015

Construction time (s)0.00200.00160.0
15 17 19 21

7 9 11 13er
1 3 5
Set numb(d) KD

Figure 6: Computation time of the space partitioning. The vertical axis shows the construction time for the spatial acceleration
data structure in seconds. The left-right and anterior-posterior axes correspond to the horizontal and vertical axes in Figure 4,
respectively. When using UG, the construction times are almost constant. The construction times for hKD are shorter than UG
if the participating media are nearly homogeneous, and become longer otherwise. The construction times for KD are much
longer than for other schemes.

using Monte Carlo path tracing (1024 samples per pixel).
We set the single scattering albedo to 0.8, and we used theisotropic phase function. When using hKD, the rendering
times for Figures 7 (a) and (b) were 26.6 minutes and 44.9
minutes on a PC with an Intel Core i7 Extreme 975 CPU.
We used only 1 core of the CPU and the image resolution
is 480 by 960. In Figure 7 (b), using UG, OC and hKD resulted in 6.72, 8.42, 13.7 times faster rendering speeds thanusing WT, respectively. The resolution of the volume data
is 2563 . The construction times of UG, OC and hKD were
0.12s, 0.70s and 2.7s, respectively, and are negligible compared to the rendering times.
7. Conclusions and Future WorkIn this paper, we have presented an estimation frameworkthat works in the 3D space for evaluating the sampling efficiency in free path sampling, given a space partitioning. We

estimated the sampling efficiency by accounting for the average of the expected numbers of iterations needed before therays encounter ‘real’ scattering events. An important contribution of this paper is the analytical formulation for estimating the average of the expected numbers of iterations.
Additionally, the formulation does not contain any heuristic user specified parameters. Then, we have shown that theestimation framework can be used to construct new automatic space partitioning schemes. Moreover, by using ourestimation framework, we are able to estimate the differencein the sampling efficiencies among different space partitioning schemes.
For the future work, we would like to extend our framework to account for 1) the occurrence of ‘real’ scattering
events inside the subspaces, and for 2) the case where therays are not uniformly distributed in the space. We believethese cases could be handled by applying some weightc 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

Y. Yue et al. / Toward Optimal Space Partitioning for Unbiased, Adaptive Free Path Sampling

Table 2: Quantitative evaluation of (a) the estimated acceleration ratios and (b) the measured acceleration ratios. In
the notation α(β, γ), α, β and γ indicate the average, minimum and maximum of the acceleration ratios for the 168
cases. For example, in (a), hKD is 1.09 times superior to
OC in average. In the worst case, the performance of hKD
is 0.68 times to OC, and in the best case, the performance of
hKD is 2.22 times to OC.
(a)
UG
OC
hKD
KD

(b)
UG
OC
hKD
KD

WT
UG
OC
hKD
3.80(1.0, 19.7)
5.47(1.0, 47.5) 1.20(1.0, 2.41)
6.35(1.0, 62.3) 1.33(0.68, 3.24) 1.09(0.68, 2.22)
6.89(1.0, 64.8) 1.47(1.0, 3.32) 1.21(1.0, 2.22) 1.13(0.84, 1.66)
WT
UG
OC
hKD
4.08(1.0, 19.7)
5.60(1.0, 39.2) 1.18(0.99, 2.30)
6.09(1.0, 48.5) 1.23(0.58, 2.61) 1.02(0.57, 1.54)
6.64(1.0, 49.9) 1.36(1.0, 2.65) 1.15(0.97, 1.57) 1.14(0.89, 1.82)(C) (23500112) and National Science Council of Taiwan(NSC99-2221-E-002-016).
References[BB09] BADAL A., BADANO A.: Accelerating Monte Carlo simulations of photon transport in a voxelized geometry using a massively parallel graphics processing unit. Medical Physics 36, 11(2009), 4878–4880. 2[BM03] B ROWN F. B., M ARTIN W. R.: Direct sampling ofMonte Carlo flight paths in media with continuously varying
cross-sections. In Proc. ANS Mathematics & Computation Topical Meeting (2003). 2[CCW72] C ARTER L., C ASHWELL E., W.M.TAYLOR: Monte
Carlo sampling with continuously varying cross sections along
flight paths. Nuclear Science and Engineering 48 (1972), 403–
411. 2[Cha50] C HANDRASEKHAR S.: Radiative Transfer. Clarendon
Press, 1950. 1[Col68] C OLEMAN W.: Mathematical verification of a certain
Monte Carlo sampling technique and applications of the technique to radiation transport problems. Nuclear Science and Engineering 32 (1968), 76–81. 2[CPP∗ 05] C EREZO E., P ÉREZ F., P UEYO X., S ERÓN F. J., S IL LION F. X.: A survey on participating media rendering techniques. The Visual Computer 21, 5 (2005), 303–328. 2[GJJD09] G UTIERREZ D., J ENSEN H. W., JAROSZ W., D ON NER C.: Scattering. In SIGGRAPH ASIA Courses (2009). 2[JC98] J ENSEN H. W., C HRISTENSEN P. H.: Efficient simulationof light transport in scenes with participating media using photon
maps. In Proc. SIGGRAPH 98 (1998), pp. 311–320. 1, 2
[Lep07] L EPPÄNEN J.: Development of a New Monte Carlo Reactor Physics Code. PhD thesis, Helsinki University of Technology,
2007. 1, 2
[LW96] L AFORTUNE E. P., W ILLEMS Y. D.: Rendering participating media with bidirectional path tracing. In Proc. EGWR 96(1996), pp. 91–100. 1[Per02] P ERLIN K.: Improving noise. In Proc. SIGGRAPH 2002
(2002), pp. 681–682. 7

(a)

(b)

Figure 7: Rendered results of a smoke under different environment light sources. The smoke (b) has 4 times higher
extinction coefficient than the smoke (a).

ing factor on the estimated values. Another interesting research direction would be to construct other space partitioning schemes based on our framework. We are seeking for a
kd-tree based space partitioning scheme that utilizes someapproximations during the partitioning so that the computation time needed for the partitioning is kept practical while
the resulting sampling efficiency is still near optimal.
Acknowledgments
We gratefully acknowledge helpful comments and thoughtful suggestions from the anonymous reviewers. This work
was supported by JSPS Grant-in-Aid for Scientific Researchc 2011 The Author(s)c 2011 The Eurographics Association and Blackwell Publishing Ltd.

[PH89] P ERLIN K., H OFFERT E. M.: Hypertexture. Computer
Graphics (Proc. SIGGRAPH 89) 23, 3 (1989), 253–262. 2[PKK00] PAULY M., KOLLIG T., K ELLER A.: Metropolis lighttransport for participating media. In Proc. EGWR 2000 (2000),
pp. 11–22. 1, 2
[RSK06] R AAB M., S EIBERT D., K ELLER A.: Unbiased global
illumination with participating media. In Proc. Monte Carlo andQuasi-Monte Carlo Methods 2006 (2006), pp. 591–605. 1, 2, 3
[SKTM11] S ZIRMAY-K ALOS L., T OTH B., M AGDIC M.: Free
path sampling in high resolution inhomogeneous participating
media. Computer Graphics Forum 30, 1 (2011), 85–97. 2[WMHL65] W OODCOCK E., M URPHY T., H EMMINGS P.,
L ONGWORTH T.: Techniques used in the GEM code for Monte
Carlo neutronics calculations in reactors and other systems ofcomplex geometry. In Proc. Conference on the Application ofComputing Methods to Reactor Problems, ANL-7050 (1965),
pp. 557–579. 1, 2
[YIC∗ 10] Y UE Y., I WASAKI K., C HEN B.-Y., D OBASHI Y.,
N ISHITA T.: Unbiased, adaptive stochastic sampling for rendering inhomogeneous participating media. ACM Transactions
on Graphics (Proc. SIGGRAPH Asia 2010) 29, 5 (2010), 177:1–
177:7. 2, 3, 6, 7