Diffusion and Multiple Anisotropic Scattering for
Global Illumination in Clouds
Nelson Max

Greg Schussman

Lawrence Livermore National Lab. Stanford Linear Accelerator Ctr.
Mail Stop L-560
Mail Stop 26
7000 East Avenue
2575 Sand Hill Road
USA 94550 Livermore CA
USA 94025 Menlo Park CA

max2@llnl.gov

schussman@SLAC.stanford.edu

Ryo Miyazaki, Kei Iwasaki,
and Tomoyuki Nishita
University of Tokyo
7-3-1 Hongo, Bunkyo-ku
Japan 113-0033 Tokyo{ryomiya, kei-i, nis}@nislab.is.s.u-tokyo.ac.jp

ABSTRACTThe diffusion method is a good approximation inside the dense core of a cloud, but not at the more tenuous
boundary regions. Also, it breaks down in regions where the density of scattering droplets is zero. We haveenhanced it by using hardware cell projection volume rendering at cloud border voxels to account for the straightline light transport across these empty regions. We have also used this hardware volume rendering at key voxels
in the low-density boundary regions to account for the multiple anisotropic scattering of the environment.

Keywords
Diffusion approximation, multiple anisotropic scattering, global illumination, participating media, clouds.

1. INTRODUCTION
The appearance of a cloud is produced by multiple
scattering of the incident illumination. The phasefunction for the result of two scattering events is the
convolution of their two phase functions, and theconvolution of any non-delta-function phase function
with itself enough times approaches a uniform
distribution, as shown in [Stam95]. In the interior of acloud, the water droplets are dense, and the mean freepath of a photon is short, so the large number ofscattering events in a voxel will cause it to radiate
scattered light diffusely. In the diffusion
approximation, explained briefly in section 2 below,
the equation for radiance transport is approximatedby an equation for diffuse light transport, which canbe more easily solved by multigrid relaxation. Kajiyaand von Herzen [Kaj84] first proposed it to thecomputer graphics community, and Stam [Stam95]
Permission to make digital or hard copies of all or part ofthis work for personal or classroom use is granted without
fee provided that copies are not made or distributed forprofit or commercial advantage and that copies bear thisnotice and the full citation on the first page. To copy
otherwise, or republish, to post on servers or to redistributeto lists, requires prior specific permission and/or a fee.
Journal of WSCG, Vol.12, No.1-3, ISSN 1213-6972
WSCG’2004, February 2-6, 2004, Plzen, Czech Republic.
Copyright UNION Agency – Science Press

gave the first practical implementation for radiancetransport in clouds. The diffusion calculations in[Stam95] assume a rectangular solid of non-zerocloud droplet density, which is unrealistic for actual
clouds. Furthermore, Mie scattering from small water
droplets is concentrated in the forward direction, andthis directional scattering is evident in the subtle
color and intensity effects at the edges of a cloud,
where the expected number of scattering events is notlarge enough to produce anisotropy. The goal of ourwork is photorealistic rendering of clouds using the
diffusion approximation on the dense core, withmodifications to account for regions of zero droplet
density, and multiple anisotropic scattering at theedges.
Older work on global illumination in participating
media is well summarized in [Per97] and thereferences cited there. The basic Monte Carlo methodof tracing random rays from the light source or the
viewpoint suffers from the problem that few paths
from the viewpoint reach the light source, and viceversa. By tracing photons from the light source,
accumulating scattering events into a voxel data
structure called a photon map, and then doing a final
gather from the viewpoint by ray tracing the photon
map, Wann Jensen and Christensen [Jen98] solve thisproblem. However, in realistically large and dense
clouds, the photons will basically make a random

