A Simple, Efficient Method
for Realistic Animation of Clouds
Yoshinori Dobashi* Kazufumi Kaneda** Hideo Yamashita** Tsuyoshi Okita* Tomoyuki Nishita******

Hiroshima City University*{doba, okita}@im.hiroshima-cu.ac.jp

Hiroshima University

**

***

University of Tokyo{kin, yama}@eml.hiroshima-u.ac.jp

Abstract
This paper proposes a simple and computationally inexpensive
method for animation of clouds. The cloud evolution is simulatedusing cellular automaton that simplifies the dynamics of cloudformation. The dynamics are expressed by several simple
transition rules and their complex motion can be simulated with asmall amount of computation. Realistic images are then createdusing one of the standard graphics APIs, OpenGL. This makes it
possible to utilize graphics hardware, resulting in fast image
generation. The proposed method can realize the realistic motionof clouds, shadows cast on the ground, and shafts of light throughclouds.
CR Categories: I.3.1 [Computer Graphics]: Hardware
Architecture; I.3.6 [Computer Graphics:] Methodology andTechniques; I.3.7 [Computer Graphics]: Three-Dimensional
Graphics and Realism; I.6.3 [Simulation and Modeling]:
Applications;

Additional Keywords: Animation, Atmospheric Effects,
Rendering, Graphics Hardware, Volume Rendering.

1. INTRODUCTION
Clouds play an important role when making images for flight
simulators or outdoor scenes. Their color and shapes change
depending on the position of the sun and the observer. This meansthat the density distribution of clouds should be defined in threedimensional space to create realistic images. Therefore, a lot of
methods have been developed to display clouds [33, 10, 5, 14, 19,
30, 25]. Using these methods, extremely realistic images can begenerated. Their main purpose is, however, to create images ofstatic clouds. Fascinating animations of clouds with changing
their shapes and color are often used, however, in movies,*3-4-1, Ozukahigashi, Asaminami-ku, Hiroshima, 731-3194 Japan(Dobashi's current address: Hokkaido University, Faculty ofEngineering, Kita 13, Nishi 8, Kita-ku, Sapporo 060-8628 Japan)**1-4-1, Kagamiyama, Higashi-hiroshima, 739-8527 Japan***7-3-1, H ongo, B unkyo-ku, Tokyo, 113-0033 Japan***nis@is.s.u-tokyo.ac.jp

commercial films and so on. They are often created by filming
them in advance and replaying the film quickly. Since generatingsuch realistic animations by computer graphics is useful, a lot of
methods have been developed [13, 29, 6, 31, 7, 22, 9, 15, 8, 32].
This paper proposes a new method for realistic animation ofclouds. Our aim is to develop a simple method that can create
realistic animation as quickly as possible, preferably in real-time.
We propose an efficient simulation method and a hardwareaccelerated rendering method. In particular, the method is suitable
for animation of cumulus-like clouds in landscape scale. Our
method has the following features.
• Our simulation method creates realistic cloud motion with asmall amount of computation. It uses cellular automaton thatcan simulate the motion just by simple Boolean operations.
• Our rendering method realizes a fast computation of photorealistic images. It can quickly calculate shadows and shafts oflight through clouds, as well as cloud color, by making the mostof graphics hardware.
A straightforward approach to creating realistic cloud motion is tosimulate the physical phenomena. That, however, is impracticalsince it is computationally expensive. Therefore, a simple andefficient method is required that maintains the visually convincing
result. As one such methods, Nagel et al. extremely simplified thecloud dynamics using cellular automaton [21]. The method cansimulate cloud formation by simple transition rules. Unfortunately,
the method is not sufficient for our purpose since their aim is not
to create realistic animation. So, we extend their work andpropose a new method to realize realistic cloud evolution.
Furthermore, we propose a method for generating realistic images
including cloud shadows and shafts of light using the standard
graphics API, OpenGL. This results in fast image generation since
we can make use of graphics hardware. Our rendering method candisplay the following three effects: 1) cloud color taking into
account the single scattering of light, 2) shadows of clouds cast on
the ground, 3) shafts of light through clouds. Most of the previous
methods render these effects by using ray-tracing, one of the mosttime-consuming methods. To overcome this problem, we proposea hardware-accelerated rendering method based on OpenGL.
Using the proposed method, a photo-realistic image can begenerated within one minute on a standard PC.
Our method is not sufficient if the user needs physically exact
cloud motion, since it is one of the numerical models that
simplifies the physical phenomena. Our method is suitable forusers who want a simple, easy-to-use, and computationallyinexpensive method that can create visually convincing results.
Furthermore, the proposed method makes as much use of graphics
hardware as possible. Since graphics hardware is becoming fasterand faster, we believe our method is one of the promisingtechniques that will realize real-time animation in the near future.

