(~~

ComputerGraphics,Volume21, Number4, July 1987

A SHADING MODEL FOR ATMOSPHERIC SCATTERING CONSIDERINGLUMINOUS INTENSITY DISTRIBUTION OF LIGHT SOURCES
Tomoyuki
Nishita
Fukuyama University
Higashimura,Fukuyama
Hiroshima 729-02, Japan

Yasuhiro Miyawaki and Eihachiro Nakamae
Hiroshima University
Saijo, Higashi-Hiroshima
Hiroshima 724, Japan
Abstract
Studio spotlights
producedazzling
shafts oflight,
while light scattered from fog illuminatedby automobile headlights renders driving difficult.
This is because the particles
in the illuminatedvolumebecome visibleby scattering
light.
A
shading model for scatteringand absorption oflight caused by particles
in the atmosphereis
proposedin this paper.
The method takes into
account luminous intensity
distributionof light
sources,shadows
due to obstacles, and density ofparticles.
The intensity
at aviewpointiscalculated
by integrationof light scattered byparticles between the viewpoint and a given pointon an object.
The regions to be treated in thismanner are localized
by considering
illumination
volumesand shadow volumes caused by obstacles in
the illumination volumes.
CR Categories
and Subject Descriptions:
1.3.7[Computer Graphics]: Three-Dimensional Graphics andRealism:
1.3.3[Computer Graphics]: Picture/Image
Generation
General Terms:Algorithm
Additional
Key
Words
and
Phrases:
atmosphericscattering,luminouscharacteristic, shadow volumes

shading,
intensity

I. INTRODUCTION
Many papers
have been published which enhance
the realism of calculated
images by considering
properties
ofobjects,such
as reflection,refraction, and transparency,
and properties
oflightsourcessuch as point sources,
linearsources, area sources,
polyhedronsources[I-3],and sky light [4].
In order to produce realistic
images, the phenomena in the space
between lightsources and objects must also be considered; that
Permi~ion ~ copy without fee aU or part of th~ martial ~ granted
prodded that the ~ p i ~ are n ~ made or dmtdbuted ~r d~e~commercial advantage, the ACM copyright noti~ and the tkle of the
pubfi~tion and ~ da~ appea~ and notice is Oven that ~pying ~ bypermi~ion of the Associ~ion ~ r Computing Machme~. To copy
~herwise, or ~ ~publish, requi~s a fee and/or specific perm~sion.

1987

ACM-0-89791-227-6/g7/007/0303

is, the effect of atmospheric particles
along thepath of the light.
Light scattered from particles
such as water and dust
in
the
air
causesilluminatedvolumes
to glow.
In particular,
understanding this effectmakes it possibletosimulate such real-world phenomena as: I) shafts oflight causedby headlights of automobiles, street
lamps, beacons, or searchlights
at night(theseeffects are important for drivers), 2) light beams
caused by spotlights in studios and on stages,
3)
shafts of light
pouring
through windows, and 4)
visibility reduction of traffic signals or traffic
signs in fog or haze.
Representativeof methods considering particles
in the air is the fog effect in day time,commonlyused in flight
simulators.
This method simply
calculates attenuation as a function of depth fromthe viewpoint.
Some light scattering models havebeen developed:
Blinn first developed
a method
adaptedto displaying
Saturn's
rings or cloudlayers[5], but this method
is limited to thinlayers.
Kajiya developed a method
taking accountof densitiesof particles [6].
Recently b~x madeit possible to display light beams passing throughgaps of clouds or trees [7].
In all these methods,
however,light sources are limited to a simple,
single parallel or isotropic source,
ignoring
thedistributionof luminous
intensityof the lightsource.
In order to display
illumination
effects
in fog,we musttake into account the luminous
intensity characteristics of actual light sources,
in particularspotlights and headlights, and also
must
consider effects
due tomultiplelightsources.
This paper presents a method considering
these factors.
In this paper we use the following
assumptions:I) light
sources treatedhere are point sources
with variable angular intensity distributions,anda parallel
source.
2) the light from multiple
scattering and the illumination onto objects due toscattering can be neglected.
3) objects
treatedhere are composed of convex polyhedra.

2. CALCULATION$00.75

FOR LIGHT SCATTERING

Illuminatedvolumesin the air glow by lightscattered from particles in the space, whilelight
traversing
the space is attenuated by absorption.
These effects
are modelledasfollows:forscattering,particles
are considered
to be aninfinitenumber of point sources; for attenuation,
the atmosphere is supposed to be a semi-transparent

303

~

SIGGRAPH'87, Anaheim,July 27-31, 1987

object.
First,
we discusstheeffectsoflightscattering
due to a point
source whose angular
range
ofilluminationislimited(e.g.
aspotlight).
As shown in Fig.
1, the intensity of lightreaching a viewpoint Pv from a point
Pi on anobjectcan be obtainedby summing the reflected
intensity, attenuated by traversing
PiPv throughthe absorbing medium, with the intensity scatteredin the direction of the viewpoint from each point
along PiPv.
The attenuation varies exponentiallywith distance (e.g., from a point on PiPv
to Pv),
as describedby Bouguer'slaw (e.g., [8]).
The
total contribution from scattering is obtainedbyintegratingover all points along PiPv .
For
example, at any point P on PiP v at distance s from
Pv, the intensity scattered in the direction of the
viewpoint
is Ip(s),and willbe attenuated bytraversing distance s to the viewpoint.
Thus, theintensity of light at the viewpoint is given by

I = Iie-~(L)+Ip(S)e-r(S)~d s

T(S)

=

f

r

Fig.

odt

~

.

PVpoint

(1)

0
where

light source

I Intensity of light reaching a viewpoint froma point Pi on an object.

(2)0o is the extinction coefficient per unit length anddepends
on particle
density(afunctionofpositionl,~(s) expresses
the optical depth ofPPv, and T(L) that of PiP v.
Consider the direct illumination due to a pointlightsource with luminous intensity I(0,~) (0 isthe angle
from illumination
axis;
~ is
the
revolution
angle;
see Appendix),and uniform
ambient light Ia, which is due to interreflectionof light
from objects'
surfaces and to multiple
scattering.
The intensity of light reaching
the
point
P from the light source is also attenuated
due to absorption by particles and subject
to theinversesquarelaw of distance,
Then Ip is givenby the following equation.

Ip(s) = ~1(0'$) w F(a) e "T(r)+ Ia
r2~

K(I+ 9coalS(a~2))-0
In applying
thisequation,
properties should be noted:

K(1+50cosd*(a/2))

thefollowing

The third
term refers to the shaft of light.
Thus, if the light ray doesn't pass throughanilluminatedvolume,the intensity is given by
the first two terms.
b) Even for uniform density, we may consider many
types
ofluminous
intensity
distributions,
therefore
numerical
integrationisusefulbecauseof the difficultyof obtainingananalytical solution.
c) Outdoors, during the day, Ia corresponds to theintensityoftheatmosphereatinfinite
distance,that is, the first
two terms ofequation(51 arc equivalent to the equation forfog effect [10].
d) For multiple light sources, the total intensity
is obtainedby summingup theintensities
yielded by equation (5) for each light source.
e) In the case of a parallel
source,~ ofequation (3) becomes:

Ip(s) = wroF(a) e-~(r)+ ia:hazy atmosphere:murky atmosphere

(51

a)

(4)

where K is a constant.

304

I=Iie-°L+Ia (1-e -aL) +jLwI~(6,¢) e -°(r+s) F(~)/r2ods(3)Here w is the reflectance,r is the distance
between the light source Q and the point P, ~(r) is
the optical depth of PQ (see equation (2)), F(a) is
the scattering
phase function, and ~ is the angle
between the forward scattering
directionand thescattered ray (see Fig. l). The phase function is afunctionof particle radius, index of refraction,
and wave length.
The simplest example of the phasefunction is that for isotropic scattering,F(a)=constant.
For extremely
small particles such as
molecules of the air~ Rayleigh scattering theory isused; F(~) : K(I+cosLa) (K is a constant).
Forrelatively
small
particles
such as fog,Miescattering theory is used.
Mie theory
is very
complicated,so this paper uses
the following
experimental
approximation
[9]
for
foggy
atmosphere,F(a ):

In equation (1), I i is the intensityof light
arriving
at point
Pi on an object, and is also
subject to attenuation due to the distance from thelight source just as in equation(3)(i.e.,
Ii
dependson
exp(-T(r))/r2).
Calculationofilluminance due to point
sourceswith variable
angular intensity distributions is performed by themethodgiven
in Reference [2].
If the density ofparticles is uniform, equation (I) can be expressed
as follows

where ~
optical(61

is the intensity of light and r(r)thedepth from the reference plane (e.g., a

(~

Q

~'

ComputerGraphics,Volume21, Number4, July 1987

light source

wi ndo w

lumination~ . / i l

~.

~

s Pv
;hadow#v o i ume

/I

Fig. 2 Calculation of scattered light using
illumination and shadow volumes.
Fig. 3 Illumination volume for a parallel source.
window planeforthrough a window).

a

shaft

of

light

pouring

3- ILLUMINATION VOLUMES AND INTEGRATION SEGMENTS
Light sources treated here are point
sourceswith variable angular intensity distributions, andparallelsources.
In order tominimizetheapplicationof equation(I) (requiring a lengthy
calculation), illumination volumes are localized by
the following principles:(I) For point sources,
directionsof radiation are limited by reflectors
and lenses.
(2)
For
axisymmetricluminous
distributions,the illumination volume is defined
as a circularcone whose vertex
is the lightcenter,and whose central axis coincides with theillumination
axis(see Fig.2).
(3) Fornonaxisymmetrieluminous intensity distributions suchas headlights,
the illumination
volumeisanellipticalcone.
(4) For a parallel light source
entering through a window, the illumination
volumeis a prism of which the base is the window and thebase is swept in the direction of light(see Fig.

3).
Shadows
on particles
in the atmosphere are asimportant as those on the surfaces of objects; when
there are objects
in an illumination
volume,non-illuminated parts arise within it.
In order to detect shadows on objects, we use
shadow volumes [11] defined by the light source andthe contour lines of each polyhedron as viewed fromthe light source [2].
These shadow volumes
arealso used for detecting shadows in the atmosphere.
The shadows on each particle
along the llne
between the viewpoint and a point on an object(hereafter called the ray) can be obtained by using
the ray tracing algorithm (in this case, the ray
tracing algorithm
is applied
for each sampling
point
on
the
ray).
However,
this is very
time-consuming and gives rise to sampling errors.
Thefollowing
methodis more efficientandaccurate.
The shadow
segments
on the ray are
obtainedby using intersection points between the

ray and both illumination
volumesandshadow
volumes,which are pre-oalculated before intensity
calculation at each pixel.
No allasingproblems
caused by shadow calculation arise because shadow
regions are accurately calculated.
The segments in which scattered
light must be
calculated(e.g.,SlS 2
and sBs 4 in Fig.2) are
extracted by the following processes.
The following preproeessing is executed:i) Extract illumination volumes for each lightsource.
ii) Extract
polyhedra intersecting illumination
volumes.
iii) Extract shadow volumes formed by the lightsources and those polyhedra.
Then the following steps are executed for eachpixel on every scan line.
iv) Extract illumination volumes intersecting
theray, and calculate the intersection points. If
there is no intersection,no calculationisrequired for scattered light.
v) Calculate the intersections between the ray andall shadow volumes withineach illumination
volume.
vi) Extract illuminated segments on the ray.
(These segments
are visible
from the lightsource.)vii) Integrate intensities of the scatteredlighton the illuminated segments of the ray.
In step ii), all polyhedra are given a sphericalbounding
volume,and the polyhedra intersecting
illumination
volumesare extractedusing thesebounding spheres.
Steps i), ii) and iii) need not
be recalculated
if the viewpoint is changed, as
they are independent of the viewpoint.
They arealso used for detection of shadows on faces.
For stepsiv) and v), in order to calculateintersections between the ray and illumination
andshadow volumes
efficiently,
the regions visible

305

~

SIGGRAPH'87, Anaheim, July 27-31, 1987

I ~ ® ~ 1Q

quantitativeinvisibility,
as inhiddenline
elimination[12].
For example, in Fig. 5, as wetraverse the ray from Pu to Pi,
the quantitativeinvisibility
is incremented
each time the rayenters a shadow volume (Sz, s3, and s e in Fig.5),
and decremented
when it leaves(s4 ands5).
Segments
for whichquantitativeinvisibility iszero are visible (sis 2 and SsS 6 in Fig.5).

light sourcei ~-~'. il l u m i n a t i o n ~vo,ume/.:;:!~
shado~., I:.':'::<'
I:.':'::<'\/.:;:!~ shado~.,

4. DENSITY DISTRIBUTION OF PARTICLES~"~

In many cases,
the effectof spotlights
instudioor stage lighting is intensified by using
dry Ace.
Dry ice, smoke, dust,
and fog or haze
exhibit non-uniform particle densities.
We discusscalculationof
scattered
light
under
such
conditions.
In a dense volume,the extinctionof light
increases due to the increased optical depths
bothfrom the lightsource to a particle and from theparticle to the viewpoint.
On the other hand,becausethe intensityof scattered
light also
increases,
dense volumes in most cases are bright.
In orderto calculate
the intensity
attheviewpoint
Pu due to scattered light from a point P
on the ray, we must know the luminous intensityoflightsource Q in the direction of P, the particle
density at P, and the optical depths of QP and PPu.
Therefore, to take into account the distribution ofdensity,
we establish
a
model
for
density
distribution
and a method for calculating opticaldepths:SC r e e i lFig. 4 Pre-caleulation of intersections between the
ray and shadow and illumination volumes on
the screen.
from the screen are pre-calculated;the contour
lines of shadow and illumination volumes viewed
from the viewpoint are projected onto the screen as
shown in Fig. 4. Since the projected contour lines
are convex polygons, they can be extracted by meansof polygon scan conversion.
Then the illuminationand shadow volumes
intersecting with the ray areeasily obtained.
In this way, the regions in which
scattered light is calculated
are localized(seedottedparts in Fig. 4). When the light source isa point source, the illumination
volumeprojects
onto the screen as a triangle whose vertex is the
projected
light center.
The projectedshadow
volumeis a polygon(see shaded area in Fig. 4)
bounded by a triangle with vertex at the projected
lightcenter.
The shadow volumes intersecting theray can be extracted efficiently by means of thesebounding
triangles.
When an illumination volume
includes the viewpoint,
the illumination
volumecan't
be projected
onto
the screen,
and theintersection test must be executedin 3-D space,even though computation time increases.
The extractionof integration segments (i.e.,the visible parts when viewed from a lightsource)in step vi)
is perfomedby using the notion ofI) Modelling of density distribution
In this paper,
we do not consider
arbitrarydensity
distributions
over the whole space, butrather the simpler case in which certain regions of
the atmosphere are occupiedby smoke or fog of
uniform density.
One researcher has represented
volume density using a 3-D grid,with numericaldensities
for each element[6], but this method
requires a great deal of memory.
We use a model inwhich
the boundaries
of layersofdifferent
densities are represented by curved surfaces.
Forafoggy medium
due to dry ice,thedifference of densities between
the fog and theatmosphereis large, and because the fog is heavy,lightsourcesource$1llumination
viewpoint

J

(,

ation

-

~

/

Fig. 5 Extraction of integration segments for
scattered light.
(a view from the light source)

306

]

•

volumeFig. 6 Calculation of optical depths from the lightsource to a particle and from the particle
to the viewpoint.

(~@ComputerGraphics,Volume21, Number4, July 1987

it lies below the atmosphere.
The boundary surface
between the layers
can be expressedby mappingheightfunctions onto a horizontal plane (see Fig.
6). Properly speaking, this boundary surface must
be calculatedby physical simulation, but here weuse a Fourier series composed of a short sum ofsine
waves,
as used to represent
cloudsinreference[13].
For smoke, we map
thicknessfunctions onto a number of ellipses.

and polka dots.
The front-left spotlight exhibitsa soft edge, while the front-right exhibits a sharpedge.
In picture (b), we see shadow effects due toobstacles; the shape of the shaft of lightisclipped.
Pictures(c) and (d) depict dry ice andsmoke, respectively.

2) Calculation of optical depth
As shown in equation (2), the optical depth canbe obtainedby integrationof densities.
We
numericallyintegrate
the density along the rayfrom the viewpoint to the object.
By integratingin
this
direction,
theopticaldepthofintermediate
pointsiseasily available
for
calculating attenuation of scattered
light.
When
density
is sampled only at uniform intervals along
the ray, aliasing problems may occur.
To overcome
this problem,
we pre-caleulate the intersections
between the ray and the boundary
surface,detectthe intervals
in which sampling errors may occur,
and add intermediate
sampling
points in thoseintervals.