walk through the dense medium, and progress astraight line distance O( n ) in n bounces, total path
length n, or computation time n. Thus illumination on
one side of a dense cloud will take a very long time toreach the other side. Recently, Premoze, Ashikmin,
and Shirley [Prem03] have applied the path integral
method from quantum mechanics to make a more
efficient approximation, but it can be biased because
it takes into account mainly the forward scattering,and does not account for all light paths.
Jensen et al. [Jens01] proposed an approximate
solution to the diffusion equations, using a dipole oftwo point sources, one above and one below a planar
object surface, and a simple analytic solution for the
diffusion equation for point sources in a
homogeneous medium. Jensen and Buehler [Jens02]
applied this to objects represented by surface point
clouds, using an octree hierarchy to cluster the effectsof irradiance at points far from the position being
rendered. Mertens et al. [Mer03] used a similar
hierarchical method on polygonal surfaces, with asemi-analytic integration for the effect of theirradiance on a triangle. Lensch et al. [Len02] used
triangle-based piecewise linear basis functions for
global transport, and texel-based piecewise constantbasis functions for local transport. Dachsbacher andStamminger [Dach03] combine this local/global
transport separation with the hierarchical evaluation
into a texture-based approximation, which can beevaluated in real time using programmable vertex andfragment shaders. The dipole approximation used in
all these methods assumes that the surface is planar,
and that the medium is uniform and optically thick.
These methods give plausible realistic renderings formore complex geometry, but cannot deal with thenon-diffuse scattering discussed above at the tenuousedges of clouds, where the mean free path of a photon
is comparable to the geometry feature size, and the
diffusion approximation is invalid.

2. DIFFUSION APPROXIMATION
The following is a brief summary of the derivation of
the diffusion approximation. For details, see [Ish78].
(Note that the derivation in [Stam95] is incorrect, andStam’s equations differ from Ishimaru’s by factors of
3 in several terms.)
The equation for light propagation in a participating
medium isdI ( X , ω )= −κ ( X ) I ( X , ω ) +ds

κ ( X ) a ⋅ p(ω , ω ′) I ( X , ω ′) dω ′,4π(1)where X is 3D position, ω and ω’ are ray directions
on the unit sphere, I(X,ω) is radiance at X flowing in
direction ω, s is length along the ray in direction ω,
κ(X) is the extinction coefficient, a is the albedo, thescattered fraction of the extinction, and p(ω,ω’) is thephase function giving the probability density thatscattered light from ω’ scatters to ω. (Note that in[Ish78], p(ω,ω’) also includes the albedo a.)
In a cloud, κ(X) = ρ(X) σt, where ρ(X) is the number
density of water droplet particles, and σt is theextinction “cross section” per particle. Also let σs = a
σt be the scattering cross section per particle, and σa= σt – σs be the absorption cross section per particle.
Assume that the phase function p(ω,ω’) depends only
on the angle θ between the unit vectors ω and ω’,and letp (ω , ω ′)(ω ⋅ ω ′) dω

p1 =
4π

be the average cosine of the scattering angle θ. For
the Henyey-Greenstein scattering phase functionp (θ ) =1− g2
34π (1 − 2 g cos θ + g 2 ) 2which we use, p1 = g.
Let Iri(X, ω) be the reduced incident radiance, which
can be computed from the sun and sky radianceIsky(ω) by the attenuation formula∞

I ri ( X , ω ) = I sky (ω ) exp − σ t ρ ( X − sω )ds .

(2)

0

Let
I d ( X , ω ) = I ( X , ω ) − I ri ( X , ω )be the “diffuse” radiance that has been scattered atleast once. Approximate Id(X, ω) with the first orderdependence on the unit direction vector ω :Id (X,ω) = I0(X) + I1(X) ω .
Substituting
I (X,ω) = Iri(X, ω) + I0(X) + I1(X) ω(3)

into equation (1), and integrating ω over the unit
sphere, we get a scalar equation. Then multiplying by
the unit vector ω and integrating over the unit sphere,
we get a vector equation. We solve this vector
equation for I1(X) in terms of I0(X) and its
derivatives, and substitute the result into the scalar
equation, which also involves the divergence of I1(X).
By doing some vector algebra, vector calculus, andmathematical manipulations on these equations, we
end up with the single “diffusion” equation

∇⋅1∇I 0 ( X ) − ρ ( X )σ a I 0 ( X ) =3ρ ( X )σ tr− Q0 ( X ) + ∇ ⋅Boundary Conditions(4)

1
Q1 ( X )4πρ ( x)σ tr

where

σ tr = σ t (1 − ap1 )is the “transport cross section”,

I d ( X , ω ) (ω n) dω = 0,

p (ω , ω ′) I ri ( X , ω ′)dω ′Q ( X , ω ) = ρ ( X )σ s

2π

4π

is the first scattered external illumination,
Q0 ( X ) =1
4π

Q( X , ω )dω

(5)

4π

is the first scattered external illumination averaged
over the unit sphere, andQ1 ( X ) =1
4π

Q( X , ω ) ω dω

For the boundary conditions at the surface of thevoxel regions containing non-zero droplet density, wereally want the inward diffuse radiance Id(X, ω) to be
exactly zero for all directions ω with ω n > 0, where
n is the inward surface normal. But this is impossible
for a function of the form I0(X,ω)+ I1(X, ω) ω unlessthe function is zero for all ω. So an approximate
condition of no total inward flow is used:

(6)

4π

is its directionally weighted average, a vector.
The incident sun and sky illumination are assigned todirection bins on the subdivided surface of a unit
cube, as shown in Figure 1b). The power in bin i is
Isky(ωi)∆ωi , where ∆ωi is the solid angle of bin i.
This power is propagated into the voxels of the cloudvolume to compute Iri(X,ωi) by tracing a dense
collection of rays in the direction ωi , starting fromthe external faces of the box enclosing the cloud.
Each ray intersects the voxels in a collection ofsegments Sj, along which the attenuation integral in(2) can be accumulated. The lengths l(Sj) of thesesegments weight the contributions of Iri(X,ωi) to theintegrals in (5) and (6) at the voxel X. So thecontribution to the integral (5) of a ray segment Sj isρ(X) σs Iri(X,ωi) ∆ωi l(Sj). The phase functionp(ω,ω’) depends only on the angle θ between ω andω’. Thus the contribution of light flowing in direction
ω’ = ωi to the integral (6) lies only along direction
ωi , since the other components cancel by symmetry.
It therefore equals ρ(X) σs Iri(X,ωi) ∆ωi l(Sj) p1 ωi .
After these contributions are summed for all rays indirection bin i, the values of the sums for (5) and (6)in each voxel are normalized by dividing by the sumof the lengths l(Sj) of the segments it contains. The
scalar Q0(X) and the vector Q1(X) are saved at eachvoxel X.
The sky and sun colors were computed taking into
account Rayleigh scattering from the molecules in theatmosphere, using the methods of [Nis96], [Dob97],
and [Pre99], and Mie scattering from aerosols, using
the methods of [Slo02].

where the integral is over the inward hemisphere.
After some mathematical manipulations, this resultsin the boundary condition:
I0 ( X ) −2n ⋅ ∇I 0 ( X ) +3 ρ ( X )σ tr

1n ⋅ Q1 = 0
2πρ ( X )σ tr(7)

For more details on the derivation of this equationand of equation (4), see [Ish78].

Finite Difference Solution
To solve the partial differential equation (4) we use
finite differences to estimate the derivatives. The
difference for the gradient of I0 is evaluated at theface centers, and for the value of ρ in the first term ofequation (4), the average of the ρ(X) for the two cells
sharing the face is used, unless one of them has ρ = 0,in which case the value at the other cell is used.
Central differences are used to evaluate the gradientfor the last term of equation (4). All this results in aset of linear equations relating each value I0(X) to thevalues at its six adjacent neighbors. The boundary
condition (7) is enforced by assuming a temporary
value of I0(Xe) at an adjacent empty cell Xe where thedroplet density is zero, set to satisfy (7), and using it
to compute the differences for equation (4). Thisaffects the coefficients of the linear equations, which
end up not involving I0(Xe), so the temporary value ofI0(Xe) can be different for different full cells adjacentto the same empty cell Xe.
The resulting set of linear equations can be solved
iteratively by the Gauss-Seidel or Jacobi method, bysetting each I0(X) to a weighted sum of its old value
and the old values of its six neighbors, plus a constantterm involving Q0 and Q1 that comes from theconstant right hand side of the linear equation. This
will take a long time for a large grid, since a brightspot caused by sunlight hitting one side of the cloud
can only propagate across one adjacent cell per
iteration, so it will take many iterations for it to affectthe opposite side of the cloud. Therefore a multigrid
method is used. This is like mip mapping. The datavolume is considered as a multi-resolution octree
hierarchy, and coarser versions of the linear equations
are written for lower levels of the hierarchy. Theirsolutions are used to correct the solutions at finerlevels of the hierarchy to allow for faster propagation,and the solution method iterates up and down thelevels in the hierarchy until convergence. For details,see [Bri00]. We used the Hypre system API and
multigrid solver from Lawrence Livermore National
Laboratory [Hypre03].

3. PROPAGATION IN CLEAR AIR

If Xe is an empty border voxel with ρ(Xe) = 0,
adjacent across face F to a full voxel X, with ρ(X) >0, as shown in Figure 1a), we place a viewing cube,
shown in Figure 1b), at the center of Xe. Then, as in
the hemicube algorithm for radiosity, we render thescene, including the clouds, onto the six faces of thisviewing cube. We used 64 by 64 resolution images
on each cube face, and then averaged the colors into
the 4 by 4 direction bins for each face, as shown in
Figure 1b), to give a total of 96 direction bins for thesix faces.

X

Xe
F

G
a)

E

E’
b)

Figure 1. a) Light propagating from an empty cellto a full cell. b) A viewing cube divided into bins.

Figure 1a) shows several rays in the direction ωi atthe center of direction bin i, crossing face F. Insteadof continuing to propagate and attenuate Iri(X,ωi) toaccumulate the values of Q0 and Q1 as describedabove, we replace the propagating value by theaverage color seen in direction bin -ωi, which isindicated in bold lines in Figure 1b), in order toinclude the inscattered light from the clouds. If F is
the only face across which Xe is adjacent to a full cell,
all directions ωi crossing F from Xe to X have flow
towards the left, so we need only render the hemicube
to the right of the bold line in Figure 1b). Similarly, if
Xe is adjacent to only two full voxels, across faces F
and G that share a common edge E, then the needed
incoming directions fill two overlapping hemicubes,
and we do not render the two half faces adjacent to
the dashed edge E’ shown in Figure 1b).
The diffusion approximation assumes that the clouds
are dense, so that enough multiple scattering takes

place to make the radiance almost independent of
direction. As described above and in [Ish78], I1(X)
can be computed from I0(X) and its derivatives.
However the linear approximation (3) for I(X) is not
valid at the edges of the cloud, where single
scattering, and orders of multiple scattering too low
to insure directional independence, are important.
Therefore the approximation (3) cannot be used torender an image of the cloud surface.
Instead, hierarchical volume rendering is used torender the volume with opacities ρ(X)σt and voxel
colors I0(X). (The solution described above is done
for the red, green, and blue wavelength bands.) An
octree hierarchy of least squares approximations tothe color and opacity values on the various level
cubes is precomputed, together with the RMS error ofeach approximation. Then, in a view dependent
manner which weights the error by the projected areaof the cube on the image, appropriate cubes are
selected, and composited in back to front order, using
the hardware cell-projection method. The details are
given in [Schu01]. For this application, the algorithm
in [Schu01] was modified to handle color, and to use
either piecewise linear or piecewise constantfunctions in the approximation hierarchy.
The revised values for Q0 and Q1 obtained bypropagating the incident radiance from this volume
rendering are used again in obtaining a new solution
to the partial differential equation (4), and thesequence of PDE solution, volume rendering, and ray
propagation is iterated until I0 stops changing
significantly. The illumination and shadowing of theground by the sky and clouds is also found using
hardware volume rendering onto hemicubes placedon a ground sample grid, and the diffusely reflecting
ground is included as a background to the volume
rendered radiance images of the cloud. The skyhemisphere, colored as described above, is also partof the background. However the radiance of the sun
exceeds the dynamic range of the fixed-point frame
buffer, so the sun is not added to the sky background.
Instead,the direct sunlight, Isun(Xe,ωsun) is
propagated and attenuated separately, and saved as
Irisun(Xe,ωsun) at every voxel. Its effects are added
separately to Q0 and Q1.

4. FINAL GATHER
The final rendering is a “final gather” ray tracing of
the volume solution to equation (4). The reduced sun
flux Irisun(Xe,ωsun) is used, together with the exact
phase function p(ωsun, ωviewing) to give an exact final
single scattering of the attenuated sunlight, similar totwo pass volume rendering for shadow effects. Thiscan give the brightening effect on the edges of acloud when the sun is directly behind it. However, itcannot account for the multiple directional scattering

of the sunlight, nor of the skylight, which contributes
to the darker and bluer color of the edges of thecloud, when seen against other non-edge cloud
regions as a background.
To include a more accurate multiple directional
scattering effect for the regions near the edges of thecloud, where the particle density is below a thresholdT necessary for the diffusion approximation to beadequate, we again use volume rendering. When
doing the final gather ray tracing per pixel, we
determine low-density ray segments Sj where the ray
stays within the edge region where 0 < ρ(X) < T. For
each such segment, we integrate ρ(X)σt along thesegment to get the optical depth Dj, put a cube at the“center” point cj where the integral reaches half thisvalue Dj, and compute the radiance on the faces of
the cube using volume rendering, as above. The x
axis of the cube is aligned with the direction of theviewing ray, as shown in Figure 2.
C

G
H

D

K
L

I

B

J
F

E

A

Figure 2. Faces and subface of the viewing cube.

Knowing Dj, the Poisson distribution [Poi97] can tellus the probability qij of scattering of scattering
exactly i times along this segment:qij =exp(− D j ) D iji!

pixelwise by the radiance images on the cube andsummed to give a convolution sample Kij. Finally the
Kij are multiplied by the probabilities qij and summed
over i to get the contribution Ej of the segment Sj to
the viewing ray through a pixel:
Ej =

(8)

This is only an approximation, because the scattering
will change the direction of the ray. However, forsmall water droplets inside clouds, the scattering ishighly peaked in the forward direction, so for the first
few orders of scattering, it is a good approximation.
To get the phase function for the result of i scattering
events, we must convolve the phase functionp(ω,ω’) with itself i - 1 times, using integration overthe unit sphere. This is particularly simple for theHenyey-Greenstein phase function, because the result
is the same formula, with g replaced by gi, as can beproved by expanding p(ω,ω’) as a Legendra series incosθ. These multiple scattering phase functionspi(ω,ω’) are precomputed, and saved for each pixelon the faces of the viewing cube shown in Figure 2.
Then for each ray segment Sj, they are multiplied(9)

If the centers cjkl of a collection of segments Sjkl for aregion of pixels (k, l) form a continuous surface, there
is a possible speed-up. The higher order (large i)scattering convolutions Kijkl of the radiance images
are blurred, and thus change slowly with the pixel
index (k, l). Therefore they can be interpolated from asparser sampling of the pixels in the region. Only theprobabilities qijkl needs to be computed at every pixel.
On the other hand, for smaller i, only the front face ofthe cube containing the viewing direction need beaccurate, because the scattering is mainly forward.
For even smaller i, only a central section of the front
face is required. This saves rendering time because
the clipping to the view frustum is performed while
traversing the octree hierarchy, so clipping eliminates
large octree cells, and their children are not evenconsidered.
We used three different viewing cube sampling
patterns. The first is for the rear five faces of theviewing cube in Figure 2: ABCD, AEFB, BFGC,
CGHD, and DHEA, which are sampled initially alongrays at every 8 output image pixels in x and y, as
indicated by the * symbols in Figure 3 below. This
figure shows the method of computing theconvolutions for these five faces, for the first lowdensity segment, if any, along the viewing rays forpart of a low resolution image.
2 1.

qij K ij .
i2 - 4 1

1 - 1 1 1 - - - 1 1

1

1 2 - 2 - 2 - - - 2 1

2 - 2 -

