Lattice-Based Volumetric Global Illumination

Feng Qiu, Fang Xu, Zhe Fan, Neophytou Neophytos,

Abstract—We describe a novel volumetric global illumination framework based on the Face-Centered Cubic (FCC) lattice. An FCC
lattice has important advantages over a Cartesian lattice. It has higher packing density in the frequency domain, which translates to
better sampling efficiency. Furthermore, it has the maximal possible kissing number (equivalent to the number of nearest neighbors of
each site), which provides optimal 3D angular discretization among all lattices. We employ a new two-pass (illumination and rendering)
global illumination scheme on an FCC lattice. This scheme exploits the angular discretization to greatly simplify the computation in
multiple scattering and to minimize illumination information storage. The GPU has been utilized to further accelerate the rendering
stage. We demonstrate our new framework with participating media and volume rendering with multiple scattering, where both are
significantly faster than traditional techniques with comparable quality.
Index Terms—Volume visualization, volume rendering, participating media, lattice, FCC lattice, sampling, multiple scattering, GPU.

✦
1 I NTRODUCTION transport in volumes, which includes multiple scattering. Blinn [2] has
Direct volume rendering has been an important technique in the graph- analytically solved the transport equation for constant density medium
ics and visualization community for many years. Usually the vol- with single scattering. Kajiya and von Herzen [17] have proposed
ume to be rendered is a function defined on 3D (or higher dimension) tracing rays in inhomogeneous volumes. To calculate multiple scatter-
space and discretized into a set of cells or points. Here, a regular ing, spherical harmonics have been used. Radiosity uses finite element
discretization has a single generating matrix, which is called a lat- methods to solve the transport equation. Rushmeier and Torrance [27]
tice [6]. Mathematically, a lattice in Rn is a discrete subgroup [9] of have exploited the zonal method for isotropic scattering. Max [20] has
Rn that can be generated from a vector basis by a linear combination extended the discrete ordinates method to capture anisotropic multiple
with integral coefficients. In other words, a lattice is the subgroup scattering. Hegeman et al. [14] have proposed a two-pass approach for
a1 v1 + a2 v2 + . . . + an vn , where v1 , v2 , . . . , vn is the vector basis, vi is strong forward scattering. Harris and Lastra [13] have used a similar
the generation vector, and ai are integers. The sites of a lattice are con- approach to render clouds. Kniss et al. [18] have introduced a volume
nected with a set of lines (or links). A site is defined as a neighbor of lighting model for GPU-accelerated volume rendering with forward
another site if they are connected by one lattice link. The lattice sites scattering using a single pass algorithm based on half angle slicing.
result from the discretization of space or space-time, and the lattice Riley et al. [26] have extended this method to render atmospheric phe-
links discretize the angular space. The most frequently used lattice for nomena. Jensen [16] has presented the photon mapping method for
volume visualization is the Cubic Cartesian (CC) lattice with gener- rendering participating media. Cerezo et al. [3] have provided a nice
ation vectors of X = (1, 0, 0), Y = (0, 1, 0), and Z = (0, 0, 1). Other survey on rendering participating media. Finally, Geist et al. [12] have
grids such as Body-Centered Cubic (BCC) and Face-Centered Cubic revised the Lattice-Boltzmann method to render participating media.
(FCC) lattices have also been investigated [22, 24, 29].
Direct volume rendering algorithms reconstruct a continuous func- In this paper, we describe a new volumetric global illumination
tion, which is projected to a 2D image. This procedure involves di- framework, which exploits both spatial and angular discretization on
mension reduction thus inevitably loses some information. To capture lattices. In computer graphics, spatial discretization has been well-
more details, many lighting and illumination methods have been de- studied to simplify the calculation of light-object interaction, but an-
veloped. Local illumination models omit sophisticated effects such as gular discretization has not been fully exploited. Specifically, we
multiple scattering and indirect illumination for the sake of rendering adopt the FCC lattice because it has better sampling efficiency com-
speed. However, they are the dominant light-object interaction for par- pared with the CC lattice and it provides optimal angular discretization
ticipating media (smoke, clouds) and many translucent materials. To among all lattices. We further describe a new two-pass algorithm to
cope with these effects, volumetric global illumination techniques are render participating media and volumes with multiple scattering. The
required in order to present important visual features [28]. idea of this algorithm is that the traced photons only move along the
Global illumination has not been widely employed in volume ren- lattice links. We call these photons “diffuse photons”. Here, the phase
dering because of the computation complexity. When a photon en- function is discretized on the lattice links to simplify the diffuse pho-
counters an object, it might be reflected, refracted, and scattered many ton tracing and radiance estimation. The storage of diffuse photons
times before finally being absorbed or exiting the scene. It is closely is minimized by storing the number of photons on lattice links. For
related to the radiative transfer problem that has been studied by physi- flexibility and extensibility, we also implemented tracing photons with
cists for decades [5]. The simulation of all kinds of interaction is accurate direction, which are called “specular photons” in this paper.
time-consuming, and many simplified models have been proposed in The O-Buffer data structure proposed by Qu et al. [25] has been ex-
computer graphics. Max [21] has evaluated several optical models for ploited to reduce the storage space of specular photons. Our volumet-
direct volume rendering and presented an integral equation for light ric global illumination framework is capable of producing high quality
images and is significantly faster than traditional methods. This gen-
eral and flexible framework can be extended to render hybrid scenes
• Feng Qiu, Fang Xu, Zhe Fan, Neophytou Neophytos, Arie Kaufman, and with both volumes and surface objects.
Klaus Mueller are with the Department of Computer Science, Stony Brook
University, Stony Brook, NY 11794-4400. E-mail: qfeng, fxu, fzhe,
nneophyt, ari, [email protected]. Our paper is structured as follows. The definition and mathematical
Manuscript received 31 March 2007; accepted 1 August 2007; posted online properties of FCC lattices are explained in Section 2. Section 3 ana-
27 October 2007. lyzes the sampling scheme on FCC lattices. Section 4 and 5 present
For information on obtaining reprints of this article, please send e-mail to: the rendering methods and implementation details, respectively. The
[email protected]. rendering results are shown in Section 6. Finally, we conclude the
paper and describe future work.
2 FCC L ATTICE
An FCC lattice consists of simple CC cubic cells with additional sam-
pling points located at the center of each cell face (see Fig. 1a). Ac-
cording to the lattice definition described in Section 1, any FCC lattice
site can be constructed via linear combination√of the three basis vectors
X = (1, 0, 0), Y = (0, 1, 0) and Z = (0.5, 0.5, 2/2), which can be ob-
tained by defining an appropriate rotated coordinate system (see Fig.
1b). This definition scheme is adopted in our paper due to its simplic-
ity, where the FCC lattice can be decomposed into two √interleaved sub
cubic lattices with a deviation vector of v = (0.5, 0.5, 2/2). This rep-
resentation provides a framework that enables quick indexing and ef- (a) (b)
ficient implementation of many basic lattice operations. For instance,
the mapping from an arbitrary FCC lattice site of index (i, j, k) to its Fig. 2. (a) The 12 neighbors (red spheres) of an FCC lattice site (green
CC coordinates (x, y, z) can be defined by the following equations: sphere) forms a cuboctahedron; (b) The voronoi cell of an FCC lattice
site (green sphere) is a rhombic dodecahedron.
x = i + (k mod 2)/2,
y = j + (k mod 2)/2,
√
z = 2k/2. (1) respectively. Thus, the FCC lattice presents a much more efficient
spacial sampling scheme over the traditional CC lattice.
Finding the nearest FCC lattice site to a sample point can be imple-
mented by first looking for two neighbors having the shortest Euclid-
ean distance to the sample point in the two sub cubic lattices using
Eqn. 1, and then selecting the closer one between them. 3 S AMPLING ON FCC