Pictures (e), (f) and (g) are examples of night
scenes outdoors.
Picture(e) shows a building
illuminatedby searchlights,
street lamps
andheadlights;
as analyticaldescriptionsof theluminous intensities of street lamps and headlights
are not available, actual numerical data are used(see
Appendix).
Note that due to the phasefunction,searchlights
appearbrighterwhen
directed toward the observerthan when directedaway.
Pictures (f) and (g) show the effect of thephase function (see equation (4)): the beams in (g)are more glaring than those in (f), where the phasefunction is neglected.
Picture (h) shows a shaftof light pouring through a window.
As shown in these examples, the proposed method
produces realistic images.

6. CONCLUSION
For the case with upper and lower layers ofdifferent
densities,
as shown in Fig. 6, the ray
can traverse at most one density boundary,so thecalculation can be simplified as follows: we assume
that the extinction
coefficientsfor upper andlower layers are q and o', and the boundary surfaceof the layer is g(x,y).
If the point Pi on an
object exists above g, equation (5) can be usedwithout any modification.
If not,the opticaldepth, • = o(r+s), of equation (5) must be replacedby• = a(r+s)-(o-a')(rl.sl)(see Fig. 6 for rl,
Sl).

5. EXAMPLES
Fig. 7
Pictureslighting:
animationshows thea special

shows examples of the proposed
method.
(a),(b) and (c) are examples for studioa scene of dancing
pencilsfrom theentitled"Feast of Light".
Picture (a)image illuminated by four spotlights
andeffects light that combines a daisy shape

This paper describes a methodof representing
illumination
effects
considering particles in theatmosphere.
The following
conclusions
canbestated from the results.
I) The proposed method could be used for lighting
designs for selecting and arranging
lights withconsiderationof atmosphericconditions, and for
estimating
illumination
effects
in studios
orstages.
2) Shafts of light due to light sources such as
spotlights
and headlights,
takingaccountofluminous
intensity
distributions, can be modelledby considering scattered light from particles
inthe atmosphere.
3) By considering illumination volumes for lightsources, rendering efficiency can be increased
bylocalizing
regions
in which light scattering iscomputed.
4) For shadows,
intersections
of
rays
with
illumination
volumes and with shadow volumes are

(b)

(a)
Fig. 7 Examples.

307

~

SIGGRAPH'87, Anaheim,July 27-31, 1987

(e)

(d)

(e)

(f)

(h)

(g)

Fig.7 Examples (continued).

308

O

~

ComputerGraphics,Volume21, Number4, July 1987

obtained efficiently by using theirprojected
contour lines. The shadow segments on the ray are
extracted using the hidden line technique, with the
ray viewed from the light source.

Acknowledgements
The authors wish to thank Bonnie Sullivan for
carefully reading and commenting on this paper.

REFERENCES

1) Warn, D.J.: Lighting Controls for Synthetic
Images, Computer Graphics, Voi. 17, No.3 (1983)
pp.13-21.
2) Nishita, T., Okamura,I. and Nakamae,E.:
Shading
Models for Point and Linear Sources, ACM Trans.
on Graphics, Vol.4, No.2 (1986) pp.124-146
3) Nishita,T.
andHakamae,E.:
Half-Tone
Representation of 3-D Objects Illuminated by
Area Sources or Polyhedron Sources, IEEE, Proc.
of COMPSAC (1983) pp.237-241.
4) Nishita, T. and Nakamae, E. : Continuous Tone
RepresentationofThree-Dimensional Objects
Illuminated by Sky Light, Computer Graphics,
Voi,19, No. 3 (1986) pp.23-30.
5) Blinn, J. F.: Light Reflection Functions for
Simulation of Clouds
and
Dusty
Surfaces,
Computer Graphics,
Voi. 16, No.3 (1982) pp.2129.
6) KaJiya, J.T. and Von Herzen, B.P.: Ray Tracing
Volume Densities,
Computer Graphicsj
Voi.18,No.S (1984) pp.165-174.
7) Max,~, L.: Atmospherical Illumination andShadows, Computer Graphics, Vol.20, No.4 (1986)

pp.117-124.
8)

IES : IES Lighting Handbook Reference Volume,(1981) p.6-6.
9) Gibbons,M.G.: Radiation Received by UncollimatedReceiver from a 4~ Source, J. Opt. Soe. ofAmerica, Vol.~8, No.8 (1958) pp.550-555.
10) Nakamae, E., Harada,K,, Ishizaki,T, Nishita,T.:A
Montage; The Overlaying of the Computer
Generated Images onto a Background Photograph,
Computer Graphics, Vol.20, No.3 (1986) pp.207214.
11) Crow, F. : Shadow Algorithms for Computer
Graphics, Computer Graphics, Voi.11, No.3 (1977)
pp.242-248.
12) Nakamae, E. and Nishita, T. : An Hidden Line
Eliminationof
Polyhedra,
Information
Processing in Japan, Voi°12 (1972)
pp.134141.
13) Gardner, Y.G.
: Visual Simulation of Clouds,
ComputerGraphics,
Voi.19,No.3(1985)
pp.297-303.

Fig. 8 Illustration of angles for non-axisymmetricluminous intensity.

APPENDIX
Calculation of Luminous Intensity Characteristics
I)

Axisymmetricintensity

distributionofluminousIn this case the intensity can be expressed by afunction of angle 0 from the illumination axis (seeFig. 8). In most cases the distribution curves of
the intensity can be expressed as a polynominal of
cose.
Here we consider the sharpness of theilluminated edge.
In ordertocontroltheillumination range, the expression l(e)= Io(cose) c
has been proposed [1] (For this expression,
theillumination
range decreases as c increases).
However, because sharp edges cannot be produced,
and because of the computational expense of the
exponent in the cosine term, this function is notwell suited to modelling spotlights. Therefore we
propose the following function for spotlights:
I(e) : I0{ (1-q)(cose-cosy)/(1-cos¥) + q }, (7)where I0 is the luminous intensity at e=0, ~ is thebeam spread, and q is a parameter which controls
the sharpness of the edge (q=1: very sharp, q=0:
very soft) (see Fig. 2 for 7).
2) Non-axisymmetric distribution of luminous
intensity
In this case the luminous intensity is expressedby I(~,~), where ¢ is the revolution angle from thereference plane which includes the illumination
axis (refer to Fig. 8). In the figure the lightcenter is Qe, the reference point is Q0, Qc is apoint on the illumination axis, and P is anarbitrary point. Then cos~ can be obtained by the
inner product of the normal vectors of triangles,
QeQcQ0 and QeP Qc.
First we discuss luminous intensity expressed byfunctions.
For example, luminous intensities forthe daisy shape and for polka dots like those
produced by a mirror ball in a disco are calculated
as follows.
In these sources, the intensity I(e,$)

309

~

S I G G R A P H '87, Anaheim, July 27-31, 1987

can be expressed as the product of I(e) of equation(7) and the following function k(¢).
k(@)= max{ p + ( 1 - p ) [ 2 m o d ( @ , A ¢ ) - n ¢ ) ] / ~ ¢ ,

0 }

For
the
polka-dot
pattern,
intensity can be obtained by adding

the luminous
the variationin the ~ d i r e c t i o n to equation (9).

That i s ,(8)I(e,¢)

where the symbol "mod" means modulo, A@ and p (p<1)
are the period of luminous intensity and the ratioof the minimum intensity to the maximum intensity,
respectively: When p is small (including negative),the radius of the polka-dots becomes small.
For
the daisy shape, the distribution of luminous
intensity is given byI(e,~)= Z(e)k(~,).

= "r(e3k(¢3k(e)v

,

where

v=mod([@/A~]+[e/A0],2) (i.e., v=O or v = l ) : v
is used to decrease the number of polka dots.

For arbitrary distributions, numerical luminous
intensity data should be used, for example, every 5
degreesone and¢.
Intervening luminous
intensities are calculated by linear interpolation.
For headlights, however, more detailed data, forexample every 2 degrees, are necessary because of
the radical change of intensities (see Fig. 93.

(9)

2~

0.
2'
4,
60
8o
. . . . . . . . . .

.~.

.,.~

o

o

,

,~

u

z

'4-

b-

~-

.!.U" 1 2 ~ 1 4 °16 ° I , 8 °FAg. 9 Distribution curve of luminous intensity for a headlight.

310

(103

20 ° 22 °