1 - - - - - - 1 - - 1 1

1 1 1 - - - -

4 - - - 4 - - - 4 - - - 4 - - - 4 - - - 4 - - - - - - - - - - - - 1 - - - - - - - - - - - - - - - - - - - - 2 - - - - - 2 - 2 - - - - - - - - - - - - - - - - - - - 1 1 1
4 - - - * - - - - - - - * - - - 4 - 2 1- - - - - - - - - - - - - - - - - - - 1
2 - 2 - - - - - - - - - - - - - - - 2 - 2
1 1 1 - - - 1 - - - - - - - - - - - - - 1 2 - - - - - - - - - - - - - 4 - - - 4
1 - - - - 1 - - - - - - - - - - - - - 1 2 - - - - - - - - - - - - - - - - - - 1 - - - - - - - - - - - - - - - - - - - 2 - * - - - - - - - * - - - - - - - * -Figure 3. Convolution computation locations.

The convolutions Kijkl are bilinearly interpolated tothe positions indicated by the symbol “–” if thevalues of the distances along the ray at which theviewing cubes were placed are close enough. If not, anew viewing cube is placed at the appropriate
distance along the ray, and new values of Kijkl arecalculated. These new calculations are done
hierarchically at the vertices of a quadtree, in order tofind values that may be satisfactory for interpolationin one or more sub-squares. Thus a 4 in Figure 3
indicates a new viewing cube along a ray at an output
image pixel with x and y coordinates divisible by 4, a
2 indicates a new viewing cube along a ray at a pixel
with x and y coordinates divisible by 2, and a 1
indicates a new viewing cube for a pixel whose x or ycoordinate is an odd number. A blank indicates thatthe viewing ray had no low-density segments.
The second sampling pattern was for the front face
EFGH of the viewing cube, which was initiallysampled every 4 pixels in x and y, with doubled
resolution, in order to give more accurate interpolated
convolutions near the peak of the forward scattering
phase function. The final pattern uses only the center
quarter IJKL of the front face of the viewing cube,
initially sampled every 2 pixels in x and y. The total
pixels were the same as on the front face in theprevious case, so the pixel spacing is half the spacing
on the whole front face, and a quarter the spacing on
the other five cube faces. This central section isdisregarded when computing the convolution for thewhole front face. Thus, every direction is counted
exactly once in the three methods, and theconvolutions which have been interpolated orcomputed for each of them can be added to get the
total convolution for each low-density ray segment.
The convolutions are computed hierarchically on 9
by 9 blocks of pixels, so that when interpolation isallowed, the convolutions to be interpolated are
already known. The blocks are computed inhorizontal sweeps across the image, so that the left
hand column of pixels is known from the previous
block, and the top row of pixels is known from theprevious sweep of blocks. Thus the convolutions
need only be saved for the right hand column of theprevious block, and for one full image row,
representing the bottom row of the previous sweep’s
blocks. Once the convolutions are known for the lowdensity segments, equations (8) and (9) are evaluatedper pixel, and the voxel segments on the viewing ray
are composited from front to back, using either the
low-density convolution solution (appropriately
weighted by the optical depth of the voxel segment),
the high-density diffusion solution, or a density
dependent interpolation of the two (in order toeliminate aliasing artifacts from switching solutions).

