University of Southern California Dissertations and Theses

Techniques for efficient cloud modeling, simulation and rendering

(USC Thesis Other)

Techniques for efficient cloud modeling, simulation and rendering

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content TECHNIQUES FOR EFFICIENT CLOUD MODELING,
SIMULATION AND RENDERING
by
Bei Wang
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Ful¯llment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(Electrical Engineering)
May 2008
Copyright 2008 Bei Wang
Dedication
This dissertation is dedicated to my family, my advisor and all friends who helped me a
lot on the tough journey...
ii
Acknowledgements
I would like to express my deepest gratitude to my advisors, Professors C.C.Jay Kuo. I
would never ¯nd better role model and advisor who provides such high quality academic
training. His guidance has been very important to me, and I feel that I have now grown
a great deal. The ¯ve and half years I spent working with him has been truly a cherished
experience I deeply appreciate in my life. It has been a great pleasure to work on this
topic,andseeitgrowsinthesepastyears. Iamveryhappytohavejoinedinthisjourney.
I am also very grateful to my qualifying examination and thesis defense committee
members, Professors Krishna Nayak, Ulrich Neumann, Richard Leahy, B.Keith Jenkins,
and Karen Liu, for their guidance on and encouragements for my research.
I am deeply thankful to my friends in Los Angeles, for their emotional support and
care. They include: Iris Shen, Kelly zhu, Qi Liu, Min Chen, Han Wang, Chunli Ren,
Chiao Wang, Chao Chen, Cong Peng, Chuping Liu, Feilong Liu, Qi Cai, and Maple Qiu.
I wish to thank my colleague and good friends, Jingliang Peng, May Kuo and Jessy
Lee who has always been there when I was stuck. Also, many thanks go to my group-
matesandfriends,whoincludes: IlyaEckstein,NacoChang,KevinChou,RonaldChang,
RogerChang,AmyLee,AthenaHuang,YoungminKwak,YuShui,AngelDai,Shih-Hung
iii
Lee,JiayingLiu,KeithZhang,DahuaXie,YuHu,QingLi,SteveZhou,AlexHsu,Gracia
Wang... I truly enjoyed the time working with these colleagues in our cozy group.
I wish to thank Dr. Chia-huang Yeh, Dr. Maverick Shih, who have been so kind,
encouraging and cooperative in working in the projects during my Ph.D. training, but
not been mentioned above. Without their their support, ¯nishing my dissertation would
be much tougher.
Lastbutnotleast,letmegivemyprofoundthankstomyhusbandYanLu,myparents,
Chuanlong Wang and Cuiyun Zhang, my parents-in-law, Bingquan Lu and Beihua Yang,
my elder brother and sister-in-law, Rui Wang and Qian Wu, my sister-in-law and her
husband Yao Lu and Yan Xia. Their constant love has always been the most important
thing to me in this work. This Dissertation is dedicated to my family. It is because of
them that I ¯nally made it through the di±cult time, and it is because of them that I
know I will never give up on my endeavors.
iv
Table of Contents
Dedication ii
Acknowledgements iii
List Of Tables viii
List Of Figures ix
Abstract xii
Chapter 1: Introduction 1
1.1 Signi¯cance of the Research . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Background of the Research . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.1 Cloud Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.2 Cloud Illumination . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Contributions of the Research . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4 Outline of the Dissertation. . . . . . . . . . . . . . . . . . . . . . . . . . . 9
Chapter 2: Research Background 10
2.1 Mathematics of Cloud Dynamics and Lighting. . . . . . . . . . . . . . . . 10
2.1.1 Classi¯cation of Clouds . . . . . . . . . . . . . . . . . . . . . . . . 10
2.1.2 Cloud Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.2.1 Ideal Gases . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.1.2.2 Parcels . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.2.3 Lapse Rate . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.2.4 Phase Transition . . . . . . . . . . . . . . . . . . . . . . . 16
2.1.2.5 Latent Heat . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.1.3 Cloud Illumination . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.1.3.1 Albedo of Cloud . . . . . . . . . . . . . . . . . . . . . . . 17
2.1.3.2 Single Scattering and Multiple Scattering . . . . . . . . . 18
2.1.3.3 Optical Distance and Phase Function . . . . . . . . . . . 19
2.1.3.4 Light Intensity Calculation . . . . . . . . . . . . . . . . . 21
2.2 Similarity Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.2.1 Objective and De¯nition . . . . . . . . . . . . . . . . . . . . . . . . 22
2.2.2 Similarity Theory and Dimensionless Analysis . . . . . . . . . . . . 24
v
2.3 Previous Work on Cloud Simulation and Rendering . . . . . . . . . . . . . 32
2.3.1 Physical-based Dynamic Simulation of Clouds . . . . . . . . . . . . 32
2.3.2 Graphic Tools for Cloud Modeling . . . . . . . . . . . . . . . . . . 36
2.3.2.1 Polygon . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.3.2.2 Procedural Noise . . . . . . . . . . . . . . . . . . . . . . . 36
2.3.2.3 Textures Sprite . . . . . . . . . . . . . . . . . . . . . . . . 37
2.3.2.4 metaball . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.3.2.5 Particle System . . . . . . . . . . . . . . . . . . . . . . . 38
2.3.2.6 Voxel Grid . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.3.3 Cloud Rendering . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.3.3.1 Discrete Ordinates . . . . . . . . . . . . . . . . . . . . . . 41
2.3.3.2 Spherical Harmonics . . . . . . . . . . . . . . . . . . . . . 41
2.3.3.3 Finite Element . . . . . . . . . . . . . . . . . . . . . . . . 43
2.3.3.4 Monte Carlo Integration. . . . . . . . . . . . . . . . . . . 45
2.3.3.5 Line Integral . . . . . . . . . . . . . . . . . . . . . . . . . 47
Chapter 3: A Similarity Approach to Particle-based Cloud Simulation 50
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.2 Thermal Simulation by Similarity Approach . . . . . . . . . . . . . . . . . 52
3.3 Thermal Movement Modeling by Particles . . . . . . . . . . . . . . . . . . 56
3.3.1 Parameters in Particle Systems . . . . . . . . . . . . . . . . . . . . 56
3.3.2 E®ect of Water Vapor Life Time . . . . . . . . . . . . . . . . . . . 57
3.3.3 Cloud Shape Control . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.3.4 Wind E®ect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.4 Phase Transition Modeling by Voxel Dynamics . . . . . . . . . . . . . . . 60
3.4.1 Water Continuity and Phase Transition . . . . . . . . . . . . . . . 61
3.4.2 Fractal Structure and Entrainment . . . . . . . . . . . . . . . . . . 62
3.4.3 E±cient Update of Water Phase Change . . . . . . . . . . . . . . . 64
3.4.4 Continuous Density Calculation . . . . . . . . . . . . . . . . . . . . 64
3.4.5 Other Implementation Details . . . . . . . . . . . . . . . . . . . . . 65
3.5 Cloud Illumination and Rendering . . . . . . . . . . . . . . . . . . . . . . 67
3.6 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
Chapter 4: A Decoupled-PDE Approach to Cloud Simulation with Volu-
metric Modeling 79
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.2 Mathematical Background . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.3 Stam's Stable Fluid Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.3.1 Helmholtz-Hodge Decomposition . . . . . . . . . . . . . . . . . . . 84
4.3.2 Navier-Stoke Solver . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.3.2.1 External Force . . . . . . . . . . . . . . . . . . . . . . . . 86
4.3.2.2 Advection . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.3.2.3 Di®usion . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.3.3 Force Extension . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
vi
4.3.3.1 Virtual Temperature Force . . . . . . . . . . . . . . . . . 89
4.3.3.2 Coriolis Force . . . . . . . . . . . . . . . . . . . . . . . . 89
4.3.3.3 Vorticity Con¯nement . . . . . . . . . . . . . . . . . . . . 90
4.4 Decoupled-PDE Cloud Simulation . . . . . . . . . . . . . . . . . . . . . . 91
4.4.1 Turbulent Plume . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.4.2 Decoupled Cloud Simulation . . . . . . . . . . . . . . . . . . . . . 94
4.4.3 Implementation Details . . . . . . . . . . . . . . . . . . . . . . . . 96
4.4.3.1 Vertical Water Vapor Distribution . . . . . . . . . . . . . 96
4.4.3.2 Coriolis Force . . . . . . . . . . . . . . . . . . . . . . . . 98
4.4.3.3 Wind E®ect . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.4.3.4 Interpolation Between Key Cloud Layers . . . . . . . . . 99
4.4.3.5 Boundary Condition . . . . . . . . . . . . . . . . . . . . . 99
4.4.4 Cloud Rendering . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
4.5 Synthesized Cloud Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
Chapter 5: Real-Time Hierarchical Cloud Simulation based on Similarity
Approach 110
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
5.2 Requirements of Real-Time Simulations . . . . . . . . . . . . . . . . . . . 111
5.3 Dynamics of Convective Thermal . . . . . . . . . . . . . . . . . . . . . . . 114
5.4 Dynamics of Thermal Inner Motion . . . . . . . . . . . . . . . . . . . . . . 117
5.4.1 Vorticity and Vortex Ring . . . . . . . . . . . . . . . . . . . . . . . 117
5.4.2 Axis-Symmetric Flow . . . . . . . . . . . . . . . . . . . . . . . . . 118
5.4.3 Relationship between Thermal and Vortex Ring. . . . . . . . . . . 121
5.5 Hierarchical Cloud Simulation . . . . . . . . . . . . . . . . . . . . . . . . . 122
5.5.1 Simulation Framework and Parameters. . . . . . . . . . . . . . . . 122
5.5.2 Velocity-Vorticity Conversion . . . . . . . . . . . . . . . . . . . . . 124
5.5.3 Rotational Velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
5.5.4 Phase Transition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
5.5.5 Thermal Interaction . . . . . . . . . . . . . . . . . . . . . . . . . . 126
5.6 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
5.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
Chapter 6: Conclusion and Future Work 134
6.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
Bibliography 139
vii
List Of Tables
4.1 Discrete approximations to several 2D gradient operators. . . . . . . . . . 82
4.2 Acceleration contributions from di®erent factors. . . . . . . . . . . . . . . 94
viii
List Of Figures
2.1 A real cloud scene to illustrate two types of light interaction. . . . . . . . 18
2.2 The angle between the incident direction and the exitance direction Á. . . 20
2.3 The light transport phenomenon in clouds. . . . . . . . . . . . . . . . . . 22
2.4 Illustration of the "top-hat" pro¯le (left) and the Gaussian pro¯le (right). 28
2.5 The incremental volume over which the vertical momentum equation is
integrated. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.1 The cloud modeling framework which consists of the two-stage cloud for-
mation process and the cloud rendering process. . . . . . . . . . . . . . . 51
3.2 The solution for radius (R), height (z), buoyancy (B) and vertical velocity
(w) of a thermal in a uniform stably strati¯ed °uid. . . . . . . . . . . . . 56
3.3 The histogram of condensation time for cloud formation. . . . . . . . . . . 59
3.4 Comparison of cloud shapes controlled by di®erent initial water vapor dis-
tributions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.5 Comparisonofcloudshapes(a)withoutand(b)withthewinde®ectduring
its formation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.6 The potential direction of the substructure. . . . . . . . . . . . . . . . . . 63
3.7 Comparison of clouds (a) with and (b) without the fractal structure. . . . 64
3.8 Illustration of the splat texture. . . . . . . . . . . . . . . . . . . . . . . . . 67
3.9 Simulation of the cloud formation process, where steps (a)-(c) show how a
cloud is formed along time. . . . . . . . . . . . . . . . . . . . . . . . . . . 71
ix
3.10 Three sampled results of a cloud formation process.. . . . . . . . . . . . . 72
3.11 Clouds generated with di®erent buoyancy frequency values: (a) N = 0:5
and (b) N =0:01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.12 Cloudsgeneratedbydi®erententrainmentcoe±cients: (a)®=0:5and(b)
®=0:03. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
3.13 Clouds rendered under di®erent lighting condition . . . . . . . . . . . . . 76
3.14 A scene of multiple clouds. . . . . . . . . . . . . . . . . . . . . . . . . . . 77
3.15 The close up of multiple clouds. . . . . . . . . . . . . . . . . . . . . . . . . 77
3.16 The same cloud scene at sunset.. . . . . . . . . . . . . . . . . . . . . . . . 78
4.1 Illustration of an uniform grid for °uid simulation. . . . . . . . . . . . . . 82
4.2 The °owchart of the Navier-Stoke solver. . . . . . . . . . . . . . . . . . . . 86
4.3 Theadvectioninsimulation:(a)thevelocity¯eld,(b)thepathofaparticle
and (c) to ¯nd the velocity ¯eld backwards from t+4t to t. . . . . . . . 87
4.4 Velocity Di®usion in °uid simulation . . . . . . . . . . . . . . . . . . . . . 88
4.5 The height dependence of several dimensionless quantities of a turbulent
plume: horizontal radius R, vertical velocity w and buoyancy B. . . . . . 92
4.6 The development of a turbulent plume: the early stage (left) and the late
stage (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.7 The proposed decoupled cloud simulation framework. . . . . . . . . . . . . 95
4.8 The plateau-like kernel for water vapor distribution. . . . . . . . . . . . . 97
4.9 A horizontal cloud layer comparison (a) with and (b) without the Coriolis
force. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.10 OpenGL implementation of modi¯ed Behren's Algorithm . . . . . . . . . 103
4.11 Cloud formation with a grid of size 64£64. . . . . . . . . . . . . . . . . . 106
4.12 A synthesized cumulus cloud with a grid of size 128£128. . . . . . . . . . 107
4.13 Examples of cumulus cloud with di®erent parameter sets. . . . . . . . . . 108
x
4.14 A scene with multiple cumulus clouds. . . . . . . . . . . . . . . . . . . . . 109
4.15 A scene with multiple cumulus clouds at sunset. . . . . . . . . . . . . . . 109
5.1 Two °uid simulation approaches: (a) the Lagrangian approach, (b) the
Eulerian approach, and (c) the 3-D Eulerian approach. . . . . . . . . . . . 113
5.2 Successive outlines of thermals traced from photograph (top) and a graph
of z
2
against t (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
5.3 Thevelocitydistributionofparticlesinsideathermal,wherethehorizontal
(left) and the vertical (right) velocities are expressed as multiples of the
vertical velocity of the front of the thermal. . . . . . . . . . . . . . . . . . 117
5.4 Illustration of a 2-D velocity ¯eld and the associated 3-D vorticity. . . . . 118
5.5 Avortexringwheretheredarrowindicatethedirectionofaxis(left)anda
two-dimensional slice of the velocity ¯eld induced by the vortex ring (right).119
5.6 Schematic sketch of a circular vortex ring. . . . . . . . . . . . . . . . . . . 120
5.7 The hierarchical cloud simulation framework. . . . . . . . . . . . . . . . . 123
5.8 The distribution of thermal velocities. . . . . . . . . . . . . . . . . . . . . 124
5.9 The rendering result of a sample thermal consisting of 150 particles. . . . 125
5.10 (a) The standard Euler forward method and (b) the rotational velocity
method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
5.11 The birth of a new thermal in the hierarchical cloud simulation. . . . . . . 127
5.12 Cloud formation based on Hierarchical structure . . . . . . . . . . . . . . 129
5.13 Statistic FPS for hierarchical cloud simulation. . . . . . . . . . . . . . . . 130
5.14 Statistic FPS comparison for hierarchical cloud simulation with and with-
out inner motion turned o®. . . . . . . . . . . . . . . . . . . . . . . . . . . 131
5.15 Cauli°ower cloud formation using hierarchical cloud simulation. . . . . . . 132
5.16 A scene with multiple cumulus clouds. . . . . . . . . . . . . . . . . . . . . 133
5.17 A scene with multiple cumulus clouds at sunset. . . . . . . . . . . . . . . 133
xi
Abstract
Realistic cloud simulation and rending ¯nd applications in many 3D applications such
as high-quality 3D ¯lm production and realistic 3D °ight simulation. However, the in¯-
nite variation in the cloud shape and appearance make cloud simulation and rendering
a challenge in terms of complexity and realisticity. The major focus of this work is on
e±cient and realistic cloud simulation. We adopt a similarity approach to describe the
characteristics of simple °uid °ow without rigorously solving the governing PDEs. Based
on the similarity solution of turbulent thermals, plumes and vortex rings in stable strati-
¯ed °uids or unstable environments, we propose three methods for cloud simulation and
modeling.
First, a two-stage cloud simulation method is proposed to model the cloud forma-
tion process. At the ¯rst stage, instead of solving PDEs for the velocity pro¯le, cloud
formation is obtained using a lower level building block; namely, particle simulation.
New attributes of particles are incorporated to meet the modeling requirement, e.g., the
condensation time and density of particles. For particle motion, the characteristics of
the convective process is captured by the thermal similarity solution consistent with the
governing Navier-Stoke equation. At the second stage, we consider the process of being
xii
condensed into the visible water droplet. The voxel is used to store the water vapor den-
sity and the water droplet density. Due to the release of latent heat, we choose to have
a small scale convective process. Finally, the cloud droplet density is rendered using a
two-pass rendering method and multiple scattering to result in the high albedo property
of the cloud.
Then, a decoupling cloud simulation method that divides three spatial dimension
into independent horizontal and vertical dimensions is proposed. With the assumption
that the radial pro¯les at all heights are similar, we choose several key frame layers only
for actual 2D calculation, interpolate other layers from them and add some turbulence.
Along the horizontal dimension, the movement is driven by heat di®erence. This virtual
temperature di®erence results in a di®erent water vapor distribution. To model the
cloud formation from a plume, we assume that the water phase change occurs at its
late development stage, where the plume spreading between levels of buoyancy vanishes
and the vertical velocity is equal to zero. At such a stage, the water vapor will be
su±ciently accumulated to experience the water phase transition. The 2D cloud layer is
then simulated by the Navier-Stoke equations. Besides the momentum equation coming
from the virtual temperature di®erence, one more force term is introduced in our model,
i.e. the Coriolis force. Finally, water droplets are rendered using a 3D texture rendering
technique. By separating the horizontal and vertical couplings in the traditional 3D
simulation, wecanreducethecomputationalcostyetcapturetheessentialcharacteristics
of the cumulus cloud.
Finally, we introduce a novel hierarchical cloud simulation method. The general
shape and movement of clouds in the simulation model are governed by the similarity
xiii
of thermals. Furthermore, their relative position to the camera is calculated. If they
are not hidden by other thermals or they are in the boundary of a cloud, we turn on
thermal'sinnermotionwhichbehaveasaweakvortexring. Wesimulatethermal'sinterior
dynamicsusingtheaxis-symmetric°ow,whichprovidesasimpleschemetodeterminethe
velocityatanyarbitrarylocationinsidethe°ow. Withtheaxis-symmetric°ow,weadopt
the Lagrangian approach for thermal's interior motion, where a bunch of water droplet
particles inside the thermal are updated at each time-step, generating the desired cloud
detailstructure. Furthermore,weconsiderthephasetransitionwithasimpli¯edstrategy.
With the entrainment mass whose coe±cients are chosen by users, we can increase the
water droplet's density, thus its size and transparency for rendering. A method for cloud
formation due to the latent heat release from the condensation is also considered.
xiv
Chapter 1
Introduction
Clouds are pictures in the sky
They stir the soul, they please the eye
They bless the thirsty earth with rain,
which nurtures life from cell to brain |
But no! They're demons, dark and dire,
hurling hail, wind, °ood, and ¯re
Killing, scarring, cruel masters
Of destruction and disasters
Clouds have such diversity |
Now blessed, now cursed,
the best, the worst
But where would life without them be?
Vollie Cotton
1
1.1 Signi¯cance of the Research
Ever since 1802, people began to know clouds better from the essay of Luke Howard, the
"GodfatheroftheCloud",ontheclassi¯cationoftheformsofclouds. Hisclassi¯cationof
clouds is still in use today, which classi¯es clouds into three main types: cumulus, stratus
andcirrus. Longbeforethatessay'scomingout,peoplealreadyviewedthecloudasahint
of weather variation, but it is Howard that had explained the formation and evolution
of clouds more clearly and scienti¯cally. Thereafter, lots of scientists have investigated
cloud dynamics and the atmospheric science. Nowadays, fairly accurate weather forecast
can be made based on extensive research on clouds so far. Realistic cloud simulation and
rendering is important to many computer applications. For instance, high-quality and
e±cient rendering of clouds is the key to realistic °ight simulation in °ight games or pilot
training softwares and computer-generated 3D clouds help make rich visual e®ects at a
relatively low budget in movies.
Cloudsimulationandrendering,albeithelpfulinmany¯elds,placesahighdemandon
computational resources due to the underlying complexity in cloud dynamics and cloud
illumination interaction. To make it practical in computer simulation, proper simpli¯ca-
tion in the dynamics and illumination of clouds is necessary. Most existent techniques
in computer simulation and rendering of clouds have to take a tradeo® between com-
putational complexity and visual quality, based on the available hardware resources. In
this research, we are concerned with simulation, modeling and rendering of clouds using
computers.
2
1.2 Background of the Research
The topic of cloud modeling and rendering has been researched for about two decades
in the ¯eld of computer graphics. With the advances of available computing resources,
researchers have kept improving the techniques in cloud simulation, modeling and ren-
dering. Some related concepts and research work are brie°y reviewed below.
1.2.1 Cloud Simulation
In general, the topic of cloud simulation involves two major parts - cloud modeling
and cloud dynamics simulation. The cloud modeling method generally determines the
methodology for cloud dynamics simulation. There are two main categories of cloud
modeling methods: the procedural-based approach and the physics-based approach.
By procedural cloud simulation, one generates clouds based on visual properties of
clouds. The looking of a cloud is typically amorphous, randomized and rich in swirling
and details, especially around cloud boundaries. Due to the fractal property of clouds,
procedural methods may use Fourier series to generate naturally-looking cloud texture
or use certain noise functions as the basis to model clouds. By utilizing a sequence of
di®erentscalesofnoiseandturbulencefunctions,clouddetailsindi®erentresolutionscan
be simulated with those procedural methods. The procedural cloud simulation approach
has several advantages as described below.
² E±cient ¯ne detail presentation
Typically, there is no limitation on the simulation resolution for a procedural cloud
simulator. In contrast, due to its computational complexity, a physics-based cloud
3
simulator cannot simulate clouds beyond a certain level of details. In a certain
sense, aproceduralcloudsimulatorhasthepotentialof\in¯nite"-detailsimulation,
with the only limitation of numerical precision provided by the underlying com-
puting devices. As a result, procedural methods are good for non-repetitive cloud
simulation in a very large area.
² Easy parameterizations
Di®erent cloud types or di®erent shapes of the same cloud type can be easily con-
trolled by parameters in the cloud simulator.
² Solid texturing
The 2D or 3D solid texturing can be used to avoid the di±culty in surface texture
mapping.
However, the procedural cloud simulation has the following disadvantages.
² Less °exibility in control
Although parameterization is an advantage of procedural methods, ¯nding the cor-
respondencebetweenmodelparametersandtheresultantvisualexpectationisoften
conducted by trial and error, and these structures are essentially random and di±-
cult to modify or control interactively in the creation of di®erent cloud types.
² Non-realistic animation
Procedural methods are built upon non-physics-based modeling so that their °exi-
bilityincloudanimationsimulationislimited. Typicalproceduralmethodsattempt
4
to perturb noise functions and/or texture coordinates, while resultant visual e®ects
are usually not satisfyingly to human eyes.
² Aliasing
Since a procedural method uses the summation of Fourier series or noise functions,
aliasing is a common problem that has to be carefully taken care of.
Since the purpose of computer graphics is to re-build the world on computers, it is
natural to follow the same rules as clouds do in the physical world. Therefore, another
importantbranchofcloudsimulationisphysics-based. Typically,physics-based°uidsim-
ulation involves partial di®erential equation (PDE) solving, which leads to high accuracy
and amorphous °uid shapes due to the non-linearity of advection. As compared with the
procedural method, the physics-based cloud simulation has the following advantages.
² Accurate and realistic animation
The way in which a physics-based simulation method generates a cloud is just like
thewayanexperimentismadeinalab. Thatis,itattemptstorepeatwhathappens
inthenaturebyfollowingcertainphysicsrules. Asaresult,theresultantanimation
with a physics-based cloud simulator closely matches what occurs in the real world.
² More °exibility in control
The type of a cloud form depends on parameters such as di®erent forces applied to
it, the density of the water vapor, etc. With these parameters directly associated
with physical events, it is much easier and °exible to simulate what happens in the
real nature. For instance, we may simulate how the cloud shape changes according
5
tothevariationinthedirectionandthemagnitudeofawindforce, whichisusually
di±cult to achieve with procedural methods.
However, a physics-based cloud simulator demands a high computational complexity.
Due to the PDE solution, the computational cost of a physics-based cloud simulator is
proportionaltothegridsizeusedinthePDEsolution. The¯nerthegridsizeis,themore
the simulation details will achieve at a higher computational cost. Although computer
graphics hardware can be exploited to accelerate the computation, this is still the major
bottleneck of physics-based cloud simulation.
1.2.2 Cloud Illumination
Cloud illumination is also a big challenge because of the complex interaction between
the light and clouds. A cloud is composed of millions of water droplets. Clouds ab-
sorb little light energy, instead, each water droplet re°ects or scatters nearly all incident
light. To avoid the complexity of cloud radiometry, approximate self-shadowing and
single-scattering have been used to approximate cloud illumination. However, accurate
and realistic generation of cloud images requires simulation of multiple light scattering
within cloud volumes. The complexity of multiple scattering makes exhaustive simula-
tion almost impossible, and simpli¯cation and approximation has been used to reduce
the computational complexity.
6
1.3 Contributions of the Research
Researchhasbeenconductedtoaddressthechallengesmentionedin thelastsection, and
several major contributions have been made in this work. They are detailed below.
² A new physics-based cloud modeling approach is proposed with improved compu-
tational e±ciency. First, we make a more realistic modeling of the cumulus cloud
formation process. With the atmosphere represented as a strati¯ed °uid, the cloud
formation process is in essence a process of turbulent plume or thermal convec-
tion. The cumulus cloud formation is closely related to an important atmospheric
property, knownastheBuoyancyfrequency, whichisameasurementofatmosphere
strati¯cation. The formation process is described below. First, an air parcel gets
heatedfromthegroundandraisedupduetobuoyancy. Asitpassestheleveloffree
convection (LFC), the vertical motion will dominate and the air parcel continues
to lift up. As it rises up, it expands and cools down, until it reaches an equilibrium
level(EL).Fromthenon,itceasesfurthermovingupandbeginstooscillatearound
the EL and move laterally, forming the cloud anvil.
² A similarity approach is proposed in cloud simulation to avoid the high compu-
tational cost associated with the PDE solution. Almost all physics-based cloud
simulators view the cloud as a simple type of °uid. For a simple °uid °ow, we are
able to describe the characteristics of °uid dynamics by inferring various constant
parameters associated with the °ow under certain conditions. Based on this obser-
vation, we propose to use the similarity approach in cloud simulation, which not
7
onlycapturesthekeycharacteristicsofthephysicalconvectiveprocessbutalsosim-
pli¯es the computation by avoiding the PDE solution as done in the physics-based
approach.
² A two-stage particle simulation is used to model the cumulus cloud formation,
whereparticledynamicsiscausedbytheturbulentthermale®ectinstablestrati¯ed
°uids. Aspeci¯cwaterphasechangeprocessis¯rstproposedtotheparticlesystem.
Then, with the propertyof a particle rendering system, we can achieve better cloud
illumination and rendering. With this proposed two-stage particle simulation, we
are able to achieve scalable solutions to meet variant quality/complexity demands.
² An approach of decoupling the cloud model on voxel grids is proposed based on the
characteristics of turbulent plumes in stable strati¯ed °uids and the physical fact
that the behavior along vertical and horizontal directions is signi¯cantly di®erent
in atmospheric science. In the proposed approach, the water vapor distribution
along the vertical dimension is determined by the plume similarity solution, while
the horizontal dimension is determined by solving the Navier-Stoke equation for
key cloud layers and interpolating other cloud layers. Finally, volume rendering
techniques are utilized to visualize resultant clouds.
² A novel real-time hierarchical cloud simulation is proposed by exploiting the fact
that a cumulus cloud is composed of a bunch of thermals that exhibit the self-
similarity property. In the proposed method, the coarse level cloud simulation is
derived by considereing thermal's movement as a whole. Then, motion dynamics
is used to generate details of cloud boundaries according to the relative position of
8
thermals and the camera. By applying the axis-symmetric °ow to thermal's inner
motion, we are able to achieve simple yet visually convincing real-time simulation.
1.4 Outline of the Dissertation
The thesis is organized as follows. In Chapter 2, the atmospheric knowledge about cloud
dynamics and illumination and the mathematics background for the similarity approach
are introduced, and related work on cloud modeling, simulation and illumination is sum-
marized and categorized. Cloud simulation based on two-stage particle modeling is de-
scribed in Chapter 3. Cloud simulationbased on a horizontal-vertical decoupling method
is discussed in Chapter 4. Real-time hierarchical cloud simulation using self-similarity
is presented in Chapter 5. Finally, the dissertation is summarized and future research
directions are pointed out in Chapter 6.
9
Chapter 2
Research Background
2.1 Mathematics of Cloud Dynamics and Lighting
Realisticvisualsimulationofcloudsmainlyfocusesontwoaspectsofcloudsinnature: dy-
namics and radiometry. Cloud dynamics simulates the process of cloud formation, cloud
movement and clouds interaction. Cloud radiometry simulates how sky light interacts
with clouds. Both these two topics correspond to highly complex processes in the real
world, thus very di±cult for e±cient simulation and presentation by computers without
simpli¯cation. In this chapter, we present the necessary mathematics and physics princi-
ples for the cloud simulation and radiometry, together with simpli¯cation techniques for
the purpose of e±cient cloud simulation.
2.1.1 Classi¯cation of Clouds
Inordertosimulatearbitrarycloudscenes,wemust¯rstbeawareofsomebasicknowledge
about clouds of di®erent types. Fortuitously, the standard classi¯cation of cloud types is
basedontheirvisualappearance[31]. Thethreefundamentalclassesofcloudsare: cirrus
10
(\curl of hair"), stratus (\layer") and cumulus (\heap"). Cirrus clouds are wispy clouds
athighaltitudes(5»13km). Stratuscloudsarelayercloudswithnodistinctdetail,lying
at low altitudes (0» 2km). Cumulus clouds are heap clouds, also lying at low altitudes
(0 » 2km). Combinations of the basic class names are used to describe clouds with
combined characteristics and altitudes of the basic types. \Cirrostratus" refers to a high
layer of cirrus clouds. \Stratocumulus",\altocumulus" and \cirrocumulus" refer to low,
intermediate, and high broken cloud layers, respectively. \Altostratus" refers to a thick
layer cloud at intermediate altitudes (2 » 8km). \Nimbostratus" refers to a dark layer
producing rain or snow, and \cumulonimus" refers to a large cumulus cloud producing a
shower.
2.1.2 Cloud Dynamics
The dynamics of cloud formation, growth, motion and dissipation are complex. In order
to e±ciently simulate these complex cloud processes with a computer, it is necessary to
¯rst understand these dynamics before further simpli¯cation and approximation can be
made.
First of all, we need to understand the concept of convective motions, by which we
mean all motions that can be attributed to the action of a steady gravitational ¯eld upon
variations of density in a °uid. Following this concept, almost all the kinetic energy
of the earth's atmosphere results from convection. In atmospheric science, however, one
generallyusesamorerestrictedde¯nitionofconvectionwhichencompassesonlyaclassof
relativelysmall-scale,thermallydirectcirculationswhichresultfromtheactionofgravity
upon an unstable vertical distribution of mass.
11
Themostcommonlyseencloudsarecumulusclouds,whicharethetypeofcloudsthat
bear most characteristics of convection. For the other types of clouds, we can also view
them as results of the convective process with some other factors involved and dominant.
In the following, the cloud formation process will be described in detail. Before that,
a commonly used concept, ideal gas, in atmospheric science is ¯rst introduced.
2.1.2.1 Ideal Gases
Air is a mixture of several gases, including nitrogen, oxygen, and several other gases.
Their mixture is relatively uniform in the atmosphere. However, it is not the case with
some gases, such as water vapor and carbon dioxide, whose density can signi¯cantly vary
over the atmosphere. Even though the amount is not too much, due to the e®ects on
radiative transfer and the fact that it is the key factor of cloud formation, water vapor
is essential in atmospheric thermodynamics [1]. Gases in the atmosphere are ideal gases
and thus follow the ideal gas law that
p=½RT (2.1)
where ½ is the gas density and R is a gas constant. In atmospheric dynamic modeling,
air is treated as a mixture of \dry air" and water vapor, known as moist air. For moist
ideal gas, Eq. (2.1) changes into
p=½R
d
T
v
; (2.2)
12
where R
d
is a dry air constant and T
v
is the virtual temperature given approximately by
T
v
=T(1+0:6q
v
); (2.3)
and where q
v
is the water vapor density and virtual temperature is de¯ned as the tem-
perature that dry air would have if its pressure and density were equal to those of a given
moist air.
2.1.2.2 Parcels
In atmospheric science, the air parcel is a conceptual tool used to study the dynamics
of atmosphere. It is a small mass of air that can be thought of as \traceable" relative
to its surroundings. Parcel approximation is very useful in developing both the math
and understanding of cloud dynamics. The lifetime of a parcel can help in a better
understanding of the cloud formation process.
As an air parcel gets heated from the ground, it begins to rise due to the density
di®erence which accounts into buoyancy, and moves through the cooler surrounding air.
When a parcel is warmed at constant pressure, it expands and its density decreases. As
long as its density is smaller than that of the surrounding air, an air parcel will continue
tomoveup. Afteritreachestheleveloffreeconvection(LFC),theliftedairwillcontinue
to lift and vertical motion dominant in this movement. It continues to move up until a
level of equilibrium level (EL) where both its temperature and its density are close as
those of the surrounding air. Then, it stops to rise further and oscillates around the EL
level. Asanairparceloscillates,itdampenswithtimeandmoisturegetaccumulatedand
13
condensed; meanwhile, it moves laterally, forming the cloud anvil. When an air parcel
changes altitude without a change in heat (not to be confused with temperature), it is
said to move adiabatically. Because air pressure (and therefore temperature) varies with
altitude, the parcel's pressure and temperature change. The idea gas law and the laws of
thermodynamics can be used to derive Poisson's equation for adiabatic processes, which
relates to the temperature and pressure of a gas under adiabatic changes [45]:
µ
T
T
0
¶
=
µ
p
p
0
¶
·
; (2.4)
whereT
0
andp
0
areinitialvaluesoftemperatureandpressure,andT andparethevalues
after an adiabatic change in altitude and
·=
R
d
c
p
=
c
p
¡c
v
c
p
¼0:286;
and where constants c
p
and c
v
are the speci¯c heat capacities of dry air at constant
pressure and volume, respectively.
In atmospheric science, there is a more convenient variable to account for adiabatic
changes of temperature and pressure, i.e., potential temperature. The potential tem-
perature µ of a parcel can be de¯ned as the ¯nal temperature that the air parcel would
have if it were moved adiabatically from pressure p and temperature T to pressure p
0
.
Mathematically, we have
µ´T
µ
p
0
p
¶
R=cp
: (2.5)
14
Thepotentialtemperatureisaconvenientconceptinatmosphericmodelingsinceitsvalue
is constant under adiabatic change of altitude.
2.1.2.3 Lapse Rate
The Earth's atmosphere is in static equilibrium. The so-called hydrostatic balance of the
opposing forces of gravity and air pressure results in an exponential decrease of pressure
with altitude [45] of the following form
p(z)=p
0
µ
1¡
z¡
T
0
¶
g=¡R
d
; (2.6)
where z is altitude in meters, g is gravitational acceleration, 9:81ms
¡2
and p
0
and T
0
are
the pressure and temperature at the base altitude. Typically, p
0
¼ 100kPa and T
0
is in
the range of 280K¡310K. The lapse rate, ¡, is the rate of decrease of temperature with
altitude. In the Earth's atmosphere, temperature decreases approximately linearly with
heightinthetroposphere, whichextendsfromthesealeveltoabout15kmup(calledthe
tropopause). Therefore, ¡ is assumed to be a constant. A typical value for ¡ is around
10Kkm
¡1
. Eqs. (2.5) and (2.6) are used to compute the environmental temperature and
pressure of the atmosphere in the absence of disturbances, which, as described in the
following subsection, are used to compare with the local temperature and pressure for
the computation of air's saturation point.
15
2.1.2.4 Phase Transition
Water has three phases: water vapor, water droplet and ice, these three phases are able
to transit from one to another. Water droplet and ice are visible. Here, we only consider
the transition between the water droplet and the water vapor phases. The process is
that the water vapor particle oscillates around EL and get accumulated. When the water
vapor density exceeds a certain amount, it gets condensed. The maximum water vapor
density can be computed via
w
max
(T)=
217:8
T
exp
µ
19:482¡
4303:4
T ¡29:5
¶
: (2.7)
As the density of water vapor particles is higher than the maximum amount given by Eq.
(2.7), water vapor transit into water droplet.
2.1.2.5 Latent Heat
The energy process is involved in each phase transition of water. For instance, water
molecules are very energetic in the vapor phase, but less so in the liquid phase. Thus,
in the process of condensation, when water transits from the vapor phase to the liquid
phase, it transits from a higher-energy phase to a lower-energy phase, part of the energy
has to be released, which is called as latent heat. The release of latent heat causes a new
heat point source, forming a small-scale convection. In our modeling, latent heat is the
source of substructure of clouds.
16
After the introduction of necessary atmospheric concepts, the basic mathematics un-
derlying the modeling of the simple convective process will be presented, which is dis-
cussed in detail in Sec. 2.2.
2.1.3 Cloud Illumination
Besides cloud dynamics, cloud illumination is another challenge in cloud rendering since
cloud's interaction with light is complicated. We call the light interaction with cloud as
radiometry. Agoodunderstandingofcloudradiometryisnecessaryingeneratingrealistic
cloud images. The visible part of the cloud is the condensed cloud droplet. The most
common light interaction that a®ects the appearance of clouds is light scattering. The
physics of cloud illumination will be reviewed in this subsection.
2.1.3.1 Albedo of Cloud
The albedo is the ratio of scattered to incident radiation power most commonly for light.
It is a unitless measure of a surface's re°ectivity. Actually, for an incident light, the
cloud has two types of interaction formats, scattering and absorption. Scattering can be
thought of as a process that an elastic collision between an object and incident photon
whichmayalsochangethedirectionofphotonwhileabsorptionreferstotheprocessthat
the photon is absorbed and the energy of the light is converted to another type through
the medium.
In general, the albedo depends on the direction and directional distribution of incom-
ing lights. No matter a photon is scattered or absorbed, it is all \gone away" in some
17
sense. Thus, in general, we call it extinction. The albedo of clouds is high so that scat-
tering is the main interaction format with light and, actually, multiple scattering cannot
be ignored. In Fig. 2.1, a real cloud scene is shown, from which we can see the visual
e®ects of two interaction types with light.
Figure 2.1: A real cloud scene to illustrate two types of light interaction.
2.1.3.2 Single Scattering and Multiple Scattering
Scattering of light by a single particle is called single scattering. Multiple scattering
is the scattering of light from multiple particles in succession along the light path. To
understandsingleandmultiplescatteringphenomenabetter,weneedtointroduceoptical
conceptsaswellasphysicalproperties,includingtheparticledensityinsideaunitvolume
d
p
, particle material properties of absorption and scattering across sections of particles,
denoted by ¿
a
and ¿
s
, respectively.
Parameter ¿
s
is the property of material measured in the unit of the inverse of length
(cm
¡1
or m
¡1
), and its reciprocal represents the distance or the free path for scattering.
If the physical extent of a medium is smaller than this distance, then on average, photon
18
passingthroughitisscatteredatmostonce,whichistypicalofsinglescattering. However,
whenthephysicalextentofamediumisbiggerthanthisdistance,multiplescatteringwill
dominate in this model. From a global view, the coe±cients of particle's absorption and
scattering are the combination of the particle density and the material property; namely,
K
a
=d
p
¿
a
;
is the absorption coe±cient while
K
s
=d
p
¿
s
is the scattering coe±cient.
2.1.3.3 Optical Distance and Phase Function
The optical distance is a dimensionless measure used to denote the opaqueness of the
mediumaslightpassingthrough,knownascoe±cient¹. Amoreintuitiveconceptsimilar
to this one is transparency or transmittancy, which can be expressed as
T =e
¡¹(s;s
0
)
; (2.8)
where ¹ is the optical depth from s to s
0
. When ¹ = 1, T =e
¡1
¼ 37%, it indicate that
there is approximate 37% light passing through the medium without being scattered or
absorbed.
Besides the optical property of the medium, there is another important optical prop-
erty that greatly e®ects the visual appearance of cloud, which is the phase function. The
19
Á Á Á
Figure 2.2: The angle between the incident direction and the exitance direction Á.
phase function is a function of direction that determines how much light from incident
direction
1
is scattered into the exitance direction
2
as shown in Fig. 2.2.
The phase function is dimensionless and normalized since it is a probability distribu-
tion that satis¯es the following relationship:
Z
4¼
P(
1
;
2
)d
2
=1
Thecharacteristicsofscatteringdependonthesizeofparticles. Scatteringoflightdue
to large particles in clouds (i.e., cloud droplets) obeys Mie scattering [33]. Simpli¯cation
has been made by Henyey and Greenstein [23], which is known as the Henyey-Greenstein
phase function in form of
P(Á)=
1
4¼
1¡g
2
(1¡2gcosÁ+g
2
)
3=2
; (2.9)
wheregistheeccentricityoftheellipsetocontroltheanisotropyofscattering. Itspositive
value indicates that most of the incident light is scattered in the forward direction while
its negative value indicates in the backward direction. When g =0, it means an isotropic
20
scattering process with no direction favored. The cloud droplet has very strong forward
scattering. The scattering of a smaller cloud particle follows the Rayleigh Phase function
of the following form:
P(Á)=
3
4
(1+cos
2
Á)
¸
4
: (2.10)
For the scattering of the cumulus cloud, since it has more humidity, the particle size is
relatively large and it is proper to adopt the forward Henyey-Greenstein phase function.
2.1.3.4 Light Intensity Calculation
Fig. 2.3alsoillustratethelighttransportinclouds. Asdescribedbefore, lightabsorption
and out-scattering are responsible for light energy loss. The next particle along the light
path receives the amount of light from in-scattering so that the total light intensity can
be written as
dL(x;
2
)
ds
=¡K(x)L(x;2)+K
s
Z
4¼
P(x;
2
;
1
)L(x;
1
)d
1
; (2.11)
where dL is the change of radiance and ds is the direction along
2
. Integration of the
equation along the ray, we get the following simpli¯ed equation [32]:
L(D;
2
)=L(0;
2
)T(0;D)+
Z
D
0
g(s)T(s;D)ds; (2.12)
where
g(s)=K
s
Z
4¼
P(x;
2
;
1
)L(x;
1
)d
1
;
21
T(s;s
0
) is the transmittance in Eq. (2.8), D is the distance along the ray passing through
thecloudfromtheedgeoftheincidentlighttotheedgeoftheexitancelight,andL(0;
2
)
is the incident radiance.
In-scattering
Out-scattering
Transmission
In-scattering
Out-scattering
Transmission
Figure 2.3: The light transport phenomenon in clouds.
2.2 Similarity Approach
Atmospheric science is very complicated. However, accurate physics is not our goal to
pursue here. Instead, we want to obtain the visually convincing result with the help of
physics. Thus,simpli¯edphysicsthatisabletocatchthebulkdynamicspropertyismore
useful to researchers in the ¯eld of computer graphics.
2.2.1 Objective and De¯nition
In the study of dynamics of turbulent °ow, discrete convective entities are de¯ned to be
compact regions of °uids that penetrate homogeneously or are stably strati¯ed. Com-
monly studied convective entities include plumes, thermals, vortices from momentary
22
sources. A large class of °uids can be built with these convective entities. In partic-
ular, thermals and plumes serve as good building blocks to provide dynamics of moist
convection and cumulus clouds.
Althoughmoderncloudresearchisdominatedby3-DPDEanalysisof°uiddynamics,
the bulk dynamic property that result in convincing visualization is most relevant to our
currentapplications. Itismeaningfultostudyresultsforplumesandthermalsconsidered
in earlier cloud models. Actually, these models still have a great impact on atmospheric
scientists today.
The following two properties are fundamental to the behavior of convective entities.
² Similarity
There exists structural similarity in the mean °ows of convective entities as a func-
tion of the distance and/or time and, consequently, it is possible to derive their
bulk solution or use the reduced form of their governing equations. It should be
emphasized that it is the similarity, rather than the approximate equations, that
simpli¯es the computation.
² Entrainment
It refers to a process that convective entities attract the ambient °uid and take the
environmental mass to go along as they penetrate into the environment.
There are four di®erent similarity types [53] for consideration to facilitate cloud sim-
ulation as described below.
² Geometric similarity
It refers to similarity in shape regardless of the actual size. As a result, all lengths
23
scaleinproportionandcanberepresentedintermsofasinglecharacteristiclength,
which can be a function of time or displacement.
² Dynamical similarity
It refers to the similarity in force. The acceleration also varies proportionally.
² Self similarity
It refers to similarity between two successive stages of a single realization, where
the magnitude may change with time or distance but in proportion. For example,
for the turbulent plume, to be discussed in 2.2.2, the average straight-sided axis-
symmetric cones with a mean radius is proportional to the axial distance. These
cones are self similar.
² Similarity between realizations
Turbulent plumes have the same asymptotic spread angle under a wide range of
source conditions, which are known as similar °ows. Generally speaking, turbulent
plumes in a stable environment spread initially as approximately straight-sided
cones until their buoyancy reverses. Then, the plume °uid overshoots, falls back
and spreads sideways. Thus, the plume is not self similar for its whole life process.
Nevertheless, a similar single dimensionless solution can still be derived in Sec.
2.2.2.
2.2.2 Similarity Theory and Dimensionless Analysis
Similarity is important because it results in simpli¯cation of °ow simulation. It occurs
wherethe°owislargelydeterminedbythegrosspropertiesofthesource,e.g.,momentum
24
andbuoyancy°uxes, withonlythetransientin°uenceofthemass°ux, andtheturbulent
intensity level in the emerging °uid. In terms of di®erential equations, similar °ows are
more likely in a marching problem where all boundary conditions are given at the source
level.
Forsomesimple°uid°ow,agreatdealofinformationaboutitscharacteristicsmaybe
inferred from various parameters associated with the °ow, without solving the governing
equations explicitly. In particular, if the °ow is steady and various dependent variables
vary in a simple way with independent variables and boundary conditions, dependent
variables may often be related to independent variables and boundary conditions with
the aid of the governing equations simply by requiring the dimension of each side of the
equalitytobealike. Theanalysisisparticularlysimpleifthedependenceofalldependent
variablesononeormoreindependentvariablesissimilar. Ifthisistrue,weareableto¯nd
a solution that is consistent with the governing equations at a much lower computational
cost.
Before ¯nding the similarity solution, it is necessary to ¯nd the important dimension-
less parameters that characterizes the underlying °ow. The foundation for dimensional
analysis is the Buckingham pi theorem as described below.
Theorem 2.2.1 If equation '(q
1
;q
2
;q
3
;:::;q
n
) = 0 gives the only relationship among
(n,q)'s and if it holds for an arbitrary choice of units in which q
1
;q
2
;q
3
;:::;q
n
are mea-
sured, then the relation '(¼
1
;¼
2
;¼
3
;:::;¼
m
) = 0 is satis¯ed, where ¼
1
;¼
2
;:::;¼
m
are in-
dependent dimensionless products of q's. Furthermore, if k is the minimum number of
primary quantities necessary to express the dimensions of q's, then m=n¡k.
25
If governing equations are known, q's can be determined. Even if the governing
equations are not known, it is still possible to guess q's sometimes. Similarity theory
and dimensional analysis were ¯rst successfully applied in the investigation of simple
convective °ows by Schmidt [51] and Batchelor [3]. Thereafter, researchers in Britain
did a succession of studies, among them the most famous were Morton [34], Taylor and
Turner [61], whose experiments support the results of analysis.
In my thesis, the concern of investigation was the convection of \plumes" and \ther-
mals" released from a point of buoyancy in an ambient °uid with simple strati¯cation
under an assumption that the ambient °uid is not a®ected by the convection. Here, a
plume is de¯ned as a buoyant jet in which the buoyancy is supplied steadily from a point
source with a continuous buoyant region while a thermal refers to a discrete buoyant
element in which the buoyancy is con¯ned to a limited volume of °uid which is similar
to the concept of parcel introduced previously.
As already described, after an air parcel reaches EL, it starts to oscillate around this
level,andtheoscillationfrequencyisrelatedtothepotentialtemperatureortheairparcel
density. Accordingly, the buoyancy frequency has the form of
N
2
=¡
g
µ
0
dµ
dz
(2.13)
where N is called the Brunt-VÄ aisÄ alÄ a or buoyancy frequency in the unit of (time)
¡1
.
In a stably strati¯ed °uid, N is the frequency at which an in¯nitesimal sample of °uid
oscillates if displaced vertically. There are two possible cases.
26
² N
2
<0.
Environmental temperature falls more rapidly with height than parcels. Thus, the
parcel always ¯nds itself `lighter' than surrounding and hence continues to rise. At
this case, N is imaginary, meaning that the displaced parcel continues to move at
an increasing speed, and this corresponds to an unstable region.
² N
2
>0.
This happens when the °uid is stably strati¯ed. One might expect that the vertical
velocity will vanish at some height above the source as the buoyancy will become
negative at some lower height due to the decrease in density of the ambient °uid
with height.
Thus, quantity N
2
is a useful measure of atmospheric strati¯cation. We might want
to ¯nd a similarity solution for turbulent convection in the stably strati¯ed °uid. when
theambient°uidisstrati¯ed,anadditionalparameterdescribingthestrati¯cationisnec-
essary to fully describe the system. It is convenient to make explicit use of the governing
equations for Momentum, Heat and continuity as given below.
Momentum:
dV
dt
=¡
1
½
0
rp+B
^
k+ºr
2
V (2.14)
Heat:
dB
dt
=·r
2
B (2.15)
Continuity:
r¢V=0 (2.16)
27
The buoyancy B in Eqs. (2.14) and (2.15) is de¯ned as
B=g
Ã
µ¡
¹
µ
¹
µ
!
(2.17)
where µ is the potential temperature of the plume and
¹
µ is the potential temperature of
the environment.
By following Morton [34], Taylor and Turner [61], we assume a particular radial
dependence of velocity and buoyancy, and integrate the governing Boussinesq equations
over a horizontal plane. The particular form of radial dependence chosen will only a®ect
the numerical value of coe±cients in the resulting relations for vertical velocity w and
buoyancy B but not their dependence on z or the boundary °ux of buoyancy. There are
two types of radial pro¯les: the \top-hat" and the \Gaussian" pro¯les as illustrated in
Fig.2.4.
Figure 2.4: Illustration of the "top-hat" pro¯le (left) and the Gaussian pro¯le (right).
The following primary assumptions are made to solve the governing equations.
1. The °ow is steady.
28
2. The radial pro¯les of the mean vertical velocity and the mean buoyancy are similar
at all heights.
3. The mean turbulent in°ow velocity is proportional to the vertical velocity.
4. The °ow is Boussinesq.
For the third assumption, we use a simple relationship u = ¡®w, where u is the
mean turbulent radial velocity and ® is a constant that is proportional to the fractional
entrainment of mass. This assumption is useful for the following integration as well as
the determination of the particle density increase rate in the proposed cloud modeling
algorithm to be presented in Chapter 3. Experimental results have demonstrated that
the water vapor content distribution inside a cloud follows a \top-hat" pro¯le [17]. Thus,
we can integrate the Boussinesq mass continuity equation in the radial coordinate over
the horizontal area of the °ume and ¯nd
Z
2¼
0
Z
R
0
1
r
@
@r
(ru)rdrdµ+
@
@z
Z
2¼
0
Z
R
0
wrdrdµ =0
Using the entrainment relation, we get
2¼®Rw =
@
@z
(¼R
2
w) (2.18)
Furthermore, we need to integrate the momentum equation. In the horizontal direc-
tion,sincethepressuretermisrelativelysmall,weneglecttheperturbationfrompressure
gradient so that the vertical momentum equation becomes
29
dw
dt
=r¢Vw =B
Thereafter,theaboveequationisintegratedovertheincrementalvolumeasillustrated
in Fig. 2.5, which can be written as
Z
2¼
0
Z
R
0
Z
z+¢z
z
r¢Vwd¿ =
Z
2¼
0
Z
R
0
Z
z+¢z
z
Bd¿
With the divergence theorem, the term on the left may be expressed as a surface
integral as
Z Z
wV¢^ ndS =
Z
2¼
0
Z
R
0
Z
z+¢z
z
Bd¿
where ^ n is the unit normal vector to the surface, S, of the volume of integration. As
w vanishes at the lateral boundaries of the plume, we get the integrated equation as
d
dz
(¼R
2
w
2
)=¼R
2
B (2.19)
For the last step, we deal with the buoyancy equation over the incremental volume as
dB
dt
=0=r¢VB
and
30
R
ΔZ
Figure 2.5: The incremental volume over which the vertical momentum equation is inte-
grated.
Z Z Z
r¢VBd¿ =
Z Z
BV¢^ ndS =0
By evaluating the surface integral over the top, bottom, and lateral surfaces in Fig.
?? and using Eq. (2.15) for the representation of buoyancy B, we have
d
dz
h
¼R
2
w(µ¡µ
0
)
i
=2¼R®w(
¹
µ¡µ
0
)
After letting N
2
´
g
µ
0
d
¹
µ
dz
, we ¯nally have the heat equation as:
d
dz
(¼R
2
wB)=¡¼R
2
wN
2
(2.20)
Theprocessofsolvingtheturbulentconvectioninastablystrati¯ed°uidasdescribed
above demonstrates a way to obtain the similarity solution. Besides, it can be used
to determine the vertical water vapor distribution in the proposed voxel grid modeling
algorithm to be discussed in Chapter 4. It is also worthwhile to point out that the
31
similarity solution obtained above matches experimental results well. Consequently, the
similarity solution can be use to describe the characteristics of a simple °uid e®ectively.
2.3 Previous Work on Cloud Simulation and Rendering
Cloud simulation and rendering o®er one of the main challenging tasks in the modeling
and rendering of natural phenomena. Many e®orts have been made in cloud modeling
andrenderinginlasttwentyyears. Twoimportantaspectsincloudmodelingandrender-
ing, i.e., dynamics and radiometry were discussed in the last section. Related previous
research on this subject will be reviewed in this section. A complete cloud simulation
system involves three main ingredients, namely, dynamic physical simulation, modeling
and rendering. They will be detailed in the following three subsections.
2.3.1 Physical-based Dynamic Simulation of Clouds
Some researchers have attempted to perform vivid animation for clouds based on proce-
dural modeling. For instance, Perlin [40] used the turbulence texture to phase shift the
domainofthe°owinthey(vertical)directiontogetcloudanimation. Joshua[24]created
clouds with an implicit function and various octave noise by translating and rotating the
texture coordinates over time so as to reach rolling and turbulent e®ects. More recently,
Wang [67] used a bunch of textured sprites to model clouds and animated them by ro-
tating or translating the representative sprites, giving the illusion that wisps are moving
and tumbling with the wind. However, the procedure approach is rather limited. For
dynamic cloud simulation, we are more interested in the physical-based cloud simulation
approach. Related work along this direction will be our focus below.
32
Physical-based cloud simulation has been adopted by meteorologists and atmospheric
scientists due to its accuracy. Generally speaking, natural phenomenon simulation and
weather forecasting are a very complex task, which involves the solution of a nonlinear
system of partial di®erential equations (PDEs). Consequently, it requires long computa-
tional time to simulate cloud development for a relative short period of time. As high
performancecomputingfacilitiesbecomemoreaccessible,physical-basedcloudsimulation
has also attracted the attention of researchers in the ¯eld of computer graphics recently.
Simple one-dimensional (1D) spatial models were used in early computer simulation
ofatmosphericscience; forinstance, theworkofSrivastava[54]. Onlytheverticalmotion
¯eld and its changes under the in°uences of condensation and precipitation were consid-
ered in Srivastava's model. Later, this model was extended to the two-dimensional (2D)
case. Since the true three dimensional spatial model was too costly at that time, the slab
symmetry or the axial symmetry [45] was used by researchers. Although the methodol-
ogy is constrained by these symmetric properties, the resultant approach can simulate
thehorizontalwind shear, which is important in cloud dynamics. Such simulation results
were reported in [59]. However, to characterize the rotational °ow such as vortices along
horizontal and vertical axes of rotation, which is common in real clouds, simulation with
a three-dimension (3D) spatial model is needed. Steiner [58] presented a full 3D model,
and 3D cloud simulation has progressed since then. For a more detailed cloud simulation
in atmospheric physics, please refer to Rogers and Yau [45] and Cotton and Anthes [69].
Thefullsimulationperformedforatmosphericphysicsisoftentoomuchforcomputer
graphics applications. For the former, computer simulation is used to understand the
underlying physical process, which includes some details that may not be visible, e.g.,
33
speci¯c tracking of water phase, water droplet size distributions, complex microphysics
and °uid dynamics at a variety of scales. For the latter, the goal is simply to create real-
istic images and animation of clouds. As a result, simpli¯ed simulation for visualization
can be used.
Visualization-driven cloud simulation was done by Kajiya and Von Herzen [25]. They
solved a simple set of PDEs to generate cloud data sets and applied their ray tracing
algorithm for rendering. The PDEs solved were the famous Navier-Stoke equations of
incompressible °uid °ow (to capture thermodynamics for temperature advection and
latent heat e®ects) and a simple water continuity equation. The simulation required
about 10 seconds per time step (which corresponds to one second of cloud evolution) to
update a 10£10£20 grid on a VAX 11/780. More recently, Overby et al. [38] presented
a similar but slightly more detailed physical model based on PDEs. Built upon the
stable °uid simulation algorithm [56], their work solved the Navier-Stoke equations of
incompressible °uid °ow. The stability of this method allows a much larger time step
so that they achieved a simulation rate of one iteration per second on 15£15£15 grid
using an 800MHz Pentium III. Harris [30] implemented faster and more realistic cloud
simulation using programmable graphics hardware with °oating points.
Miyazaki et al. [48] used a coupled map lattice (CML), an extension of the coupled
map lattice model from the physics literature [26], [71], for cloud simulation. CML is an
extension of cellular automaton (CA) with continuous state values (rather than discrete
values) at each cell. Harris [21] implemented the CML method on graphics hardware.
Rulesfromatmospheric°uiddynamics,suchasapproximateincompressibility,advection,
vapor and temperature di®usion, buoyancy and phase changes are incorporated in CML.
34
Itwasreportedin[48]thatthesimulationofatimestepof3-5secondsona256£256£40
lattices demands about 10 seconds on a 1 GHz Pentium III.
Besides mathematical equations, some researchers tried simpler but less realistic rule-
basedsimulationtechniques. Neyret[36]usedananimatedparticlesystemtomodelcloud
behaviorwithasetofheuristics toapproximatethe rollingbehaviorof convectiveclouds.
Dobashi et al. [72] adopted a simple CA model introduced by Nagel and Raschke [35] to
animateclouds. TheoriginalCAin[35]includedrulesforthespreadofhumiditybetween
neighboring cells for the formation of clouds in humid cells, but had no mechanism for
evaporation. Dobashi et al. [72] added a stochastic rule for evaporation so that clouds
can grow and dissipate. Their model achieved a simulation time of about 0.5 seconds on
a 256£256£20 volume using a system of dual 500MHz Pentium III.
As compared with previous work in cloud formation for computer graphics applica-
tions, our work to be presented in this proposal has several unique features. Research in
Chapter 3 is similar to Neyret's work [36] in terms of modeling tools. That is, both use
particles to describe the cloud formation process. However, the cloud simulation process
in Chapter 3 is not based on heuristic approximation, where the similarity solution is
adopted to model the convective process and thus the cloud behavior. The similarity
solution for cloud simulation as presented in Chapter 4 is similar to that of Harris [30],
Kajiya and Von Herzen [25] and Overby et al. [38]. However, our work focuses more on
simulation simpli¯cation with the similarity approach.
35
2.3.2 Graphic Tools for Cloud Modeling
There exist many graphic tools for cloud modeling. An explicit representation for each
visible cloud droplet requires too much computation and storage. Thus, coarser graphic
models are often proposed. In this subsection, six commonly used graphic tools for cloud
modeling are discussed. They are polygon, procedural noise, textures sprite, metaball,
particle system and voxel volumes. However, these models are not mutually exclusive
and a proper combination of these tools may generate better results.
2.3.2.1 Polygon
For cloud modeling in the context of TV weather forecasting rendering, geometry prim-
itive polygons is suitable as reported in [2]. The polygon is generated from the original
data grid with the marching cubes algorithm. According to the di®erent viewing angle
and distance, the technique of level of details (LOD) can be used to save computation.
Afterwards, the triangle geometry is then rendered with semi-transparent texturing. The
hardware was used in [2] to compute cloud opacity.
2.3.2.2 Procedural Noise
Procedural noise was used to model clouds at the very beginning of the cloud modeling
history, e.g. [29,39], due to its simplicity and capability for ¯ne detail representation. It
usesanoisefunctionasthebasisandperformsacombinationofmultiplenoisefunctionsto
generatea3Drandom¯eldwithacontinuousdensityto¯llthecloudvolume. Ebert[12{
14] made signi¯cant contributions to this topic by considering both o®-line computation
of realistic cloud images and real-time computation of clouds, smokes and steams with
36
graphichardware. Ebert[13]usedanunionofimplicitfunctionstocontrolcloud'sgeneral
shape, perturbed the solid space de¯ned by implicit functions with procedural noise to
createdetail,andthenusedascanlinerendererforrendering. Joshua et al.[24]extended
Ebert's work for real-time animation and rendering with the aid of GPU. Hugo [15]
extended the procedural noise method to render the cloud-related natural scenery.
2.3.2.3 Textures Sprite
Researchers used textures to represent the cloud surface (instead of its volume). Gardner
[20] applied fractal texturing to the surface of ellipsoids to create the cloud appearance.
With the combination of multiple textured and self-shaded ellipsoids, he synthesized
convincing cloudy scenes that is reminiscent of landscape paintings of John Constable.
Lewis [29] used ellipsoids to model clouds together with procedural noise to generate ¯ne
detail. Elinas and StÄ urzlinger [16] used a variation of Gardner's method to render clouds
composed of multiple ellipsoids on hardware interactively. Wang [67] modeled each cloud
pu® by a bunch of textured sprites, where the general cloud shape is manually designed
by an artist in 3DS MAX. This technique was used in Microsoft's Flight Simulator 2004.
Harris and Lastra [22] used the impostor technique to accelerate routines for real-time
rendering of distant clouds.
2.3.2.4 metaball
By the metaball technique [8], a volume is represented as a superposition of potential
¯elds from a set of sources, each of which is de¯ned by a center, radius and strength,
which is also called a 'blob'. Simple logical operations such as AND/OR can be de¯ned
37
to manipulate these metaballs. Due to its simple mathematical expression, the ¯eld can
be visualized by ray tracing relatively quickly. However, it is still too slow for interactive
display. Polygonized routines were developed for speed-up. That is, a complete metaball
area is divided into a 3D grid and a calculation is performed for each edge in the grid.
A vertex is created at the location where the formula has a turning point. However,
the isosurface extracted may not be suitable for cloud rendering. Nishita et al. [60] used
metaballs to model clouds by hand-placing several metaballs manually to create a rough
shape and then added detail with a fractal method to generate new metaballs on the
surface of existing ones. Dabashi et al. [72] used metballs to model clouds from satellite
images and other natural atmospheric scene. Furthermore, they simulated clouds on a
voxel grid with logic binary operations and then converted the grid representation into
metaballs for rendering with splatting.
2.3.2.5 Particle System
Another popular method for natural scene modeling is the use of particles [41], especially
for the modeling of objects/sceneries with an amorphous shape or ¯ne detail, which do
not have well-de¯ned boundaries. A particle is a simple computer graphics primitive
with a single 3D position and a small number of attributes such as the radius, color,
texture,lifetime,etc. Inaparticlesystem,objectsaremodeledbyacollectionofparticles.
Reeves [41] applied the particle system to the modeling of clouds and other \fuzzy"
objects/phenomena. Later,ReevesandBlau[42]proposedastructureparticlesystemand
developed a probabilistic algorithm for particle shading and rendering. Particles can be
created manually using some modeling tool, procedurally generated, or the combination
38
of these two. There are many ways for particle rendering. The static cloud was modeled
by particles and rendered as a small, textured sprites (or \splat" [68]) by Harris [22].
This will be used to model and render clouds in Chapter 3.
The particle system is attractive since its codes are simple to maintain and it can
render amorphous objects e®ectively. Since a particle represents a spherical volume im-
plicitly, a coarse cloud modeled by particles requires much less storage than that built
with other methods. However, this advantage may diminish when the detail is increased,
since many tiny particles are needed to achieve the desired quality. Under this case,
some alternative techniques might be used or combined. Another modeling tool, which
is similar to the particle system, was proposed by Stamminger [57], where point-based
modeling and rendering were considered for trees and clouds. The point sample in his
model has similar attributes as the particle in a particle system.
2.3.2.6 Voxel Grid
Avoxelisthe3Danalogofapixel. Itisasingleunitofaregulargridsubdivisionofarect-
angular volume. A voxel can be used to model and render volumetric objects, especially
those commonly used in physical-based modeling that involves the solution of PDEs.
Examples of voxel-based cloud modeling and rendering can be found in [9,27,28,68,70].
Kajiva and Von Herzen [25] performed a simple physical cloud simulation, recorded the
result in the voxel volume, and then applied ray-tracing for cloud illumination. Dobashi
et al.[72]simulatedcloudsonvoxelgridusingCAwithbinarylogicoperations,converted
the result into the metaball representation and rendered it using splatting. Miyazaki et
al. [48] used the extension of CA, called the coupled map lattice (CML), to simulate the
39
cloud and then rendered it in the same way as [72]. Overby et al. [38] solved PDEs to
generate various cloud types with the voxel grid and then rendered them using Harris's
SkyWorks rendering engine [22]. Later, Harris [30] solved PDEs with computer graphics
hardwaretoachievereal-timecloudsimulationonthevoxelgrid. Thevoxelgridapproach
plays a critical role in the PDE-based physical simulation system, which creates realistic
clouds. However, due to the high computational cost, the voxel grid cannot be too large
which may sacri¯ce the cloud detail.
2.3.3 Cloud Rendering
Cloud rendering plays a critical role in the visualization of a generated cloud model.
Early work on light scattering for the computer graphics application was conducted by
Blinns [7]. To render the rings of Saturn, Blinns [7] proposed an approximate method
to compute the appearance of cloudy or dusty surfaces via statistical simulation of light-
matter interaction under the assumption that the primary e®ect of light scattering is
resulted by the re°ection from a single particle in the medium while multiple re°ections
canbeignored. Thissinglescatteringassumptioniscommonlyusedincomputergraphics.
However, the single scattering assumption is only valid for media with particles of a low
single scattering albedo.
Accuratecomputationoflightscatteringinmediawithhighsinglescatteringalbedois
expensivesinceitrequirestheevaluationofadoubleintegral. Incomputergraphics,either
simplifying assumptions are adopted to reduce the complexity or o®-line computation is
performed. Work done on this topic can be classi¯ed into ¯ve categories: 1) discrete
40
ordinates, 2) spherical harmonics, 3) ¯nite element, 4) Monte Carlo integration and 5)
line integral methods. They are reviewed below.
2.3.3.1 Discrete Ordinates
Indiscreteordinates,acollectionofM discretedirectionsisassignedtostoretheradiosity
exiting each volume element, where the intensity along each direction is assumed to be
constant. This method can yield anisotropic scattering. If the interaction between a pair
of elements can be represented by only one direction, the number of non-zero elements
in a matrix of linear coe±cients is MN
2
, where N = n
3
is the number of elements in
the volume [32]. Chandresekhar [49] used discrete ordinates in the simulation of radia-
tive transfer. Later, Max [32] pointed out sampling artifacts in Chandresekhar's method
due to a ¯nite partition of the spatial angle into discrete directions. To improve the
basic method of discrete ordinates, Max [32] spread the shot radiosity over the entire
direction e±ciently, rather than along discrete directions. The method achieves a sig-
ni¯cant speedup by dealing a whole plane of source elements simultaneously, where the
computation time is reduced to O(MNlogN)+M
2
N.
2.3.3.2 Spherical Harmonics
ThesphericalharmonicsY
m
l
(µ;Á)aretheangularsolutionofLaplace'sequationinspheri-
calcoordinates. Thesphericalharmonicsformacompleteorthonormalbasis. Thismeans
41
that an arbitrary function f(µ;Á) can be represented by an in¯nite series expansion in
terms of spherical harmonics, i.e.
f(µ;Á)=
1
X
l=0
l
X
m=0
A
m
l
Y
m
l
(µ;Á):
Givenanumberofsamples(valuesoff),aseriesoflinearequationscanbeformulatedand
solved for coe±cients A
m
l
. The spherical harmonics method has been used by Bhate and
Tokuta [5], Kajiya and Von Herzon [25] and Stam [55] to implement multiple scattering.
Kajiya and Von Herzen [25] presented a ray tracing technique to render an arbitrary
volumeofscatteringmedia. Inadditiontoasimplesinglescatteringmodel,theydescribed
a method for multiple scattering that uses spherical harmonics. Their single scattering
simulation method stores the cloud density and illumination data on voxel grids. Their
algorithmisatwo-passprocess. Inthe¯rstpass,scatteringandabsorptionareintegrated
along paths from the light source through the cloud to each voxel where the resulting
intensitiesarestored. Inthesecondpass,raysaretracedthroughthevolumeofintensities
and light scattering to eyes was computed, resulting in a cloud image. For multiple
scattering, they derived a discrete spherical harmonics approximation to the multiple
scattering equation, and solved the resulting matrix of PDEs using relaxation. This
matrix replaces the ¯rst integration pass of the single scattering algorithm, which is
known as the P
N
-method in the transport theory literature, where N is the degree of the
highest harmonic in the spherical harmonic expansion.
While Kajiya and Von Herzen [25] derived a general N-term expression for multi-
ple scattering using a spherical harmonics expansion, Stam [55] truncated the expansion
42
after the ¯rst term to produce results of a di®usion type for the scattered portion of
the illumination ¯eld. In \optically thick" media, where scattering events are frequent,
multiple scattering can be approximated as di®usion of the light energy. That is, at any
location in the medium, photons can be found traveling in arbitrary directions. The light
is said to be di®used. Stam [55] studied this di®usion approximation in more detail.
Like [25], Stam represented the scattering medium on a voxel grid and described two
ways to solve for the intensity. In the ¯rst method, he discretized the di®usion approx-
imation on the grid to formulate a system of linear equations, and then solved it using
the multigrid method. The second method is a ¯nite element method, in which he used a
series expansion of basis functions, speci¯cally the Gaussian kernel of the distance. This
expansionleadstoamatrixsystemwhichissolvedusingtheLUdecomposition. Inspired
by [25] and [55], two-pass algorithms for computing light scattering in volumetric media
are common nowadays.
2.3.3.3 Finite Element
The ¯nite element method provides another technique for the solution of integral equa-
tions, which has been used for light transport calculation. In the ¯nite element method,
anunknownfunctionisapproximatedbydividingthedomainofafunctionintoanumber
ofsmallpiecescalledelements,overwhichthefunctioncanbeapproximatedusingsimple
basis functions. Then, the unknown function can be represented with a ¯nite number of
unknowns and solved numerically.
One application of ¯nite elements in computer graphics is the radiosity method in
computingdi®usere°ectionamongsurfaces. Inthismethod, surfacesofascenerepresent
43
the domain of a radiosity function and an integral equation characterizes the intensity
of light re°ected from these surfaces. To solve the radiosity equation, surfaces are ¯rst
subdivided into a number of small elements and then radiosity is represented by a sum
of weighted basis functions (or elements). This formulation results in a system of linear
equationsthatcanbesolvedfortheseweights. Thecoe±cientsofthissystemareobtained
by integrals over surfaces. Intuitively, light incident on an arbitrary point in the scene
can be re°ected to any other point so that the coe±cients are integrals over the entire
scene. This is called the form factor by Cohen [11], which has to be evaluated for every
pair of elements in the scene.
There are a couple of methods originally developed for radiant heat transfer analysis.
Oneistheradiositymethodasdescribedaboveandtheotheristhezonalmethodproposed
by Rushmeier and Torrance in [47], which is an extension of the radiosity method by
including radiative transfer in volumes of participating media. The zonal method divides
the volume of a participating medium into ¯nite elements which are assumed to have
constant radiosity. Similar to the radiosity method, form factors are computed for every
pair combination of surface elements in the scene and every pair of volume elements
and all surface-volume pairs. Like the radiosity method, a system of simultaneous linear
radiosity equations is generated based on these form factors. The solution of this system
is the steady-state di®use radiosity at each element of the environment, including the
scattering and absorption e®ects of the participating medium. Rushmeier and Torrance's
zonal method is restricted to isotropic scattering media.
Nishita et al. [60] introduced approximations and a rendering technique for global
illumination of clouds, accounting for multiple anisotropic scattering and skylight. This
44
method can be viewed as a ¯nite element method since the volume is divided into voxels
and radiative transfer between voxels is computed. They made two simplifying obser-
vations to reduce the computational cost. The ¯rst one was that the phase function
of cloud water droplets is highly anisotropic, favoring forward direction scattering. Its
consequence is that not all directions contribute strongly to the illumination of a given
volume element. Therefore, one can compute a \reference pattern" of voxels that con-
tributes signi¯cantly to a given point. This pattern is constant at every position in the
volume since the sun is located very far away. Then, the same sampling pattern can be
used to update the illumination of each voxel. The second observation was that only
the ¯rst few orders of scattering contribute strongly to the illumination of a given voxel.
Therefore, only up to the third order of scattering was computed in [60]. Harris [22]
extended the result in [60] by exploiting the anisotropy of scattering for cloud droplets in
cloud illumination.
2.3.3.4 Monte Carlo Integration
Monte Carlo integration is a statistical method that uses sequences of random numbers
tosolveintegralequations. Incomplexproblemslikelighttransport,wherecomputingall
possible light-matter interactions would be impossible, Monte Carlo integration reduces
the complexity by randomly sampling the integration domain. With enough samples,
whicharechosenintelligentlybasedonimportance,anaccuratesolutioncanbefoundwith
much less computation than that would be required by a complete model. The technique
of intelligently choosing samples is called importance sampling, and the speci¯c method
depends on the problem being solved. A common application of Monte Carlo methods in
45
computergraphicsisMonteCarloraytracing[63]. Wheneveralightraytraversesascene
interactswithmatter(asolidsurfaceoraparticipatingmedium), statisticalmethodscan
be used to determine whether the light is absorbed or scattered. If the light is scattered,
the scattered ray direction is also chosen statistically. Importance sampling is used to
determine the direction by evaluating a probability function.
A technique to render an arbitrary volume of participating media using Monte Carlo
ray tracing was proposed by Blasi et al. [6]. They placed no restriction on the medium,
allowingarbitrarydistributionsofthedensityandthephasefunctionsandaccountingfor
multiple scattering. They demonstrated an importance sampling technique that uses the
phase function as a probability function to determine the outgoing direction of scattered
rays. In this manner, the in-scattering integral does not have to be evaluated over the
entire sphere of incoming directions and a large amount of computation can be saved.
Using the phase function for importance sampling ensures that the most signi¯cant con-
tributions to scattering are used to determine the intensity.
Jensen [64] proposed a photon mapping technique, which is a variation of pure Monte
Carloraytracing, wherephotons(particlesofradiantenergy)aretracedthroughascene.
Starting at the light sources, many photons are traced through the scene. Whenever a
photon lands on a non-specular surface, it is stored in a photon map, which is a data
structure to store the position, the incoming direction, and radiance of each photon hit.
The radiance on a surface can be estimated at any point from photons closest to that
point. Photon mapping demands two passes: to build the photon map in the ¯rst pass
and to generate an image from the photon map in the second pass. Image generation
is typically performed using ray tracing from eyes. The photon map has the °exibility
46
of Monte Carlo ray tracing, but avoids the grainy noise. Jensn and Christensen [65]
extended the basic photon map to incorporate the e®ects of participating media. To do
so, they introduced a volume photon map to store photons within participating media
and derived a formula to estimate the radiance in the media using this map. Their
techniquesenable thesimulationof multiplescattering, volumecaustics(focusing of light
onto participating media caused by specular re°ection or refraction), and color transfer
between surfaces and volumes of participating media.
2.3.3.5 Line Integral
Ascomputergamesbecomepopular,interestsinsimulatinglightscatteringforinteractive
applicationshavegrown. Forview-dependente®ectsanddynamicphenomena,techniques
described in previous sections are not practical. While those techniques accurately por-
tray the e®ects of multiple scattering, they require a large amount of computation. For
interactive applications, simpli¯cation should be made.
First, we may ignore volumetric scattering altogether to simplify the computation.
With or without scattering, visualization of shadowing e®ects due to absorption by the
medium is desirable. This requires at least one pass through the volume (along the
directions of light propagation) to integrate the intensity of transmitted light. Since the
simpli¯ed method performs the intensity integration along lines from the light source
through the volume, Harris [22] called it the line integral method. Kajiya and Von
Herzon's original single scattering algorithm is a line integral method. Intuitively, line
integral method is limited to single scattering since they cannot propagate light back to
points already traversed. In real time cloud rendering, Harris [22] demonstrated that the
47
line integral method can be used to compute multiple scattering in the forward direction
so that it can be used for interactive cloud rendering.
Dobashi et al. [72] also described a simple line integral technique to compute the
cloudilluminationusingthestandardblendingoperationsprovidedbycomputergraphics
APIs such as OpenGL . They represented clouds as collections of large \particles" in
form of textured billboards. To compute illumination, they rendered particles in an
increasing distance from an initially white frame bu®er to the sun. They con¯gured
OpenGL blending operations so that each pixel covered by a particle was darkened by
an amount proportional to the attenuation of the particle. After rendering a particle,
they read the color of the pixel at the center of projection of the particle from the frame
bu®er and stored this value as the intensity of incident light that reached the particle
through the cloud. Traversal of the particles in the order of increasing distance from the
light source evaluates the line integral of extinction through each pixel. Because pixels
are darkened by every particle that overlaps them, this method computes accurate self-
shadowing of the cloud. After this ¯rst pass, they rendered particles from back to front
with respect to the viewing point, using the intensities computed in the ¯rst pass. They
con¯guredblendingtointegrateabsorptionandsinglescatteringalonglinesthrougheach
pixel of the image, resulting in a realistic image of the cloud. They further enhanced this
realism by computing the shadowing of the terrain by clouds and shafts of light between
clouds and the ground.
The volume must be traversed at least once to integrate attenuation due to absorp-
tion. During this traversal, it makes sense to integrate scattering along the direction of
traversal. Thus, Harris extended Dobashi's work to compute multiple forward scattering.
48
In the cloud rendering system given in Chapter 4, we follow Harris' cloud illumination
model with some modi¯cation.
A similar line integral approach was proposed by Kniss, et al. [27] for absorption and
multiple forward scattering in the context of direct volume rendering. They rendered
volumes of translucent media from 3D textures by rendering slices of the volume oriented
to the face along the halfway vector between the light and view directions. This \half
angle slice" technique allows to interleave light transport integration with the volume
display. The method traverses the volume slices in the order of increasing distance from
the light source, performing alternate display and illumination passes. Three bu®ers are
maintained: twoforthecomputationoftheilluminationofthevolume(currentandnext)
and one (typically the frame bu®er) for display of the volume. In the display pass, the
current slice is rendered from observer's point of view. The slice is textured with the 3D
volumetextureblendedwiththecurrentilluminationbu®er. Thisresultinself-shadowing
of the volume as well as incorporation of scattering during the illumination pass as done
in [22]. During the illumination pass, the slice is rendered into the next illumination
bu®er from light's point of view, and blended with the current illumination bu®er to
computethenextstepinthelineintegralofextinctionandforwardin-scattering. During
thisblending, thecurrentbu®erissampledmultipletimesatjitteredlocations, andthese
samples are averaged. This accomplishes a blurring of forward-scattered light and ad hoc
approximationofmultiplescatteringoverasmallslideangelaroundtheforwarddirection.
Even though this method is ad hoc, it is physically meaningful since multiple scattering
in media with a high single scattering albedo result in \blurring" of light intensity.
49
Chapter 3
A Similarity Approach to Particle-based Cloud Simulation
Actually, in a turbulent °ow, the path of individual particles is
highly irregular and only a statistical average over time or space
reveals a systematic velocity pro¯le.
|Atmosphere convection
3.1 Introduction
A new methodology for cloud simulationand modeling is proposed in this chapter, where
a similarity approach is adopted for cloud simulation and a particle system is used for
cloud modeling. Our methodology is motivated by the following observations. Cloud
simulation usually involves the solution of partial di®erential equations (PDEs) and the
complexity is very high. For a simple °uid °ow, it is possible to use a similarity approach
to simplify the computational cost. We show that this is a feasible solution for cloud
simulation. Furthermore, we know that the air parcel's life is from water vapor to water
50
droplet in cloud dynamics. Since a commonly used primitive in computer graphics, i.e.
particles, bear similarity with parcel, we choose particles for cloud modeling.
Wecandivide cloud dynamicsintotwostages. The ¯rst stage involveswater vapor's
movement and the wind e®ect on the vapor. Before water vapor changes its phase, its
velocity follows the pro¯le obtained by the similarity approach. At the second stage, the
water vapor changes its phase and substructure due to latent heat release. A general
framework of cloud modeling is shown in Fig. 3.1 to capture the water vapor behavior in
cloud formation at these two stages.
Based on the physics obtained by a similarity solution of turbulent thermal in the
stable strati¯ed °uid presented in Section 3.2, a particle system is used to model the
cloud formation in Section 3.3.
Two-stage Simulation
Particle-based
simulation stage
Water vapor movement
Wind effect
Voxel stage
Phase transition
Substructure
Cloud Rendering
Two pass rendering
(Multiple scattering)
Figure3.1: Thecloudmodelingframeworkwhichconsistsofthetwo-stagecloudformation
process and the cloud rendering process.
Cloudsynthesisischallengingnotonlyonitsmodelingandsimulationwhichgenerate
amorphous shape, but also on illumination. Thus, after cloud formation, we consider
cloud illumination and rendering in Section 3.4. Cloud is a high albedo scene. Multiple
scattering should be considered in the lighting model if we want to get a realistic cloud
image. Multiple scattering can be achieved by a two-pass rendering process. The light
51
intensity incident to each particle is calculated in the ¯rst pass while the light incident
to viewer's eye is determined in the second pass.
Synthesized clouds are shown in Section 3.5 to demonstrate the performance of our
proposed cloud simulation, modeling and rendering techniques.
3.2 Thermal Simulation by Similarity Approach
Thermalisusedtorepresenttheairparcelinthischapter. Solvingthesimilaritysolution
for turbulent thermal is similar to the solution for convective plumes as described in the
lastchapter. Ifwetreattime(ratherthanheight)astheimportantindependentvariable,
many assumptions made concerning the behavior of plumes may also apply to thermal.
They are stated below.
1. The radial pro¯les of the mean vertical velocity and the mean buoyancy are similar
at all time.
2. The mean entrainment velocity is proportional to the mean vertical velocity.
3. Thedensityperturbationinthethermalissmallascomparedwiththemeandensity
(Boussinesq approximation).
Byemployingproceduressimilartothoseusedtodrivetheradiallyintegratedplume
equations in Chapter 2, we may derive the following conservation equations integrated
over the volume of the thermal.
Mass:
d
dt
(
4
3
¼R
3
)=4¼R
2
®w: (3.1)
52
Momentum:
d
dt
(
4
3
¼R
3
w)=
4
3
¼R
3
®B: (3.2)
Heat:
d
dt
(
4
3
¼R
3
B)=¡
4
3
¼R
3
®wN
2
: (3.3)
The reason to use the similarity approach is that we want to use a limited number
of parameters to characterize a simple °uid. The turbulent thermal is a good example
for dimensional analysis from a maintained source of buoyancy in a semi-in¯nite, homo-
geneous °uid (above the ground). We assume the °uid is fully turbulent which is true in
atmosphere. Then, it should be independent of the magnitude of the molecular di®usive-
ness. If the Boussinesq approximation is applicable, then the only relevant dimensional
parameter in the problem is rate F at which buoyancy is supplied by the point source
and a parameter N
2
introduced from the governing equation.
In the case of constant stable strati¯cation (where N
2
is a positive constant), we are
abletosolvealloftheaboveequationsanalytically. Byrephrasingtheequationsinterms
of new dependent variables; namely,
M ´R
3
w; V ´R
3
; F ´R
3
B;
it is easier to solve them analytically.
Ifthecloudhaszeroradiusandnomomentumatthetimeofitsrelease,theboundary
conditions at t=0 will be
M =V =0; and F =F
0
:
53
The resulting equations are conveniently rendered dimensionless by the follow scaling:
F
¤
= F
0
f;
M
¤
=
3
4¼
F
0
N
¡1
m;
V
¤
=
µ
3
¼
¶
3=4
®
3=4
F
3=4
0
N
¡3=2
v;
t
¤
= N
¡1
t;
where asterisks denote the dimensional values. If we de¯ne the dimensionless forms of
radius, vertical velocity, and buoyancy as R, w and B, respectively, then the solution of
the ¯rst complete oscillation is
² 0·t·¼ :
R = (1¡cost)
1=4
;
w =
sint
(1¡cost)
3=4
;
B =
cost
(1¡cost)
3=4
;
z = 4(1¡cost)
1=4
:
² pi·t·2¼ :
R = (3+cost)
1=2
;
w =
sint
(3+cost)
3=4
;
54
B =
cost
(3+cost)
3=4
;
z = 2
13=4
¡4(3+cost)
1=4
::::
After getting dimensionless result, we can recover dimensional results easily. For
example, we can get height z
¤
from z via
z
¤
=
1
4
µ
3
¼
¶
1=4
®
1=4
F
1=4
0
N
¡1=2
z; (3.4)
R
¤
from R via
R
¤
=
Ã
µ
3
¼
¶
3=4
®
3=4
F
3=4
0
N
¡3=2
R
!
1=3
; (3.5)
where ®, N and F
0
are all constant parameters.
Ten oscillations were solved from Eqs. (3.1)-(3.3) for the turbulent thermal. We
show the damped oscillatory behavior of the thermal with 10 oscillations in Fig. 3.2.
We see from Fig. 3.2 the characteristics of thermal's movement, where the x-axis is
timeandthey-axisisthedimensionlessheight. Individualparticlesofthe°uidovershoot
thelevelatwhichtheirbuoyancy¯rstvanishesandthendeceleratetozerovelocity. Then,
a thermal oscillates around its equilibrium level and gets dampened. In the meanwhile,
it continues to expand and its radius increases.
55
5pi 10pi 15pi 20pi 25pi 30pi
-1
0
1
2
3
4
5
Time(t)
Dimensionless Height
R’-t W’-t B’-t Z’-t
Figure 3.2: The solution for radius (R), height (z), buoyancy (B) and vertical velocity
(w) of a thermal in a uniform stably strati¯ed °uid.
3.3 Thermal Movement Modeling by Particles
3.3.1 Parameters in Particle Systems
We use an individual particle to model thermal's lifetime. As it moves, it follows the mo-
tion described by the similarity solution. One new attribute, i.e. the density, is assigned
toaparticle, whichwillbeusedinthecloudformationprocess. Anothernewattributeof
aparticleisitscondensationtime. Duetorandomnessof theinitial momentum, di®erent
particles have a distinct ability to oscillate before it get condensed. Thus, we give each
of them a condensation time, which is actually the lifetime of the water vapor particle.
When the condensation time reaches, this particle is condensed at the voxel and it ¯nally
dies.
56
At each time step, particles are generated at the ground level and the water vapor
densityisupdatedinvoxels. Boththenumberofwatervaporparticlesandtheintervalto
release water vapor particles can be °exibly selected by users. Thus, our implementation
is not constrained by the volume resolution, which is the biggest limitation on the PDE-
based simulation method despite their high accuracy.
Although the physics process is not the greatest motivation while the visual quality
is in computer graphics, the more close to real physical process, the more °exibility
extent we can make in further work. For physics-based cloud formation, the Navier
Stoke equation is only the governing equation for density movement. However, the real
formation involves more. For example, cloud anvil caused by thermal's oscillation and
cloud accumulation cannot be well represented by the governing equation. Furthermore,
the cloud vertical extent is limited by Buoyancy frequency. The similarity solution helps
avoid an explicit solution of the governing PDE and deal with these facts that go beyond
the density-governing equation.
Even though the velocity similarity for water vapor particles is considered by our
model, a real cloud may have an amorphous shape, which implies that there are other
factors to result in these di®erent shapes. They are identi¯ed in following sections.
3.3.2 E®ect of Water Vapor Life Time
Duetotherandomnessoftheinitialmomentum, di®erentparticlesexhibitdi®erentoscil-
lating patterns. We assign each particle a condensation time, which is the lifetime of the
water vapor. When the condensation time reaches, the particle is dead and condensed
57
into a voxel. As discussed in Section 3.2, a thermal (or a water vapor particle) may oscil-
latearoundtheEL.Duringthisprocess,iteitheraccumulatesandcondensestothewater
droplet or continues to oscillate, and moves laterally causing the cloud anvil. Although
the amount of heat that each water vapor particle gets from the ground area is similar,
it still varies and this results in di®erent heights as given in Eq. (3.4).
As a consequence of the water vapor transition process, the higher the altitude
is, the less likely the water vapor will remain. Thus, a water vapor particle with a
higher momentum tends to reach its condensation level faster and its condensation time
is shorter. In contrast, for water vapor particles heated by mean F
0
at the ground, they
tend to oscillate more around the EL level. During the oscillation, it moves laterally to
form the cloud anvil. Since it has a random initial horizontal velocity, the result also
varies.
Since the condensation time is related to the heat source, we can use the heat source
to control the condensation time. Note that, typically, condensation does not happen
during the down-draft motion of the water vapor particle. When the water vapor moves
downward, it gets compressed and a higher pressure makes the phase change less likely.
To exploit this observation for model simpli¯cation, condensation only happens in the
up-draft motion in our model. Fig. ?? shows a statistical distribution of condensation
time in the cloud formation, and the corresponding cloud shape is given in Fig. 3.4(a).
3.3.3 Cloud Shape Control
A cloud shape is generally related to the water vapor source on the ground. Depending
on vegetation (e.g., barren/grass) and the soil type, a water vapor can be generated
58
2 4 6 8 10
0
200
400
600
800
Condensation time (oscillation)
Particle number
Figure 3.3: The histogram of condensation time for cloud formation.
di®erently. By the atmospheric measurement, the horizontal water content concentration
ofacumuluscloudwithoutthewinde®ectisGaussiandistributed[17]. Thus,aGaussian-
shaped vapor source is typically chosen on the ground. This model controls the water
vapor momentum as well as its emission frequency. The higher emission frequency a
water vapor source has, the more water vapor particles to be released to the space. We
canspecifythemeanvalueofthewatervaporsourceasafunctionoftimeandspace. The
actualwatersourcecanbeperturbedbyPerlinnoisetogeneratethee®ectofrandomness.
Fig. 3.4 shows a pair of cloud formation results with di®erent water vapor initial
distributions. Fig. 3.4(a) has quite a few particles with short condensation time so that
the cloud has an obvious top. On the other hand, if fewer particles have high heat °ux,
the resultant cloud does not have a obvious top as shown in Fig. 3.4(b).
59
(a) (b)
Figure 3.4: Comparison of cloud shapes controlled by di®erent initial water vapor distri-
butions.
3.3.4 Wind E®ect
When a wind is blown to a cloud, the e®ect is more obvious at the top (rather than the
bottom) of a cloud [17]. Thus, we add one more attribute to the particle, i.e. the vertical
distancefromtheapproximatecloudbase,whichisdeterminedbythedi®erencefromthe
initial heat °ux of the current particle to the minimum heat °ux of water vapor particles
in our model. The degree of horizontal movement due to the wind can be controlled by
this distance. In Fig. 3.5, we compare a pair of clouds with and without the wind e®ect
using this model.
3.4 Phase Transition Modeling by Voxel Dynamics
Besides thermal's movement, another important factor in cloud's formation is the phase
transition. Water vapor particles oscillate until condensation in the cloud formation
process. The latter is easy to implement using the voxel model with a grid structure. As
60
(a) (b)
Figure 3.5: Comparison of cloud shapes (a) without and (b) with the wind e®ect during
its formation.
stated in Sec. 3.1, we adopt a particle system in the ¯rst stage and the voxel dynamics
in the second stage for cloud synthesis. Users can choose the grid size of a bounding
box °exibly in the vorxel model. Its size is decided by the initial parameter range. The
larger the resolution, the ¯ner detail for observation. The water vapor/droplet density
in a voxel is updated at each time step. When the water phase transition occurs, the
released latent heat could be the new source of water vapor generation, which leads to
the fractal structure of clouds.
3.4.1 Water Continuity and Phase Transition
Cloud's EL is determined by the user. According to the saturation of the water vapor
density given in Eq. (2.7), one can generate a saturation density look-up table at the
initialization stage. The water vapor density in a voxel is compared with the saturation
densitytodecidewhetherthewaterphasetransitionoccursateachtimestep. Besidesthe
saturationpoint,thereisanotherfactorthata®ectsthewaterdropletdensity;namely,the
61
cloud condensation nuclei (CCN), which are also known as cloud seeds. They are small
particles about which cloud droplets coalesce. Water requires a non-gaseous surface to
makethetransitionfromavaportoaliquid. Intheatmosphere,thissurfacepresentsitself
astinysolidorliquidparticlescalledCCNs. WhennoCCNsarepresent, thewatervapor
can be supercooled below 0
o
C before droplets spontaneously form. To take this factor
into account, we assign a continuous random number to each voxel at the initialization
stage to indicate the distinct condensation capability.
We see from Fig. 3.2 that radius R is expanding until it ¯nally reaches a stable
value during thermal's life time. This process can be explained as follows. A thermal is
heated and lifted. Then, it gets expanded and cooled. For the same thermal, the bigger
it can reach, the less pressure di®erence between its inside and outside. Thus, it is more
probable for the water phase transition to occur, which corresponds to the death of a
water vapor particle. The water phase transition is primarily determined by radius R.
3.4.2 Fractal Structure and Entrainment
A procedural cloud yields good-looking image with a di®erent level of details since its
structureisfractal-like. Byde¯nition,fractalshaveaself-similarstructureatallscales. It
isassumedbymanyresearchersthatthefractal-likecloudappearanceisduetothefractal
nature of turbulence. Furthermore, the water vapor is not isolated in the atmosphere.
It actually attracts ambient mass and mixes with the environmental air. Then, the
environmental air becomes part of the cloud, which is called entrainment. It is observed
experimentally that entrainment happens mostly at the cloud boundary, where there is
more turbulence. This o®ers physical theory on the birth of a new thermal due to latent
62
heatreleaseatthecloudboundary. Thus,weimplementfractalsonlyatcloudboundaries
by assigning local water vapor particles a short lifetime.
Since latent heat is released in the second stage, the substructure is considered in
this stage. The latent heat released from the water phase change provides a local heat
source for buoyancy, which is governed by the similarity solution with the heat source
parameter from the phase change. The substructure is only generated at the surface by
the evidence of atmospheric science and the consideration of computational e±ciency.
For each voxel, we use two FLAGs to denote whether it contains the water vapor and the
water droplet. When the phase change occurs, we check positions (neighbors above the
current voxel) marked with a star in Fig. 3.6 with the boolean operation. We compare a
pair of clouds with and without the fractal structure in Fig. 3.7 while other parameters
are kept the same. It is obvious that the cloud with the fractal structure looks more
realistic to human eyes.
Figure 3.6: The potential direction of the substructure.
63
(a) (b)
Figure 3.7: Comparison of clouds (a) with and (b) without the fractal structure.
3.4.3 E±cient Update of Water Phase Change
The computation consists with two parts in our model. The ¯rst part is the update
of particles at every time step in a particle system. The second part is the update of
water phase change in the voxel stage. If we divide the bounding box into N
3
voxels in
the water phase change stage, the complexity is O(N
3
) as we need to traverse all voxels
to check whether the water vapor/droplet has changed. We can improve the update
procedure by only calculating voxels that have new water particles that die from the last
time step. By updating a Boolean check-in table for water phase change at every time
step,ourcomputationalcomplexitycanimprove10¡20%accordingtothedi®erentvoxel
resolution.
3.4.4 Continuous Density Calculation
Even by releasing water vapor particles and updating the phase transition at every time
step, and letting the initial horizontal velocity direction of each particle be uniformly
64
distributed in all directions, the condensed water droplet density in voxels may still not
be a smooth function. The discontinuity may result in some visual artifact. Thus, in the
¯nalstepofvoxelmodeling, weapplyasimplespatialsmooth¯lterofthefollowingform:
d sm(i;j;k) =
1
12
(d(i¡1;j;k)+d(i+1;j;k)
+d(i;j¡1;k)+d(i;j+1;k)
+d(i;j;k¡1)+d(i;j;k+1)
+6¤(d(i;j;k))); (3.6)
where d(i;j;k) is the original water droplet value and d sm(i;j;k) is the smoothed water
droplet value.
After applying the ¯lter, we have a smooth water droplet density ¯eld ready for
rendering. A two-pass illumination and rendering method is used in our work, which is
similartotheworkofDobashi[72]andHarris[22]. Atthecenterofeachvoxel,weplacea
particle whose splatting texture is selected according to the corresponding water droplet
density. Based on the possible density value, we pre-calculate 32 textures for look-up
during the rendering process.
3.4.5 Other Implementation Details
In addition to the above main implementation technique, there are two more implemen-
tation details to be discussed below.
² Billboarding
65
Billboarding is a technique that adjusts an object's orientation so that it faces
some target, which is usually the camera. This technique is popular in games and
applications that require a large number of primitives. Regardless of the speed of
a graphics card, it appears to be never su±cient. Billboarding can be used to cut
down the number of primitives by replacing geometry with an imposter texture to
model a scene. Thus, billboarding is used to provide an illusion of a 3D object in
cloud modeling in our proposed particle system.
² Splat Texture
Another technique is to use the pre-computed splat texture to account for the light
intensity loss as the light passes through the volume of cloud particles. As less light
passes through the center than edges, the texture exhibits a fall-o® from the center
to the edges. Here, we use an interpolating polynomial that varies from v
1
to 0 in
a smooth way of the following form
v =v
1
£(2f
3
¡3f
2
+1); (3.7)
where v
1
is a parameter and 0 · f · 1 is the distance between the target pixel
and the center normalized by the texture radius. Since the water droplet's density
is continuous, we may pre-compute the splat texture for each one. However, this
involves too much computation. Here, we quantize the value into N
q
levels, and
generate N
q
splat textures in advance. 32 splat textures are pre-computed in our
implementation. Fig. 3.8 gives the splat texture we use in our simulation model.
66
Figure 3.8: Illustration of the splat texture.
3.5 Cloud Illumination and Rendering
The two-pass rendering method is adopted for cloud illumination and rendering. The
particle model is used to compute complicate lighting interaction to achieve multiple
scattering illumination. First, clouds are stored as an array of particles to represent
water droplets. They have attributes including position, density, size, albedo (the ratio
of re°ected light to the incoming light) and extinction (reduction of light intensity). For
cloud lighting, we approximate light scattered in the light's direction (from the cloud
center to the sun). To do this, we sort particles away from the light position so that
the closest particle will be rendered ¯rst. This is obvious since the particle which is not
occluded in any way should receive the largest amount of light. After sorting particles,
we set up the camera to be placed in the sun's position, viewing the cloud center. The
projection should map the cloud onto the whole viewing port. We use the orthographic
projection since it will not deform far away particles.
67
The lighting equation as proposed by Harris [22] is given by
I
k
=
8
>
>
>
<
>
>
>
:
g
k¡1
+T
k¡1
£I
k¡1
(k >0)
I
0
(k =0)
(3.8)
This equation relates the intensity of the current particle I
k
with the intensity of the
previous particle I
k¡1
and the transparency of the previous particle T
k¡1
. The other
term is g
k¡1
, where
g
k
=a
k
£¿
k
£p(!;!
0
)£I
k
£°=4¼:
Since g
k¡1
is related to the intensity of the k¡1 fragment, we should read this value
from the frame bu®er. All other values in the above equation, such as albedo, extinction,
solid angle, are constants and p(!;!
0
) is a phase function which will be discussed later.
The original intensity I
0
is fully bright so that the bu®er is ¯rst set to white.
Eq. 3.8canbeencodedforgraphicshardwarethroughblending. Blendingcalculates
the sum of the incoming fragment multiplied by the source factor and the existing frag-
ment in the bu®er multiplied by the destination factor. In our case, the source factor is 1
(as g
k¡1
does not have any coe±cient) and the destination factor is T
k¡1
, the transmit-
tance of the fragment in the splat texture. The result of the blend operation is a color
which is then read back for the next particle. Since opacity is stored in the splat texture,
transmittance T
k¡1
is equal to one minus opacity. This provides the blend function with
proper parameters.
Thelightingprocesscanberoughlydescribedasfollows. Westartfromafullybright
bu®er and use the splat textures to decrease the luminance, thus \darkening" the color
68
of the particles as they get farther from the light. Then, the calculated color will be
stored for later rendering. When we splat particles for lighting, we need to include the
phase function. This is equal to 1.5 in the direction of the light so that we scale up the
color by this value before rendering this particle as a billboard. To simulate multiple
scattering in the eye direction, a phase function is used. The phase function allows the
calculationofthedistributionoflightscatteringalongthedirectionofincidentlight. The
phase function takes parameters from two directions; namely, the incident light direction
and the direction of light arriving at the viewer in our current case. The incident light
directionistheoppositedirectionofthelight. Thedirectionoflightarrivingattheviewer
is the vector between the particle position and the camera viewing point.
The following simple Rayleigh scattering function is used in the our lighting model
p(!;!
0
)=0:75£(1+(cosµ)
2
); (3.9)
where µ is the angle between ! and !
0
. When ! is in the direction of !
0
, the function is
equal to 1.5 giving the value used in lighting.
3.6 Simulation Results
In the cloud synthesis, we do not solve PDEs directly to save the computational cost.
The main computational cost comes from two parts:
² update of water vapor particles; where a water vapor particle may die (i.e. being
condensed into voxels) and it can be re-used afterwards; and
69
² water phase change, which depends on the voxel size.
The formation process of a cumulus is shown in Fig. 3.9. In this formation, we
control the real vertical range to a reasonable value so that the cloud anvil is very clear
to see. Also, from the cloud formation process, we see that the basic cloud shape is
kept during the formation since the source releasing water vapor is the same. This result
matches the experiment result in the lab from Scorer in 1957 [52]. Another formation
process of a cumulus cloud is illustrated in Fig. 3.10, where a cumulus cloud gradually
develops from Fig. 3.10(a) to Fig. 3.10(b) and ¯nally to Fig. 3.10(c). For this example,
a grid size of 64£64£32 is used in the voxel stage for phase transition modeling. Both
particle modeling and voxel modeling do not involve heavy computations. The cloud
synthesis speed is very fast. It takes around 9.4 seconds to get a mature piece cloud as
shown in Fig. 3.10(c) in a PC with Pentium 4 3.0GHz CPU & 1G RAM. To increase the
resolution and quality of rendered results, we may decrease the water vapor density for
each particle or increase the voxel grid number. Consequently, there are more particles
to render in the rendering pass. The tradeo® is a higher computational cost.
We show experimental results with di®erent parameter settings in Figs. 3.11 and
3.11. We present a pair of cloud results with a di®erent Buoyancy frequency value in Fig.
3.11. The buoyancy frequency values in Figs. 3.11(a) and (b) are equal to 0:5 and 0:01,
respectively. By the de¯nition of buoyancy frequency, it is the temperature di®erence for
a unit step vertically. A larger value in N means the temperature di®erence is bigger.
Then, there exists more mixing during water vapor's lifting so that the height a particle
can reach is lower. Since the buoyancy frequency in Fig. 3.11(a) is larger than that in
70
(a) (b)
(c)
Figure 3.9: Simulation of the cloud formation process, where steps (a)-(c) show how a
cloud is formed along time.
Fig. 3.11(b), the vertical extent of the resultant cumulus cloud becomes smaller while
other parameters are kept the same.
We show a pair of clouds with di®erent entrainment coe±cient values in Fig. 3.12.
ThemassentrainmentvaluesinFigs. 3.12(a)and(b)aresetto0.5and0.03,respectively.
The mass entrainment coe±cient represents the ability to interact with the environment
air. For a large ® value, there is more mixing with surrounding air, which in return
71
(a) (b)
(c)
Figure 3.10: Three sampled results of a cloud formation process.
limits the height a particle can reach. The ® value in Fig. 3.12(a) is larger than that in
Fig. 3.12(b). Thus, the vertical extent of Fig. 3.12(a) is smaller. Its oscillation height
is smaller and its ability to hold the water vapor is higher, which means the saturation
point for the water phase transition is higher than that in Fig. 3.10(b). As a result,
the cloud in Fig. 3.12(a) is not as °u®y as that in Fig. 3.12(b). Furthermore, the grid
number in the water phase transition stage increases to 80£80£64 in Fig. 3.12. As
compared with Fig. 3.11, there are more particles to render and more details to observe
72
in Fig. 3.12. The cloud synthesis time for the result shown in Fig. 3.12 is 27 sec. This
time can be further reduced to around 24 sec (about 11% improvement) using a more
e±cient water phase update as presented in Sec. 3.4.3.
We show a pair of cloud result under di®erent lighting condition in Fig. 3.13. Since
cloud in Fig. 3.13(a)is close to the Sun, the anisotropic lighting e®ect is very obvious.
The cumuluscloud under a di®erent lightingcondition is givenin Fig. 3.13(b), where the
shadow is more clear to view.
A scene containing several synthesized clouds is shown in Fig. 3.14. A close-up of
the scene in Fig. 3.14 is given in Fig. 3.15. Finally, the same cloud scene rendered at
sunset is shown in Fig. 3.16.
3.7 Conclusion
A new cloud simulation model based on the similarity approach was presented in this
chapter. An user can play the water vapor source on the ground freely to control the
generalcloudshape. Atwo-passrenderingalgorithmwasimplementedtoachievemultiple
scattering. Without solving the PDEs directly, we can represent the convective cloud
with a set of parameters to allow e±cient computation while catching the property of the
convective cloud realistically. In rendering, some modi¯cation has been done to improve
Harris' illumination model by using multiple spatting textures to represent the density
at di®erent levels. A couple of rendered cloud examples were shown to demonstrate
the visual e®ect of the proposed new methodology of cloud simulation, modeling and
rendering.
73
(a)
(b)
Figure3.11: Cloudsgeneratedwithdi®erentbuoyancyfrequencyvalues: (a) N =0:5and
(b) N =0:01.
74
(a)
(b)
Figure 3.12: Clouds generated by di®erent entrainment coe±cients: (a) ® = 0:5 and (b)
®=0:03.
75
(a)
(b)
Figure 3.13: Clouds rendered under di®erent lighting condition
76
Figure 3.14: A scene of multiple clouds.
Figure 3.15: The close up of multiple clouds.
77
Figure 3.16: The same cloud scene at sunset.
78
Chapter 4
A Decoupled-PDE Approach to Cloud Simulation with
Volumetric Modeling
Time-lapse photography reveals that cumulus clouds develop by a
series of rising plumes or turrets.
|Atmosphere convection
4.1 Introduction
Due to the contributions of Euler, Navier and Stokes from 1750's to 1850's, people have
a good understanding of the cloud phenomenon nowadays. The cloud is a mixture of
air, water vapor and condensed water droplet. The system of Navier-Stokes equations
provide an accurate mathematical model for °uid °ows in the nature with the cloud as a
special example. However, these equations can only be solved analytically under a very
simple setting, and most practical solutions are obtained numerically. In the context of
computergraphics, itisworthwhiletotradeaccuracyforcomputationale±ciencyaslong
as the visual appearance of the rendered cloud is convincing to human eyes.
79
FosterandMetaxas[19]presentedahotgassimulationwithexplicitEulerianschemes,
where they introduced solutions from numerical °uid dynamics to computer graphics.
however, the time step used in their simulation was restricted to a small value due to the
stability consideration. Later, Stam [56] proposed a stable °uid solver, which can choose
a much larger time step in solving these PDEs. Stam's stable °uid solver has a great
impact on physics-based computer graphics simulation.
In Chapter 3, we introduced a graphic approach for cloud rendering that begin with
ananalyticsolution of cloudformation withseveralcontrolparameters. Suchananalytic
solution serves as a basis for the simulation, modeling and rendering of clouds. Di®erent
cumulus cloud shapes can be synthesized by varying these parameters. In this chapter,
weconsider another graphicapproach for cumulus cloud synthesiswhich relies on the nu-
merical solution of the governing partial di®erential equations (PDEs). However, we take
another short cut to reduce the computational cost. That is, we consider the decoupling
of the horizontal- and vertical- dimensions in the °uid solver for cloud simulation. In this
decoupled scheme, layers for cumulus clouds are generated along the horizontal direction
while water vapor input is synthesized to account for the plume characteristics of clouds
along the vertical direction. Each cloud layer requires the simulation of a continuous
volume of °uid on a two-dimensional grid.
4.2 Mathematical Background
Mathematically, the state of a °uid at a given time is modeled as a velocity vector ¯eld
that associates a velocity vector to a point in the 3D space. The system of Navier-Stokes
80
equations provides a precise description of how a velocity ¯eld evolves along time. In the
model,arectangular3Dspaceisuniformlypartitionedbyagrid(i.e. theCartesiangrid),
each of them has a 2D vector u(t) = (u(t);v(t)) to represent the velocity at a discrete
point x=(x;y).
To examine the °uid motion by the Lagrangian formulation, i.e., we apply Newton's
Law to a moving particle using a ¯xed coordinate system. The rate of velocity change
comesfromtwoparts: oneduetothechangeintimewhiletheotherduetothechangein
space. Once, the velocity ¯eld is obtained, we are able to determine other quantities such
as heat, densities, cloud components, etc. Thus, the velocity ¯eld is the key component
in cloud simulation. The governing Navier-Stoke equation for the velocity ¯eld can be
expressed as
@u
@t
=¡(u¢r)u+ºr
2
u¡
1
½
rp+f; (4.1)
with the following mass continuity constraint
r¢u=0; (4.2)
where º is the kinematic viscosity of the °uid, ½ is the °uid density and f is the external
force. The ¯rst three items on the right-hand-side of Eq. (4.1) are the acceleration due
to advection, the di®usion and the pressure gradient, respectively.
Aftergettingthevelocity¯eldfromtheNavier-Stokeequation,wecandeterminethe
density and other scalar quantities via
@d
@t
=¡(u¢r)d+·r
2
u+Sr¡®d; (4.3)
81
( , ) 0 0 ( , ) 0 1 ( , ) 0 2 ( , ) 0 3
( , ) 1 0 ( , ) 1 1 ( , ) 1 2 ( , ) 1 3
( , ) 2 0
( , ) 3 0
( , ) 2 3 ( , ) 2 2 ( , ) 2 1
( , ) 3 3 ( , ) 3 2 ( , ) 3 1
Δx
Δy
Figure 4.1: Illustration of an uniform grid for °uid simulation.
where · is the di®usion constant, Sr is a source term and ® is the dissipation rate.
Before discussing each term in the Navier-Stoke equation, we have a brief review of
vector calculus below. In Eqs. (4.1) and (4.2), symbol \¢" (dot) denotes a dot product
between vectors while symbol r (nabla) is the vector of spatial partial derivatives, i.e.
r=
³
@
@x
;
@
@y
´
for the 2D case. Discrete approximations to several 2D gradient operators
derivedfromthemaregiveninTable4.1,wheresubscripts(i;j)denoteadiscreteposition
in Fig. 4.1, and 4x and 4y are the spatial spacings along the x- and the y-direction,
respectively.
Operator Abbrev. De¯nition Finite Di®erence Form
Gradient grad f rf =
³
@f
@x
;
@f
@y
´
rf =
³
@f
i+1;j
¡f
i¡1;j
24x
;
@f
i;j+1
¡f
i;j¡1
24y
´
Divergence div F r¢F=
@u
@x
+
@v
@y
r¢F=
u
i+1;j
¡u
i¡1;j
24x
;
v
i;j+1
¡v
i;j¡1
24y
Curl N/A curlF=r£F curlF=
v
i+1;j
¡v
i¡1;j
24x
;
u
i;j+1
¡u
i;j¡1
24y
Laplace N/A r
2
f =
@
2
f
@x
2
+
@
2
f
@y
2
r
2
f =
f
i+1;j
¡2f
i;j
+f
i¡1;j
(4x)
2
+
f
i;j+1
¡2f
i;j
+f
i;j¡1
(4y)
2
Table 4.1: Discrete approximations to several 2D gradient operators.
82
Thegradientofascalarisavectorfunctionofapositionobtainedbyapplyingthedel
operator to a scalar function. The divergence is the scalar (or dot) product between the
del operator and F, which has an important physics meaning, i.e., the extent to which
the vector ¯eld °ow behaves like a source or a sink at a given point. The divergence
operator is applied to the velocity ¯eld of a °uid °ow to measure the net change of
the °uid amount around a small spatial volume in Navier-Stoke equations. Under the
incompressibility condition, we obtain the mass continuity equation in Eq. (4.2), which
ensures that the °uid always has zero divergence. This is equivalent to that the point of
concern is neither a source nor a sink.
After introducing the two basic vector operators (i.e. gradient and divergence), we
can combine them to get the divergence of a gradient byr¢r=r
2
, which yields the fa-
moussecond-orderLaplacianoperator. TheLaplacianoperatoriscommonlyencountered
in physics e.g. the di®usion and wave equations. The Laplacian operator captures the
`smoothness' of a function. It measures the di®erence between the value of f at a point
anditsmeanvalueatsurroundingpoints. Inthesteadystate,boththedi®usionequation
and the wave equation are reduced to a simple form as r
2
f = 0, which is known as the
Laplace equation. When the right-hand-side this equation is non-zero, i.e. r
2
f = a, it
is the Poisson equation.
4.3 Stam's Stable Fluid Solver
Stam's stable °uid solver [56] is presented in this section.
83
4.3.1 Helmholtz-Hodge Decomposition
A smooth vector ¯eld w can be uniquely decomposed into a set of vector by [10]
w=rÁ+r£Ã+h; (4.4)
where Á is a scalar potential ¯eld satisfying r£(rÁ) = 0, and Ã is a vector potential
¯eld satisfying r¢(r£Á) = 0 and h is a harmonic vector ¯eld which satisfy r¢h = 0
andr£h=0 with appropriate boundary conditions. However, for incompressible °ows,
only two ¯elds are su±cient to represent the vector ¯eld; namely,
w=rÁ+h: (4.5)
In other words, under the incompressible assumption, the vector ¯eld is decomposed into
a divergence free vector and a gradient of a scalar ¯eld. We rede¯ne these two vectors by
u and p as
w=u+rp: (4.6)
The mass continuity assumption in Eq. (4.2) implies that there is no sink or source
¯eld in the ¯nal velocity ¯eld at each time step. Equivalently, we need a divergence-
free velocity ¯eld. By combining this with Eq. (4.2), which assumes a zero-divergence
¯eld r¢u = 0, Stam de¯ned an operator P that projects the vector ¯eld w to the
divergence-free part u=Pw. Following Eq. (4.6), we have
Pw=Pu+P(rp):
84
Furthermore, the projection of a vector on itself is itself so that Pu = u. Then, we
obtain P(rp) = 0. The operator P helps eliminate a factor in Eq. (4.6). Multiplying
both sides in Eq. (4.6) by operatorr, we have
r¢w=r
2
(p): (4.7)
The resulting equation is the Poisson equation for the scalar ¯eld p with the Newmann
boundary condition
@p
@n
=0 or, equivalently, u¢^ n=0 on @D.
After solving Eq. (4.7), we can get the divergence-free velocity ¯eld via
u=w¡rp: (4.8)
By applying the projection operator to Eq. (4.1), we obtain a single equation for the
velocity
P
@u
@t
=P
µ
¡(u¢r)u+ºr
2
u+f¡
1
½
rp
¶
: (4.9)
Together with terms P(rp)=0 and Pu = u, Eq. (4.9) can be rewritten as
@u
@t
=P
³
¡(u¢r)u+ºr
2
u+f
´
: (4.10)
4.3.2 Navier-Stoke Solver
Eq. (4.10) leads to a straightforward way to solve the Navier-Stoke equation. The
corresponding°owchartisshowninFig. 4.2,whichbeginswithinitialvelocityW0. After
85
W0 W1 W2 W3 W4
Force
Application
Advection Diffusion Projection
Figure 4.2: The °owchart of the Navier-Stoke solver.
force application, advection, di®usion and projection, we get a divergence-free velocity
¯eld. These procedures are detailed below.
4.3.2.1 External Force
Force accelerates the °uid. During the simulation, we need °exibility to control the °uid
movement, and force can be applied for this purpose. Given force f, we can solve for the
velocity ¯eld using
@u
@t
=f: (4.11)
The numerical solution to (4.11) can be achieved using the ¯nite di®erence method:
u
t+4t
=u
t
+4t¢f:
Temperature di®erence, buoyancy, gravity and etc. can all be sources of external force
that a®ect °uid movement.
4.3.2.2 Advection
Advectionreferstothetransportofsomepropertyand/orsubstanceofa°uid. Advection
iscalculatedbyasemi-LagrangianmethodinStam's°uidsolver. Insteadofkeepingtrack
86
U t ( )
X
U t Δt ( )
Figure 4.3: The advection in simulation:(a) the velocity ¯eld, (b) the path of a particle
and (c) to ¯nd the velocity ¯eld backwards from t+4t to t.
of the movement of individual °uid particles, we divide a rectangular space uniformly
into voxels, and assign a particle to each voxel center. At each time step, we update the
velocity and acceleration for each grid point as shown in Fig. 4.3 (a) & (b). This can be
achieved by Euler's method as
X
t+4t
=X
t
+u
t
4t: (4.12)
Thismethod issimple forexplicit(orforward)integrationofdi®erentialequations. How-
ever, there is a problem with this method. That is, when the time step is large, the
explicit method will yield an unstable result.
Instead of traversing the velocity ¯eld in the forward direction, an alternative yet
stable solution was proposed by Stam [56], which traces the trajectory of a particle at
each grid point backwards as shown in Fig. 4.3 (c). Suppose that the desired location
is X. Then, by interpolating the actual value of four grid point around position X at
previous time step, we get the particle velocity at a grid point updated via
X
t+4t
(~ x)=X
t
[~ x¡u
t
4t] (4.13)
87
Figure 4.4: Velocity Di®usion in °uid simulation
Eq. (4.13) is the ¯rst-order Euler method. There are more accurate methods, such as
the second-order Runge Kutta method and the midpoint method, etc.
4.3.2.3 Di®usion
Viscous di®usion characterize °uid's resistance to °ow. In Fig. 4.4, we show a velocity
di®usion process in two time steps. Velocity di®usion is often modeled by the following
heat equation
@u
@t
=ºr
2
u: (4.14)
There are several methods to solve the di®usion equation. However, to achieve stability,
we should follow stam's implicit method and solve the problem backwards in time as:
(I¡º4tr
2
)u
t+4t
(~ x)=u
t
(~ x); (4.15)
where I is the identity operator. When the di®usion operator is discretized, we are led
to a sparse linear system of unknowns u
t+4t
. There is however an e±cient way to solve
this system with the Helmholtz-Hodge decomposition as described before.
88
4.3.3 Force Extension
Several force terms are considered in our simulation.
4.3.3.1 Virtual Temperature Force
The primary convection process is driven by the buoyancy force from the ground upward
alongtheverticaldirection. Inthedecoupledmodel,wenolongusethebuoyancyforceas
a direct driving force, since its e®ect has already been incorporated in the plume model.
However, along the horizontal dimension, the heat di®usion and virtual temperature
di®erence will result in the movement of water vapor particles, where the heat di®erence
is due to the di®erence in water vapor distribution. A simple way to incorporate this
e®ect is to add a new force term along the direction with the largest local gradient; i.e.,
F
vt
=¿(T
v
¡T
vi
); (4.16)
where T
v
is the virtual temperature of the current voxel in a grid, and T
vi
is the virtual
temperatureofitsneighborthathastheminimumvirtualtemperature,and¿ isaconstant
scale factor.
4.3.3.2 Coriolis Force
The coriolis force is one of the main forces to be considered in the horizontal direction as
shown in Table 4.2. In the decoupled cloud simulation, the coriolis force is comparable
89
with other contributors like the pressure gradient etc. Thus, we need to incorporate this
factor. The Coriolis force can be expressed as
F
cor
=¡fk£u; (4.17)
where f =2sin' is the coriolis parameter and u is the horizontal velocity vector.
4.3.3.3 Vorticity Con¯nement
Fluids such as smoke and convective clouds typically contain rotational °ows at a variety
of scales. As explained by Fedkiw et al. [18], numerical dissipation in simulation on a
coarse grid damps out these interesting features. Vorticity con¯nement is used to restore
these motions. It works by computing vorticity ~ ! =r£~ u ¯rst. From the vorticity, we
obtain a normalized vorticity gradient ¯eld
~
ª=
~ ´
j~ ´j
; where ~ ´ =r~ !:
Vectors in this vector ¯eld point from areas of lower vorticity to areas of higher one. The
vorticity con¯nement is computed from equations as
F
vc
=¾(
~
ª£~ !); (4.18)
where ¾ is a user-control scale parameter.
90
4.4 Decoupled-PDE Cloud Simulation
4.4.1 Turbulent Plume
The characteristics of the convective plume were described in Chapter 2. To get the sim-
ilarity solution of a turbulent plume in a stable strati¯ed °uid, we assume that the radial
pro¯le of velocity and buoyancy are geometrically similar at all heights. Mathematically,
this is equivalent to saying that the vertical and radial dependencies are separable, and
that the radial dependencies of each variable are the same.
Following the turbulent plume assumption in Chapter 2, the set of equations to be
solved is Eqs. 2.18, 2.19 and 2.20. Similarly to the solution to turbulent thermal, we can
rephrase variables to get dimensionless ones as
V ´Rw; U ´R
2
w; F ´R
2
wB:
At the source, the width and the momentum of the plume are assumed to vanish while
the buoyancy °ux is equal to that of the source; namely,
V =U =0 and F =
2
¼
F
0
at z =0:
There is no closed-form solution to this set equation. With the numerical solution
for dimensionless variables, we can get the vertical velocity w, radius R and buoyancy
B as shown in Fig. 4.5. The vertical velocity vanishes at a dimensionless height of 2.8
whilethebuoyancy¯rstvanisheswhen z =2:125. Theindividualparticlesof°uidwithin
the plume overshoot the level at which their buoyancy ¯rst vanishes, then decelerate to
91
the zero velocity while presumably spreading away from the central axis. The bulk of
the spreading should occur between the levels of vanishing buoyancy and zero vertical
velocity.
1 2 3 4 5
0
0.5
1
1.5
2
2.5
3
R/W/B
Height(z)

R−z
w−z
B−z
Figure 4.5: The height dependence of several dimensionless quantities of a turbulent
plume: horizontal radius R, vertical velocity w and buoyancy B.
The resulting solution leads to a successful representation of the convective process
that matches the experiment data. The photographs of time exposures of a plume in a
stable strati¯ed °uid at the early and the late stages in its development from Morton,
Taylor and Turner [34] are shown in Fig. 4.6.
Thedimensionlesssolutiondepictsthecharacteristicsofaturbulentplume. Aplume
is determined by three independent parameters N, F
0
and ® as
z
¤
=2
¡7=8
¼
¡1=4
®
¡1=2
F
1=4
0
N
¡3=4
z; (4.19)
92
Figure 4.6: The development of a turbulent plume: the early stage (left) and the late
stage (right).
where z is dimensionless height, N is atmosphere's strati¯cation, F
0
is the heat source,
and ® is a constant which is proportional to the entrainment of mass same as that in
Chapter3.
There are physical justi¯cation to decouple the horizontal and vertical phenomena
and treat them independently for simpli¯cation. Table 4.2 gives the acceleration con-
tributions from di®erent factors based on the following dynamic governing equation in
atmospheric science
dv
dt
=¡fk£v¡rÁ¡
1
½
a
rp
a
+
´
a
½
a
r
2
v+
1
½
a
(r¢½
a
K
m
r)v: (4.20)
Weseefromthistablethatthescalesoftheverticaldirectionandthehorizontaldirection
behave quite di®erent in cloud formation. The scale di®erence justi¯es their di®erent
treatment in computer graphics.
93
Term Acceleration or Force/ Horizontal Accel. Vertical Accel.
Mass Expression (ms
¡2
) (ms
¡2
)
Local ¹ a
l
=
d
¹
v
dt
=
@
¹
v
@t
+(¹ v¢r)¹ v 10
¡4
10
¡7
»1
acceleration
Coriolis force
F
c
M
a
=¡fk£v 10
¡3
0
per unit mass
E®ective gravitational
F
g
Ma
=¡rÁ 0 10
force per unit mass
Pressure gradient
F
p
Ma
=
1
½a
rp
a
10
¡3
10
force per unit mass
Viscous force
F
v
M
a
=
´
a
½
a
r
2
v 10
¡12
»10
¡3
10
¡15
»10
¡5
per unit mass
Turbulent °ux divergence
F
t
M
a
=
1
½
a
(r¢½
a
K
m
r)v 0»0:005 0»1
of momentum
Table 4.2: Acceleration contributions from di®erent factors.
4.4.2 Decoupled Cloud Simulation
Based on the discussion in the previous subsection, we propose a decoupled cloud sim-
ulation framework as shown in Fig. 4.7. The cumulus cloud is formed by quite a few
plumes, and we can develop several turbulent plumes together. Due to the similarity
assumption in turbulent plumes (i.e. the radial pro¯les are similar at all heights), we
select several key frames for the horizontal development. Over the horizontal domain,
we solve the Navier-stoke equation driven by the virtual temperature di®erence. The
vertical dimension relationship is linked by convective plumes at di®erent heights. Hori-
zontal heat di®erence is caused by the virtual temperature di®erence of a distinct water
vapor source distribution with Eq.2.3. Finally, we render the cloud using the 3D texture
rendering method with volume shadowing.
For initialization, a user can choose a set of independent variables freely. This pa-
rametersetincludesatmosphere'sstrati¯cationN, heatsourceF
0
, andmassentrainment
94
constant®. IfF
0
isbigger, itmeanstheplumehasmoremomentumsothatithastheca-
pabilitytoreachahigherelevation, leadingtoacumuluscloudofalargerverticalextent.
For a larger value of ®, it means that more entrainment mass is involved and there is
moremixingwiththeenvironmentalatmosphere. Then, ithasalowerlikelihoodtoreach
a high vertical extent. The parameter N is a measurement of atmosphere strati¯cation.
If N is large, the potential temperature di®erence is big for the same height step and it
will result in more heat exchange, more mixing, and a lower vertical extent. Through
this initialization, a user can control the general shape of a cloud with great °exibility.
Key fr ame
User choic e of general s hape c ontrol
Horiz ont al
N avier-St oke Equati on
W ater phas e c hang e
V er t ic al
W ater v apor di st r i bu t ion
and i np ut
Int erp ol a t ion of ot her layer s
an d add perturbati o n
Rend ering
Figure 4.7: The proposed decoupled cloud simulation framework.
Although the 3D simulation is decoupled into horizontal and vertical directions, the
simulation along the horizontal direction still follows the dynamic °uid °ow PDEs. Our
horizontal simulation tool is built on top of a stable °uid simulator, which is very similar
95
to that of Stam [56]. It solves the Navier-Stoke equation with the following two-step
technique.
1. The semi-Lagrangian advection technique proposed by Stam is used to compute
the intermediate velocity ¯eld ~ u, external force including virtual temperature force
by Eq. (4.16), coriolis force by Eq. (4.17) and vorticity con¯nement force by Eq.
(4.18). Then, this step solves Eq. (4.1) without the pressure term.
2. The intermediate ¯eld ~ u is made incompressible (so that it satis¯es Eq. 4.2) using
the Helmholtz-Hodge decomposition. The water vapor is injected to the key frame
layer at each time step. In the meanwhile, the water vapor moves along the veloc-
ity ¯eld, coagulates to a certain level decided by the environmental pressure and
temperature with Eq. (2.7), and condenses to the water droplet. Afterwards, water
droplet and water vapor continue to move along the velocity ¯eld.
4.4.3 Implementation Details
4.4.3.1 Vertical Water Vapor Distribution
The plume has a di®erent developing stage, e.g. the early stage and the later spreading
stage. Cloud formation is a physical process involving both dynamic water vapor move-
ment and water phase change. When the water vapor density coagulates to a certain
level, the water phase change begins to happen. Since we need more water vapor sources
toenablewaterphasechange,wemaysaythatcloudformationhappensatthelateplume
spreading stage. As shown in Fig. 4.6, the spreading stage happens between levels of
96
vanishing buoyancy and zero vertical velocity, which are determined by the user through
variables N and F
0
.
Another clue that helps determine the water vapor distribution is the atmospheric
fact. The horizontal water content distribution is of either the \Gaussian" or the \top-
hat" pro¯le. We can use the plateau-like kernel to give these e®ects [43]. The kernel
equation to decide the water vapor distribution at each layer is given by
'(r)=e
¡a¢
(
r
b
)
n
1¡
(
r
b
)
n
; for 0·r·b; (4.21)
where b is the radius at the height where buoyancy begin to vanish, a is the width of the
kernel, and n controls the slope of drop-o®. As shown in Fig. 4.8, a larger a value makes
a narrower kernel while a larger n value generates a plateau-like kernel.
n=2
a=1
n=2
a=8
n=3
a=5
n=3
a=20
n=5
a=4
n=5
a=50
Figure 4.8: The plateau-like kernel for water vapor distribution.
97
Figure 4.9: A horizontal cloud layer comparison (a) with and (b) without the Coriolis
force.
4.4.3.2 Coriolis Force
The Coriolis force, which is computed by Eq. (4.17), is introduced since it cannot be
neglected along the horizontal dimension. Fig. 4.9 shows a comparison with and with-
out this term for a horizontal layer. The incorporation of the Coriolis force makes the
rotational feature in the horizontal layer clearer.
4.4.3.3 Wind E®ect
As mentioned in Chapter 3, the top of a cumulus cloud has a larger turbulent e®ect than
that of the cloud base. This characteristics is easy to generate in our decoupled model.
At a di®erent cloud layer, we can set a di®erent air resistance parameter to achieve this
wind e®ect; namely, the cloud top has a smaller air resistance parameter while the cloud
base has a bigger one. The wind e®ect (from left to right) can be observed in Figs. 4.11
(e) and (f).
98
4.4.3.4 Interpolation Between Key Cloud Layers
Under the assumption of geometrical similarity at all heights, we do not have to perform
the horizontal computation at each height. Instead, we choose several cloud layers along
the vertical axis, which we call `key cloud layers', for the horizontal simulations. The
density values of the remaining cloud layers are obtained through interpolation with
noise perturbation. The number of key cloud layers is a user-de¯nable parameter.
4.4.3.5 Boundary Condition
Decoupling simulation has no vertical boundary limits since plume has vertical extents.
In the horizontal direction, we let the velocity vanish at boundary.
4.4.4 Cloud Rendering
The water droplet density is stored in a grid based volume, which can be viewed as the
3D texture so that the 3D texture rendering technique can be easily applied. The simple
3D illumination can be done by casting shadows inside a volume. Behrens [4] proposed
a method to calculate the shadow of 3D texture, which takes into account of some key
aspectsofshadowswithrespectofsemi-transparentmaterial. Generallyspeaking,opaque
objects cast stronger shadow than semi-transparent ones, since one would expect more
light to pass through a semi-transparent object than an opaque subject.
The algorithm for calculating shadows work with 3D texture represented as object
oriented slices. The slices are processed in order starting with the one that is closest to
light. The light is assumed to be directional so the light direction is the same at every
point of a surface. Behrens' method is summarized below.
99
1. Simply draw a single slice p
i
into the bu®er.
2. Get s
i
by blending p
i¡1
into s
i¡1
via
s
i
=p
i¡1
¯s
i¡1
:
3. Get p
¤
i
by blending s
i
into p
i
via
p
¤
i
=s
i
¯p
i
:
4. Read the resulting slice p
¤
i
out of the bu®er.
In this algorithm, p
¤
i
is the ¯nal shadowed version of slice p
i
, and \x¯y" should be
interpreted as \blending x into y". The ¯rst blend adds the shadow cast by slice p
i¡1
to
the accumulated shadow cast by all previous slices. The second blend function applies
the new calculated shadow to the current slice to produce the shadowed version of the
slice. The results are independent of the viewing direction but dependent on the light
direction. When the light direction changes, the shadow has to be calculated again.
There are two blending functions in Behrens' method, where the alpha blending
technique is used. Strong shadow and no shadow e®ects are resulted with ® = 1 and
®=0, respectively. The blending function in Step 2 is achieved by
®=1¡(1¡®
s
)(1¡®
d
); (4.22)
100
where ®
d
is the alpha value in s
i¡1
and ®
s
is the alpha value of slice p
i¡1
. Note that
®
s
=1 and 0 denote that the slice is opaque and completely transparent, respectively.
The blending function in Step 3 is
c=(1¡®
d
®
s
)c
d
; and ®=®
d
; (4.23)
where ®
d
is the alpha value of slice p
i
and ®
s
is the alpha value of shadow slice s
i
.
Although the basic idea of Behrens' method is simple, its implementation consists
of many details. For the two blending functions in Behrens' method, there is no simple
OpenGLblendingfunctiontoimplementthemdirectly. Thus,wehavetomodifyBehrens'
algorithm to make it easy to implement with OpenGL. The modi¯ed Behrens' algorithm
consists of the following six steps, and Fig. 4.10 summarized the algorithm:
² For all slices from p
1
to p
n
,
1. Draw image slice p
i
into bu®er1;
2. Get s
¤
i
by blending p
i
into s
i
in bu®er3 with 4:25;
3. Get p
¤
i
by blending s
¤
i
into p
i
in bu®er1 with 4:26;
4. Get shadow s
i+1
by blending p
i
to shadow s
i
in bu®er2 with 4:27;
5. Make a copy of new shadow s
i+1
in bu®er3;
6. Read shadowed slice p
¤
i
out of bu®er1.
Take a close view at Eq. (4.23) where treats ®
s
= 1 as a strong shadow in the
original algorithm, in order to implement the algorithm in OpenGL, we have to invert
101
the meaning of source alpha value, making that ®
s
= 1 means no shadow [62]. Now Eq.
(4.23) becomes as:
c = (1¡®
d
®
s
)c
d
= (1¡®
d
(1¡®
s
))c
d
= (1¡(®
d
¡®
s
®
d
))c
d
: (4.24)
Step 3 in the modi¯ed algorithm is to compute
® = ®
s
¡®
s
®
d
= (1¡®
d
)®
s
; (4.25)
and it is implemented in OpenGL with the blend function
glBlendFunc(GL ONE MINUS DST ALPHA, GL ZERO).
Step 4 in the modi¯ed algorithm is to compute
c = (1¡®
s
)c
d
® = ®
d
; (4.26)
and it is implemented in OpenGL with the blend function glBlendFunc(GL ZERO,
GL ONE MINUS SRC ALPHA).
102
Finally, Step 5 in the modi¯ed algorithm is to compute
® = (1¡®
s
)®
d
(4.27)
and it is implemented in OpenGL with the blend function
glBlendFunc(GL ZERO,GL ONE MINUS SRC ALPHA).
As suggested by Fig. 4.10, all of those operations are carried out in the frame bu®er.
…
…
VolumeV
p1
pi
pn
Buffer1 Buffer2 Buffer3
pi si* si
VolumeV*
P1*
Pi*
Pn*
…
…
1
2
3
5
4
6
…
…
VolumeV
p1
pi
pn
…
…
…
…
VolumeV
p1
pi
pn
Buffer1 Buffer2 Buffer3
pi si* si
Buffer1 Buffer2 Buffer3
pi si* si
VolumeV*
P1*
Pi*
Pn*
…
…
VolumeV*
P1*
Pi*
Pn*
…
…
…
…
1 1 1
2 2 2
3 3 3
5 5 5
4 4 4
6 6 6
Figure 4.10: OpenGL implementation of modi¯ed Behren's Algorithm
4.5 Synthesized Cloud Results
The formation process of a cumulus is shown in Fig. 4.11. Fig. 4.11 (a)-(d) show the
cumulus formation process, and Fig. 4.11 (e)-(f) shows the cumulus cloud under a wind
direction from left to right. In this implementation, the horizontal size is 64£64, 10 key
framelayersischosen. Thecomputationalcostismuchlowerwhilethecharacteristicsfor
thecumuluscloudiswellcaptured. ThecomputationisconductedonaPentium43.0GHz
103
machine. Since we decompose the 3D computation to the computation in the 1D vertical
and the 2D horizontal directions, the complexity is reduced signi¯cantly. For a grid of
size 128£128£32, we need an average of 3:01 seconds for each iteration while Harris'
method demands 4:16 seconds on a comparable machine. For a grid of size 64£64£64,
our algorithm demands 1:13 seconds for each iteration while Harris' algorithm needs 1:75
seconds [30].
Furthermore, in our cloud simulation, each vertical key frame layer does not have to
`talk' to each other, which is controlled by the plume similarity solution independently.
Thus, parallel computing can improve the speed furthermore. Meanwhile, the fast 2D
°uid solver on GPU can help reduce the computational cost. The grid of size 128£128
is used in Fig. 4.12, where more details is observed.
To illustrate the e®ect of constant parameters on the cloud shape, we show in Fig.
4.13 three cumulus clouds synthesized with di®erent sets of parameters. In Fig. 4.13,
each cumulus cloud is composed of three plumes. For the example in Fig. 4.13(a), the
values of heat °ux for the three plumes are 2:3, 5:5 and 2:5, respectively; the buoyancy
frequency is 0.01; the entrainment coe±cient is 0.3; the width and the slope drop-o® of
the plateau kernel are 1 and 2, respectively. For the example in Fig. 4.13(b), we keep
the plateau kernel parameters the same as those in Fig. 4.13(a), and only change the
heat °ux values for the three plumes to 1:6, 9:8 and 2:0. As compared with Fig. 4.13(a),
a bigger vertical extent of the cloud is observed in Fig. 4.13(b) due to a bigger heat
°ux value. In Fig. 4.13(c), all atmospheric parameters F, N , ® are kept the same as
those in Fig. 4.13(a), while the width and the slope drop-o® of the plateau kernel are
set to 50 and 5, respectively. We observe that the cloud in Fig. 4.13(c) is more °u®y
104
than that in Fig. 4.13(a) since the width of the plateau kernel is closer to `top-hat' for
the cloud in Fig. 4.13(c). It is demonstrated in Fig. 4.13 that, with the proposed cloud
simulator, the general shape of a cumulus cloud can be °exibly controlled by adjusting
the parameters [66].
AskyscenewithmultiplecloudsispresentedinFig. 4.14,wherecloudsaregenerated
using four di®erent sets of parameters. The same scene under a di®erent illumination is
rendered in Fig. 4.15. The visual results as shown in Fig. 4.14 and 4.15 demonstrate
the potential of our algorithm for realistic cloud modeling and illumination. The cloud
imagesin Fig. 4.14 and4.15 are rendered with Behrens's method [4]. The rendering time
is less than ¯ve seconds for each scene.
4.6 Conclusion
A new cloud simulation process based on horizontal- and vertical-decoupling was pro-
posed to reduce the computational cost. The synthesized cloud was rendered using the
3D texture rendering technique with volume shadowing. The proposed decoupled-PDE
approach to cloud simulation allows a user the °exibility to control the general shape by
choosinginitializationparametersandthehorizontalwatervapordistribution. Then, the
2D PDE of key frame layers is solved by the stable °uid solver at each time step and the
plateau-like kernel water vapor distribution source is injected at each vertical cloud layer
according to the turbulent plume similarity solution. The wind e®ect to the cumulus
cloud is also easy to synthesize by choosing the air-resistance parameter in our proposed
method.
105
(a) (b)
(c) (d)
(e) (f)
Figure 4.11: Cloud formation with a grid of size 64£64.
106
Figure 4.12: A synthesized cumulus cloud with a grid of size 128£128.
107
(a) (b)
(c)
Figure 4.13: Examples of cumulus cloud with di®erent parameter sets.
a)F=f2.3, 5.5, 2.5g;N=0.01, ®=0.3;n=2, a=1;
b)F=f1.6, 9.8, 2.0g;N=0.01, ®=0.3;n=2, a=1;
c)F=f2.3, 5.5, 2.5g;N=0.01, ®=0.3;n=5, a=50.
108
Figure 4.14: A scene with multiple cumulus clouds.
Figure 4.15: A scene with multiple cumulus clouds at sunset.
109
Chapter 5
Real-Time Hierarchical Cloud Simulation based on
Similarity Approach
From the behavior of cumulus cloud, it consists largely of isolated
massesofbuoyantairrisingintoandmixingwiththeirsurroundings
from intermittent sources on the ground or within clouds
|Atmosphere convection
5.1 Introduction
A novel real-time hierarchical cloud simulation is proposed in this chapter. Our moti-
vation comes from the early study of the cloud model in atmospheric science. Although
atmosphericcloudstudyisdominatedbythree-dimensional(3D)numericalanalysisnow,
the early thermal model still provides an useful and intuitive explanation of the cloud
formingprocess. Beingabasicconvectiveentity,thermalstudybasedonthesimilarityap-
proach catches the characteristic behavior of clouds. As a cloud moves up, self-similarity
keeps its rough shape while its front edge has various appearance caused by the inner
110
convective circulation of the thermal. It is straightforward to derive a hierarchical model
for cloud simulation based on the above description. That is, the general shape and
movement of a cloud can be governed by the similarity approach as described in sec. 2.2
while its detail structure caused by thermal inner dynamic motion can be modeled as a
weak vortex ring.
The rest of this chapter is organized as follows. Requirements for real-time °uid
simulation is stated in Sec. 5.2. Thermal's similarity analysis is presented in Sec.5.3.
The behavior of the convective entity modeled by the vortex ring is described in Sec.5.4.
Thehierarchicalcloudsimulationframeworkanditsimplementationdetailsareexplained
in Sec.5.5. Finally, simulation results are given in Sec.5.6.
5.2 Requirements of Real-Time Simulations
With the increasing computational power of tomorrow's game consoles and PC's, more
and more physical e®ects (so far only seen in feature ¯lms) will be adopted in real-time
computergames. However, thereisstillagapbetweenpre-computingaphysicale®ecton
a workstation for hours and simulating the same e®ect in real-time, which is typically at
40-60framespersecondonasingleconsole. Thereareseveralrequirementsonsimulation
algorithms for interactive applications.
² Cheap to compute
A 3D computer game typically runs at a constant rate of 40-60 frames per second.
Among the total computational time allocated to a single frame, physically-based
simulation only gets a fraction due to the need of other tasks such as rendering
111
and/or character animation. Recently, physics simulation has been moved from the
CPU to the GPU (NVidia) or a dedicated Physics Processing Unit (PPU). The
faster the algorithm the more the number of objects being simulated in real time.
Thus, low complexity remains to be one of the most important criteria.
² Low memory consumption
While one may add as much memory as needed to simulate a desired e®ect in o®-
line computation, the memory on a game console or a single PC is limited and it
has to be shared among di®erent tasks that have comparable computational load.
Thus, memory e±ciency of the simulation algorithm has to be carefully considered,
too.
² Stability
Adaptivetime-steppingcanbeusedino®-linesimulationwhenthestabilityproblem
arises. In contrast, since a game runs at a ¯xed frame rate, the simulation method
has to be stable for a chosen time step. On the other hand, the value of external
forces and object movements can get very high. For example, an object may move
at a speed of 100mph in certain games. A real-time simulation method needs to
cope with these situations while remaining to be stable.
² Plausibility
By reducing the computational and spatial complexity so as to preserve stability,
the visual quality tends to degrade. What we require from a real-time simulation
method is visual plausibility (rather than scenes that are undistinguishable from
the real world).
112
In the context of °uid simulation, there are two commonly used approaches; namely,
Lagrangian and Eulerian. The position of parcels (or simulated quantities) of a °uid is
updatedintheLagrangianapproachasshowninFig. 5.1(a). Incontrast,quantitiesgoing
through¯xedpositionsareupdatedintheEulerianapproachasshowninFig.5.1(b). The
Eulerian approach is quite popular in computer graphics. Speci¯cally, °uid simulation is
achieved by animating the velocity in a grid as shown in Fig.5.1(c). The governing rule
for °uid simulation is the Navier-Stokes equations, which were described in Sec. 4.3.2.
There several drawbacks to the Eulerian velocity grid approach:
1. The amount of data is huge.
2. There is a lot of dissipation as opposed to the Lagrangian approach.
3. The simulation changes if the grid size is changed.
Duetothereal-timeconstraintanddrawbacksoftheEulerianmethodin°uidsimulation,
we choose the Lagrangian method over the Eulerian method in our simulation model.
(a) (b)
ijk
v
(c)
Figure 5.1: Two °uid simulation approaches: (a) the Lagrangian approach, (b) the Eu-
lerian approach, and (c) the 3-D Eulerian approach.
113
5.3 Dynamics of Convective Thermal
Thermal's dynamic behavior in a stable strati¯ed environment was explained in Sec.
2.2. In this section, we examine thermal's similarity analysis in an unstable atmospheric
environment.
In terms of thermal's dynamics, the velocity is determined by the buoyancy force.
The local buoyancy force on a °uid of unit mass can be described as
g
½
0
¡½
i
½
0
=gB; (5.1)
where g is the gravity acceleration parameter, and ½
i
and ½
0
are interior and exterior
densities of the thermal, respectively. The ratio B as de¯ned in (5.1) varies within the
thermal because of varying dilution. When the distribution of B is unknown, the mean
value B is used instead in the following description.
Under the assumption that the approximate horizontal section of a thermal is a
circle, we use r to denote the radius of the largest horizontal section. Then, the vertical
velocity of the frontal of a thermal follows the form
w
2
=C
2
gBr; (5.2)
where C is a Froude number, which is a dimensionless quantity denoting the ratio of
inertia and the buoyancy force. Note that the Froude number can be used to quantify
the resistance of an object moving through the °uid with di®erent sizes and speeds [46].
114
Let z be the height of the foremost part of the thermal. If the two velocities, dz=dt
and dr=dt, are proportional to each other, we can derive the following relationship of
thermal's radius r and height z:
z =nr; (5.3)
where n is a coe±cent indicating the widening rate. With the support of cumulus
observational study [50], we choose broadening coe±cient n from 3:1 » 3:7 in our sim-
ulation. The Froude number, C, whose physical meaning is the resistance of an object,
has a simple relationship with broadening coe±cient n as
C¼
1
2n
:
Since only one representative velocity exists, the motion is similar in all stages in the life
of a thermal and all thermals are similar.
ThevolumeV ofathermalthatpossessesaclearlyde¯nedoutlineintheexperiment
is given by
V =mr
3
; (5.4)
wheremisaparameterdeterminedbytheexperimentwhosetypicalvalueis3. Thetotal
buoyancy of the thermal is constant at all stage of a thermal. By using subscript 0 to
denote values at the moment of release, we have
gBmr
3
=gB
0
V
0
(5.5)
115
By replacing w = dz=dt, integrating Eq. 5.2 with the help of Eq. 5.5 and choosing the
origin of t appropriately, we obtain
kz
2
=t; (5.6)
where
k =m
1=2
=2nC(gB
0
V
0
)
1=2
: (5.7)
Typical thermal outlines during its lifetime are shown in Fig. 5.2. It is obvious that,
as the thermal moves up, its radius is proportional to its reached height.
Z
2
Time
Second
12
8
4
100 200 300
Figure 5.2: Successive outlines of thermals traced from photograph (top) and a graph of
z
2
against t (bottom).
116
5.4 Dynamics of Thermal Inner Motion
We see from Fig. 5.2 that, as the thermal arises, the front edge of its successive outlines
varies along time although its general shape is preserved. To observe the interior motion
of a thermal, Woodward performed an experiment by releasing a transparent thermal
containing white pellets with a diameter of 5mm [46]. The isopleths of horizontal and
vertical velocity components observed in the experiment are shown in Fig. 5.3. Such a
velocity distribution can be related to vorticity and vortex ring as discussed below.
Figure 5.3: The velocity distribution of particles inside a thermal, where the horizontal
(left) and the vertical (right) velocities are expressed as multiples of the vertical velocity
of the front of the thermal.
5.4.1 Vorticity and Vortex Ring
Vorticity ! is written as curl(u) or ! =r£u. Under the incompressibility assumption,
mass conservation can be expressed as div(u) = 0. In the 2D case, the Lagrangian
vorticity form of the Navier-Stokes equation for inviscid °uids is simply
d!
dt
=0. It means
that, oncecreated, vorticityneverdiesbutitissimplyadvectedalongthe¯eld. Handling
the 2-D case is simple since it only requires vortices placed at isolated particles.
117
We show a 2-D velocity plane which contains rotational °ow in Fig. 5.4, where the
vorticityisrepresentedasanisolatedparticleatthecenteroftherotation°ow. Extended
to the 3-D velocity space, the vorticity ¯eld can expressed as a curve connected by those
isolated vorticity particles(i.e. the orange curves in Fig. 5.4). 3-D vorticity consists of
several types: isolated voriticty particles and curves. When the curve is in form of a
closed circle, it is called a vortex ring. The vortex ring structure is shown in Fig. 5.5.
Figure 5.4: Illustration of a 2-D velocity ¯eld and the associated 3-D vorticity.
5.4.2 Axis-Symmetric Flow
By using the vorticity form, °ow's representation can be greatly simpli¯ed. Speci¯cally,
we can avoid the drawback of numerical dissipation resulted from the Eulerian °uid
simulation method. We are interested in how the mark density (the cloud droplet in our
case) move in the vorticity-represented °ow.
The dynamics of thermal interior motion is circularly symmetric. When the °uid
°ow is symmetrical with respect to the z-axis, it can be described by a function of r and
118
Figure 5.5: A vortex ring where the red arrow indicate the direction of axis (left) and a
two-dimensional slice of the velocity ¯eld induced by the vortex ring (right).
z [37]. If the azimuthal component of the velocity is zero, the vorticity ¯eld consists of a
set of circular vortex ¯laments or vortex rings. The governing equations can be written
as
@»
@t
+u
r
@»
@r
+u
z
@»
@z
=0; (5.8)
and
@
2
Ã
@z
2
+
@
2
Ã
@r
2
¡
1
r
@Ã
@r
=¡r!; (5.9)
whereu
r
andu
z
aretheradialandsteam-wise velocitycomponents, respectively, ! is the
vorticity (»´!=r), and Ã is the stream distribution. These quantities are related via
u
z
=
1
r
@Ã
@r
; (5.10)
and
u
r
=¡
1
r
@Ã
@z
: (5.11)
119
Eq. 5.8 is a transport equation for » derived from the vorticity transport equa-
tion. For inviscid °ows, it represents a pure convective equation similar to the vorticity
transport equation for 2-D °ows. Although » is no longer conserved for viscous °ow,
it underlies most of the existing vortex particle method for axis-symmetric °ows. Since
there is no dependence on the µ coordinate, the computational domain reduces to an
(r;z) plane as shown in Fig. 5.6. With this representation, vortex rings are described by
a single point which is their intersection with a selected meridional plane. For inviscid
°ows, the»-equation can be interpretedbythe Kevin and Helmholtz theorems. Since the
circulation is constant along any vortex ¯lament, any increase in the vortex ring radius
must be accompanied by an increased vorticity.
r
0
z
0
¡
z
r
r
0
r
0
r
0
z
0
z
0
z
0
¡
z
r
¡ ¡
z
r
z
r
z
r
Figure 5.6: Schematic sketch of a circular vortex ring.
An integral relation between the velocity ¯eld and the vorticity ¯eld can also be de-
rivedfromthevelocityatagivenpoint,(r;z), inducedbythevortexringwithcirculation
¡
0
located at point (r
0
;z
0
) [44]; namely,
u
z
(r;z;r
0
;z
0
)=
¡
0
2¼
r
³
(r+r
0
)
2
+(z¡z
0
)
2
´
"
K(k)¡
(z¡z
0
)
2
+r
2
¡r
2
0
(r¡r
0
)
2
+(z¡z
0
)
2
"(k)
#
120
u
r
(r;z;r
0
;z
0
)=
¡¡
0
(z¡z
0
)
2¼r
r
³
(r+r
0
)
2
+(z¡z
0
)
2
´
"
K(k)¡
(z¡z
0
)
2
+r
2
+r
2
0
(r¡r
0
)
2
+(z¡z
0
)
2
"(k)
#
;(5.12)
whereK and " are the ¯rst and the second order elliptic functions, respectively, and
k´
4rr
0
(r+r
0
)
2
+(z¡z
0
)
2
: (5.13)
Since the velocity at position (r;z), vortex ring location (r
0
;z
0
) and scaler circulation
¡
0
can be calculated, it makes the Lagrangian approach a natural choice to simulate the
°uid.
5.4.3 Relationship between Thermal and Vortex Ring
The thermal does not possess any circulation at the moment of release. Actually, circula-
tionisgeneratedgraduallyintheearlystageofitsexistence[46]. Intheearlystage,there
is a buoyancy force acting on a column of buoyant °uid and increasing the circulation. In
thelaterstage,muchofthecolumnisreplacedbytheexterior°uid. Whenthecirculation
around a thermal reaches its ¯nal con¯guration, it is proportional to wz.
Theconstantcharacteristicsofathermalincludeitstotalweight, gBV,andthe°uid
density, ½. It can be shown by dimensional analysis that circulation is equal to
¡
0
=const(gBV=½)
1=2
: (5.14)
121
Since thermal's interior circulation is bounded by its parameters, a thermal behaves like
a weak vortex ring at its mature stage. This relationship provides a link between the
thermal's movement as a whole and thermal's inner motion in our model.
5.5 Hierarchical Cloud Simulation
5.5.1 Simulation Framework and Parameters
Based on the previous discussion, we propose a hierarchical cloud simulation framework,
whose block-diagram is shown in Fig. 5.7. The user can set some initial variables for the
whole cloud simulation °exibly, for example, the initial thermal number and the approx-
imate location range. At every time step, thermal's position and velocity are updated by
thermal's similarity property. After that, we check whether the inner dynamics of the
thermal is needed according to thermal's position with respect to the camera position. If
the thermal is hidden by others from the current viewpoint, inner motion computation is
not needed. On the other hand, if a thermal is close to the camera and at the boundary
of the cloud, the time-variant shape caused by inner motion will a®ect the visual result.
Thus, inner particles inside the thermal has to be updated accordingly. By using the
hierarchical cloud simulation framework, we can perform the cloud simulation in real
time.
Oneadvantageofthesimilarityapproachisthatthemodelcatchesthebulkdynamics
using a characteristic length. The user can change these parameters °exibly to meet the
desired cloud simulation.
122
Figure 5.7: The hierarchical cloud simulation framework.
² In a stable strati¯ed environment, a user can specify thermal's initial parameter,
including heat °ux (F
0
), buoyancy frequency (N) and entrainment involving coef-
¯cient (alpha).
² In an unstable environment, a user can specify thermal's physical parameters such
as the entrainment involving coe±cient (alpha), the widening rate (n) and the
volume increase rate (m).
Besides, the number of particles inside a thermal can be speci¯ed. The number of
inner particles may change according to thermal's radius as it moves in the atmosphere.
123
The user can set the initial number and the maximum number of particles of a ther-
mal according to the relative position to the camera. Randomness is introduced at the
initialization stage, particles inside the thermal have di®erent initial positions, sizes and
densities.
5.5.2 Velocity-Vorticity Conversion
Weareinterestedinparticle'svelocityinsideathermal. Asliceofthevelocitydistribution
insideathermalisgiveninFig. 5.8, wheretheredcirclegivestherangeofathermal, the
blue point is the location where the vortex ring intersects with the slice and white lines
indicate the velocity ¯eld. Fig. 5.9 gives a rendering result for a sample thermal. The
vortex ring is a connected vorticity circle. Furthermore, we extend the axis-symmetric
°ow to the 3-D space by treating each vertical slice in the thermal separately.
Figure 5.8: The distribution of thermal velocities.
124
Figure 5.9: The rendering result of a sample thermal consisting of 150 particles.
5.5.3 Rotational Velocity
If we follow the standard Euler forward method to update inner particle's velocity and
position, their values could deviate from the original rotation path as shown in Fig.
5.10(a). When the time step for the simulation is too big, the simulation can become
unstable. Throughrotation,particle'spositionchangesontheorbit. Sincethetranslation
is actually due to rotation, we convert the distance along the tangential direction to the
arc length at every time step, thus getting the rotation angle. Then, we can move the
particle based on this rotation angle as illustrated in Fig. 5.10(b). With this method, the
time step in our simulation could be larger to meet the real-time simulation requirement.
(A) (B)
Figure 5.10: (a) The standard Euler forward method and (b) the rotational velocity
method.
125
5.5.4 Phase Transition
Another important physical process of cloud simulation is the phase transition; namely,
the water vapor transits to the visible water droplet. As thermal moves up, we need to
consider two processes. First, as the altitude increases, the environmental temperature
and the pressure decrease, the ability to hold the water vapor decreases. Second, during
thermal's lifetime, it continues to attract the environmental mass, which becomes part
of the thermal. The water droplet density increases due to these two processes. In our
simpli¯ed model, phase transition is achieved by increasing particle's transparent value
which is a function of particle's elevation and lifetime.
5.5.5 Thermal Interaction
An individual thermal moves by following the dynamics derived from the similarity ap-
proach. During the process, a thermal may collide into another, which enables the water
vapor to get accumulated and condensed. The latent heat releases as a result of the
condensation from the water vapor to the water droplet. The released latent heat can
be viewed as a new heat °ux source, which could result in a new thermal. Although
individual thermal's condensation can generate a new thermal theoretically, thermals'
collision leads to more condensation and more latent heat release. For simplicity, our
model does not allow a new thermal generated from the condensation of an isolated
individual thermal.
At each time step, thermal's relative distance is computed and compared with a
threshold to check whether a new thermal is born. New thermal's grow direction is
determined based on the combined local temperature gradient carried by each thermal.
126
(a) (b)
(c) (d)
Figure 5.11: The birth of a new thermal in the hierarchical cloud simulation.
Figs. 5.11 (a)-(d) show a new thermal born process. It is shown in particle form so that
the dynamics process can be observed clearly. Two thermals collide in Fig. 5.11 (a) and
a new thermal is born in Fig. 5.11 (c). Afterwards, the new thermal moves on its own as
shown in Fig. 5.11(d).
5.6 Simulation Results
Simulation results are given in this section to demonstrate the performance of the algo-
rithm presented in previous sections.
127
A whole cloud scene is shown in Fig.5.12, which is initialized with ¯ve thermals,
each of which follows the atmospheric movement and formation according to thermal's
similarity approach. The inner motion of the thermal is computed by the similar axis-
symmetric °ow of a weak vortex ring. As a thermal moves and grows furthermore, it
collides with others and creates a new thermal. The new thermal will inherit some of
parental thermal's property and move on its own. All these thermals have a similar dy-
namic property but di®er in the radius broadening rate. This process can be speci¯ed
by several coe±cients such as broadening coe±cient, n, and entrainment mass involving
coe±cient, ®. Since the inner motion introduces enough details, we can provide satisfac-
tory visual quality with a small number of thermals in our simulation. For example, the
total number of thermals in Fig. 5.12 (d) is only 13, where each thermal has less then
400 particles. The computational cost is very low, which enables real-time simulation.
In Fig. 5.13, we show the FPS statistic analysis of the scene presented in Fig. 5.12.
The left yaxis (shown in blue) records the value of FPS while the right yaxis (shown in
red)denotesthetotalnumberofparticlesforthewholescene. FPSdecreasesdramatically
as the total number of particles increases. With an acceptable number of thermals, we
can achieve a frame rate higher than 50, which meets the real-time requirement. This
simulation is performed by a PC with 2GHz CPU.
In Fig. 5.14, we show the FPS statistics for the same cloud simulation as given in
Fig. 5.12 by adaptively turning on the inner motion computation according to thermal's
relativepositionwithrespecttothecamera. Inthe¯gure,thebluelinecorrespondstothe
cloudsimulationalwayswiththeinnermotioncomputationwhiletheredlinecorresponds
tothecasewheretheinnermotioncomputationisturnedo®ifthecorrespondingthermal
128
(a) (b)
(c) (d)
Figure 5.12: Cloud formation based on Hierarchical structure
is hidden behind other thermals from the viewpoint. Since the bulk computation of this
simulationistheinnerdynamicsofallthermals,weobservethattheFPSofthesimulation
improves as the number of thermals increases by enabling the turn-o® feature.
In Fig. 5.14, a cauli°ower cloud formation and evolution with a di®erent initial
thermal setting is shown. The illumination and rendering method we used here is similar
to that in Chap.3 Sec.3.5, the scene is under a strong sun light at noon, thus the top
cloud which is closer to sun has very bright color compare to the cloud bottom.
Fig. 5.16 and Fig. 5.17 shows another cauli°ower cloud scene with di®erent initial
settings under two di®erent lighting conditions. In Fig. 5.16, the sun is locating behind
129
50
70
90
110
130
150
Thermal Number
FPS
FPS Statistics of Hierachical Cloud Simulation
7 8 9 10 11 12 13 14 15 16
0
400
800
1200
1600
2000
Particle Number
Figure 5.13: Statistic FPS for hierarchical cloud simulation.
thecloud, thusthesilverliningfeaturesformultiplescatteringisobvioustoobserve. Fig.
5.17 gives the cloud image at sunset. With the multiple scattering illumination apply on
the cloud scene, the simulation result looks more realistic.
5.7 Conclusion
A novel real-time hierarchical cloud simulation was proposed based the similarity ap-
proach. Thecoarselevelcloudsimulationwasderivedbythermal'smovementasawhole.
Then, according to the relative position between the thermal and the camera, the corre-
sponding interior dynamics motion can be performed to generate a lot of cloud detail at
the cloud boundaries. By apply the axis-symmetric °ow into the thermal's inner motion,
we can achieve simple while visually convincing real-time simulation results.
130
7 8 9 10 11 12 13 14 15 16
50
60
70
80
90
100
110
Thermal Number
FPS
FPS Statistics Comparison between Two Simulations with/o Inner Motion Turning off

simulation without turn off
simulation with turning off
Figure5.14: StatisticFPScomparisonforhierarchicalcloudsimulationwithandwithout
inner motion turned o®.
131
(a) (b)
(c) (d)
(e) (f)
Figure 5.15: Cauli°ower cloud formation using hierarchical cloud simulation.
132
Figure 5.16: A scene with multiple cumulus clouds.
Figure 5.17: A scene with multiple cumulus clouds at sunset.
133
Chapter 6
Conclusion and Future Work
6.1 Conclusion
In this thesis, the physics of cumulus cloud formation was reviewed and used as the basis
for cloud modeling and simulation. We also proposed to leverage similarity approach to
describe the characteristics of simple °uid °ow without rigorously solving the governing
PDEs. Based on the similarity solution of turbulent thermals, plumes and vortex rings
in stable strati¯ed °uids or unstable environments, we proposed three methods for cloud
simulationmodeling. Themainresults,aspresentedinChapter3,Chapter4andChapter
5, are summarized below.
In Chapter 3, a two-stage cloud simulation method was proposed to model the cloud
formationprocess. Atthe¯rststage,insteadofsolvingPDEsforthevelocitypro¯le,cloud
formation was obtained using a lower level building block; namely, particle simulation.
Sincetheindividualparticlepathisirregularintherealworld,onlythestatisticalaverage
in time and/or space provides a systematic velocity pro¯le. New attributes of particles
were incorporated to meet the modeling requirement, e.g., the condensation time and
134
densityofparticles. Inourproposedmodelingmethod,particlesactasanimportantunit,
i.e.,theairparcelinatmosphereandtheyareusedtomodelanairparcel'swholelifetime
before it condenses into a visible cloud droplet. For particle motion, the characteristics
of the convective process is captured by the thermal similarity solution consistent with
the governing Navier-Stoke equation. At the second stage, we consider the process of
beingcondensedintothevisiblewaterdroplet. Thevoxelisusedtostorethewatervapor
density and the water droplet density. Due to the release of latent heat, we choose to
have a small scale convective process. However, we only simulate this sub-scale process
at the cloud surface area for computational e±ciency. Finally, the cloud droplet density
is rendered using a two-pass rendering method and multiple scattering to result in the
high albedo property of the cloud. To conclude, without solving the governing equations
rigorously, we are still able to capture the naturalism of cumulus cloud formation at a
low computational cost.
InChapter4,Adecouplingcloudsimulationmethodthatdividesthreespatialdimen-
sion into independent horizontal and vertical dimensions was proposed. This approach
is motivated by the atmospheric observation that cumulus clouds develop by a series of
rising plumes. With the assumption that the radial pro¯les at all heights are similar, we
choose several key frame layers only for actual 2D calculation, interpolate other layers
from them and add some turbulence. Along the horizontal dimension, the movement is
driven by heat di®erence. This virtual temperature di®erence results in a di®erent water
vapordistribution. Tomodelthecloudformationfromaplume,weassumethatthewater
phase change occurs at its late development stage, where the plume spreading between
levels of buoyancy vanishes and the vertical velocity is equal to zero. At such a stage,
135
thewatervaporwillbesu±cientlyaccumulatedtoexperiencethewaterphasetransition.
The 2D cloud layer is then simulated by the Navier-Stoke equations. Besides the mo-
mentum equation coming from the virtual temperature di®erence, one more force term is
introducedinourmodel, i.e. theCoriolisforce. Thisforcetermiscomparablewithother
horizontal acceleration factor in atmosphere. The link between layers is their location at
the turbulent plume. Each key frame layer gets the water vapor input at each time step,
where the input horizontal water content distribution is assumed to be either \Gaussian"
or \top-hat". A plateau-like kernel is adopted to generate the water vapor distribution
°exibly. Finally, water droplets are rendered using a 3D texture rendering technique. By
separating the horizontal and vertical couplings in the traditional 3D simulation, we can
reduce the computational cost yet capture the essential characteristics of the cumulus
cloud. Further computational speed-up can be achieved by parallel computing and GPU
processing, which are research topics in the future.
In Chapter 5, we introduced a novel hierarchical cloud simulation method. The gen-
eralshapeandmovementofcloudsinthesimulationmodelaregovernedbythesimilarity
of thermals. Users can initialize the thermal number and the rough initial position at
every time step, and thermals update their locations by a similarity approach. In the
meanwhile, their relative position to the camera is calculated. If they are not hidden by
other thermals or they are in the boundary of a cloud, we turn on thermal's inner motion
which behave as a weak vortex ring. We simulate thermal's interior dynamics using the
axis-symmetric °ow, which provides a simple scheme to determine the velocity at any ar-
bitrary location inside the °ow. With the axis-symmetric °ow, we adopt the Lagrangian
approach for thermal's interior motion, where a bunch of water droplet particles inside
136
the thermal are updated at each time-step, generating the desired cloud detail structure.
Furthermore, we consider the phase transition with a simpli¯ed strategy. With the en-
trainmentmasswhosecoe±cientsarechosenbyusers,wecanincreasethewaterdroplet's
density, thus its size and transparency for rendering. A method for cloud formation due
to the latent heat release from the condensation is also considered. Finally, we adopt a
two-pass cloud illumination model which is similar to that in Chapter 3 to render the
cloud. Due to its low complexity, the hierarchical cloud simulation can be conducted in
real time for gaming and °ight simulation applications.
6.2 Future Work
To make the current research work more complete and suitable for the °ight simulation
in a gaming environment, it is worthwhile to get e®ective modeling of more di®erent
cloud types, further speed-up in rendering, and more functionalities such as °y-in/out
interactivity with the cloud. Some possible future research topics are stated below.
² Hardware Application
Simpli¯edcloudsimulationmodelshavebeenproposedtoreducethecomputational
cost. Further computational complexity reduction can be achieved by considering
hardware speed-up. With parallel/grid computing and GPU implementation, the
simulation and the rendering speed should improve a lot. This will make real-time
high quality cloud simulation possible.
² Stratus and Cirrus Cloud Modeling
137
As described in Chapter 2, almost all cloud types can be viewed as the result of a
convective process. However, only the cumulus cloud bears the extreme property of
theconvectivecloudwhileothersaredominatedbyotherphysicalprocessesaswell.
Thus, we need to ¯nd the approximation to other physical processes, and extend
our current convective cloud model so as to generate di®erent cloud types.
² Cloud Illumination
Currently, the computation cost is still quite high for accurate illumination. It is
interesting to ¯nd e±cient yet accurate cloud illumination model for our simulation
model.
² Cloud Interactivity
Flight simulation is one of the most important ¯elds that need realistic cloud sim-
ulation and rendering. For this application, new functionalities such as °y-in and
°y-outinteractivitywiththecloudisingreatdemand. E±cientrenderingisneeded
for this type of application.
138
Bibliography
[1] D.G.Andrews,An Introduction to Atomospheric Physics,1sted. Cambridge,2000.
[2] A. B. Andrzej Trembilski, \Transparency for polygon based cloud rendering," Sym-
posium on Applied Computing, Feb 2002.
[3] G.K.Batchelor,\Heatconvectionandbuoyante®ectsin°uid,"inQuarterlyJournal
of the Royal Meteorological Society, vol. 80, 1954.
[4] U. Behrens and R. Ratering, \Adding shadows to a texture-based volume renderer,"
in IEEE Symposium on Volume Visualization, 1998, pp. 39{46.
[5] N. Bhate and A. Tokuta, \Photorealistic volume rendering of media with directional
scattering," in In Proceedings of Eurographics, 1992, pp. 227{245.
[6] L. S. B. Blasi, P. and C. Schlick, \Arendering algorithm for discrete volume density
objects," in In Proceedings of Eurographics, 1993, pp. 201{210.
[7] J. Blinn, \Light re°ection functions for simulaiton of clouds and dusty surfaces," in
SIGGRAPH'82, 1982, pp. 21{29.
[8] J. F. Blinn, \A generalization of algebraic surface drawing," ACM Trans. Graph.,
vol. 1, no. 3, pp. 235{256, 1982.
[9] B. Cabral, N. Cam, and J. Foran, \Accelerated volume rendering and tomographic
reconstruction using texture mapping hardware," in VVS '94: Proceedings of the
1994 symposium on Volume visualization. New York, NY, USA: ACM Press, 1994,
pp. 91{98.
[10] A.ChorinandJ.Marsden, A Mathematical Introduction to Fluid Mechanics, 3rded.
Springer, New York, 1993.
[11] M. F. Cohen and J. R. Wallace, Radiosity and Realistic Image Synthe-
sis. Boston, MA: Academic Press Professional, 1993. [Online]. Available:
citeseer.ist.psu.edu/cohen93radiosity.html
[12] D. P. K. P. D.S Ebert, F.K. Musgrave and S. Worley, Textureing & Modeling: A
Procedural Approach, 3rd ed. Morgan Kaufman, 2002.
[13] D. Ebert and E. Bedwell, \Implicit modeling with procedural techniques," in Pro-
ceedings Implicit Surfaces, 1998.
139
[14] D. S. Ebert, \Volumetric modeling with implicit functions: A cloud is born," in
ACM SIGGRAPH, 1997.
[15] H. elias, \Cloud," http://freespace.virgin.net/hugo.elias/models/m clouds.htm.
[16] P. Elinas and W. StÄ urzlinger, \Real-time rendering of 3d clouds," J. Graph. Tools,
vol. 5, no. 4, pp. 33{45, 2000.
[17] K. A. Emanuel, Atmospheric Convection, 1st ed. Oxford University Press, 1994.
[18] R. Fedkiw, J. Stam, and H. W. Jensen, \Visual simulation of smoke," in ACM
SIGGRAPH, 2001, pp. 15{22.
[19] N. Foster and D. Metaxas, \Realistic animation of liquids," in Graphical Models and
Image Processing, no. 58, 1996.
[20] G. Y. Gardner, \Visual simulation of clouds," in ACM SIGGRAPH, 1985.
[21] M.J.Harris, \Implementationofacmlboilingsimulationusinggraphicshardware,"
University of North Carolina," Technical Report TR02-015.
[22] M. J. Harris and A. Lastra, \Real-time cloud rendering," in EG 2001 Proceedings,
A. Chalmers and T.-M. Rhyne, Eds. Blackwell Publishing, 2001, vol. 20(3), pp.
76{84.
[23] L. Henyey and J. Greenstein, \Di®use radiation in the galaxy," The Astrophysical
Journal, vol. 90, pp. 70{83, 1941.
[24] D. S. E. Joshua Schopik, Joseph Simons and C. Hansen, \A Real-Time cloud mod-
eling, rendering, and animations system," in Eurographics/SIGGRAPH Symposium
on Computer Animation, 2003.
[25] J. T. Kajiya and B. P. V. Herzen, \Ray tracing volume densities," in SIGGRAPH
'84: Proceedings of the 11th annual conference on Computer graphics and interactive
techniques. New York, NY, USA: ACM Press, 1984, pp. 165{174.
[26] e.Kaneko,K.,\Theoryandapplicationsofcoupledmaplattices,"NonlinearScience:
theory and applications. Wiley, pp. 169{189, 1993.
[27] J. Kniss, S. Premoze, C. Hansen, and D. Ebert, \Interactive translucent volume
rendering and procedural modeling," in VIS '02: Proceedings of the conference on
Visualization '02. Washington, DC, USA: IEEE Computer Society, 2002.
[28] M. Levoy, \Display of surfaces from volume data," IEEE Comput. Graph. Appl.,
vol. 8, no. 3, pp. 29{37, 1988.
[29] J. P. Lewis, \Algorithms for solid noise synthesis," in ACM SIGGRAPH, 1989, pp.
263{270.
[30] T. S. A. L. Mark J. Harris, William V. Baxter III, \Simulation of cloud dynamics
on graphics hardware," in EUROGRAPHICS, 2003.
140
[31] B. Mason, The Physics of Clouds. Oxford University Press, London, 1957.
[32] N. L. Max, \E±cient Light Propagation for Multiple Anisotropic Volume Scatter-
ing," in Fifth Eurographics Workshop on Rendering, 1994, pp. 87{104.
[33] G. Mie, \Bietage zur optik truver medien speziell kolloidaler metallosungen," An-
nallen der Physik, vol. 25, no. 3, p. 377, 1908.
[34] B. R. Morton, \Buoyant plumes in a moist atmosphere," no. 2, 1957.
[35] K. Nagel and E. Raschke, \Self-organizing criticality in cloud formation?" Physica
A Statistical Mechanics and its Applications, vol. 182, pp. 519{531, Apr. 1992.
[36] F. Neyret, \Qualitative simulation of convective cloud formation and evolution," in
In Eurographics Workshop on Animation and Simulation, 1997, pp. 113{124.
[37] A. Ogawa, Vortex Flow. CRC Press, 1993.
[38] D.Overby, Z.Melek, andJ.Keyser, \Interactivephysically-basedcloudsimulation,"
in PG '02: Proceedings of the 10th Paci¯c Conference on Computer Graphics and
Applications. IEEE Computer Society, 2002, p. 469.
[39] Perlin, \An image synthesizer," in ACM SIGGRAPH, 1985, pp. 287{296.
[40] K. Perlin, \Making noise."
[41] W. T. Reeves, \Particle systems technique for modeling a class of fuzzy objects,"
ACM Trans. Graph., vol. 2, no. 2, pp. 91{108, 1983.
[42] W. T. Reeves and R. Blau, \Approximate and probabilistic algorithms for shading
andrenderingstructuredparticlesystems,"inACMSIGGRAPH,1985,pp.313{322.
[43] P. G. Renato Pajarola, Miguel Sainz, \Object-space blending and splatting of
points," UCI-ICS Technical Report," Technical Report, 2003.
[44] E. Rivoalen and S. Huberson, \Numerical simulation of axisymmetric viscous °ows
by means of a particle method," Journal of Computational Physics, 1998.
[45] R. Rogers and M. Yau, A Short Course in CLoud Physics, 3rd ed. International
Series in Natural Philosophy. Butterworth Heinemann, Burlington, MA, 1989.
[46] R.S.Scorer,\Experimentsonconvectionofisolatedmassesofbuoyant°uid," Journal
of Fluid Mechanics, vol. 2, pp. 583{594, 1957.
[47] H. E. Rushmeier and K. E. Torrance, \The zonal method for calculating light inten-
sitiesinthepresenceofaparticipatingmedium," in SIGGRAPH '87: Proceedings of
the 14th annual conference on Computer graphics and interactive techniques. New
York, NY, USA: ACM Press, 1987, pp. 293{302.
[48] Y. D. T. N. Ryo Miyazaki, Satoru Yoshida, \A method for modeling clouds based
on atmospheric °uid dynamics," in Paci¯c Graphics, 2001.
141
[49] C. S., Radiative Transfer. Dover, New York, 1961.
[50] P. M. Saunders, \An observational study of cumulus," Journal of the Atmospheric
Sciences, pp. 451{467, Aug 1961.
[51] W. Schmidt, \Turbulent propagation of a stream of heated air," in Z. angew. Math.
Mech, no. 21, 1941.
[52] R. Scorer, \Experiments on convection of isolated masses of buoyant °uid," Journal
Fluid Mech., vol. 2, 1957.
[53] R. K. Smith, The Physcis and Parameteriazation of Moist Atmospheric Convection.
Philadelphia, PA, USA: Society for Industrial and Applied Mathematics, 1997.
[54] R. Srivastava, \A study of the e®ect of precipitatoin on cumulus dynamics," Journal
of Atmospheric Science, vol. 24, pp. 36{45, 1967.
[55] J. Stam, \Multiple Scattering as a Di®usion Process," in Rendering Techniques '95
(Proceedings of the Sixth Eurographics Workshop on Rendering), P. M. Hanrahan
and W. Purgathofer, Eds. New York, NY: Springer-Verlag, 1995, pp. 41{50.
[Online]. Available: citeseer.ist.psu.edu/stam95multiple.html
[56] J. Stame, \Stable °uids," in ACM SIGGRAPH, Aug. 1999.
[57] M. Stamminger and G. Drettakis, \Interactive sampling and rendering for
complex and procedural geometry," in Rendering Techniques 2001 (Proceedings
of the Eurographics Workshop on Rendering 01), K. Myskowski and S. Gortler,
Eds., Eurographics. Springer Verlag, 2001. [Online]. Available: http://www-
sop.inria.fr/reves/publications/data/2001/SD01
[58] J. Steiner, \A three-dimensional model of cumulus cloud development," Journal of
Atmospheric Science, vol. 30, pp. 414{435, 1973.
[59] T. Takeda, \Numerical simulation of a precipitation convective cloud: the formation
of a "long-lasting" cloud," Journal of Atmospheric Science, vol. 28, pp. 350{376,
1971.
[60] E. N. Tomoyuki Nishita, Yoshinori Dbashi, \Display of clouds taking into account
multiple anisotropic scattering and sky light," in ACM SIGGRAPH, July 1996.
[61] J. S. Turner, \The starting plume in neutral surroundings," J. Fluid Mech., no. 13,
1962.
[62] E. Vane, \Cloud rendering using 3d textures," University of Waterloo," Technical
Report, 2004.
[63] E. Veach, \Robust monte carlo methods for light transport simulation," Ph.D. dis-
sertation, Stanford University, 1997.
[64] J. H. W., \Global illumination using photon maps," in In Proceedings of the 7th
Eurographics Workshop on Rendering, 1996, pp. 21{30.
142
[65] J. H. W. and P. Christensen, \E±cient simulation of light transport in scenes with
participating media using photon maps," in ACM SIGGRAPH, 1998, pp. 311{320.
[66] B. Wang, \E±cient and realistic cumulus cloud simulation based on similarity ap-
proach," in ISVC: International Symposium on Visual Computing, 2007.
[67] N. Wang, \Realistic and fast cloud rendering," in ACM SIGGRAPH, 2003.
[68] L. Westover, \Footprint evaluation for volume rendering," in SIGGRAPH '90: Pro-
ceedings of the 17th annual conference on Computer graphics and interactive tech-
niques. New York, NY, USA: ACM Press, 1990, pp. 367{376.
[69] R.A.A.WilliamR.Cotton, Storm and Cloud Dynamics. NewYork: McGraw-Hill,
1982.
[70] O. Wilson, A. VanGelder, and J. Wilhelms, \Direct volume rendering
via 3d textures, Tech. Rep. UCSC-CRL-94-19, 1994. [Online]. Available:
citeseer.ist.psu.edu/wilson94direct.html
[71] T. Yanagita and K. Kaneko, \Modeling and characterization of cloud dynamics," in
Physical Review Letters, vol. 78, no. 22, 1997.
[72] T. O. T. N. Yoshinori Dbashi, Kazufumi Kaneda, \A simple, e±cient method for
realistic animation of clouds," in ACM SIGGRAPH, July 2000.
143

Abstract (if available)

Abstract Realistic cloud simulation and rending find applications in many 3D applications such as high-quality 3D film production and realistic 3D flight simulation. However, the infinite variation in the cloud shape and appearance make cloud simulation and rendering a challenge in termsof complexity and realisticity. The major focus of this work is on efficient and realistic cloud simulation. We adopt a similarity approach to describe the characteristics of simple fluid flow without rigorously solving the governing PDEs. Based on the similarity solution of turbulent thermals, plumes and vortex rings in stable stratified fluids or unstable environments, we propose three methods for cloud simulation and modeling.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Advanced liquid simulation techniques for computer graphics applications

PDF

City-scale aerial LiDAR point cloud visualization

PDF

3D deep learning for perception and modeling

PDF

Object detection and recognition from 3D point clouds

PDF

Interactive rapid part-based 3d modeling from a single image and its applications

PDF

Variational techniques for cardiac image analysis: algorithms and applications

PDF

Three-dimensional kinetic simulations of whistler turbulence in solar wind on parallel supercomputers

PDF

Focus mismatch compensation and complexity reduction techniques for multiview video coding

PDF

3D object detection in industrial site point clouds

PDF

Semantic modeling of outdoor scenes for the creation of virtual environments and simulations

PDF

Green learning for 3D point cloud data processing

PDF

Kinetic studies of collisionless mesothermal plasma flow dynamics

PDF

Explainable and green solutions to point cloud classification and segmentation

PDF

Point-based representations for 3D perception and reconstruction

PDF

Plant substructuring and real-time simulation using model reduction

PDF

Efficient coding techniques for high definition video

PDF

Advanced techniques for human action classification and text localization

PDF

Hybrid methods for robust image matching and its application in augmented reality

PDF

Labeling cost reduction techniques for deep learning: methodologies and applications

PDF

Grid-based Vlasov method for kinetic plasma simulations

Asset Metadata

Creator Wang, Bei (author)

Core Title Techniques for efficient cloud modeling, simulation and rendering

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering (Multimedia and Creative Technology)

Publication Date 02/01/2010

Defense Date 12/04/2007

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag Cloud,Cumulus,OAI-PMH Harvest,similarity approach,simulation

Language English

Advisor Kuo, C. C. Jay (committee chair), Nayak, Krishna (committee member), Neumann, Ulrich (committee member)

Creator Email beiwang@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m1003

Unique identifier UC1230327

Identifier etd-Wang-20080201 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-40172 (legacy record id),usctheses-m1003 (legacy record id)

Legacy Identifier etd-Wang-20080201.pdf

Dmrecord 40172

Document Type Dissertation

Rights Wang, Bei

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu