Graphics Hardware (2002), pp. 1–10
Thomas Ertl, Wolfgang Heidrich, and Michael Doggett (Editors)

Interactive Rendering of Atmospheric Scattering Effects
Using Graphics Hardware
Yoshinori Dobashi† ,Tsuyoshi Yamamoto† ,† Hokkaido University

Tomoyuki Nishita†††† The University of Tokyo

Abstract
To create realistic images using computer graphics, an important element to consider is atmospheric scattering,that is, the phenomenon by which light is scattered by small particles in the air. This effect is the cause of the lightbeams produced by spotlights, shafts of light, foggy scenes, the bluish appearance of the earth’s atmosphere, andso on. This paper proposes a fast method for rendering the atmospheric scattering effects based on actual physicalphenomena. In the proposed method, look-up tables are prepared to store the intensities of the scattered light, andthese are then used as textures. Realistic images are then created at interactive rates by making use of graphics
hardware.
Categories and Subject Descriptors (according to ACM CCS): I.3.1 [Computer Graphics]: Hardware Architecture—
Graphics processors; I.3.3 [Computer Graphics]: Picture/Image Generation—Bitmap and frame buffer operations;I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism—Color, Shading, Shadowing and Texture

1. Introduction
In the field of computer graphics, realistic image synthesis
is one of the most important research subjects. When there
are small particles such as dust in the air, light is scatteredand absorbed by them. This is called atmospheric scattering. Several different methods have been developed to simulate this. These methods make it possible to display various effects, such as the light beams caused by spotlights,
shafts of light through clouds, foggy scenes, the sky, and theearth viewed from space1, 2, 3, 4, 5 . These research results are
now widely used for a variety of applications, such as stage
lighting designs, driving simulations, and space/flight simulations. However, the computational cost is very expensive,
since the intensities of the scattered light must be integrated
along with each viewing ray. A fast and precise rendering
method must be developed, especially for drive simulatorsand space/flight simulators.
This paper proposes a method for rendering these atmospheric scattering effects at interactive rates. Our method canhandle a wide range of atmospheric scattering effects, including the scattering observed in a smoky room, scattering
due to fog, and scattering due to the earth’s atmosphere. In
the case of a smoky room, the light beams can be producedc The Eurographics Association 2002.


by point light sources and/or sunlight, and our technique cantake this into account. The scattering due to the earth’s atmosphere is indispensable for rendering the sky or the earth
viewed from space, and we have to take into account the following phenomena. That is, the sky is bluish at daytime andreddish in the evening, and the earth as viewed from spacebecomes bluish. For the earth viewed from space, clouds
also play an important role in creating realistic images. Our
method realizes a fast display of these effects. For the earth’s
atmosphere, we propose a method for displaying the atmospheric effects viewed not only from the ground, but also as
seen from space. Note that the rendering of the atmosphereviewed from the ground is the same as rendering the sky.
Furthermore, a method for rendering clouds as viewed fromspace is developed to create realistic images of the earth.
The rendering processes are accelerated by using graphics hardware. The intensities of the scattered light arecomputed in a pre-process and stored in look-up tables.
These tables are used as textures. Next, the images are created by hardware-accelerated volume rendering techniques6 .
Dobashi et al. have already proposed a method for rendering
light beams by applying volume rendering techniques5, 7, 8 .
This paper discusses the problems of the approach taken byDobashi, Yamamoto, Nishita / Interactive Rendering of Atmospheric Effects

these previous methods and addresses them. The differences
between this paper and those previous ones are: (1) the accurate sampling of shadows and the use of texture for efficientevaluation of the intensity of scattered light and (2) an extension to the rendering of the sky and the earth as viewed
from space. The proposed method improves the image quality significantly without increasing the rendering time. Our
method realizes a fast rendering of the light beams producedby point and/or infinite light sources. The proposed method
also realizes fast rendering of the earth’s atmosphere.

sampling plane I light sourcelscreenviewpointts
objectIv
Ib
α P
scattering lightT

Figure 1: Calculation of atmospheric scattering of light.
3. Overview of Our Method

2. Previous Work
Many methods have been developed for generating imagesby taking into account the scattering and absorption of light
due to small particles in the air1, 2, 3, 4, 9, 10, 11 . Although thesemethods can create realistic images, most of these employ a
ray-tracing algorithm. Therefore, it takes several minutes tocreate a single image.
To address the problem, hardware-accelerated methods
have been developed. Since the rendering of atmospheric effects can be considered as a kind of volume rendering, it isnatural to apply the techniques developed for volume rendering. Many of them use hardware texture mapping functions to render volumes defined on voxels12, 13, 14, 15, 16, 17 . In
these methods, the particle density is stored at each voxel.
Then the volume data is transferred to texture memories
as a 2D/3D texture. These methods have already been applied to rendering smoke and clouds. The volume data,
however, consume a large amount of texture memory. Furthermore, for rendering light beams, sharp edges exist in
the atmosphere due to the luminous intensity of the lightsource and shadows of objects. It is very difficult to capture the edges by using a voxel-based approach. To address
the problem, methods for rendering shafts of light have been
proposed5, 7, 8, 18, 19 . However, these still have problems insampling shadows, so artifacts appear in the resulting images. Keller et al. presented a general and elegant solution to
reduce the artifacts due to the sampling errors20 . However,
our method achieves a faster image generation for the specific problem treated in this paper, that is, the atmosphericscattering. Mech proposed an altenative method for displaying shafts of light by rendering boundaries of the light beams
represented by polygons21 . However, the method cannot display the shadows of objects in the atmosphere. For the purpose of scientific visualization, hardware-accelerated methods for rendering the earth as viewed from space have been
proposed22, 23, 16 . In these methods, however, realistic images
are not generated since the physical properties of the lightscattering due to the earth’s atmosphere are not taken into
account. Our method takes these into account and accelerates the computation of the scattering of light by using twodimensional textures mapping functions.

For rendering light beams, we assume the density of atmospheric particles is constant. For the earth’s atmosphere, theshape of the earth is assumed to be a sphere and the density of particles decreases exponentially according to theirheight above the ground. Multiple scattering is important foratmospheric scattering but is time-consuming. Our method
approximates multiple scattering as a constant ambient term.
Our method utilizes the hardware-accelerated volume rendering technique5, 6, 7, 8 . As shown in Fig. 1, planes are placedin front of the screen, and they are subdivided into a mesh.
We call the planes sampling planes. The atmospheric scattering effects are rendered by drawing the sampling planes.
The intensities of the light scattered at points on the planes
are obtained by using the look-up tables computed in a preprocess. The look-up tables are used as textures and aremapped onto the sampling planes. The intensity of scattered
light at a point depends on the light source position, theviewing direction, and so on. Therefore high (three- to five) dimensional look-up tables are required, and it requires alarge amount of memory. To address this, we propose themethod to store the intensity in two-dimensional textures.
Furthermore, sub-planes are introduced for the precise sampling of object shadows and the luminous intensity distribution of the light source.
The above method is extended to the rendering of the sky.
To represent the earth as viewed from space, however, it is
not appropriate to use the method described in Fig. 1. In
this case, the atmosphere is considered as a very thin spherical layer, and therefore a large number of planes are required for accurate sampling of the intensities of the scattered light, which results in increased computation time.
Therefore, sample spheres are generated around the earth.
The centers of the sample spheres coincide with the earth’s
center, and we then use the sample spheres instead of thesampling planes. Furthermore, clouds viewed from space arerendered to enhance the reality.
4. Rendering of Light Beams
Applying the method in Fig. 1 to rendering light beams in astraightforward way causes serious problems due to the limitations of the available graphics hardware. In this section,c The Eurographics Association 2002.
Dobashi, Yamamoto, Nishita / Interactive Rendering of Atmospheric Effects

the problems are firstly described in detail and then solutions to them are proposed. In the following, let us assume
that the number of light sources is one for the sake of simplicity. When there are multiple light sources, the intensityof the light reaching the viewpoint is obtained by summing
the contribution from each light source.
4.1. Problems in Rendering Light Beams
As shown in Fig. 1, the intensity of light reaching the viewpoint, Iv , is computed by summing the intensity of light reflected from an object and the intensities of light scatteredby atmospheric particles. Note that Fig. 1 assumes a point
light source. Since we assume the density of the particles isa constant ρc , Iv is given by:
Iv (λ) = Ib (λ)g(T, λ) + Is (λ) + Ia (λ),
Is (λ) =T0Il (t, λ)H(t)g(t, λ)dt,g(t, λ) = ρc F(α, λ) exp(−βρc (s + t))/s2 ,

(1)
(2)
(3)

where Ib is the intensity of the object, λ is a wavelength, T
is the distance between the viewpoint and the object, Il is theintensity of light emitted from the light source toward point
P in Fig. 1 (Il is obtained by using the luminous intensity
distribution of the light source, or light map25 ). Without lossof generality, we can assume Il is always less than 1.0. s andt are the distances from point P to the light source and theviewpoint respectively, Ia is the ambient term that approximates the multiple scattering, F is a phase function of theparticles, and α is the phase angle (see Fig. 1). H is a visibility function that returns the answer 1 if the light sourceis visible from point P, and otherwise returns 0 if it is not.
g is a function of the attenuation ratio of the intensity oflight reaching the viewpoint after it is scattered at point P,
and β is the extinction coefficient. Note that g is given byg(t, λ) = γ exp(−βρct) for an infinite light source, where γ
is a constant that represents the attenuation ratio of light. In
Eq. 1, the second term, Is , is directly related to the renderingof the light beams. This subsection discusses the problems incomputing Is . One simple approach to rendering light beams
is to apply the traditional texture-based volume rendering
technique5, 6, 7, 8, 17 . First, the intensities of the scattered lightare computed at the vertices of the mesh representing eachsampling plane and stored as colors of the vertices. Next, theplanes are drawn and the intensities are accumulated in theframe buffer by using color-blending functions. This simple
approach, however, causes the following problems.
In general, the values of the intensity distribution of thelight source and the visibility function H change rapidly.
Therefore, a lot of planes are required to achieve accurate
sampling of these functions. However, using many planes results in poor images due to color quantization errors. Colors
assigned at the vertices are stored with a limited precision (8bit or 12-bit on high-end machines) in most available graphics hardware. Therefore, quantization errors are accumulated
c The Eurographics Association 2002.


in proportion to the number of the sapling planes used. In our
experiment, we found that the attenuation functions, g, areprone to be influenced by the quantization errors. The reason for this is that these functions are very small and quantized to zero in most cases. This problem may be resolvedin the future when the graphics hardware with much higher
bit depth becomes available. Even so, increasing the number of the planes is not desirable, since the rendering timeis proportional to the total number. In addition, the sampling
planes have to be subdivided into a finer mesh for the accurate evaluation of the intensity of the scattered light. This
also results in significant increase in the rendering time.
4.2. High Quality Rendering
There are two key ideas in our approach. One is to use texture to evaluate the intensity of the scattered light at points
on the sampling planes, and the other is to use sub-planes forthe accurate sampling of shadows and the luminous intensity
distribution of the light source. We first investigate Eq. 2 indetail before explaining the proposed method.
Is is computed numerically by subdividing the integral of
Eq. 2 into n sub-segments. The sub-segments are of equal
length, that is, ∆t. The kth sub-segment is expressed as[tk ,tk + ∆t]. Eq. 2 is rewritten as follows.
Is (λ) =∆Is (tk , λ) =n∑ ∆Is (tk , λ),k=1
tk +∆ttk

Il (t, λ)H(t)g(t, λ)dt.

(4)
(5)

The functions on the right-hand side of Eq. 5 are classified
into two groups. One group consists of functions that requirea high sampling ratio (i.e., Il and H). The rest of them are included in another group. Eq. 5 is approximated by the product of these groups, that is,∆Is (tk , λ) ≈ fh (tk , λ) × fl (tk , λ),(6)fh (tk , λ) = ttkk +∆t Il (t, λ)H(t)dt.
 tk +∆tfl (tk , λ) = tkg(t, λ)dt(7)

Eq. 5 is approximated very well by Eq. 6 since fl changes
smoothly. Moreover, fh does not include g, which is sensitive to quantization errors. In particular, the errors do notoccur in evaluating the visibility function H since it is a twovalued function (0 or 1). So, fh can be sampled at a shorter
interval than ∆t. On the other hand, fl includes g. fl is evaluated precisely by using texture mapping. This prevents thequantization errors from accumulating. In the following, themethods for evaluating fl , fh are proposed.
4.2.1. Textures for Intensity of Scattering
This subsection describes the creation of textures for fl . Let
us explain our algorithm for the point light source only. ItDobashi, Yamamoto, Nishita / Interactive Rendering of Atmospheric Effects

bounding volume
light sourceof light beam

sampling
sampling
plane kplane k+1
light source
∆t
t
U
s

viewpointtk

αP(u,v)

objectshadow
sampling planeP’( u’,v) V

d0

viewpointQviewing ray

sub-plane

Figure 2: Creating the textures for light beams.
is straightforward to extend the method for an infinite lightsource. fl is given by the following equation.
fl (tk , λ) =tk +∆ttkexp(−βρc (s + t))
ρc F(α, λ)
dt.
s2

(8)

The basic idea of our method is to store the integrated values expressed by the above equation in a texture. This idea issimilar to the pre-integrated volume rendering method proposed by Engel et al24 . To store the integrated value in Eq.
8, high-dimensional tables are required since fl is a functionof the light source position, the viewpoint, and the directionof the viewing ray. However, in the case where there is a
uniform density, a simpler expression exists. As shown in
Fig. 2, let us assume a local coordinate system, UV , with
its origin at the position of the light source. The U axis isparallel to the viewing ray and it passes through the lightsource position. The V axis is perpendicular to the U axis,
and it passes through a point Q on the viewing ray that isclosest to the light source. Using this coordinate system, let
us denote the coordinate of point P on sampling plane k as(u, v). The distance from the viewpoint to P is tk (see Fig. 2).
The coordinate of point P is denoted by (u , v). t is the distance between the viewpoint and point P . Then, the following relationships are obtained: t = u −u+tk , cos α = −u /s,
s2 = u2 + v2k . The exponential term in Eq. (8) is rewritten as:exp(−βρc (s + t)) = exp(−βρc ( u2 + v2 + u − u + tk ))= exp(−βρc ( u2 + v2 + u − u))× exp(−βρctk ).
As a result, fl is expressed by:
fl (u, v, λ) = ck q(u, v, λ),(9)where,
ck = exp(−βρctk ),q(u, v, λ) =ρc F(cos−1 (−u / u2 + v2 ), λ)u√exp(−βρc ( u2 + v2 + u − u)) du×u2 + v2u+∆t

ck depends only on the distance from the viewpoint to thepoints on the sampling plane. ck is evaluated by assigning it
as the color of vertices of the mesh representing the plane.

Figure 3: Rendering light beams using sub-planes.
q(u, v, λ) is computed in a pre-process step for different values of u and v. The sampling interval ∆t is determined by
the intersection of the ray with the sampling planes, and so∆t is different from pixel to pixel theoretically. In our implementation, however, ∆t is set to a constant, that is, the interval of the sampling planes. This makes it possible to store
q(u, v, λ) in a two-dimensional look-up table. One advantageof this expression is that fl is evaluated precisely, since thetable stores the integrated values. The tables are used as the
texture. λ is sampled at wavelengths that correspond to the
RGB components. The remaining issues are how to determine the value of ∆t (i.e. the interval of the sampling planes)
and the number of the planes.
In most graphics hardware, values stored in the texture
are clamped from zero to a certain value Imax (usually, Imax
is 1.0). So, we determine ∆t so that the values in the table ofq do not exceed Imax . In addition, a maximum limit, ∆max ,is specified since the use of a too large value of ∆t makes
the approximation of Eq. 6 poor. First, ∆t is set to a smallvalue, ∆min , and the table of q is computed. We extract the
maximum value from the table. Then, ∆t is incremented by∆min , and the maximum value is computed again. This process is repeated until ∆t reaches ∆max or the maximum value
exceeds Imax . The resulting value is used as ∆t.
The number of the sampling planes is automatically computed at every frame. We first define a bounding volume of
the light beam (see Fig. 3). We use a pyramid as the bounding volume. The apex of the pyramid is at the position of
the light source and the bottom of the pyramid is placed at
a distance d0 from the light source. d0 is specified by theuser. Then, the nearest and farthest points of the pyramidfrom the viewpoint are computed. The planes are generated
between these two. When the distances from the viewpoint
to these points are tmin and tmax , the number of the planes is(tmax − tmin )/∆t + 1.
4.2.2. Precise Sampling using Sub-planes
This subsection describes the method for computing fh . As
shown in Fig. 3, the basic concept is to insert sub-planes between the sampling planes (or the main planes). Each subplane is represented by a single quadrilateral. The sub-planesc The Eurographics Association 2002.
Dobashi, Yamamoto, Nishita / Interactive Rendering of Atmospheric Effects
fl0.05

number of
sub-planes
12

0.04

10

0.03

8
6

0.02

4

0.01
27.1

24.7

22.3

19.9

17.5

15.1

12.7

7.9

10.3

5.5

3.1

2
0.7

0

0

distance from viewpoint [m]

(a)

(b)

(c)

(d)

Figure 4: Examples of artifacts. Bottom column shows theclose-up around the sphere. (a) shows sampling errors of
shadows, (b) quantization errors of intensity of scattered
light, (c) artifacts are removed by our method, (d) a raytraced solution.
do not need to be represented by the mesh. As shown in
Fig. 3, the appropriate number of the sub-planes is determined adaptively depending on the intensity of the scattered
light. Let us assume that mk sub-planes are generated between plane k and k + 1. Then, fh is given by:fh (tk , λ) =1 mk −1∑ H(tk + j∆r)Il (tk + j∆r, λ), (10)mk j=0(∆r = ∆t/(mk − 1)).
The computation of Eq. 10 is accelerated by the graphics
hardware. First, the visibility term H is obtained by computing the shadows on the sub-planes by a hardware shadow
mapping technique25 . Next, Il is obtained by projecting alight map texture onto the sub-planes by using a projective
texture mapping technique25 . After these processes, the subplanes are rendered and their colors are accumulated into theframe buffer. Note that only a single quadrilateral is rendered
for each sub-plane. The pixel values of the resulting image
are equal to fh . This image is stored as a texture. When creating the images, this texture and the texture for fl describedin Section 4.2.2 are mapped onto the sampling planes using
the multiplicative texture mapping function.
The number of the sub-planes, mk , is determined as follows. From Eqs. (6) and (10), the contribution, ξ j , of subplane j to the pixel intensity is estimated from the following
equation.
ξ j = H(tk + j∆r)Il (tk + j∆r) fl (tk , λ)/mk ≤ fl (tk , λ)/mk .
The inequality in the above equation is due to the fact thatH ≤ 1 and Il < 1. If fl (tk , λ)/mk is less than a small value
ε, we cannot expect sub-plane j to significantly improve thequality of the resulting image. Therefore, it is natural to select the number mk so that fl (tk , λ)/mk is equal to or greaterthan ε. To implement this idea, we first shoot a ray toward
the center of the screen and compute the distance, dk , fromc The Eurographics Association 2002.


Figure 5: The distribution of fl and the number of subplanes. The curved line shows the distribution of fl and thebar chart shows the number of sub-planes for each subsegment.
the viewpoint to the intersection of the ray with each sampling plane. Then, the minimum number that satisfies thefollowing condition is used as mk .
max fl (dk , λ)/mk ≥ ε,λwhere ε is a threshold specified by the user. fl (dk , λ) is obtained by the table for fl described in the previous section.
4.3. Experimental Result
In this section, the validity of the proposed method is verified by applying it to a simple example where a sphere isilluminated by a single spotlight as shown in Fig. 4. In Fig.
4, the images in the bottom column show the close-up viewsof the images in the top column. The spotlight is placed at
20 [m] high. Figs. 4(a) and (b) are generated using the previous method7 . We used a desktop PC (Athlon 1.7 GHz) witha NVIDIA GeForce3. The size of each image is 300×450.
In Fig. 4(a), the number of sampling planes is 40. It took
0.13 [sec.] to render this image. It is obvious that the shadows are not sampled precisely. In Fig. 4(b), the number ofthe planes is increased, 160. However, pseudo-contours now
appear due to the quantization errors. The rendering time isalso increased to 0.50 [sec]. In these images, each plane isrepresented by a 60×90 mesh to sample the intensity of thescattered light accurately. This also increases the rendering
time. Fig. 4(c) is generated by the proposed method. Clearly,
the quality of the image is improved. The rendering time isstill fast, 0.08 [sec.]. In this case, the 10×15 mesh is sufficient for each sampling plane since the intensity of the scattered light is evaluated accurately on a pixel basis by a hardware texture mapping function. For comparison, Fig. 4(d) is
generated by using the ray tracing method1 . It took 29 seconds. The average error in intensities between Figs. 4(c) and(d) is 1.5 [%]. Fig. 5 shows the distribution of the function
fl on the ray passing through the center of the screen. The
horizontal axis indicates the distance from the viewpoint.
The sampling planes are placed at an interval of 1.2[m]. The

Dobashi, Yamamoto, Nishita / Interactive Rendering of Atmospheric Effects

threshold ε to insert sub-planes (see Section 4.2.2) is set to
1 [%] of the maximum value Imax that can be handled by thegraphics hardware (in this experiment, ε = 0.01 since Imax =1.0). The bar chart in Fig. 5 shows the number of sub-planes.
In total, the numbers of the main planes and sub-planes are
23 and 174, respectively. As shown in Fig. 5, the number ofthe sub-planes is increased where the value of fl is big, thatis, the intensity of the scattered light is strong.
Finally, we compare our method with the method proposed by Keller et al20 . This method can also reduce the artifacts due to the sampling error. The basic idea of Keller’s
method is to create multiple images with different positionsof the planes and accumulate them into the final image. We
accumulate four images with forty planes to render the same
scene as shown in Fig. 4. These numbers are chosen so thatthe sampling interval is approximately the same as the oneused by our method. Note that the sizes of the four images
are same as the one used by our method (i.e. 300×450),
to make the number of samples per each pixel equal. Eachplane is represented by the 15×10 mesh. Although Keller’s
method created the image with the same quality, it took 0.12
[sec.]. Our method is 50 % faster. The reason is as follows.
Keller’s method needs to draw 40 planes represented by themesh for each of four images. In total, 160 planes represented by the mesh are rendered. In our method, 23 planes
represented by the mesh and 174 sub-planes represented
by the quadrilateral are rendered. Therefore, the number ofpolygons rendered is smaller in our method. Furthermore,
our method does not generate sub-planes where the intensityof the scattered light is small. This also results in decrease in
the rendering time.
The rendering speed of our method depends on the ratio ofgeometric throughput and pixel fill rate. Especially, for creating high resolution images, our method becomes slow. Thisis because the fill rate becomes the bottleneck. This problem is resolved by combining our method with the Keller’s
method. That is, the high resolution image is created by compositing multiple low resolution images with the interleaved
sampling technique20 .
5. Rendering the Earth’s AtmosphereThis section describes the method for rendering the earth’s
atmosphere. For the earth’s atmosphere, however, we can ignore the shadows of objects. Therefore, we focus on how tocreate textures for the intensity of the scattered light. In thefollowing, we will first describe the principle of the intensity calculation for the atmospheric scattering, and then our
methods for rendering the earth’s atmosphere are proposed.
5.1. Intensity Calculation for Atmospheric Scattering
The earth’s atmosphere consists of air molecules andaerosols. Aerosols are particles of dusts floating in the
air. The scattering properties of the air molecules andsun light IsunPsθθv sun s
P’∆tt
top ofP
atmosphere Thviewpoint
earth

Figure 6: Rendering of atmosphere viewed from ground.
sun light Isunsample sphereviewpoint
θv
Pk+1

n

T
θsun

θg

Pk

hk

k+1

top ofatmosphere

Ie

k

earth

2

S

Pg

1

Figure 7: Rendering of atmosphere viewed from space.
the aerosols obey Rayleigh and Mie scattering theories,
respectively4 . Figs. 6 and 7 show the principles of the intensity calculation. The viewpoint is on the ground in Fig. 6
or is in space in Fig. 7. The intensity of light reaching the
viewpoint Iv is given by Eq. 1110 .
Iv (λ) = Ie (λ)gv (T, λ)+ Isun (λ)

T
0R(α,t, λ)gl (s, λ)gv (t, λ)dt, (11)where T is the distance to the top of the atmosphere (see Fig.
6) or to the ground (see Fig. 7). In Eq. 11, subscripts "r" and"m" represent the air molecules (Rayleigh scattering) and theaerosols (Mie scattering), respectively. Ie , R, gv , and gl are
given by the following equations.
Ie (λ) = Isun (λ)r(λ)gl (S, λ) cosθg ,R(α,t, λ) = Kr (λ)ρr (t)Fr (α) + Km (λ)ρm (t)Fm (α),
gv (t, λ) = exp(−τ(t, λ)), gl (s, λ) = exp(−τ(s, λ)),where Ie is the intensity of light reflected on the ground, Isunis the intensity of the sunlight, r is the reflectivity at point Pg
on the ground, θg is the angle between the normal vector at
Pg and the light vector, and S is the distance between Pg andthe top of the atmosphere (see Fig. 7). Ie is 0 if the viewing
ray does not hit the earth. gv is a function of the attenuation
ratio between the viewpoint and point P, and gl is a functionof the attenuation ratio of the sun light reaching point P. ρrand ρm are the densities of the atmospheric particles, and t is
the optical length. Fr and Fm are phase functions. For moredetail, see4 .
c The Eurographics Association 2002.
Dobashi, Yamamoto, Nishita / Interactive Rendering of Atmospheric Effects
sun light
Isun

5.2. Rendering Sky
In a similar way to the method for the light beams, the integral term in Eq. 11 is evaluated by subdividing it into nsub-segments. The length of each sub-segment is ∆t.
First, let us discuss the attenuation ratio of the sunlightreaching point P on the viewing ray (see Fig. 6). Let us
assume that the distance between the viewpoint and point
P is t, the height of point P is h, and the angle between
the sunlight direction and the vertical direction at P isθsun . Then, the attenuation ratio is expressed as a function of the height h and the angle θsun since the densityof the atmospheric particles depends only on their heightfrom the ground (see Appendix A.1). That is, gl (s, λ) =gl (h, θsun , λ)(= exp(−t(PPs , λ))). Note that gl (h, θsun , λ) =0 if the sun as viewed from point P is behind the earth.
gl (h, θsun , λ) is computed in a pre-process and stored in a
texture.
Next, let us consider the ray that passes through subsegment k, [tk ,tk + ∆t]. Eq. 11 is rewritten as follows.
n−1

Iv (λ) = Isun (λ)gl (0, θg , λ) cosθg × ∏ ∆gv (t j , λ)
j=0

n−1

+ Isun (λ) ∑ gl (hk , θk , λ)[Fr (α, λ)∆Ir (tk , λ)
k=0


k−1
+ Fm (α, λ)∆Im (tk , λ)] ∏ ∆gv (t j , λ) ,

(12)

j=0

where ∆Ir (tk , λ) and ∆Im (tk , λ) are the intensities of scattering and ∆gv (t j , λ) is the attenuation ratio of the segment[t j ,t j + ∆t]. These are given by the following equations.
∆Ir (tk , λ) =∆Im (tk , λ) =tk +∆ttk
tk +∆ttk

Kr (λ)ρr (t) exp(−τ(t − tk , λ))dt,
Km (λ)ρm (t) exp(−τ(t − tk , λ))dt,∆gv (tk , λ) = exp(−τ(tk → tk + ∆t, λ)),where τ(tk → tk + ∆t, λ) indicates the optical length of thesegment [tk ,tk + ∆t]. ∆Ir , ∆Im , and ∆gv can also be stored in
the look-up tables. Let us consider a ray that passes throughsegment PP as shown in Fig. 6. The angle between the vertical direction at P and the segment PP is θv . Then, in asimilar way to gl , ∆Ir , ∆Im , and ∆gv are expressed as functions of the height h and the angle θv (see appendix A.2),
that is, ∆Ir (t, λ) = ∆Ir (h, θv , λ), ∆Im (t, λ) = ∆Im (h, θv , λ),∆gv = ∆gv (h, θv , λ). These functions are also pre-computed.
The interval ∆t is determined in a similar way to the method
for the light beams (see Section 4.2.1).
The look-up tables for gl (h, θsun , λ), ∆Ir (h, θv , λ),∆Im (h, θv , λ), and ∆gv (h, θv , λ) are used as textures. The
wavelength λ is sampled at RGB. These textures are mapped
onto the sampling planes. The phase functions are computedat the vertices of the planes and interpolated by the Gouroudc The Eurographics Association 2002.


cloudscloud sphere
m
m-1
m-2

earth
1

Figure 8: Rendering of clouds viewed from space.
shading function. Then the planes are rendered in back tofront order. Only the planes that intersect with the atmosphere are generated at the interval ∆t. The images are generated by computing Eq. 12 by hardware color blending functions.
5.3. Rendering the Atmosphere as Viewed from SpaceAs shown in Fig. 7, we use sample spheres for rendering
the atmosphere as viewed from space. The integral term in
Eq. 11 is evaluated by using the sample spheres. The heightof sample sphere k above the ground is hk and the viewing
ray hits sample sphere k at Pk (see Fig. 7). The intensity of
the light reaching Pk is obtained by multiplying the intensityof the sunlight and the attenuation function gl (hk , θsun,k , λ)as described in Section 5.2. In the following, we will discussand ∆Im,kand the attenuation ratio, ∆gv,k ,the scattering, ∆Ir,kdue to particles in Pk Pk+1 .
Let us assume that the angle between the vertical direction at Pk and the direction of Pk Pk+1 is θv (see Fig. 7).
, ∆Im,k, and ∆gv,k are expressed as functions ofThen, ∆Ir,kθv (see appendix A.2). Therefore, these are stored as a onedimensional texture for each sample sphere. Eq. 12 is eval(θv , λ), ∆Im,k(θv , λ), anduated by using the textures for ∆Ir,k
∆gv,k (θv , λ) instead of ∆Ir (θv , λ), ∆Im (θv , λ), and ∆gv (θv , λ).
The images are generated by rendering the sample spheres
mapped with these textures. We specify the texture coordinates on a per-vertex basis. Westermann et al. proposed
an accurate method16 for mapping textures onto spheres byspecifying the texture coordinates on a per-pixel basis withthe pixel texgen OpenGL extension. In our case, however,assigning the texture coordinates on a per-vertex basis issufficient since the intensity of the scattered light is verysmooth. The phase functions are computed at the verticesof the mesh of each sphere.
5.4. Rendering of Clouds
For rendering the earth viewed from space, threedimensional clouds are generated from satellite images.
Methods developed in the field of remote sensing enable us
to determine the height of a cloud top26 . Since any estimation of the density inside the clouds and the height of thecloud base are difficult problems to solve, the density andthe height are simply assumed to be a constant and specified
by the user. Then, as shown in Fig. 8, m cloud spheres areDobashi, Yamamoto, Nishita / Interactive Rendering of Atmospheric Effects

generated. For each cloud sphere, a texture Ti (i = 1, · · · , m)
is mapped. Let us call Ti cloud textures in the following.
The cloud texture represents the density distribution on eachcloud sphere. Each element of the cloud texture Ti stores thedensity obtained from the satellite image. The alpha channels of the cloud texture store the attenuation ratio of thelight that passes through the clouds as represented by cloudsphere i.
The rendering method for clouds is as follows. In the following, let us assume for the sake of simplicity that there
is no atmosphere. The attenuation due to the atmosphereis computed by using the texture for the attenuation function gl (h, θsun , λ). The rendering of clouds consists of two
steps. The first step is to compute the attenuation ratio of
the sunlight reaching points on the cloud spheres. To do this,the camera position is placed at the position of the sun andcloud spheres are rendered in descending order based on
their distances from the center of the earth. After drawing
cloud sphere i, the pixel values in the frame buffer represent
the attenuation ratio between the sun and points on the cloudsphere. This image is stored as the attenuation texture Gi .
These processes are repeated for all cloud spheres.
Next, in the second step, images are generated by using
the cloud texture and the attenuation texture. First, cloud texture Ti and attenuation texture Gi are mapped onto cloudsphere i by using the multiple texture mapping function.
Next, the cloud spheres are rendered in back to front order,that is, the inmost sphere is rendered first. The phase function of the cloud particles is evaluated at the vertices of themesh of the cloud spheres.
6. Results
In this section, several images generated by our method aredemonstrated† . Fig. 9 shows examples of light beams rendered by the method described in Section 4. In Fig. 9(a),the proposed method is applied to the lighting design of astudio. Fig. 9(b) shows typical shafts of light caused by sunlight. The sunlight passes through the stained glass windows
of a church. Next, Figs. 10 and 11 show examples of rendering the earth’s atmosphere. Fig. 10(a) is an example ofrendering the sky. This image also demonstrates shafts oflight due to mountains that shut out parts of the sun’s rays.
Figs. 10(b), (c), and 11 show the earth viewed from space.
In Fig. 10(b), realistic clouds as well as the atmosphere arerendered. In Fig. 10(c), a thin atmospheric layer covers theearth. Fig. 11 shows a sequence of images of the sunrise
viewed from space. The atmosphere exhibits beautiful color
variations from red to blue. In Fig. 10(b), the method described in Section 5.3 is used, since the viewpoint is far fromthe earth. In Figs. 10(c) and 11, the method described in Section 5.2 is used, since the viewpoint is near to the surface.
† see color section for Figs. 9(b), 10(c), and 11

Table 1: Computation time.
time [sec.]

No. of planes[sub-planes]

Fig. 9(a)0.32

9[309], 13[174],
16[162]

Fig. 9(b)0.1230[263]Fig. 10(a)0.16

100(sky),30(shaft)[135]Fig. 10(b)0.06

10(sample spheres)Figs. 10(c) & 11

0.10

300

Some of software packages approximate the atmosphere byusing spheres with a simple hue shift through the atmosphere(e.g. 27 ). However, it is difficult to simulate the color variations shown in Figs. 10 and 11 by using a simple hue shift.
Table 1 shows computation times. We use the same computer as the one described in Section 4.3. The size of each
image is 720×480. The numbers of the main planes and subplanes are also shown in the table. For Fig. 9(a), the numbers
for each of three spotlights are shown. For Fig. 10 (c), thenumber of the sample sphere is shown. In Figs. 9 and 10(a),the threshold to insert the sub-planes is 0.01 (see Section
4.2). In Fig. 10(a), 100 planes are used for computing thesky colors and 30 planes are used for rendering the shafts oflight due to the mountains. In Fig. 10(c), the number of cloudspheres used is 5. As shown in Table 1, our method can create the images at interactive rates. The sub-planes and 3D
clouds can be omitted for real-time previews. For example,
Fig. 9(b) is rendered in 90 fps (frame per second) without
the sub-planes and Fig. 10(c) is rendered in 30 fps without
3D clouds.

7. Conclusion
We have proposed a fast rendering method for atmosphericscattering effects. In the method, two-dimensional texturesthat store the intensities of the scattered light are preprepared. The luminous intensity distribution of the lightsource and the object shadows are sampled precisely by using the sub-planes. We have also realized the efficient rendering of the earth’s atmosphere. The sky and the earth asviewed from space are rendered almost in real-time. We havedeveloped a method of rendering realistic clouds by using
satellite images.
In future work, we need to discuss the problem posed by
the multiple scattering of light. For participating media witha high albedo, the realism of images might be enhancedby taking into account the multiple scattering. An efficientmethod for the multiple scattering has to be developed.
c The Eurographics Association 2002.
Dobashi, Yamamoto, Nishita / Interactive Rendering of Atmospheric Effects

Acknowledgements
The authors would like to acknowledge Miss Keiko Sotome(Institute of Advanced Media Arts and Sciences) for her help
in making the geometric data for the church and the opera
house in Fig. 9. This research was partially supported by the
Ministry of Education, Japan, Grant-in-Aid for Encouragement of Young Scientists (B) (2) 14750271, 2002.
References
1.

T. Nishita, Y. Miyawaki, E. Nakamae, "A Shading
Model for Atmospheric Scattering Considering Distribution of Light Sources," Computer Graphics,
21(4):303–310, 1987.

2.

N. Max, "Atmospheric Illumination and Shadows,"
Computer Graphics, 20(4):117–124, 1986.

3.

A. Preetham, P. Shirley, B. Smits, "A Practical Analytic
Model for Daylight," Proc. SIGGRAPH’99, 91–199,
1999.

4.

5.

T. Nishita, T. Shirai, K. Tadamura, E. Nakamae, "Display of The Earth Taking into account Atmospheric
Scattering," Proc. SIGGRAPH’93, 175–182, 1993.
Y. Dobashi, K. Kaneda, H. Yamashita, T. Okita,
T. Nishita, "A Simple, Efficient Method for Realistic
Animation of Clouds," Proc. SIGGRAPH 2000, 19–
28, 2000.

12. U. Behrens, R. Ratering, "Adding Shadows to a
Texture-based Volume Renderer," Proc. 1998 Simposium on Volume Visualization, 39–46, 1998.
13. E. Lamar, B. Hamann, K. I. Joy, "Multiresolution Techniques for Interactive Texture-Based Volume Visualization," Proc. IEEE Visualization’99, 355–362, 1999.
14. J. Stam "Stable Fluids," Proc. SIGGRAPH’99, 121–
128, 1999.
15. M. Nulkar, K. Mueller, "Splatting with Sadows," Proc.
Volume Graphics 2001, 35–49, 2001.
16. R. Westermann, T. Ertl "Efficiently Using Graphics
Hardware in Volume Rendering Applications," Proc.
SIGGRAPH’98, 169–177, 1998.
17. T. J. Cullip, U. Neumann "Accelerating Volume Reconstruction with 3D Texture Hardware," Technical Report
TR93-027, University of North Carolina, Chapel Hill
1993.
18. M. Kawara, http://www.daionet.gr.jp/ masa/ ishadowmap/index.html
19. C. Everitt, http://www.r3.nu/ cass/shadowsandstuff/
20. A. Keller, W. Heidrich, "Interleaved Sampling," Proc.
Rendering Technique 2001, 269–276, 2001.
21. R. Mech "Hardware-Accelerated Real-Time Renderingof Gaseous Phenomena," Journal of Graphics Tools,
6(3):1–16, 2001.

6.

B. Cabral, N. Cam, J. Foran, "Accelerated Volume Rendering and Tomographic Reconstruction Using Texture
Mapping Hardware," Proc. 1994 Simposium on Volume
Visualization, 91–98, 1994.

22. C. A. Lee, C. Kesselman, "Near-Real-Time Satellite
Image Processing: Metacomputing in CC++," IEEE
Computer Graphics and Applications, 16(4):79–84,
1996.

7.

Y. Dobashi, T. Yamamot, T. Nishita, "Interactive Rendering Method for Displaying Shafts of Light," Proc.
Pacific Graphics 2000, 31–37, 2000.

23. N. Max, R. Crawfix, D. Williams, "Visualization for
Climate Modeling," IEEE Computer Graphics and Applications, 13(4):34–40, 1993.

8.

Y. Dobashi, T. Yamamot, T. Nishita, "An Accurate,
Fast Method Using Graphics Hardware for Rendering
Shafts of Light," The Journal of The Institute of Image
Information and Television Engineers, 55(7):362–370,
2001.

24. K. Engel, M. Klaus, T. Ertl, "High-Quality PreIntegrated Volume Rendering Using HardwareAccelerated Pixel Shading," Proc. Graphics Hardware
2001, 9–16, 2001.

9.

H. W. Jansen, P. H. Christensen "Efficient Simulationof Light Transport in Scenes with Participating Media
using Photon Maps," Proc. SIGGRAPH’98, 311–320,
1998.

10. T. Nishita, Y. Dobashi, K. Kaneda, H. Yamashita, "Display Method of the Sky Color Taking into AccountMultiple Scattering," Proc. Pacific Graphics’96, 117–
132, 1996.
11. H. E .Rushmeier, K. E. Torrance, "The Zonal Method
for Calculating Light Intensities in the Presenceof a Participating Medium," Computer Graphics,
21(4):293–302, 1987.
c The Eurographics Association 2002.


25. M. Segal, C. Koroblkin, R. V. Widenfelt, J. Foran,
P. E. Haeberli, "Fast Shadows and Lighting Effects Using Texture Mapping," Computer Graphics,
26(2):249–252, 1992.
26. C. D. Rodgers, "Retrieval of Atmospheric Temperature and Composition from Remote Measurements ofThermal Radiation," Review of Geophysics and Space
Physics, 14:609–624, 1976.
27. http://ennui.shatters.net/celestia
A. Appendix
Let us assume a coordinate system with its origin at the center of the earth. In this system, the coordinates of point P in

Dobashi, Yamamoto, Nishita / Interactive Rendering of Atmospheric Effects

Fig. 6 are assumed to be (0, 0, Re + h), where Re is the radiusof the earth.

A.1 Attenuation of Sunlight
The optical length τ is expressed by the following equation.
τ(t, λ) =t0βr (λ)ρr (l) + βm (λ)ρm (l)dl,(A.1)

where βr and βm are the extinction coefficients, and ρrand ρm are the densities of air molecules and aerosols, respectively. These are given by ρr = exp(−h/Hr ) and ρm =exp(−h/Hm ), where h is the height from the ground, and Hrand Hm are constants (Hr = 7994[m], Hm = 1.2[km]).
Assuming the height of the top of the atmosphere is
H. A vector in the direction of the sun is given as(0, sin θsun , cosθsun ). The length s of PPs in Fig. 6 is:
s = −(Re + h) cos θsun + (Re + H)2 − (Re + h)2 sin2 θsun .
(A.2)
A coordinate of a point on PPs is given by (0, s sin θsun , Re +h + s cos θsun ), where s is the length between point P and
the point on PPs . The height of this point is given byhs (s ) = s2 sin2 θsun + (Re + h + s cos θsun )2 − Re .
(A.3)
Eqs. (A.1), (A.2), and (A.3) indicate that the optical lengthτ(s, λ) between PPs is expressed as a function of h and θsun .
So, gl is also expressed as a function of h and θsun . That is,
gl (h, θsun , λ) = exp(−s0θv , that is,τ(PP , λ) = −∆t0βr (λ)ρr (h(t  )) + βm ρm (h(t  ))dt  .

(A.6)
In the case of the sample spheres, by replacing P and P with


∆Im,k
, and ∆gv,kPk and Pk+1 in Fig. 7, it is obvious that ∆Ir,kare expressed as functions of h and θv . The height from the, ∆Im,k, andground is fixed at hk for sample sphere k. So ∆Ir,k
∆gv,k are functions of θk .

(a) Light beams caused by spotlights.

Figure 9: Examples of rendering light beams.

βr (λ)ρr (hs (s ))+βm ρm (hs (s ))ds ).
(A.4)

A.2 Intensity of Light Scattering and Attenuation Ratio
In Fig. 6, let us consider a point on PP , denoted by Q.
The coordinate of point Q is given by (0,t  sin θv , Re + h +t  cos θv ), where t  is the length of PQ. The height of point
Q, hv (t  ), is given by replacing s and θsun in Eq. (A.2) witht  and θv , respectively. Therefore, hv (t  ) is a function of t  , h,and θv . ∆Ir is equivalent to the integrated value of the scattering intensity between PP in Fig. 6. The density at point
Q and the optical length between PQ are computed by using
hv (t  ). Therefore, ∆Ir is expressed as a function of h and θv .
That is,∆Ir (h, θv , λ) =∆t
0(a) Sky colors and shafts of light.

Kr (λ)ρ(hv (t  )) exp(−τ(hv (t  ), λ))dt  .

(A.5)
Similarly, ∆Im is expressed as a function of h and θv . The attenuation ratio ∆gv is obviously a function of h and θv since
the optical length between PP is a function of h and θv ,that is,Similarly, ∆Im is expressed as a function of h and θv .
The attenuation ratio ∆gv is obviously a function of h and θv
since the optical length between PP is a function of h and(b) The earth covered by the bluish atmosphere.

Figure 10: Examples of rendering earth’s atmosphere.
c The Eurographics Association 2002.
Dobashi, Yamamoto, Nishita / Interactive Rendering of Atmospheric Effects(b) Light beams through stained glass windows.
Figure 9: Examples of rendering light beams.

(c) The earth with the atmosphere.
Figure 10: Examples of rendering earth’s atmosphere.

Figure 11: The sunrise viewed from space.

c The Eurographics Association 2002.