2. PREVIOUS WORKIn this section, previous works are briefly reviewed. Methods
related to the simulation and methods to rendering are separately
reviewed.
SIMULATION: In computer graphics, there are two categories
to simulate the gaseous motion like clouds. One is to simulate thephysical process of fluid dynamics [13, 29, 31, 9, 32]. The other isa heuristic approach [10, 6, 7, 22, 15, 8]. Most of the methods inthe former category need a large amount of computation time.
Stam, however, developed a fast simulation method bysimplifying fluid dynamics [32]. He demonstrated a real-timeanimation of smoke on a high-end workstation. However, since
our purpose is to simulate clouds covering a large area in the sky
using a standard PC, it is sill time-consuming. Furthermore, thephase transition effects from vapor to water should beincorporated to simulate cloud formation. To our knowledge, themethod developed by Kajiya et al. is the only one to include thephenomena [13]. It is, however, very complex and timeconsuming. Controlling cloud shapes is also difficult by using themethods in this category.
On the other hand, the latter approaches, such as procedural
modeling, are computationally inexpensive and much easier toimplement. The disadvantage in these methods, however, is thatthe user has to search parameters by trial and error to create
realistic animation. Creating realistic-looking motion using themdoes not seem to be an easy task.
Our method lies in an intermediate position between these twocategories. The method reflects the physical formation process ofclouds in part, and it is computationally more efficient and easier
to implement than are the previous physical simulation methods.
Our method can create more realistic cloudy scenes than thephysical based methods do. Controlling the shapes and theirmotion to create realistic animation is an easy task using our
method. Dobashi et al. tried to develop such a method using
cellular automaton [3]. Unnatural animation is created by that,
however, since the formation and extinction of clouds are
frequently repeated. The method in this paper is an extension totheir method.
RENDERING: One of the simplest ways to display clouds is touse mapping techniques, such as the method developed by
Gardner [10]. In order to display photo-realistic images, however,
it is desirable to use the physical model, taking into accountscattering/absorption due to particles. Many such methods havetherefore been developed [13, 5, 14, 28, 29, 19, 31, 4, 32]. Someof them take into account multiple scattering of light [13, 19, 31,25]. Additionally, Nishita et al. take into account the effect ofskylight on the cloud color [25]. Multiple scattering and skylightare important for realistic image synthesis but is time-consuming.
Our method approximates them as a constant ambient term. Oneof the major approaches to rendering the volume density similarly
to clouds is to use 3D textures. Stam used 3D hardware texture
mapping to display gases [32]. With the help of the high-end
graphical workstation, the method can generate realistic images inreal-time by combining 3D textures and advecting cloud texturesdeveloped by Max et al. [18]. Unfortunately, his method is not
sufficient for our purpose since the method does not include theatmospheric effects such as shafts of light, one of the essential
factors in generating realistic images of outdoor scenes. Our
method can handle these effects. Although using 3D textures is
simpler and efficient, the 3D texture mapping hardware is still
expensive and not universally available as are 2D textures. Since
one of our objectives is to create realistic animation on a standardzxy

tim e t i+1

tim e t iact = 0act = 1

hum = 0hum = 1

cld = 0cld = 1(a) Sim ulation process.
sm oothing

volum e rendering(b) Rendering process.
Figure 1: O verview of our m ethod.

PC, we propose a method using 2D textures. The idea of themethod is of course applicable to 3D textures.
There are also a lot of methods for generating shafts of light [23,11, 12, 16, 17, 27, 24]; most of them use the ray-tracing algorithm
and require dozens of minutes. Although Nishita and Nakamae
developed a faster method using a scanline accumulation buffer[24], it still requires several minutes. We propose a much fastermethod using hardware color blending and texture mapping
functions. To our knowledge, no other methods are available
making use of graphics hardware to render shafts of light.

3. BASIC IDEA
Fig. 1 shows an overview of our method. Our method consists of
two processes, simulation and rendering. As shown in Fig. 1(a),the simulation space is divided into voxels. The voxels correspondto cells used in the cellular automaton. At each cell, three logicalvariables, vapor/humidity (hum), clouds (cld), and phasetransition (or activation) factors (act) are assigned. The state ofeach variable is either 0 or 1. Cloud evolution is simulated byapplying simple transition rules at each time step. The transition
rules represent formation, extinction, and advection by winds.
Since the state is either 0 or 1, the rules can be expressed by
Boolean operations. Therefore, each variable can be stored in one
bit to save the memory cost and the simulation process is
accelerated by using bit field manipulation functions.
Images are generated in the rendering process by making use of
the simulation results (Fig. 1(b)). As described above, what we
can obtain from the simulation is no more than there are clouds(cld = 1) or, there are not-clouds (cld = 0) at each voxel.
Therefore, a density at each point is calculated by smoothing thebinary distribution as shown in Fig. 1(b). The clouds are thenrendered using volume rendering techniques. The rendering
process consists of two steps. The first step calculates the intensityof light reaching the center of each voxel. Cloud shadows are also
calculated in this step. The shadows are obtained as a texture.
Then, in the second step, images are generated. Clouds are
rendered by using a splatting method [1]. To render shafts of light,
we consider multi spherical shells with their center at the
viewpoint (see Fig. 7). The shells are then drawn from back tofront using the hardware alpha-blending function. Shafts of lightare rendered by mapping the shadow texture on the shells (seesection 5.2).

tim e ti

tim e t i+1

4. SIMULATION METHOD
We extend the following four points to Nagel’s method [21]:• Extinction of clouds• Wind effects• Speeding up of the simulation• Controlling cloud motion
Details of the above extensions are explained in sections 4.2
through 4.5 after the brief description of Nagel’s method in thenext section.

4.1 Growth Simulation

z
x

z
y

x

y

hum
cld
act = 0act = 1

hum = 0hum = 1

cld = 0cld = 1Figure2: B asic transition rules.

In this section, Nagel’s method to simulate the cloud formation
process is described briefly. The physical processes of cloudformation are outlined as follows. Clouds are formed as a bubbleof air is heated by underlying terrain heat, causing the bubble tobecome less dense, and to rise into regions of lower pressure inwhich the bubble expands. Expansion cools the bubble, increasing
the relative humidity inside the bubble. The phenomenon calledphase transition then occurs, that is, water vapor in the bubble
becomes water droplets, or clouds. Nagel et al. used a cellular
automaton [35] to simulate these processes in the following way.
For simplicity, the simulation space is aligned parallel to xyz axes
and the number of cells is assumed to be nx × n y × nz . As
mentioned before, three logical variables, hum, act, and cld, areassigned at each cell (see Fig. 1(a)). Each represents vapor, phasetransition factor, and clouds. The state of each variable is either 0
or 1. hum=1 means there is enough vapor to form clouds, act=1
means the phase transition from vapor to water (clouds) is readyto occur, and cld=1 means there are clouds. In the following,A ∧ B and A ∨ B indicate conjunction and disjunctinon between
A and B, respectively, and ¬A indicates negation of A. Theirtransition rules are given as follows.
hum(i, j , k , t i +1 ) = hum(i, j , k , t i ) ∧ ¬act (i, j , k , t i )cld (i, j , k , t i +1 ) = cld (i, j , k , t i ) ∨ act (i, j, k , t i )

act

,

,

act (i, j , k , ti +1) = ¬act (i, j , k , ti ) ∧ hum(i, j , k , ti ) ∧ f act (i, j , k ) ,

(1)
(2)
(3)

where fact(i,j,k) is a Boolean function and its value is calculated by
the status of act around the cell. The following function is used bytaking into account the fact that clouds grow upward andhorizontally.
f act (i, j , k ) = act (i + 1, j , k , ti ) ∨ act (i, j + 1, k , ti )∨ act (i, j, k + 1, ti ) ∨ act (i − 1, j , k , ti ) ∨ act (i, j − 1, k , ti )∨ act (i, j, k − 1, ti ) ∨ act (i − 2, j , k , ti ) ∨ act (i + 2, j, k , ti )determined randomly and cld is set to zero. For more details,
please refer to [21].

4.2 Cloud Extinction
One of the disadvantages of Nagel’s method is that cloud
extinction never occurs since cld, after it has become 1, remains 1
forever. Dobashi et al. introduced a new state variable to solve
this problem [3]. In that method, however, formation andextinction are repeated frequently, resulting in unnatural
animation. In the real world, cloud extinction is caused by gradualtransition of water droplets to vapor. Our method simulates thecloud extinction as follows. First, the animator specifies cloud
extinction probability, pext. Next, at each cell whose cld is 1, a
random number, rnd (0 ≤ rnd ≤ 1), is generated and cld ischanged to 0 if rnd < pext. By changing the probability at each cellat different times, the animator can specify regions where cloud
extinction occurs frequently. Although this realizes the cloud
extinction, there remains another problem. Clouds are nevergenerated after the extinction at the cell. To solve this, vapor(hum) and phase transition factors (act) are supplied at specified
time intervals. Similar to extinction, vapor probability, phum, andphase transition probability, pact, are used to set them randomly.
That is, hum is changed to 1 if rnd < phum and act is changed to 1
if rnd < pact. Cloud motion can be controlled by controlling theprobabilities, phum, pact, and pext at each cell at each time step. The
methods described in this section are summarized by thefollowing three transition rules.
cld (i, j , k , t i +1 ) = cld (i, j, k , t i ) ∧ IS(rnd > p ext (i, j , k , t i )) ,hum(i, j , k , t i +1 ) = hum(i, j , k , t i ) ∨ IS(rnd < p hum (i, j , k , t i )) , (6)act (i, j , k , t i +1 ) = act (i, j , k , t i ) ∨ IS(rnd < p act (i, j , k , t i )) ,(4)

∨ act (i, j − 2, k , ti ) ∨ act (i, j + 2, k , ti ) ∨ act (i, j, k − 2, ti ).
Of course, there are variations of the above function. We have
tried some of them and we couldn’t see significant differences ofthe simulation. So, we use the function in the original paper. The
rules are summarized in Fig. 2. As shown in the top column of Fig.
2, act becomes 1 if hum is 1 and the state of act of one of theshaded cells around the center cell (i, j, k) is 1. Then, hum
becomes 0 as shown in the middle. Finally, as shown in thebottom of Fig. 2, cld becomes 1. As a boundary condition, theirstates are assumed to be 0 outside the simulation space. Beginningfrom initial random status (all 3 status are set randomly), cloud
growth is simulated by updating the state of each variable using
Eqs. 1 through 4. As for the initialization, hum and act are

(5)
(7)

where rnd is a uniform random number, IS(e) is a Booleanfunction that returns 1 if the expression e is true, otherwise returns
0.

4.3 Advection by Wind
We can observe clouds moving in one direction, blown by winds.
New transition rules are introduced to include the wind effect. The
idea is simply to shift all the variables toward the wind direction.
We assume, for simplicity, the wind blows toward the direction of
x-axis. Other cases can be handled by rotating the simulation
space according to the wind direction. Furthermore, it is wellknown that the wind velocity is different depending on the heightfrom the ground. The wind velocity, v(zk), is therefore specified asa function of z-coordinate of each cell (i,j,k). To implement thewind effect in the context of the cellular automaton, the function,

v(zk), is assumed to return integer values. The transition rules are
as follows.
sunlight

hum(i − v( zk ), j , k , ti ), i − v( zk ) > 0
hum(i, j , k , ti +1) = ,otherwise
 0,

(8)

cld (i − v( zk ), j , k , ti ),cld (i, j, k , ti +1 ) = 
 0,

(9)

i − v( zk ) > 0,
otherwisei − v ( zk ) > 0
act (i − v( zk ), j , k , ti ),act (i, j , k , ti +1 ) = ,(10)otherwise 0,
In this paper, the velocity function v(zk) is specified as a piecewise
linear function.

view point

4.5 Controlling Cloud Motion Using
Ellipsoids
As mentioned in section 4.2, the animator can design the cloud
motion by controlling vapor probability, phase transition
probability, and cloud extinction probability. Ellipsoids are used
to do this in this paper. When wet air parcels move upward andreach the height of the dew point, clouds are gradually formed.
Ellipsoids are used to simulate the air parcels. The vapor
probability and phase transition probability are assumed to behigher at their centers than at their edges. Inversely, the cloud
extinction probability is assumed to be lower at the center since
the extinction hardly ever occurs at the center of the air parcel.
Ellipsoids also move in the direction of the wind. By controlling
ellipsoid parameters, such as sizes and positions, different kindsof clouds can be simulated. The animator specifies the regions for
ellipsoids to be generated. In our experiment, even the ellipsoids
generated using uniform random numbers result in a realistic
animation as shown in section 6.

5. RENDERING METHOD
Methods for generating realistic images are proposed in thissection. First, the density distribution of clouds is calculated bymaking use of the results of the simulation. Images are thenrendered using OpenGL. Details of the methods are described in
the following sections.

5.1 Continuous Density Distribution
Calculation
The density distribution of clouds in the real world is continuousfrom 0 to 1. The distribution obtained from the simulation,

scattered light
Figure 3: C alculation of cloud color.

4.4 Fast Simulation Using Bit Field
Manipulation Functions
Each variable can be stored in one bit since its state is either 0 or
1. This means that simulations with large numbers of cells can be
executed in a small amount of memory. The computation time isalso reduced because of the following reasons. Let us assume all
the variables are stored in an array of unsigned integers. Let m be
the bit length of the unsigned integer variable. By making use of
bit field manipulation functions of higher level language, such as
C language, transitions of m cells can be computed at the sametime. This realizes fast simulation. Most difficulties in
implementing this idea lie in transition rules concerning cloud
extinction, expressed by Eqs. 5 through 7, since random numbers
have to be generated at each cell, i.e., each bit field. This may
result in increasing the computation time. We used a look-up tablethat stores random bit sequences to save the computational cost.
See Appendix A for more details.

light frombackgroundcum ulative density
attenuation ratio

store

billboardlightm etaballFigure 4: Billboard and its texture.

however, has only two values, that is, 0 or 1. Therefore, theproposed method calculates continuous distribution by smoothing
the binary distribution, or two-valued distribution. First, thedensity, q(i,j,k,ti), is calculated at each cell (i,j,k) at each time stepti using the following equation.
q (i , j , k , t i ) =
k0

j0

1
(2t0 + 1)(2k0 + 1)(2 j0 + 1)(2i0 + 1)

t0∑t '= − t0

(11)

i0

∑ ∑ ∑ w(i', j', k ', t' )cld (i + i', j + j' , k + k ', t + t ' ),ik '= − k 0 j '= − j0 i '= −i0

where w is a weighting function and i0, j0, k0, t0 are sizes for thesmoothing. As expressed by Eq. 11, we include time as well as
space for the smoothing since the distribution is discrete in space
and time. The density at an arbitrary point, x, is then obtained as aweighted sum of a simple basis function, f. Gaussians are often
used for the basis function [29, 31, 34]. In this paper, however,
we use a field function of metaballs proposed by Wyvill et al [36].
The reason for this is as follows. A metaball has a parameter, aneffective radius, which represents its size. This means that it ismuch easier to specify the domain of influence than Gaussians
that have an infinite domain. Furthermore, the shape of Wyvill’sfield function is very similar to the Gaussians [36]. As a result, thedensity at point x is given by the following equation.
Nρ ( x, t i ) =∑ q(i, j, k , t ) f (| x − xii , j , k |) ,(12)

i , j , k∈Ω ( x , R )where R is the effective radius, Ω(x,R) is a set of cells thosecenters are within the distance R from the point x, N is the numberof elements of Ω(x,R), and xi,j,k is the coordinate corresponding to
the center of the cell (i,j,k). For the field function, f, see Appendix
B. As shown in Eq. 12, the continuous density distribution isexpressed by a set of metaballs. The user specifies the effective

fram e bu ffer billboards(initializ ed as 1.0)
A

C

D

sun
sun

read intensity

B

C

E

sun

B

A

D

C

E

projection

m etaballs

B

A

D

E

shadow texture

m ultiply attenuation
(a)

(c)

(b)

Figure 5 : Algorithm for calculating the intensity of light reaching the center of m etaballs. (a) B illboards are placed at thecenters of m etaballs and sorted based on their distances from the sun. T he fram e buffer is initialized as 1.0. (b) Billboards
are projected onto the im age plane. The colors in the fram e buffer are m ultiplied by attenuation stored in billboard tex tures.
(c) Shadow texture is obtained in the fram e buffer. Each elem ent stores the attenuation of light passing through clouds.

radius, R. The time step does not correspond to animation frames.
Therefore, the density ρ(x,t) at time t that corresponds to eachanimation frame is calculated by the linear interpolation ofdensities, ρ(x,tk) and ρ(x,tk+1), where tk < t < tk+1.

5.2 Hardware-accelerated Rendering
Using OpenGL
This section describes methods for displaying clouds as well as
shafts of light using graphics hardware. Note that the algorithm isillustrated by a pseudo-code in appendix C.

5.2.1 Rendering Clouds
Rendering of clouds is based on the splatting algorithm using
billboards. Details of the splatting method are well described in[34, 1, 20]. Therefore, let us omit the details. The basic idea for
applying it to cloud display is described here.
Fig. 3 shows the idea of calculating the color of clouds taking into
account the single scattering of light. First, the sum of thescattered light reaching from the sun on the viewing ray is
calculated. The attenuated light reaching from behind the clouds
is also calculated. The light reaching the viewpoint is the sum of
those two. Therefore, the color of a voxel depends on the scatteredcolor of the sun, the transmitted color of the sky, and theattenuation due to cloud particles. Calculation of cloud color
using splatting is as follows. First, as shown in Fig. 4, textures forbillboards are precalculated. Each element of the texture stores theattenuation ratio and cumulative density of the light passing
through the metaball (see Fig. 4). Since the attenuation is notproportional to it, the texture has to be prepared for all meatballs
when their center densities are different. However, this requires alarge amount of memory. So, the density is discretized into nq
levels and nq textures are prepared. In this paper, nq is 64. The
texture corresponding to the nearest density of each metaball ismapped onto the corresponding billboard. An image is calculated
in two steps using the texture-mapped billboards. In the first step,
the intensity of the light is calculated reaching from the sun ateach metaball. The shadows of the clouds are also calculated inthis step. In the second step, the image viewed from the viewpointis generated. The two steps are as follows.
Fig. 5 shows the idea of the first step. The basic idea is tocalculate an image viewed from the sun direction to obtain theintensity of light reaching each metaball. First, the viewpoint isplaced at the sun position and the parallel projection is assumed.
The frame buffer is initialized as 1.0. Then the billboards areplaced at the center of each metaball with their normals oriented
to the sun direction as shown in Fig. 5(a). Next, attenuation ratio
between the center of each metaball and the sun is calculated. For
example, the attenuation ratio between the metaball C and the sun
is obtained by multiplying the attenuation ratio of metaballs A, B,
and D (see Fig. 5(a)). To do this for all metaballs, the billboards
are sorted in ascending order using the distance from the sun (theorder is B-E-A-D-C in Fig. 5). Then, beginning from metaball B,
they are projected onto the image plane. The values in the frame
buffer are multiplied by their attenuation ratios that are stored in
the billboard texture (Fig. 5(b)). This can be easily done by using
blending functions of OpenGL. Then the pixel value
corresponding to the center of the metaball is read from the frame
buffer. The value obtained is the attenuation ratio between the sun
and the metaball. The color of the metaball is obtained bymultiplying the pixel value by the sunlight color. These processes
are repeated for all metaballs. After all the metaballs are processed,
the image in the frame buffer stores the attenuation ratio of thesunlight passing through the clouds (Fig. 5(c)). The image isfram e bu fferbillboardBArender objects
other than clouds

view pointEDCm etaball(a)
BA

blending
CD

E
projection(b)
Figure 6: Algorithm for generating im ages. (a) Billboards areoriented to the view point and sorted based on their distancesfrom the view point. (b) B illboards are projected onto the im age
plane. The colors in the fram e buffer are attenuated and blended w ith the colors in the billborad textures.

stored as a light map texture [1] to cast shadows on the ground.
shell n s

In the second step, the image is generated by using the color of
the metaball obtained in the first step. First, all the objects except
clouds are rendered. Next, as shown in Fig. 6(a), the billboards
are faced perpendicularly to the viewpoint and sorted indescending order based on distances from the viewpoint (theorder is E-B-D-A-C). Then they are projected onto the image
plane in back-to-front order (Fig. 6(b)). The color in the frame
buffer is blended with that of the billboard texture. The blending
process is the same as the one used in the splatting method (see[1]). That is, the colors in the frame buffer are multiplied by theattenuation ratio of the billboard texture and then the colors in the
texture are added. The process is repeated for all metaballs.

view ing
pyram idspherical shell
shell kcloudsshell 3shell 2Isshell 1

5.2.2 Rendering Shafts of LightsFig. 7 shows the idea of calculating the shafts of light. Shafts oflight are caused by particles in the atmosphere. The sunlight
passing through gaps in clouds is scattered by the particles at P in
Fig. 7. The scattered light, Is, reaching the viewpoint is recognized
as shafts of light. The scattering/absorption due to the atmosphericparticles must therefore be taken into account. The intensity oflight reaching the viewpoint is obtained by the following equation.
T∫

I = I c β (T ) + γ ( s ) I s ( s ) β ( s )ds ,

(13)

0

where Ic is the cloud color, β(s) is the attenuation ratio from theviewpoint to P due to atmospheric particles, γ(s) is the attenuation
ratio due to cloud particles from the sun to P, and Is(s) is theintensity of the light scattered at P due to atmospheric particles.
The first term in the right hand side of Eq. 13 indicates theattenuation ratio of the intensity of light from clouds. The second
term is related directly to the shafts of light. Preetham et al. take
into account the scattering of sky light as well as the sunlight torender the aerial perspective [26]. In this paper, however, we
ignore the scattering of sky light since it has little effect on shaftsof light. We also assume that the density of the atmosphericparticles decreases exponentially to the height from the ground.
Under these assumptions, the attenuation ratio, β(s), and thescattered light, Is(s), can be calculated analytically based on thepositions of the viewpoint and P [14]. Furthermore, we assume
the shafts of light are only visible under the cloud bottom, zc. Thatis, the attenuation due to clouds, γ(s), is 1.0 if zp>zc, where zp is
the z coordinate of P. In this case, γ(s) is the attenuation ratio of
the sunlight passing through the clouds. This means it has been
stored in the shadow texture obtained in the first step described in
the previous section. The shafts of light are rendered by
calculating Eq. 13 as follows.
Eq. 13 is discretized as the following equation.

∑ γ (k∆s)I (k∆s)β (k∆s)∆s ,
s

Ic
P
∆s

T

zc

view point
groundm apping shadow texture to shells

shadow texture

Figure 7: Rendering shafts of light.

Is(k∆s)β(k∆s) ∆s is stored as the colors of vertices of all thepolygons in the viewing pyramid. Finally, the second term iscomputed by rendering the shells with OpenGL’s additive
blending function. To render the shafts of light, the colors of thepolygons have to be multiplied by attenuation ratio due toclouds, γ(s). This can be easily achieved just by mapping the
shadow texture onto the polygons using OpenGL’s texture
mapping function (see Fig. 7). The function can map the texture,
multiplying the polygon’s color by the values stored in the
shadow texture. Since we have assumed γ(s) is 1 above the cloudbottom, zc, only the polygons under the cloud bottom need thismapping process.
The second step in the previous section has to be modified torender the shafts of light together with clouds. Since both of thebillboards and the shells are transparent objects, they have to berendered in back-to-front order. The procedure is as follows.
1.
2.

Calculate the colors of vertices of the polygons of the shells
in the viewing pyramid.
Repeat the following steps for k = ns, ns -1, …, 1.
2.1 Render the shell k with additive blending function. Map
the shadow texture for polygons under the cloud bottom.
2.2 Render billboards for displaying clouds between the shell
k-1 and the shell k.

6. RESULTS

ns

I = I c β (T ) +

sun

(14)

k =0

where ns and ∆s are the number of samples and the samplinginterval for the integral in Eq. 13, respectively. The attenuation,β(T), in the first term is calculated analytically by the positions of
the viewpoint and each metaball. The color of each metaball isthen attenuated by multiplying it. To calculate the second term,
spherical shells are considered as shown in Fig. 7. Their centers
are placed at viewpoint and their radii are determined so that theintervals of shells coincide to ∆s. The shells are approximated bya set of polygons to render them using OpenGL. Polygons outside
the viewing pyramid are discarded. Next, the intensity of the lightscattered at each vertex and the attenuation ratio of the path
between the viewpoint and the vertex are calculated. Then

We have made an animation to demonstrate the usefulness of our
method. Figs. 8 through 11 show sequences of images from theanimation.
Figs. 8 and 9 show the simulation of cloud formation. The numberof cells is 256x128x20. The memory for store the binary
distribution is 80 KB. In Fig. 8, clouds are formed above themountains. Ellipsoids described in section 4.5 are randomly
generated. Probabilities for vapor, phase transition and extinction
are set to 0.1, 0.001, and 0.1, respectively, at the centers of
ellipsoids. Images at every 50 steps are shown. The viewpoint isplaced above clouds. Fig. 9 shows clouds formed around themountain. In this example, ellipsoids are manually placed aroundthe top of the mountain as an initial state. Vapor probability,
phase transition probability, and extinction probability are thesame as those of Fig. 8 although they are forced to be 0 inside themountain. These examples show the formation of clouds taking
into account obstacles.
Figs. 10 and 11 show the effect of shafts of light. The cloud
evolution is simulated on 256x256x20 cells. The memory forstoring the binary distribution is 160 KB. Fig. 10 shows examples
in daytime, and Fig. 11 is in the evening. The color of the sky is
calculated by using the method proposed by Dobashi et al. [2]. In
Fig. 10, shadows on the ground and mountains are also visible.
The number of the spherical shells to calculate shafts of light is 40.
As shown in these examples, the shafts of light have a strong
visual impact to enhance the reality. In particular, the cloud color
turns red in the evening and this results in a fascinating animation
when combined with their movement.
The calculation was done on Intergraph TDZ 2000 GX1(PentiumIII 500MHz Dual). In Figs. 8 and 9, the computation
time for the simulation takes 0.3 seconds per one time step on
average. In Figs. 10 and 11, the simulation takes 0.5 seconds per
one time step. Images are rendered using the same machine. The
image sizes are 640x480. The computation time for each image in
Fig. 8 was less than 10 seconds. For Fig. 9, it took about 20
seconds. For images in Figs. 10 and 11, the computation time
ranged from 20 to 30 seconds. These results indicate that the
proposed method realizes the interactive simulation and fast
generation of realistic images.

7. CONCLUSION
In this paper, we have proposed a realistic animation method forclouds. The cloud motion is simulated using the cellular
automation. Realistic images including the shadows and the shaftsof light are generated using OpenGL. Our method has thefollowing advantages.
(1) Simulation of the cloud evolution requires only a small
amount of computation since it is executed by Boolean
operations.
(2) The memory requirement of the simulation is also small.
(3) Images can be rendered quickly by making use of graphics
hardware.
(4) Shadows of clouds and shafts of light can also be rendered.
There remain a few things to be done in the future. First, thesimulation should take into account the effects of terrain underclouds. In this case, the wind no longer blows in one direction;
hence, our simulation method must be capable of handling
multiple wind directions, or velocity fields. For faster image
generation, our rendering method should include the idea of thelevel of detail. One possible approach is to represent the voxels
hierarchically and to use coarser voxels in the distant regions fromthe viewer. This could reduce the computation time of the
splatting process since the number of metaballs could bedecreased.

APPENDIX
A. Implementation of Transition Rules
Let us assume the simulation space is nx × n y × nz = (m × n) ×
n y × nz for simplicity. n is the bit length and m is an integer value.
The necessary size of an array to store state variables is
m × n y × n z . Let us denote the arrays for hum, act, cld at time
step ti as h[m][ny][nz][ti], a[m][ny][nz][ti], c[m][ny][nz][ti],respectively. Rules expressed by Eqs. 1 and 2 can be executed in astraightforward way:h[i][j][k][ti+1]=h[i][j][k][ti] ∧ ¬ a[i][j][k][ti],(A.1)
c[i][j][k][ti+1] = c[i][j][k][ti] ∨ a[i][j][k][ti],(A.2)(i=0, …, m; j=0, …, ny; k=0, …, nz).
For the rule of Eq. 3, arrays that store bit patterns of act shifted in
one and two bit are prepared. Let us denote them as a_right1,
a_right2, a_left1, a_left2, respectively. Then Eq. 3 is calculated by
the following equation.
a[i][j][k][ti+1] = ¬ a[i][j][k][t i ] ∧ h[i][j][k][ti] ∧(a_left1[i][j][k] ∨ a_left2[i][j][k] ∨ a_right1[i][j][k] ∨a_right2[i][j][k] ∨ a[i][j+1][k][ti] ∨ a[i][j+2][k][ti] ∨a[i][j-1][k][ti] ∨ a[i][j-2][k][ti] ∨ a[i][j][k+1][ti] ∨a[i][j][k-1][ti] ∨ a[i][j][k-2][ti]),(A.3)(i=0, …, s; j=0, …, ny; k=0, …, nz).
To execute the rules of Eqs. 5 throught 7, nt bit-sequences areprepared by setting each bit to 0 or 1 using random numbers thatobey the probability i/np (i=0,..,np). The sequences are stored in
the look-up table, p[i][j] (i=0,…,np; j=1,…,nt). Using the table,
the rules are executed as follows:c[i][j][k][ti+1]= c[i][j][k][ti] ∧ p[int(pextxnp)][int(rndxnt)], (A.4)
h[i][j][k][ti+1]= c[i][j][k][ti] ∧ p[int(phumxnp)][int(rndzxnt)], (A.5)
a[i][j][k][ti+1]= c[i][j][k][ti] ∧ p[int(pactxnp)][int(rndxnt)], (A.6)
where int(x) indicates integer parts of x.