5. RESULTS
We used the final time step of a cloud simulation
produced by a hierarchical octree enhancement of thesimulation method of [Miy02], using the stable semiLagrangian advection scheme of [Stam99] and thevorticity confinement method described in [Fed01].
The simulation was at resolution 200 x 160 x 120,and was averaged down to 100 x 80 x 60 cells for
rendering. A rectangular solid enclosing the cells with
non-zero droplet density had size 91 x 51 x 25.
Figure 4 shows a cloud with only single scattering of
the sun illumination. Since the phase function for the
droplets has p1 = g = 0.9, and the scattering isstrongly peaked in the forward direction, with little
scattering towards the viewpoint, the cloud is very
dark and has been artificially brightened. Figure 5
show the final gather of the diffusion approximation
solution with both sun and sky illumination, but nocorrection for the low-density boundary regions, andno propagation across the clear air. Figure 6 shows
the added effects of the volume rendering method of
accounting for the low-density regions at the edge ofthe cloud, using only the initial sun, sky and groundillumination in the diffusion solution. Figure 7 shows
the result of one iteration of the clear air propagation,using volume rendering at the border voxels to image
the clouds. If one flips between Figures 6 and 7 in ascreen image viewer, one can see that some regionsof the cloud look brighter in Figure 7. An image from
two iterations of the clear air propagation isindistinguishable from Figure 7.
Figure 7 took 21 minutes and 33 seconds to compute,on a 1.7 GHz Pentium 4 with nVidia GeForce3
graphics. Of this time, 46 seconds was spent setting
up and solving the linear equations for diffusing theincident sunlight, 2 seconds for setting up the octree
hierarchy, 227 seconds to render the 12460 viewing
cubes at the border cells around the clouds, at
resolution 64 by 64 per face, and the 32 by 32 arrayof hemicubes on the ground sample grid, 94 seconds
to propagate the flux from the ground, sky, and cloudin 96 directions though the cloud, another 46 seconds
to solve the revised linear equations taking into
account this extra flux, and 854 seconds to render theimage, including the time to render the viewing cubesat the low-density segments, compute theconvolutions, interpolate and weight them, and ray
trace the 480 by 360 image.
Figure 8, at a time near sunset, took 24 minutes and
45 seconds. The final rendering time took longer,1034 seconds, because there were more cloud pixels
in this image, due to different viewing parameters.

resolution to the final gather pass, by ray tracing theoriginal volume instead of the averaged version, andinterpolating the radiance solution. We also hope to
add a subvoxel-resolution Perlin noise texture and toadvect the texture coordinates using the velocity fromthe weather simulation.

Figure 4. Single scattering only.

Figure 7. Diffusion solution with correction for
low-density regions, after one iteration of
propagation through clear air.

Figure 5. Initial diffusion solution.

Figure 8. Image at sunset.

Figure 6. Initial diffusion solution with correction
for low-density regions.

6. DISCUSSION AND FUTURE WORKThe hardware rendering component of our work
suffered from the dynamic range limitations of the 8
bit fixed point format of our graphics hardware
pipeline. We hope to produce more accurate results
using floating-point graphics pipelines and frame
buffers. The computing time and memory sizes limit
the data resolution for our diffusion and lightpropagation passes, but we hope to add extra

The artifacts near the bottom of Figures 6 and 7 comefrom the transitions where two low-density segments
join into one. The probabilities qij for larger numbersof scattering events are higher for the larger opticaldepth of the single low-density segment, so that more
illumination is included from the sides, and less fromthe front. We should add some kind of transition to
alleviate the edges caused by this abrupt switch.
Currently for the volume rendering in the final gather,
all cells except the one containing the viewpoint cj
are rendered. We should also remove any cellintersecting the low-density segment Sj, in order toavoid double counting the scattering along thissegment.

7. ACKNOWLEDGMENTS
This work was performed under the auspices of the
U.S. Department of Energy a) by the University ofCalifornia, Lawrence Livermore National Laboratory
under contract number W-7405-ENG-48, and b) at
the University of California, Davis under
Memorandum Agreement No. DE-FC02-01ER41202
through the SciDAC Program. We thank the WSCG
reviewers, Kwan-Liu Ma, and Fabrice Neyret forsuggestions that improved the exposition, Fabrice
Neyret for the conversations that inspired this work,
Jos Stam and Phil Collela for advice, Jim Jones and
Rob Falgout for help with Hypre, Jos Stam for his
multigrid code which we used before switching toHypre, and Robin Bing-Yu Chen, Henry Johan, andKenichi Amou for system help at the University ofTokyo, where the first author did much of this work
as a visiting professor.

