Modeling Skylight and Aerial Perspective

Arcot J. Preetham
ATI Research
preetham@[Link]

Theory of Scattering in Atmosphere

Atmosphere
Earth is enveloped by a vast amount of air called the atmosphere. It is relative thin
compared to the size of the earth and fades away with increasing distance from the earth's
surface. The earth's atmosphere is categorized into 4 layers - troposphere, stratosphere,
mesosphere and thermosphere. Troposphere where all the weather takes place and
stratosphere, the two closest layers to earth's surface constitute more than 99% of the
atmosphere. The troposphere and stratosphere extend up to 12 km and 53 km respectively
from the earth's surface. The earth's atmosphere is composed of many gases. Nitrogen
from decay of biological products constitutes 98% and Oxygen from photosynthesis
constitutes 21% of the atmosphere. Gravity holds the atmosphere close to the earth's
surface and explains why the density of the atmosphere decreases with altitude.

The density and pressure of the atmosphere vary with altitude and depends on solar
heating and geomagnetic activity. The simplest is an exponential fall off model where
pressure and density decrease exponentially with altitude. In 1976, US Standard
Atmosphere Model was adopted by COESA (Committee on Extension to the Standard
Atmosphere) which describes the earth's atmosphere as composed of 7 layers up to 86 km
and pressure, temperature and density are specified for each layer [McCartney1976]. The
density is calculated using a perfect gas relationship.

In addition to the various gases, atmosphere also contains water vapor, dust particles, etc.
The molecules and particles absorb energy at discrete wavelengths, which are determined
by their internal properties. For example, molecular oxygen and ozone absorbs light in
the ultra violet spectrum. Water vapor, methane, nitrous oxide, ozone, and CO2 absorb
light in the infrared range.

In addition to absorption, molecules and particles also scatter energy out from its original
direction. Sun’s white light is scattered once (primary scattering) or multiple times
(secondary scattering) into the viewing ray as shown in Figure 1. The scattered light is
received at the earth’s surface from all directions as diffuse skylight or daylight.
 secondary primary

Atmosphere Earth

Figure 1: Scattering of sunlight in atmosphere.

Earth’s surface is not flat and this plays a very important role in atmospheric optics.
Optical mass of a path is defined as the mass of the medium in that path of unit cross-
section and is given by
s
m = ∫ ρ ( x)dx ,
0

where ρ(x) is the density of medium. Optical length for a path is defined as the optical
mass divided by the density at the earth’s surface ρ0 and is given by
1 s
l= ∫ ρ ( x)dx
ρ0 0

Optical length has dimensions of length. The optical length in the zenith direction for
molecules is 8.4 km and for aerosols is 1.25 km. Figure 2 shows the optical lengths for
different directions in the atmosphere.

The relative optical length is defined as the ratio of the optical length of any path to
optical length at zenith direction and is given by following approximation [Iqbal1983]

l (θ s ) 1
= ,
l zenith cosθ s + 0.15(93.885 − θ s ) −1.253
where θs is the angle from zenith in degrees.
 8.4 km, 1.25 km
s 11.9 km, 1.8 km
16.9 km, 2.54 km
o
45
o
60 32.6 km, 4.9 km
o
75
θ

Earth Atmosphere

Figure 2: Optical lengths for molecules and aerosols for different paths in the atmosphere.

The sun's light travels through a much larger atmosphere when it is close to the horizon.
Therefore, a larger amount of blue light is scattered away causing the sun to appear
orange-red. Increased atmosphere presence in the horizon direction makes the stars
appear brighter in the zenith direction compared to the horizon direction.

Scattering
Scattering is a process by which a particle redistributes a fraction of the incident energy
into a total solid angle. The scattering properties depend on the refractive index and size
of the particles.

It is common for smaller particles to scatter uniformly in the forward and backward
directions and for larger particles to scatter strongly in the forward directions. Scattering
by one particle is independent of the other as long as the distance between the particles
are greater than the particle size. This is known as independent scattering. A scattering
event for the first time is known as first order scattering. Scattered light may be scattered
again by another particle and is said to undergo second order scattering. In reality, sun’s
light is scattered multiple times in the atmosphere.