In any lattice used for volume visualization, the lattice sites are dis-
cretized or sampled from Rn , and an efficient, preferably optimal, sam-
pling scheme is paramount. An optimal lattice structure captures in-
formation in the hyper-volume Rn using the least number of sampling
points. Assuming an isotropic and band-limited sampling function, the
resulting frequency support of a sampling point is a hyper-sphere, sur-
rounded by a set of alias replicas. Hence, the most efficient sampling
scheme arranges the replicated (hyper-spherical) frequency response
y y as densely as possible in the frequency domain to avoid overlapping of
the aliased spectra. As demonstrated in multi-dimensional signal the-
z x z
x ory [8] an optimal sampling scheme is obtained when the frequency
response of the sampling lattice is an optimal sphere packing lattice
(a) (b) [6]. Optimal sampling lattices can achieve up to 13.4%, 29.3%, and
50% of savings in 2, 3 and 4 dimensions, and they have been used in
Fig. 1. An FCC lattice can be constructed by (a) adding (blue) sites at volume visualization [10, 22, 29] with high quality image results.
the centers of the faces of cubic cells; (b) interleaving 2D square grids
In three dimensions there are infinite optimal sphere packings in-
(red and green semi-transparent slices), where slice 2i+1 is shifted from
cluding the FCC lattice and the HCP (hexagonal closed packing). The
slice 2i by v. The second form (b) can be constructed by rotating the first
one (a) π4 about the z axis.
spatial equivalent of the FCC lattice in the frequency domain is the
BCC lattice, which is the inverse Fourier transform of the FCC and
vice-versa. The FCC lattice in the spatial domain corresponds to the
The geometric layout of the FCC lattice also gives rise to its higher BCC lattice in the frequency domain which is not an optimal sphere
angular discretization granularity than both CC and body-centered cu- packing, and the FCC lattice achieves about 23% of savings over the
bic (BCC) lattices, which is important for our rendering framework. CC lattice in terms of sampling efficiency. It was chosen for our global
Each site in the FCC lattice has direct links to a total of 12 nearest illumination framework because it is the lattice that maximizes uni-
neighbors, in contrast to 8 and 6 in the BCC and CC lattices, respec- form angular discretization with its kissing number of 12. The HCP
tively. This is the best angular discretization rate that any 3D regular is another candidate with an optimal kissing number of 12. In strict
lattice can achieve, since in R3 the maximum number of sphere of ra- mathematical definition, HCP is not a lattice but can be defined
dius 1 that can simultaneously touch the unit sphere (i.e., the kissing √ as the
union of the lattice L with generation vectors (1, 0,
√ p 0), (1/2, 3/2, 0)
number) is 12 [6]. This unique feature has important implications for p
and (0, 0, 8/3) and the translate L + (1/2, 1/ 12, 2/3). How-
sampling and interpolation, as we will discuss further below. For ex-
ever, it has a bias towards the z-direction. For any link with direc-
ample, when a particle is scattered at an FCC lattice site, it has 12
tion d = (dx , dy , dz ) such that dz 6= 0, the link of direction −d does
possible moving directions, which is 50% more than a BCC lattice
not exist. When a ray or photon moves from one of such links and is
and 100% more than a CC lattice. In addition, the 12 links of an FCC
not absorbed or scattered at the lattice site, it cannot continue straight
lattice site are symmetric under rotation and reflection, which supports
along its original direction. Therefore, HCP is not symmetric and thus
a relatively simple computation framework. Fig. 2a shows the cuboc-
unsuitable for use in our framework.
tahedron defined by the 12 closest neighbors of a FCC site, while Fig.
2b shows the Voronoi cell of an FCC lattice site, which is essentially a Given the initial assumption that the represented function is hyper-
rhombic dodecahedron. spherically band-limited, the ideal choice for the reconstruction func-
Finally, the reciprocal lattice of FCC yields its representation in the tion is also a radially symmetric kernel. We have studied a set
frequency domain, which is essentially a BCC lattice [29]. As it will of Gaussian reconstruction filters and have found that the relatively
2
be described in the next section, this property gives rise to its near- narrow Gaussian f (r) = e−2r , offers good frequency behavior and
optimal sampling behavior that is capable of reconstructing original reasonable overlap between neighboring sites. We have also tested
signals with a minimum number of samples. It has been shown in this experimentally by reconstructing a constant cubic volume (where
the literature [29] that the FCC sampling scheme requires 13.4% and f (x, y, z) = 1) sampled on an FCC lattice and measuring the mean
23% fewer samples in R2 and R3 domain compared to CC lattices, squared error for varying reconstruction kernel extents.
4 R ENDERING M ETHODS ray originated from site xi with direction ω intersects the cuboctahe-
There are three major processes for the interaction between light and dron at point v on one of its 14 faces. Denoting the vertices of the face
volumetric objects: emission, absorption and scattering. The differen- containing v as v0 , v1 , · · · , vm (m = 2, 3) in counterclockwise order and
tial equation [21] of light propagation in volumetric objects is letting ωk = vk − xi , the probability of the photon direction changing
to ωk is defined as the barycentric coordinates [11]:
dI(x, ω ) wk
= −σt (x, ω )I(x, ω ) + σe (x, ω ) + pk = , wk = ∏ A(v, vt , vt+1 ) (5)
ds Z ∑k wk t6=k−1,k
′ ′ ′
σs (x, ω ) f (x, ω , ω )I(x, ω )d ω (2)
4π
where A(v, vt , vt+1 ) is the area of triangle v, vt , vt+1 . The barycentric
where I(x, ω ) is the light intensity at position x in direction ω , σt is coordinates of v are well defined with Eqn. 5 because all the faces of
the extinction or attenuation coefficient (combination of absorption the cuboctahedron are regular polygons.
and out-scattering, σt = σa + σs ), σs is the scattering coefficient or When arriving at a lattice site xi from direction ω j the photon might
the scattering probability per unit distance. f (x, ω , ω ′ ) is the phase be absorbed, be scattered or continue moving along the extended link
function representing the conditional probability of a photon from di- of ω j . We use the absorption and the scattering coefficients, σa (xi , ω j )
rection ω to be scattered ′ and σs (xi , ω j ), to represent the probabilities of such events. (Please
R in direction ω assuming that the photon is
scattered, and it obeys 4π f (x, ω , ω ′ )d ω = 4π . In most naturally oc- note that this approximation is only correct when the link length l is
curring media, σt and σs are independent of direction ω . Also, most small compared to 1/σt . The actual absorption probability should be
exp(− 0l σa (xi , ω j ))ds. Because the links on an FCC lattice are of
R
media is isotropic and the phase functions only depend on the cosine
of the angle θ between ω and ω ′ . The solution of Eqn. 2 is: identical length, we use σa , σs and σt here for convenience.) The
Z D Z D  Russian roulette technique [23] is used to determine whether the pho-
ton is absorbed, scattered or transmitted. In detail, the program gener-
I(x, ω ) = R(x − sω , ω )exp σt (x − t ω )dt ds, (3)
0 0 ates a random number ξ ∈ [0, 1) and the photon:
Z
R(x, ω ) = σs (x, ω ) f (x, ω , ω )I(x, ω ′ )d ω ′
′
(4)
4π