8. REFERENCES[Bri00] Briggs, W., Henson, V. E., and Steve
McCormick, S., “A multigrid tutorial: Second
edition”, SIAM, Philadelphia, 2000.
[Dach03] Dachsbacher, C., and Stamminger, M.,“Translucent shadow maps”, Eurographics
Symposium on Rendering 2003, pp. 197 – 201.
[Dob97] Dobashi, Y., Nishita T., Kaneda, K., andYamashita, H., “A fast display method of skycolor using basis functions”, The journal ofvisualization and computer animation, Vol. 8 No.
2, pp 115 – 127, 1997.
[Fed01] Fedkiw, R., Stam, J., and Jensen, H. V.,“Visual simulation of smoke”, Proceedings ofSiggraph ’01, pp. 15 – 22, 200.
[Ish78] Ishimaru, A., “Wave propagation andscattering in random media. Volume 1: Single
scattering and transport theory”, Academic Press,
New York, 1978.
[Jen98] Jensen, H. W., and Christensen, P., “Efficientsimulation of light transport in scenes withparticipating media using photon maps”,
Proceedings of Siggraph 1998, pp. 311 – 329.
[Jen01] Jensen, H. W., Marchner, S., R., Levoy, M.,
and Hanrahan, P., “A practical model forsubsurface light transport”, Proceedings ofSiggraph 2001, pp. 511 – 518.
[Jen02] Jensen, H. W. and Buhler, J., “A rapid
hierarchical rendering technique for translucent
materials, ACM Transactions on Graphics Vol.
21 No. 3 pp. 576 – 581, 2002.

[Kaj84] Kajiya, J., and von Herzen, B., “Ray tracing
volume densities”, ACM Computer Graphics(Proceedings of Siggraph 1984) Vol.18 No. 3,
pp. 165 – 174.
[Hyp03] http://www.llnl.gov/CASC/hypre[Len02] Lensch, H., Goesele, M., Bekaert, Ph.,
Kautz, J., Magnor, M., Lang, J., and Seidel, HP., “Interactive rendering of translucent objects”,
Proceedings of Pacific Graphics 2002, pp. 214 –224.
[Mer03] Mertens, T., Kautz, J., Bekaert, Ph., Seidel,
H-P., and Reeth, F. V., “Interactive rendering oftranslucent deformable objects”, Eurographics
Symposium on Rendering 2003, pp. 130 – 140.
[Miy02] Miyazaki, R., Dobashi, Y., and Nishita, T.,“Simulation of cumuliform clouds based oncomputational fluid dynamics”, Proceedings ofEurographics 2002 Short Presentations, pp. 405410, 2002.
[Nis96] Nishita T., Dobashi, Y., Kaneda, K., andYamashita, H., “Display method of sky color
taking into account multiple scattering”, PacificGraphics ’96, pp. 117 – 132, 1996.
[Per97] Perez, F., Pueyo, X., and Sillion, F., “Global
illumination techniques for the simulation ofparticipating media”, in “Rendering Techniques’97”, Dorsey and Slusallek, editors, Springer,
Vienna, pp. 309 – 320, 1997.
[Poi97]http://info.bio.cmu.edu/Courses/03438/PBC97
Poisson/PoissonPage.html[Pre99] Preetham, A. J., Shirley, P., and Smits, B.,“A practical analytic model for daylight”,
Proceedings of Siggraph 1999, pp. 91 – 100.
[Prem03] Premoze, S., Ashikmin, M., and Shirley,
P., “Path integration for light transport involumes”, Proceedings of the Eurographics
Symposium on Rendering, pp. 52 – 63, 2003.
[Schu01] Schussman, G., and Max, N., “Hierarchicalperspective volume rendering using triangle fans”,
Volume Graphics 2001, Stony Brook, NY, pp.
195 – 200.
[Slo02] Sloup, J., “A survey of the modelling andrendering of the earth’s atmosphere”, SCCG
2002, pp. 141 – 149, 2002 .
[Stam95] Stam, J., “Multiple scattering as a diffusion
process”, in “Rendering Techniques ’95”,
Hanrahan and Purgathofer, editors, pp. 41 – 50,
Springer, Vienna, 1995.
[Stam99] Stam, J., “Stable fluids”, Proceedings ofSiggraph ’99, pp. 121 – 128, 1999.