The amount of scattering is linearly proportional to the density of the atmosphere, which
varies with altitude. A choice of an analytic exponential model or a lookup table based on
US Standard Atmospheres for density is available for use.

Rayleigh scattering
Particles smaller than the wavelength and usually less than 0.1 times the wavelength of
light exhibit Rayleigh scattering [Rayleigh1871]. Discovered by the Nobel Prize winner
Lord Rayleigh, Rayleigh scattering is observed by molecules in the earth's atmosphere.
The amount of scattering for such particles is inversely proportional to the 4th power of
the wavelength. These particles scatter equally in the forward and backward directions.
The total scattering coefficient β and the angular scattering coefficient β(θ) are given by
 8π 3 (n 2 − 1) 2 6 + 3 pn
β= ( ),
3 Nλ4 6 − 7 pn
π 2 (n 2 − 1) 2 6 + 3 pn
β (θ ) = ( )(1 + cos 2 θ ) ,
2 Nλ 4
6 − 7 pn

where n is the refractive index of air and is 1.0003 in the visible spectrum, N is the
number of molecules per unit volume and is 2.545x1025 for air at standard temperature
and pressure, pn is the depolarization factor with a value of 0.0035 standard for air. The
total scattering coefficient for blue light (400 nm) is 2.44x10-5m-1, for green light (530
nm) is 1.18x10-5m-1 and for red light (700 nm) is 6.95x10-6m-1. This means that the blue
light is scattered more than the red light, which explains the blue color of the sky and the
red color of the sun at low altitudes.

Angular scattering coefficient is also equivalent to the total scattering coefficient times
the phase function. The phase function for Rayleigh scattering fair(θ) is given by
3
f air (θ ) = (1 + cos 2 θ )
16π

Mie Scattering
Larger particles scatter strongly in the forward direction and this scattering phenomenon
is called Mie scattering named after Gustav Mie. The scattering is inversely proportional
to the second order of the size of the particles and is independent of wavelength. The
phase function for angular scattering was approximated by Henyey-Greenstein and is
given by the following equation [Henyey1941].
1 1− g 2
f HG (θ ) = .
(
4π 1 − 2 g cosθ + g 2 3 / 2 )
Positive values of g represent forward scattering and negative values of g represent
backward scattering. The total scattering coefficient is given by
2π
β = 0.434cπ ( )v − 2 K ,
λ
where c is the concentration factor and varies around 6x10-17 to 25x10-17 as the turbidity
increases from clear to overcast, v is the Junge’s exponent and a value of 4 is standard for
sky model, and K varies from 0.656 for 400 nm to 0.69 for 770 nm. For more details, see
[Bullrich1964].
Skylight Models
The atmosphere scatters sun's light multiple times and scattered light is received at the
earth's surface in all directions collectively known as skylight.

Is
Q
S
ωs
R Ii
ω’
R
Atmosphere
ω x
I i+1

P
Earth

Figure 3: The variables used in skylight computation.

In Figure 3, we would like to compute the skylight in the viewing direction (ω) PS. R is a
variable point on PS and RQ is the direction (ωs) of sunlight. Points Q and S are at the top
of the atmosphere. Light from the sun is attenuated as it travels the distance QR, is
scattered by particles at R in the direction RP, and then attenuated as it travels a distance
RP. Let lAB denote the optical length from A to B. The differential amount of light
reaching point P through the path QRP, scattered by a differential volume R is given by
− βl
dI = Es e QR β (ω , ω s )e − βlRP dx ,
where Es is the irradiance of the sun outside the earth’s atmosphere, β is the total
scattering coefficient, β(ω,ωs) is the angular scattering coefficient between directions ω
and ωs and dx is the differential optical length. The total light received at P from all
points R on PS is given by
S − βlQR
I 1 (ω ) = ∫ E s e β (ω , ω s )e − βl dx .
RP

Superscript “1” on I(ω) is used to indicate first order scattered light. Light from the sun is
scattered more than once before reaching the earth’s surface. Let Ii(ω’) be the light
reaching point R from direction ω’ after being scattered i times. Light scattered into
viewing direction ω from all directions ω’ at point R is denoted by S(ω,x) and is given by
the integral
4π
S (ω , x) = ∫ I i (ω ' ) β (ω , ω ' )dω '
0
The light scattered in direction ω at R is attenuated as it travels the distance RP through
the atmosphere. Light received at earth’s atmosphere after i+1 scattering events is given
by
S
I i +1 (ω ) = ∫ e − βlRP S (ω , x)dx
P
S 4π
I i +1
(ω ) = ∫ e − βl RP
∫ I i (ω ' ) β (ω , ω ' )dω 'dx
P 0

Starting from the equations for I1(ω) and Ii+1(ω) in terms of Ii(ω), one can compute I2(ω),
I3(ω) and so on, where I2(ω) and I3(ω) is the second and third order scattered light into
viewing direction ω. Total light scattered into the viewing ray I(ω) is the sum of first
order, second order and higher orders of scattered light and is given by
I (ω ) = I 1 (ω ) + I 2 (ω ) + I 3 (ω ) + L
The above equations assume only one kind of particles in the atmosphere. Real
atmosphere consists a number of different kinds of particles, of most importance being
molecules and aerosols. The above equations can be extended to multiple particles and is
left as an exercise for the reader.

Simulation Based Methods

All simulation methods are based on the general atmosphere equations presented above.
The variations arise from use of different models for atmosphere density, different
atmosphere composition, scattering coefficients etc.

One of the first models for atmospheric scattering was presented by Klassen
[Klassen1987] which is a must read for anyone interested in modeling atmospheric
scattering. He used a simple constant density atmosphere model on a flat earth surface.
Flat earth model performs poorly for skylight computations on viewing rays close to the
horizon. For aerial perspective discussed later, where distances viewed are usually of the
order of a few tens of kilometers, the flat earth model is a good approximation.

Kaneda et al employed similar concepts to Klassen and used a more realistic atmospheric
model for his simulations [Kaneda1991]. He modeled a spherical earth with an
exponential decay density distribution.

Nishita et al [Nishita1996] take a step closer to reality and model higher orders of
scattering, which is responsible for the whitening effect of the sky. He proposes a fast
method for single scattering computations.

Analytic Models and Approximations

Simulation based methods are computationally very expensive making it unusable for all
practical purposes. Analytic or approximate models for skies that represent simulated
data or real world atmospheric data collected over the years is preferable for ease of use
in computer graphics. Of all real skies, the two extremes are clear and overcast skies.

CIE Overcast sky luminance model

An overcast sky is equivalent to a dark sky with plenty of clouds. Moon and Spencer
proposed a luminance model for overcast skies. This was later simplified and adopted by
the International Commission on Illumination - abbreviated as CIE from its French title
Commission Internationale de l'Eclairage [CIE1994]. The overcast sky luminance YOC is
given by
1 + 2 cosθ
YOC = Yz ,
3
where θ is the angle from the zenith and Yz is the zenith for overcast skies and can be
obtained from analytic formulas adopted by CIE [CIE1994].

CIE Clear sky luminance model

Pokrowski proposed a sky luminance model from theory and sky measurements for clear
sky. Kittler improved this clear sky model, which was later adopted by the CIE for the
clear sky model in 1973 [CIE1994]. The clear sky luminance YC is given by
(0.91 + 10e −3γ + 0.45 cos 2 γ )(1 − e −0.32 cosθ )
YC = Yz ,
(0.91 + 10e − 3θ s + 0.45 cos 2 θ s )(1 − e − 0.32 )
where Yz is the zenith luminance and the angles are defined in Figure 4.
The CIE luminance models provide us with relative luminance and not absolute
luminance.

ASRC-CIE model luminance model

ASRC-CIE model is a linear combination of four luminance models - the standard CIE
cloudless sky, a high turbidity formulation of the latter, a realistic formulation for
intermediate skies proposed by Nakamura and the standard CIE overcast sky. The sky
clearness and sky brightness factors are used to determine the weights for these four skies
[Littlefair1994].

θ
θs v

W N
γ

φs φ
S
E

Figure 4: Angles and directions on the sky dome [Preetham1999].