B. Field Function of Metaball as Basis
Function
The field function proposed by Wyvill et al. is used [36]. The
function is given by: 4 6 17 4 22 2
− a + a −a + 1, (r ≤ R)h( r ) =  9
,
(B.1)
9
9
0,
(r > R)
where a = r/R, r is the distance from the center of a metaball to a
calculation point, and R the effective radius of the metaball. To
use the function as basis function, we normalize it by its total
density. That is, the basis function is given by f = h(r)/c, where cis the normalizing factor. c is given by the following equation.
R748
c = 4π h(r )dr =
(B.2)
πR .
405
0

∫

C. Rendering Algorithm
Pseudo-code for the rendering clouds including the shafts of lightis given here. In each element of the billboard texture, RGB
components store the cumulative density and A component stores
the attenuation ratio of the light passing through the metaball.
PROCEDURE DisplayImage()
Place the camera at the sun position.
Set the parallel projection.
Clear screen with RGBA = (1.0, 1.0, 1.0, 1.0).
ShadeClouds()
Place the camera at the viewpoint.
Set the perspective projection.
Clear screen with background color.
RenderObject() /* Rendering objects except clouds */RenderClouds()
END PROCEDURE
PROCEDURE ShadeClouds()
Sort metaballs in ascending order from the sun.
glDisable(GL_DEPTH_TEST)
glBlendFunc(GL_ZERO, GL_SRC_ALPHA)
glEnable(GL_BLEND)
FOR k = each metaball DO

Place the billboard at the center of metaball k.
Rotate the billboard so that its normal is oriented to the sun.
Set the billboard color as RGBA = (1.0, 1.0, 1.0, 1.0).
Map the billboard texture with GL_MODULATE.
Render the billboard.
Read the pixel value corresponding to the center of metaball k.
Multiply the pixel value by the sunlight color.
Store the color into an array C[k] as the color of the billboard.
END FORStore the image (T) in the frame buffer as a light map texture.
END PROCEDURE
PROCEDURE RenderClouds()
Sort metaballs in descending order from the viewpoint.
Store distances between metaballls and the viewpoint in an array
D.
glDisable(GL_DEPTH_TEST)
glEnable(GL_BLEND)
k = number of shells
WHILE D[0] < distance to shell k DO
RenderShell( k )k=k–1
END WHILE
glBlendFunc(GL_ONE, GL_SRC_ALPHA)
FOR n = each metaball DO
IF D[n] < distance to shell k DO
Render shell(k)k=k–1
glBlendFunc(GL_ONE, GL_SRC_ALPHA)
END IF
Place the billboard at the center of the corresponding metaball n.
Rotate the billboard so that its normal is oriented to theviewpoint.
Set the billboard color as C[n].
Map the billboard texture.
Render the billboard with the blending function.
END FOR
WHILE k > 0 DO
RenderShell(k)k=k–1
END WHILE
END PROCEDURE
PROCEDURE RenderShell(k)Calculate colors of vertices of shell k.
Map the light map texture T with GL_MODULATE.
Render shell k.
glBlendFunc(GL_ONE, GL_ONE)END PROCEDURE