is absorbed
 if ξ ∈ [0, σa (xi , ω j )),
where D is the distance from x to the view point and R(x, ω ) is the is scattered if ξ ∈ [σa (xi , ω j ), σa (xi , ω j ) + σs (xi , ω j )), (6)
radiance at position x from direction ω . if ξ ∈ [σa (xi , ω j ) + σs (xi , ω j ), 1).

is transmitted
We describe a new two-pass rendering algorithm on FCC lattices.
In the first pass, photons are emitted from light sources and the photon In contrast to the deterministic procedures of absorption and trans-
energy is distributed in the scene, illuminating the media. In the sec- mission, stochastic scattering requires further processing. In tradi-
ond pass, a ray tracing method is used to generate the final image. At tional methods, a phase function f (x, ω , ω ′ ) is utilized to describe
each sampling point x on the ray of direction −ω , the radiance R(x, ω ) the probability of a photon being scattered at location x with an in-
is estimated by the photons in a small neighborhood of x for shading. put direction ω and an output direction ω ′ . The computation of the
new direction is performed using importance sampling. In the widely
4.1 Diffuse Photon Tracing
used models such as the Schlick model [1], the probability depends on
For volumetric objects where the dominant effect is diffusion, we have cos(θ ) = ω · ω ′ only, thus the importance sampling returns the value
developed a new algorithm to trace photons on the lattice links. The of cos(θ ). In order to compute ω ′ , a local coordinate system at the
volumetric objects are sampled in a finite region of the FCC lattice. scattering position has to be constructed to convert spherical angles to
As shown in Fig. 3 on a 2D FCC lattice, when a photon is emitted direction vectors, which is computationally expensive.
from a light source, the nearest lattice site to the first hit point on the In our lattice illumination framework, the computation of the
lattice boundary is calculated and the moving direction of the photon is scattering process is greatly simplified because photons only move
discretized to one of the link directions. The photon will be traced on along discretized lattice links. Here, the continuous phase function
the links between lattice sites and its path is composed of line segments f (x, ω , ω ′ ) is discretized to f (i, j, k), which represents the probability
of lattice links. of a photon located at site xi being scattered from the input link ω j
to the new output link ωk . In practice, this discrete phase function is
constructed as a 2D table of n × n resolution to represent all possible
combinations of input/output links on a lattice site (n = 12 for FCC
lattices). The generation of the discrete phase function via sampling
and normalizing the continuous phase function is described by the fol-
lowing equation:

f (xi , ω j , ωk )
f (i, j, k) = n−1
. (7)
∑t=0 f (xi , ω j , ωt )
(a) (b) A more accurate method is to calculate the integral over the angular
space Ω defined in Eqn. 8, which can be solved numerically for arbi-
Fig. 3. Illustration of tracing photons on a 2D FCC lattice. (a) The photon trary continuous phase function:
path in green color in the real medium. (b) The photon path on a 2D FCC Z
lattice. The circles represent lattice sites and the black line segments f (i, j, k) = closest(ω , ωk ) f (xi , ω j , ω )d ω
represent lattice links. The red circle is the nearest lattice site to the
photon intersection point with the medium boundary. The photon path (Ω
is composed of the green lattice links. 1 if the closest link to d~ is ωk ,
closest(ω , ωk ) = (8)
0 otherwise.
A photon emitted from the light sources has an arbitrary direction
ω . We first convert the photon direction stochastically to one of the Importance sampling of discrete phase functions is simply a binary
search for a given random number ξ such that ∑t=0 k f (i, ω j , ωt ) ≤ ξ
lattice links. An FCC lattice site xi has 12 closest neighboring sites of
k+1
equal distance, which forms a cuboctahedron as shown in Fig. 2a. The and ∑t=0 f (i, ω j , ωt ) > ξ . The complexity of this binary search is
O(log(n)), which is very efficient because n is usually very small (n = The O-Buffer data structure [25] is used to store the photon infor-
12 for FCC lattices). mation compactly, where each lattice site stores a sequence of photons
After determining the photon behavior of each encounter, we record in its Voronoi cell (Fig. 2b). For each stored photon at position x, the
their activity information on the lattice sites. Basically, a stored photon nearest lattice site xi is computed and only the offset o from x to xi is
represents a possible light path from the light sources to its location. stored. The offset o is quantized into 256 levels in each axis so that
This information is used in the following rendering pass, where the ir- o can be compactly represented in 3 bytes. This 3-byte representation
radiance of sampling positions is estimated with the stored light paths increases the photon position accuracy by 256 times of the lattice res-
within a small neighborhood. The lattice-based framework enables us olution, while it only needs 25% of the storage space of the traditional
to store photons in a 3D array and the position of the photons is im- floating point representation. Because the search radius for the radi-
plicitly defined by their index in the array. Moreover, the quantized ance estimation is usually much larger than the link length, it is good
directions are encoded in a byte using the link index. Since the link enough for most rendering applications. The Voronoi cell (Fig. 2b) of
vectors are known a-priori, we adopt an optimized solution for photon an FCC lattice site is a rhombic dodecahedron. Assuming a unit dis-
direction storage, where photons with the same direction are grouped tance between neighboring sites, √the distance from the lattice site to
together. Here, a 1D array E(ω j ) of n elements is used for a lattice the vertices of its Voronoi cell is 22 . The maximum error introduced
site xi , such that E(ω j ) is the total energy of the photons at xi with √ √
direction ω j . Due to the employed Russian roulette technique [23], by this offset quantization scheme is 23 · 255
1
· 22 .
the photon energy does not change until it is absorbed. Therefore, the For photon direction encoding, vectors ω are converted to spherical
stored photons all have the same energy and only the integral num- coordinates and represented with 2 bytes [15]. Because the photon en-
ber of photons need to be stored at link (i, j). Given the maximum ergy does not change, only one byte is used to record the color channel
possible number of photons stored on the links, the unsigned byte or of the photon. In total, the storage space of one photon is just 6 bytes.
short format can be used to represent the actual photon counts instead In the rendering pass, the radiance at each sampling point is esti-
of storing individual floating point energy values. In other words, we mated with the photons stored in the lattice photon O-Buffer. For a
store all photons in a 4D integer array E(i, ω j ) with three dimensions spherical search neighborhood
√ S with radius r, the lattice sites in the
of site position and one dimension for link direction. sphere S′ of radius r + 22 are retrieved because the maximum distance
After the recording of diffuse photons, rays are traced from the im- √
age plane into the FCC lattice to collect irradiance values. At each from a lattice site to the photons stored in it is 22 . Then, the photons
sampling point x, the radiance is estimated by the photons inside a stored in these lattice sites are visited. If the distance to the sampling
small spherical region centered at x. With our 4D array of photon point x is larger than r, the photon is discarded. The radiance at x is
counts, the radiance in Eqn. 4 is calculated with the following simpli- summed over all the photons inside S with
fied formula:
R(x, ω ) = E ∑ σs (ω p , x) f (x, ω , ω p )g(x p − x) (10)
R(x, ω ) = ∑ ∑ σs (x, ω j ) f (x, ω , ω j )g(x′ − x)E(x′ , ω j )dx′ (9) p
X j
where x p and ω p are the photon position and direction, respectively.
where ω is the reverse ray direction, f is the continuous phase function The term E is the energy of the photon, which is the same for all
and X is the set of lattice sites in the search region. g is a normalized photons since the Russian Roulette scheme is used.
smoothing filtering function used for removing high-frequency noise.
Our diffuse photon tracing algorithm is capable of tracing multi-
Because the lattice sites are positioned regularly, searching for the lat-
million photons in seconds. With the FCC lattice, the photons move
tice sites in a sphere is simple and efficient. In our experiments, the
on the lattice links and are scattered only on the lattice sites. There-
medium is isotropic and the phase function only depends on the angle
fore, the most time consuming steps in traditional methods such as
between the ray direction and the lattice link ω j . The dot product of ω
sampling the lighting properties, calculating scattered directions with
and ω j can be pre-computed and reused for all the sampling points on
phase functions are greatly simplified. Its major disadvantage is the
one ray in the rendering process, which yields an efficient implemen-
ray effect, which cannot be neglected in certain cases, for example,
tation of the radiance estimation.
in reflection, refraction or scattering on specular surfaces, and hard
Our new algorithm greatly simplifies the computation of photon-
shadows. The specular photon tracing solves this problem but is much
volume interaction and photon storage. Therefore, our program can
slower.
trace millions of photons in a short time, which is good for improving
image quality by removing the stochastic noise and variance without
excessive smoothing in the radiance estimation. Moreover, our method 5 I MPLEMENTATION
is general and can use arbitrary phase functions including those of
strong backward scattering. We have implemented the algorithms described in Section 4 to render
participating media and volumetric datasets. Since the reconstruction
4.2 Specular Photon Tracing process in the current scanning modalities such as MRI and CT only
The method described in Section 4.1 is capable of calculating multiple produce rectilinear data, the FCC data we used are generated by re-
scattering events for participating media and volumes where diffusion sampling existing rectilinear volumes or voxelizing geometric objects.
is the dominant effect. However, specular reflection and refraction In the resampling process, a windowed sinc filter has been used to
may exhibit ray effects when discretized rays hit smooth specular sur- obtain optimal signal reconstruction quality.
faces. This effect has been discovered by the radiative heat transfer We use the incremental triangle voxelization method [7] to vox-
community [4] and found to be caused by the discretization of scatter- elize polygonal models to the FCC lattice. The original method was
ing directions when accurate directions are needed for specularity. proposed for volumes of CC lattices, but the employed distance-based
To mitigate this ray effect, we describe an enhanced algorithm method enables its direct application to the FCC lattice.
called specular photon tracing, where every photon is associated with For surface voxelization, the volume density of a lattice site p is
its accurate direction ω and start position x. In each time step, the pho- determined by the distance between p and its closest triangle. Each
ton moves on the FCC lattice and the new sampling position is calcu- triangle is processed in the following manner. For each lattice site p
lated by x = d × ω where d is the step size. The lighting properties σa in the neighborhood of the triangle, the distance d between p and the
and σs are sampled to decide whether the photon is absorbed, scattered triangle is computed. The distance d is positive if p is in the normal
or transmitted at x (Eqn. 6). If the photon is scattered, the new direc- direction of the triangle, or negative if p is in the reverse direction. If
tion ω is computed with the continuous phase function. The Russian |d| is smaller than the previously stored absolute value of the distance,
roulette technique is again used to avoid photon energy change. |d| replaces the previous stored value and the density of p is updated
with the following equation: 6 R ESULTS
 We have implemented our algorithms on a 3.4GHz PC with 3GB mem-
 1 if d < −W, ory and a Geforce 8800 GTX graphics card. All the resulting images
density = 0 if d > W, (11) are of 512 × 512 resolution and cropped in Figs. 4-8.
 0.5 × (1 − d ) otherwise,
W Fig. 4 demonstrates the rendering results of participating media
with our algorithms. We use a single light of white color placed on top
where W is the width of the oriented box filter. √ While for a CC lat- of the smoke dataset. Table 1 gives the times of different algorithms
tice the optimized width was estimated to be 2 3 voxel units [7], we used to render the corresponding images.
estimate that for an FCC lattice, the optimized width of the filter is In Fig. 4b, the eccentricity coefficient k of the Schlick phase
decreased to 2 lattice units. Based on the surface voxelization result, function is set to 0.2. The absorption and scattering coefficients are
the interior of the solid is voxelized using seed growing. σa = 0.08 and σs = 0.2, respectively. Fig. 4c uses the same coef-
For inhomogeneous participating media, such as clouds and smoke, ficients except that the eccentricity k is −0.5 for strongly backward
the absorption coefficient σa (x) and scattering coefficient σs (x) typi- scattering. The photon tracing of 1 million photons and subsequent
cally do not rely on the light direction, although our algorithm is ca- rendering passes cost about 11 and 48 seconds on the CPU, respec-
pable of rendering anisotropic media. The volume data of the medium tively. With GPU acceleration, the time of the rendering stage is re-
usually defines the density field ρ (x) of particles. We assume that the duced to 2.4s. This performance is significantly faster than the orig-
σa (x) and σs (x) are proportional to the local density ρ (x). For FCC inal photon mapping method [16], where the tracing 10K photons in
lattices where the lattice links have unit lengths, it implies: a cloud model of similar size takes 8 seconds on an HP computer of
16 180MHz PA-8000 processors, while rendering a 1024 pixel wide
σa (x, ω j ) = σa ρ (x), and σs (x, ω j ) = σs ρ (x). (12) image takes 92 seconds. Note that although a higher resolution is used
in Jensen and Christensen’s method [16], the complexity of the al-
gorithms is mainly determined by the number of photons, of which
where σa and σs are user-defined scaling coefficients. In real world our example generates 100 times more. A major further enhancement
phenomena, most practical participating media are isotropic and the of our implementation can be achieved by incorporating empty space
phase function f (x, ω , ω ′ ) does not vary upon position x. In our imple- skipping or adaptive sampling techniques that are used by Jensen and
mentation, f (x, ω , ω ′ ) = f (ω · ω ′ ) is described with the Schlick model Christensen [16].
[1]. The participating medium is represented with an FCC lattice of Fig. 4d has been rendered using specular photon tracing with the
densities and all the coefficients are calculated by scaling ρ (x). same medium properties as Fig. 4b, and the rendering time is sim-
Our algorithms support chromatic participating media, where the ilar to traditional photon mapping methods [16]. With our compact
coefficients σa (x, λ ) and σs (x, λ ) are wavelength-dependent. We im- FCC O-Buffer data structure, 5.8 million stored photons only consume
plement it by defining scaling factors of absorption and scattering in 35MB of memory space. The search radius for radiance estimation is
RGB channels. The photons emitted from the light sources can be changed from 2.0 to 3.0 to remove the stochastic variance. Given the
red, green or blue randomly. The photon tracing program calculates same number of photons, diffuse photon tracing is approximately 21
σa (x, λ ) and σs (x, λ ) with proper scaling factors in the color channel times faster than specular photon tracing and the corresponding ren-
of the traced photons. In the rendering pass, the opacity value α of a dering pass is 15 times faster. We observed that the image produced
sampling point is the average extinction of those in three channels: from our lattice-based framework using diffuse photon tracing (Fig.
4b) is comparable to the image using specular photon tracing (Fig.
1 4d), and has a similar quality and appearance of those computed with
α = (σt (x, λr )d + σt (x, λg )d + σt (x, λb )d ) (13)
3 traditional photon mapping methods (such as Figure 12.10 of [23]).
However, our framework has much better performance.
where d is the step size. In Fig. 4e, the absorption coefficient has been made wavelength-
For general volume datasets, transfer functions have been exploited dependent and the values in RGB channels are σa = (0.08, 0.15, 0.3),
to define lighting properties from the density field. Some 2D transfer while the scattering coefficient is the same as in Fig. 4b. The time of
functions might also use the gradient information. Our framework is the rendering pass increases to 52.7 seconds mainly because the radi-
general and capable of integrating any transfer function as long as the ance estimation is performed in 3 channels. In Fig. 4f, the scattering
input data (density, gradient or any other data) of the transfer function coefficient in the blue channel has been changed to 0.4 and the eccen-
is defined on the lattice. In our current implementation, a 1D transfer tricity coefficient k has been changed to 0.5. In Figs. 4e and 4f, the
function is defined for σa and σs in each RGB channel. Eqn. 13 can upper part of the smoke is grey but the lower part under the shadow of
be used to compute the opacity α or a separate transfer function can the upper part becomes orange because of the different absorption and
be defined for α . scattering coefficients in the RGB channels.
Following the photon tracing computation, a single-pass ray-casting Fig. 5 demonstrates another example of participating media. The
approach is employed on the GPU to render the diffuse photon tracing eccentricity coefficient k of the Schlick phase function is 0.2. The
results. We used an algorithm similar to [19], except that the sampled absorption and scattering coefficients are σa = 0.05 and σs = 0.1, re-
volume density is used for obtaining the scattered coefficient as well spectively. One million photons have been traced in 12.7 seconds and
as transparency values to composite the final radiance results. Here, in the rendering pass amounted to 98.5 seconds on the CPU and 5.5 sec-
addition to the density volume, an additional photon storage volume onds on the GPU.
that records the diffuse photon distribution on each lattice site is sam- Fig. 6 displays the foot of the visible human CT data. Fig. 6a
pled to composite ray values. Both the density and the diffuse photon is rendered using a ray casting method with local illumination. Fig .
volumes are FCC lattices, stored and indexed as described in Section 6b has been rendered using our diffuse photon tracing algorithm. The
2. The diffuse photon storage table is essentially a 12-element array, bone appears semi-translucent and brighter and has soft self-shadows
each of which records the number of photons stored along 12 different due to multiple scattering. In Fig. 6c, the muscle and soft tissue are
lattice links. A two-byte unsigned short is allocated for each such node displayed with red color with the absorption coefficient similar to 6a.
to provide a count of up to 65536. To compute the radiance estimation Fig. 7 and Fig. 8 are the rendering results of the engine and lobster
at each sampling point, dot-products of each lattice link with the cur- data, respectively. Fig. 7a has been rendered using ray casting with
rent viewing ray are used to weight individual diffuse radiance value, local illumination and the other two images have been rendered with
which is given by indexing the photon distribution previously com- our framework and indeed the objects appear substantially different.
puted on the lattice. The 12 weighted values are then summed up to In Fig. 7b, the material absorbs green and blue photons faster and
yield the final contribution. Filtering on all lattices uses the Gaussian becomes more red gradually through the light direction. In Fig. 7c,
kernel of size 2 and the sum is normalized at the end. the high density region appears saturated red for emphasis and the sur-
 (a) (b) (c)

(d) (e) (f)

Fig. 4. Inhomogeneous smoke rendered with global illumination (multiple scattering) and an anisotropic phase function. The original data is
100 × 100 × 40 and we sampled it into a 108 × 100 × 56 FCC lattice with a windowed sinc filter. The algorithms used are (a) ray casting; (b) our diffuse
photon tracing; (c) our diffuse photon tracing with strong backward scattering; (d) our specular photon tracing; (e) and (f) our multi-channel diffuse
photon tracing.

Table 1. Times used to render the smoke in Fig. 4. From left to right, the columns represent the algorithm used, the number of photon traced (in
millions), photon tracing time (in seconds), and rendering time (in seconds) on the CPU and the GPU.
Algorithm Photon count (M) Photon tracing (s) CPU rendering (s) GPU rendering (s)
Fig. 4a ray casting N/A N/A 20.6 1.0
Fig. 4b diffuse photon tracing 1.0 11.3 48.1 2.4
Fig. 4c diffuse photon tracing 1.0 11.7 48.1 2.4
Fig. 4d specular photon tracing 0.1 27.0 729.0 N/A
Fig. 4e multi-channel diffuse photon tracing 3.0 48.5 52.7 2.6
Fig. 4f multi-channel diffuse photon tracing 3.0 48.6 52.7 2.6

7 D ISCUSSION
In our framework, a photon is saved at every step of the first pass,
regardless of whether it is scattered or not. A stored photon repre-
sents a possible path from the light sources. In the ray tracing pass,
the radiance estimation actually calculates the density of photon paths
at the sampling positions. In contrast, the traditional photon mapping
method uses the probabilistic sampling technique to calculate the step
size, and the expected step size is 1/σt . Usually σt is small and trac-
ing a photon can generate many more stored photons in our algorithms
than in photon mapping. Moreover, the FCC lattice provides a more
efficient data structure to store photons. In diffuse photon tracing, each
lattice site stores multiple photons. In radiance estimation, the contri-
bution of multiple photons on a lattice link is computed jointly using
Eqn. 9. Suppose the radiance estimation searches photons in the spher-
ical neighborhood S, and there are k0 sites and k1 photons inside S. It
Fig. 5. Cloud rendered with our diffuse photon tracing. The resolution costs O(k0 ) time for diffuse photon tracing and O(k0 + k1 ) for specu-
of the data is 96 × 74 × 143. lar photon tracing. However, it costs O(k1 + log n) time with the k-d
tree data structure in photon mapping, where n is the total number of
photons. Also, the k-d tree need to be built before the rendering pass,
which costs O(n log n) time. However, our algorithms do not need
such a preprocessing step.
rounding region is less saturated red due to color bleeding. The mater- There are other simplified lighting models for participating media
ial of the objects is easily controlled and changed using user specified and volumes [13, 14, 18, 26]. For all these methods, only forward
transfer functions. In Fig. 8b, the lobster shell absorbs green and blue scattering is considered and lighting values are propagated from slice
photons faster than red ones and scatters red photons more than green to slice. The value of each pixel is calculated by sampling and at-
and blue ones. The shell appears red and the muscle casts soft shadows tenuating neighboring pixels (up to 4) on the previous slice. In other
onto itself. Fig. 8c displays the data from a different camera position. words, forward scattering is calculated in several directions in the 2π
 (a) (b) (c)

Fig. 6. Global illumination of a CT scan of the visible human foot. The original data is 1283 and the sampled FCC lattice is 128 × 128 × 180. The
images are rendered with (a) ray casting with local illumination; (b) our diffuse photon tracing; (c) our multi-channel diffuse photon tracing.

(a) (b) (c)

Fig. 7. Global illumination of an industrial CT scan of an engine. The original data is 128 × 128 × 64 and the sampled FCC lattice is 136 × 136 × 98.
The images are rendered with (a) ray casting with local illumination, (b) and (c) our multi-channel diffuse photon tracing.

(a) (b) (c)

Fig. 8. Global illumination of a CT scan of a lobster. The original data is 128 × 127 × 28 and the sampled FCC lattice is 114 × 113 × 35. The images
are rendered with (a) ray casting with local illumination; (b) and (c) our multi-channel diffuse photon tracing.

Table 2. Rendering of the foot, engine and lobster data in Figs. 6, 7 and 8. The meaning of the columns are the same as in Table 1.
Algorithm Photon count (M) Photon tracing (s) CPU rendering (s) GPU rendering (s)
Fig. 6b diffuse photon tracing 1.0 27.0 98.0 4.7
Fig. 6c multi-channel diffuse photon tracing 3.0 86.1 118.6 5.8
Fig. 7b multi-channel diffuse photon tracing 9.0 136.4 110.6 5.4
Fig. 7c multi-channel diffuse photon tracing 9.0 147.9 110.7 5.4
Fig. 8b multi-channel diffuse photon tracing 3.0 26.9 40.7 2.0
Fig. 8c multi-channel diffuse photon tracing 6.0 53.9 37.9 1.9
solid angle. Our algorithms can handle scattering within the entire 4π [3] E. Cerezo, F. Pérez, X. Pueyo, F. Seron, and F. Sillion. A survey on par-
solid angle. Also, our methods can store photons from multiple light ticipating media rendering techniques. The Visual Computer, 21(5):303–
sources in one volume, while previous methods do not support multi- 328, June 2005.
ple light sources. The method of Max [20, 21] is more accurate than [4] J. C. Chai, H. S. Lee, and S. V. Patankar. Ray effect and false scattering in
ours but runs slowly. the discrete ordinates method. Numerical Heat Transfer Part B, 24:373–
The idea of tracing photons on the lattice links is general and might 389, 1993.
be applied to other lattices such as the CC lattice. However, a CC lat- [5] S. Chandrasekhar. Radiative Transfer. Dover Publications, 1960.
[6] J. H. Conway, N. J. A. Sloane, and E. Bannai. Sphere-packings, Lattices,
tice site only has 6 nearest neighbors, which is definitely not enough
and Groups. Springer-Verlag, New York, NY, USA, 1987.
for discretizing the phase function. Consider a strongly forward scat-
[7] F. Dachille and A. Kaufman. Incremental triangle voxelization. Graphics
tering phase function, a photon arriving from a link can only move Interface, pages 205–212, May 2000.
forward on its original direction or be scattered to 4 directions perpen- [8] D. E. Dudgeon and R. M. Mersereau. Multidimensional Digital Signal
dicular to its incoming direction. Although we can add links between Processing. Prentice Hall Professional Technical Reference, 1990.
the secondary or tertiary neighbors, this solution needs to calculate [9] D. S. Dummit and R. M. Foote. Abstract Algebra. John Wiley and Sons,
and store the absorption
√ √and scattering probabilities of links with vary- second edition, 1999.
ing length (1, 2 and 3), thus making the photon tracing algorithm [10] A. Entezari, R. Dyer, and T. Moller. Linear and cubic box splines for the
more complicated and time consuming. Instead, a FCC lattice site has body centered cubic lattice. Proceedings of IEEE Visualization, pages
12 nearest neighbors, which is sufficient in photon tracing (in previous 11–18, 2004.
systems, only 4 or 5 directions are used for forward scattering). The [11] M. S. Floater, K. Hormann, and G. Kós. A general construction of
link length in the FCC construction is uniform, which greatly simpli- barycentric coordinates over convex polygons. Advances in Computa-
fies the computation and storage: the absorption and scattering proba- tional Mathematics, 24(1–4):311–331, Jan. 2006.
bilities on 12 links of a site are the same and stored only once on the [12] R. Geist, K. Rasche, J. Westall, and R. J. Schalkoff. Lattice-boltzmann
site. Moreover, the FCC lattice has better sampling efficiency than the lighting. Proceedings of the Eurographics Workshop on Rendering Tech-
CC lattice. With the same number of sites, FCC captures 23% more niques, pages 355–362, 2004.
[13] M. J. Harris and A. Lastra. Real-time cloud rendering. Computer Graph-
information. In addition,
√ the maximum distance
√ of an arbitrary point
ics Forum, 20(3):76–84, 2001.
to its nearest site is 22 in FCC instead of 23 in CC, which means [14] K. Hegeman, M. Ashikhmin, and S. Premože. A lighting model for gen-
18.4% less quantization error of photon positions with the O-Buffer eral participating media. Proceedings of Symposium on Interactive 3D
data structure. Graphics and Games, pages 117–124, 2005.
[15] H. W. Jensen. Global illumination using photon maps. Proceedings of
8 C ONCLUSIONS AND O NGOING W ORK Eurographics Workshop on Rendering Techniques, pages 21–30, 1996.
[16] H. W. Jensen and P. H. Christensen. Efficient simulation of light transport
We have presented a novel framework for volumetric global illumina- in scenes with participating media using photon maps. SIGGRAPH, pages
tion based on FCC lattices. Benefitting from the FCC lattice unique 311–320, 1998.
geometric and sampling properties, we were able to develop algo- [17] J. T. Kajiya and B. P. V. Herzen. Ray tracing volume densities. SIG-
rithms that can render participating media and volumes with multi- GRAPH, pages 165–174, 1984.
ple scattering effects. Our diffuse photon tracing algorithm renders [18] J. Kniss, S. Premoze, C. Hansen, P. Shirley, and A. McPherson. A model
high quality images at a speed significantly faster than conventional for volume lighting and modeling. IEEE Transactions on Visualization
methods. We have also described a memory efficient specular photon and Computer Graphics, 9(2):150–162, April-June 2003.
tracing algorithm. In future work, our plan is to extend our framework [19] J. Kruger and R. Westermann. Acceleration techniques for GPU-based
to render hybrid scenes of volumetric datasets and surface objects with volume rendering. Proceedings of IEEE Visualization, pages 287–292,
specular reflection and refraction. 2003.
[20] N. L. Max. Efficient light propagation for multiple anisotropic volume
Our current implementation only uses 12 links to the nearest sites on
scattering. Eurographics Workshop on Rendering, pages 87–104, 1994.
the FCC lattice. For future efforts, we consider incorporating the links
[21] N. L. Max. Optical models for direct volume rendering. IEEE Transac-
connecting the secondary and tertiary neighbors to increase the angu- tions on Visualization and Computer Graphics, 1(2):99–108, 1995.
lar discretization granularity. With these additional links, the phase [22] N. Neophytou and K. Mueller. Space-time points: 4D splatting on ef-
function can be discretized with even higher precision, thus more ac- ficient grids. Proceedings of Symposium on Volume Visualization and
curate radiance estimation will be √ obtained. In doing so, a total of 42 Graphics, pages 97–106, 2002.
neighbors with distance less than 3 will have to be considered and [23] M. Pharr and G. Humphreys. Physically Based Rendering: From Theory
different link lengths will participate in the computation of the absorp- to Implementation. Morgan Kaufmann Publishers Inc., San Francisco,
tion and scattering coefficients. Hence, a more efficient data structure CA, USA, 2004.
will be required for diffuse photon storage. We also consider applying [24] W. Qiao, D. S. Ebert, A. Entezari, M. Korkusinski, and G. Klimeck.
our general framework to other types of particles or waves, such as VolQD: Direct volume rendering of multi-million atom quantum dot sim-
phonons and Huygens’ wavelets, to solve their energy and wave prop- ulations. Proceedings of IEEE Visualization, pages 319–326, 2005.
agation problems. Last, we would like to accelerate the photon tracing [25] H. Qu and A. Kaufman. O-buffer: A framework for sample-based
pass on the GPU to improve the overall framework performance. graphics. IEEE Transactions on Visualization and Computer Graphics,
10(4):410–421, July-August 2004.
[26] K. Riley, D. S. Ebert, M. Kraus, J. Tessendorf, and C. D. Hansen. Effi-
ACKNOWLEDGEMENTS cient rendering of atmospheric phenomena. Proceedings of Eurographics
This work was supported, in part, by NSF CCF-0702699, NSF CA- Symposium on Rendering, pages 374–386, 2004.
REER grant ACI-0093157, NIH grant CA082402, and NIH grant [27] H. E. Rushmeier and K. E. Torrance. The zonal method for calculating
5R21EB004099-02. The smoke data is courtesy of Dunc Nguyen and light intensities in the presence of a participating medium. SIGGRAPH,
Ron Fedkiw. The engine data is courtesy of GE. The visible human pages 293–302, 1987.
data is courtesy of National Library of Medicine. [28] L. M. Sobierajski and A. E. Kaufman. Volumetric ray tracing. Symposium
on Volume visualization, pages 11–18, 1994.
[29] T. Theußl, T. Möller, and M. E. Gröller. Optimal regular volume sam-
R EFERENCES pling. Proceedings of IEEE Visualization, pages 91–98, 2001.
[1] P. Blasi, B. L. Saëc, and C. Schlick. A rendering algorithm for dis-
crete volume density objects. Computer Graphics Forum, 12(3):201–210,
1993.
[2] J. F. Blinn. Light reflection functions for simulation of clouds and dusty
surfaces. SIGGRAPH, pages 21–29, 1982.