© Copyright 1999 by ACM, Inc.
Perez all weather luminance model
Perez et al proposed an all weather sky luminance model [Perez1993]. This model was
based on 5 different parameters, which related to darkening or brightening of the horizon,
luminance gradient near the horizon, relative intensity of the circumsolar region, width of
the circumsolar region and relative backscattered light. The Perez model is given by
F (θ , γ ) = (1 + Ae B cos θ )( 1 + Ce D γ + E cos 2 γ ) ,
where A, B, C, D and E are the distribution coefficients and γ and θ are the angles shown
in Figure 4. The Perez model sky luminance YP is given by
F (θ , γ )
YP = Yz .
F (0,θ s )
Perez model is similar to CIE model and has been found to be more accurate and has also
been validated by Ineichen [Ineichen1994].

Preetham et al spectral radiance model

Similar to the Perez's all weather luminance model, Preetham et al proposed an analytical
model for spectral radiance of sky [Preetham1999]. This model shows variation from
clear to overcast sky through a parameter called turbidity. Analytic model was arrived
from atmospheric simulations using US Standard Atmospheres, and modeling skylight up
to second order scattering. Luminance Y and chromaticity values x and y are given by
F (θ , γ )
Y = Yz ,
F (0,θ s )
F (θ , γ )
x = xz , and
F (0,θ s )
F (θ , γ )
y = yz ,
F (0,θ s )
where F(θ, γ ) is from Perez’s model with different values of A, B, C, D and E for Y, x
and y. The distribution coefficients for luminance and chromaticity values x and y,
absolute values for zenith luminance Yz, zenith chromaticity xz and yz are all given in
[Preetham1999]. Conversion of chromaticity values into spectral radiance is given by the
following standard method for D-illuminants [Wyszecki1982]. The relative spectral
radiant power SD(λ) is given by
S D (λ ) = S0 (λ ) + M 1S1 (λ ) + M 2 S 2 (λ ) ,
where S0(λ) is the mean spectral radiant power and S1(λ) and S2(λ) are the first two eigen
vector functions used in daylight illuminants. M1 and M2 are functions of chromaticity
values x and y and are given by
− 1.3515 − 1.7703x + 5.9114 y
M1 = ,
0.0241 + 0.2562 x − 0.7341y
0.0300 − 31.4424 x + 30.0717 y
M2 = .
0.0241 + 0.2562 x − 0.7341y

A scene rendered using this model is shown in Figure 5.

 Figure 5: Spectral radiance model for sky [Preetham1999].
© Copyright 1999 by ACM, Inc.

Aerial Perspective Model

Many artists have recorded the effect of distant mountains fading away in history. In the
presence of a medium such as atmosphere distant objects appear blue. This shift in color
helps us in perceiving depth and is known as aerial perspective.

Distant objects appear hazier and this is attributed to scattering and absorption along the
viewing ray as light travels from the source to the viewer. The light from the source loses
intensity and undergoes a spectral shift as scattering and absorption depend on
wavelength. In addition to this loss of light, light from other sources like sun, sky and
ground are scattered into the viewing ray.
 L0 Ls
S
ω inscatter
x R

Ls
s P

Figure 6: The variables in aerial perspective computation.

In Figure 6, L0 is the radiance of the distant hill and Ls is the radiance of the ray at the
viewer. If f is the extinction factor and Lin is the in-scattered light as L0 travels a distance
s to the eye, then
Ls = fL0 + Lin
The extinction factor f for light traveling the path PS with optical length lPS is given by
f = e − βl PS
Let Ls(ω) denote the spectral radiance of the sun and sky in the direction ω. Let S(ω, x) be
the term that denotes the light scattered from direction ω into the viewing direction at
point R. Therefore,
4π
S (ω , x) = ∫ Ls (ω ' ) β (ω ,ω ' )dω '
0
The total light scattered into viewing direction at point R is
S
Lin = ∫ e − βl RP S (ω , x)dx ,
P
S 4π
Lin = ∫ e − βl RP ∫ Ls (ω ' ) β (ω ,ω ' )dω 'dx .
P 0

For landscape scenes that focus on aerial perspective, the viewing rays are close to the
earth’s surface and it can safely be assumed that the density of the medium is a constant
and is equal to that at the earth’s surface. For such rays, the optical length lAB is equal to
the distance AB. Therefore, lPS is equal to s, lRP is equal to (s-x) and the above equations
simplify to
f = e − βs and
s 4π
Lin = ∫ e − β ( s − x ) ∫ Ls (ω ' ) β (ω ,ω ' )dω 'dx .
0 0
The primary contribution to in-scattered light is the direct light from the sun. Therefore,
we can safely ignore second order scattering i.e. the light from the sky without loss of
quality. Therefore, if Es is the irradiance from the sun at the earth’s surface and ωs is the
sun’s direction, then
4π
∫0
Ls (ω ' ) β (ω ,ω ' )dω ' = E s β (ω , ω s ) .
And the total in-scattered light simplifies to
s
Lin = ∫ e − β ( s − x ) E s β (ω , ω s )dx ,
0

β (ω ,ω s )
Lin = E s (1 − e − βs ) .
β
While the above equations for aerial perspective are valid for scattering by one kind of
particles, they can easily be extended to many kinds of particles, and is left as an exercise
for the reader.

Accelerated Techniques for Modeling Sky and Aerial

Perspective
Accurate calculation of the sky color is expensive because of numerical integration of
scattered light along the viewing ray. Modeling higher orders of scattering results in
double and triple integration and the computational cost increases exponentially. Dobashi
et al discuss various ways to represent and evaluate the sky intensity across the
hemispherical domain for various altitudes of the sun [Dobashi1995]. A simple and naive
approach would be to represent the sky dome by a mesh grid of sample points. Skylight
intensity at any point on the hemisphere can be obtained by linearly interpolating
between these sampled points. Any function can be represented by the summation of a set
of basis functions and another approach to representing the sky function over the
hemisphere is to use spherical harmonics. The weights to these basis functions can be
calculated in a preprocess step.

Dobashi et al proposed another set of basis functions called the cosine functions, which
use less memory than the spherical basis functions and is faster to evaluate. With recent
advances in the area of programmable graphics processors, these cosine functions can be
evaluated in real time in shaders.

Researchers have constantly looked into using graphics hardware to accelerate modeling
of sky and aerial perspective [Dobashi2000][Dobashi2002][Hoffman2002].

Rendering shafts of light through clouds has been done in CPU by numerical integration
of in-scattered light along the viewing ray taking into account the visibility of the sample
points. The numerical integration can be approximated by a summation of terms along
the viewing ray. Dobashi et al achieved this summation by rendering many virtual planes
along the viewing ray and accumulating the various terms using blending [Dobashi2000].
The visibility information at any virtual plane is calculated using the standard shadow
map technique [Williams1978]. The light information at the virtual plane is evaluated
using projective light textures technique. By careful choice of the number of virtual
planes and the resolution of the mesh for the virtual plane, photo realistic scenes with
atmospheric scattering can be generated at interactive frame rates.
Dobashi et al accelerated their previous techniques of rendering a sequence of virtual or
sampling planes and accumulating the terms using blending [Dobashi2002]. The novelty
in aerial perspective is that the summation is written as a product of a high frequency
term and a low frequency term. The low frequency term can be evaluated accurately in a
preprocess step and is written as a product of two terms; one term is stored in the vertices
of the mesh grid of the sampling plane and the other is stored as a texture map. A
sequence of planes is drawn with blending to get the final image. A rendering from their
publication is shown in Figure 7.

Figure 7: Interactive atmospheric scattering [Dobashi2002].

© Copyright 2002 by Eurographics. Included by permission.

Modern graphics hardware (e.g. Radeon 8500, Radeon 9700, GeForce 4, GeForce FX)
allows the user to specify a vertex shader for transformation and lighting. Hoffman et al
presented a new technique to render sky color and aerial perspective in real time using
programmable graphics hardware [Hoffman2002]. The two terms for aerial perspective
are extinction factor f and in-scattered light Lin. These are extended to two kinds of
particles - molecules and aerosols. f and Lin are calculated in the vertex shader on a per
vertex basis and is passed on to the fragment program through the color or texture
registers. The final color is computed in the fragment program using the equation
Ls = fL0 + Lin .

With the latest generation of graphics hardware boasting of a full floating point pipeline
(e.g. Radeon 9700, Geforce FX), the entire computation of aerial perspective terms (f and
Lin) can be done at a pixel level rather than a vertex level as was presented. The
advantage of a per pixel computation is that it does not require very fine tessellation of
the geometry to avoid artifacts due to linear interpolation.

The image in Figure 8 was rendered on 600MHz Pentium III with a Radeon 8500 at
about 60 frames per second.

Figure 8: Real time atmospheric scattering [Hoffman2002].

References

[Bullrich1964] Bullrich, K. Scattered radiation in the atmosphere. In Advances in

Geophysics, 10, Academic Press, 1964.

[CIE1994] CIE-110-1994. Spatial distribution of daylight – luminance distributions of

various reference skies. Tech. Rep., International Commission on Illumination, 1994.

[Coulson1975] Coulson, K.L. Solar and terrestrial radiation. Academic Press, 1975.

[Dobashi1995] Dobashi, Y., Nishita, T., Kaneda, K., Yamashita, H. Fast display method
of sky color using basis functions. In Pacific Graphics,1995.

[Dobashi2000] Dobashi, Y., Yamamoto, T., Nishita, T. Interactive rendering method for
displaying shafts of light, In Pacific Graphics, 2000.

[Dobashi2002] Dobashi, Y., Yamamoto, T., Nishita, T. Interactive rendering of

atmospheric scattering effects using graphics hardware. Proceedings of the conference
on Graphics hardware, 2002.

[Hoffman2002] Hoffman, N., Preetham, A.J. Rendering outdoor light scattering in real
time. Game Developers Conference, 2002.

[Henyey1941] Henyey, L.G., and Greenstein, J.L. Diffuse reflection in the Galaxy.
Astrophysical Journal 93:70, 1941.

[Ineichen1994] Ineichen, P., Molineaux, B., and Perez, R. Sky luminance data validation:
comparison of seven models with four data banks. Solar Energy 52, 4, 337-346, 1994.

[Iqbal1983] Iqbal, M. An introduction to Solar Radiation. Academic Press, 1983.

[Kaneda1991] Kaneda, K., Okamoto, T., Nakame, E., and Nishita, T. Photorealistic
image synthesis for outdoor scenery under various atmospheric conditions. The Visual
Computer 7, 5 and 6, 247-258, 1991.

[Klassen1987] Klassen, R.V. Modeling the effect of the atmosphere on light. ACM
Transactions on Graphics 6, 3, 215-237, 1987.

[Littlefair1994] Littlefair, P.J. A comparison of sky luminance models with measure data
from Garston, United Kingdom. Solar Energy, 53, 4, 315-322, 1994.

[Lynch1995] Lynch, D.K., and Livingston, W. Color and light in nature, Cambridge
University Press, 1995.

[McCartney1976] McCartney, E.J. Optics of the atmosphere. Wiley publication, 1976.

[Minnaert1954] Minnaert, M. Light and color in the open air. Dover, 1954.

[Nishita1996] Nishita, T., Dobashi, Y., Kaneda, K., and Yamashita, H. Display method
of the sky color taking into account multiple scattering. In Pacific Graphics, pp. 117-132,
1996.

[Preetham1999] Preetham, A.J., Shirley, P., Smits, B.E. A practical analytical model for
daylight. In Proceedings of Siggraph, 91-100, 1999.

[Robinson1966] Robinson, N. Solar radiation. Elsevier Publishing Company, 1966.

[Rayleigh1871] Rayleigh, L. On the scattering of light by small particles. Philosophical

Magazine 41, 447-451, 1871.

[VanDeHulst1981] van De Hulst, H.C. Light scattering by small particles. Dover

Publications, 1981.

[Williams1978] Williams, L. Casting curved shadows on curved surfaces. In Proceedings

of Siggraph, 270-274, 1978.

[Wyszecki1982] Wyszecki, G., and Stiles, W.S. Color Science. Wiley-Interscience

publication, 1982.