REFERENCES[1] D. Blythe, “Advanced Graphics Programming Techniques
Using OpenGL,” Course Note #29 of SIGGRAPH 99, 1999.
[2] Y. Dobashi, T. Nishita, K. Kaneda, H. Yamashita, “A Fast
Display Method of Sky Color Using Basis Functions,” The
Journal of Visualization and Computer Graphics, Vol. 8, No.
2, 1997, pp. 115-127.
[3] Y. Dobashi, T. Nishita, T. Okita, “Animation of Clouds
Using Cellular Automaton,” Proc. of Computer Graphics
and Imaging’98, 1998, pp. 251-256.
[4] Y. Dobashi, T. Nishita, H. Yamashita, T. Okita, “Using
Metaballs to Modeling and Animate Clouds from Satellite
Images,” The Visual Computer, Vol. 15, No. 9, 1998, pp.
471-482.
[5] D. S. Ebert, R. E. Parent, “Rendering and Animation of
Gaseous Phenomena by Combining Fast Volume andScanline A-Buffer Techniques,” Computer Graphics, Vol. 24,
No. 4, 1990, pp. 357-366.
[6] D. S. Ebert, W. E. Carlson, R. E. Parent, “Solid Spaces andInverse Particle Systems for Controlling the Animation ofGases and Fluids,” The Visual Computer, 10, 1990, pp. 471483.
[7] D. S. Ebert, “Volumetric Modeling with Implicit Functions:
A Cloud is Born,” Visual Proc. of SIGGRAPH’97, 1997, pp.
147.
[8] D. S. Ebert, “Simulating Nature: From Theory to
Application,” Course Note #26 of SIGGRAPH 99, 1999, pp.
5.1-5.52.
[9] N. Foster, D. Metaxas, “Modeling the Motion of a Hot,
Turbulent Gas,” Proc. of SIGGRAPH’97, 1997, pp. 181-188.
[10] G.Y. Gardner, “Visual Simulation of Clouds,” Computer
Graphics, Vol.19, No. 3, 1985, pp. 279-303.
[11] M. Inakage, “Volume Tracing of Atmospheric
Environments, ” The Visual Computer, 7, 1991, pp. 104-113.
[12] H. W. Jansen, P. H. Christensen, “Efficient Simulation ofLight Transport in Scenes with Participating Media using
Photon Maps,” Proc. of SIGGRAPH’98, 1998, pp. 311-320.
[13] J. T. Kajiya, B. P. V. Herzen, “Ray Tracing VolumeDensities,” Computer Graphics, 1984, Vol. 18, No. 3, pp.
165-174.
[14] K. Kaneda, T. Okamoto, E. Nakamae, T. Nishita,“Photorealistic Image Synthesis for Outdoor Scenery under
Various Atmospheric Conditions,” The Visual Computer,7(5&6), 1991, pp. 247-258.
[15] T. Kikuchi, K. Muraoka, and N. Chiba, “Visual Simulationof Cumulonimbus Clouds,” The Journal of The Institute of
Image Electronics and Electronics Engineers of Japan, Vol.
27, No. 4, 1998, pp. 317-326 (in Japanese).
[16] N. Max, “Light Diffusion through Clouds and Haze,”
Graphics and Image Processing, Vol. 13, No. 3, 1986, pp.
280-292.
[17] N. Max, “Atmospheric Illumination and Shadows,”
Computer Graphics, Vol. 20, No. 4, 1986, pp. 117-124.
[18] N. Max, R. Crawfis, D. Williams, “Visualizing Wind
Velocities by Advecting Cloud Textures,” Proc. ofVisualization’92, 1992, pp. 179-183.
[19] N. Max, “Efficient Light Propagation for Multiple
Anisotropic Volume Scattering,” Proc. of the Fifth
Eurographics Workshop on Rendering, 1994, pp. 87-104.
[20] K. Meuller, N. Shareef, J. Huang, R. Crawfis, “Hight-Quality
Splatting on Rectilinear Grids with Efficient Culling ofOccluded Voxels,” IEEE Trans. on Visualization andComputer Graphics, Vol. 5, No. 2, 1999, pp. 116-134.
[21] K. Nagel, E. Raschke, “Self-Organizing Criticality in Cloud
Formation?,” Physica A, 182, 1992, pp. 519-531.
[22] F. Neyret, “Qualitative Simulation of Convective Clouds
Formation and Evolution,” Proc of Eurographics Computer
Animation and Simulation Workshop’97, 1997, pp. 113-124.
[23] T. Nishita, Y. Miyawaki, E. Nakamae, “A Shading Model for
Atmospheric Scattering Considering Distribution of Light
Sources,” Computer Graphics, Vol. 21, No. 4, 1987, pp.
303-310.
[24] T. Nishita, E. Nakamae, “Method of Displaying Optical
Effects within Water using Accumulation Buffer,” Proc. of
SIGGRAPH’94, 1994, pp. 373-379.
[25] T. Nishita, Y. Dobashi, E. Nakamae, “Display of Clouds
Taking into Account Multiple Anisotropic Scattering and
Sky Light,” Proc. of SIGGRAPH’96, 1996, pp. 379-386.
[26] A. J. Preetham, P. Shirley, B. Smits, “A Practical Analytic
Model for Daylight,” Proc. of SIGGRAPH’99, 1999, pp. 91100.

[27] H. E. Rushmeier, K. E. Torrance, “The Zonal Method for
Calculating Light Intensities in The Presence of a
Participating Medium,” Computer Graphics, Vol. 21, No. 4,
1987, pp. 293-302.
[28] G. Sakas, M. Gerth, “Sampling and Anti-Aliasing of Discrete
3-D
VolumeDensity
Textures,”
Proc.
of
EUROGRAPHICS'91, 1991, pp. 87-102.
[29] J. Stam, E. Fiume, “Turbulent Wind Fields for Gaseous
Phenomena,” Proc. of SIGGRAPH'93, 1993, pp. 369-376.
[30] J. Stam, “Stochastic Rendering of Density Fields,” Proc. ofGraphics Interface’94, 1994, pp. 51-58.
[31] J. Stam, E. Fiume, “Dipicting Fire and Other Gaseous
Phenomena Using Diffusion Processes,” Proc. of
SIGGRAPH'95, 1995, pp. 129-136.

(a) ti = 0(b) ti = 50

[32] J. Stam, “Stable Fluids,” Proc. of SIGGRAPH'99, 1999, pp.
121-128.
[33] R. Voss, “Fourier Synthesis of Gaussian Fractals: 1/f noises,
landscapes, and flakes,” SIGGRAPH’83: Tutorial on State of
the Art Image Synthesis, 10, 1983.
[34] L. Westover, “Footprint Evaluation for Volume Rendering,”
Computer Graphics, Vol. 24, No. 4, 1990, pp. 367-376.
[35] S. Wolfram, “Cellular automata as models of complexity,”
Nature, Vol. 311, No. 4, 1984, pp. 419-424.
[36] G. Wyvill, A. Trotman, “Ray-Tracing Soft Objects,” Proc. of
CG International, 1990, pp. 439-475.

(c) ti = 100

Figure 8: Simulation of cloud formation.

(a) ti = 0(b) ti = 200(c) ti = 400(d) ti = 600
Figure 9: Cloud formation around mountains.

(d) ti = 150

(a) ti = 100(b) ti = 200(c) ti = 300(d) ti = 400
Figure 10: Examples of shafts of light (daytime).

(a) ti = 100(b) ti = 150(c) ti = 200(d) ti = 250
Figure 11: Examples of shafts of light (evening).