UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The computer modelling of fallen snow Fearing, Paul 2000

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_2000-565408.pdf [ 53.75MB ]
JSON: 831-1.0051217.json
JSON-LD: 831-1.0051217-ld.json
RDF/XML (Pretty): 831-1.0051217-rdf.xml
RDF/JSON: 831-1.0051217-rdf.json
Turtle: 831-1.0051217-turtle.txt
N-Triples: 831-1.0051217-rdf-ntriples.txt
Original Record: 831-1.0051217-source.json
Full Text

Full Text

THE C O M P U T E R M O D E L L I N G OF F A L L E N SNOW By Paul Fearing B . Sc. (Computer Science), Queen's University M . Sc. (Computer Science), University of British Columbia A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F THE REQUIREMENTS FOR T H E DEGREE OF DOCTOR OF PHILOSOPHY  in T H E F A C U L T Y O F G R A D U A T E STUDIES COMPUTER SCIENCE  We accept this thesis as conforming to the required standard  T H E  U N I V E R S I T Y  OF J U L Y  ©  B R I T I S H  2000  PAUL FEARING,  2000  C O L U M B I A  I n p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of the r e q u i r e m e n t s f o r an advanced degree a t the U n i v e r s i t y of B r i t i s h Columbia, I agree t h a t the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s t u d y . I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g of t h i s t h e s i s f o r s c h o l a r l y purposes may be g r a n t e d by the head of my department or by h i s or her r e p r e s e n t a t i v e s . I t i s u n d e r s t o o d t h a t c o p y i n g or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n .  Department of  Science  The U n i v e r s i t y of B r i t i s h Columbia Vancouver, Canada  The Computer Modelling of Fallen S n o w  Paul Fearing k  Abstract  One of nature's greatest beauties is the way fresh snow covers the world in a perfect blanket of crystalline white. Snow replaces sharp angles with gentle curves, and clings to surfaces to form ghostly silhouettes. It is said the Inuit have 50 different words for snow, yet even they can be left speechless, as snow is one of the most complex natural materials in existence. This research presents a new model of visual snow accumulation for computer graphics. We are primarily concerned with creating and simulating fallen snow (not falling snow - an important distinction), with our ultimate goal to produce view-independent, static, 3D snow surface models that can be used in artistic and scientific visualisation, film, and advertising. Our contribution is divided into two major components: snow placement and snow stability. Each are essential for modelling the appearance of a thick layer of snowfall on the ground.  Snow placement requires us to determine how much snow falls upon the scene, and where it accumulates. We simulate this with an adaptive particle/surface hybrid system that allows for such phenomena as flake flutter, flake dusting and wind-blown snow. We compute snow accumulation by shooting particles upwards towards the sky, giving each source surface independent control over its own sampling density, accuracy and computation time.  Importance ordering minimises sampling effort while max-  imising visual information, generating smoothly improving global results that can be interrupted at any point. Once snow lands on the ground, our stability model moves material away from physically unstable areas in a series of small, simultaneous avalanches. We use a simple local stability test that handles very steep surfaces, obstacles, edges, and snow transit due to wind. Our stability algorithm is flexible enough to simulate other materials, such as flour, sand, and flowing water. We show physical plausibility by comparing various aspects of our approach with real snow images. As proof that our algorithm is flexible and usable, we provide several examples of snow on complex models containing hundreds of thousands of polygons. The completed 3D snow surface model can be easily imported into commercial modelling and rendering software, allowing users to convert existing animations to a brand new season.  ii  Table of Contents  Abstract  n  List of Tables  xi  List of Figures  xii  Acknowledgements  xix  Copyrights A n d Credits  •  Additional Material  xxi xxii  1  Introduction  1  2  Snow  4  2.1  Snow in the Air  5  2.1.1  Precipitation  5  2.1.2  Formation of Ice Crystals  7  2.1.3  Ice Crystal Growth  7  Crystal Growth Due to Riming  9  Crystal Growth Due to Collision  9  2.2  Snow on the Ground 2.2.1  9  Initial Conditions of Fallen Snow  11  Density  11  Grain Shape . .  12  Grain Size  12  Liquid Water Content  13  Snow Strength  14  Impurities  14  in  2.2.2  2.3  2.4  2.5  2.6  3  Snow Hardness  14  Snow Temperature  15  Initial Rounding and Settling  16  16  Metamorphism due to Implicit Crystal Shape  2.2.3  Sintering  18  2.2.4  Crystal Growth and Recrystalization  18  2.2.5  Wet Snow Metamorphism  20  2.2.6  Wind and Snow  21  2.2.7  Additional Snow Forms  23  Sun Crust  23  Perforated Crust  23  Sun Cups  24  Snow-like Effects  24  2.3.1  Rime  24  2.3.2  Surface Hoar  25  Snow Strength  26  2.4.1  Compressive Stresses  26  2.4.2  Shear Stresses  27  2.4.3  Tensile Stresses  27  2.4.4  Angle of Repose  27  Avalanches  29  2.5.1  Loose Snow Avalanches  29  2.5.2  Slab Avalanches  30  Illumination Properties  31  2.6.1  Reflectance  31  2.6.2  Extinction  33  Previous Work  35  3.1  Geospecific Rendering of Alpine Terrain  36  3.2  Rendering Snow With Metaballs  38  3.2.1  41  Other Snow Rendering Work  iv  3.3  3.4  3.5 4  Particle Systems 3.3.1  Particle Animation and Rendering Using Parallel Computation  41  3.3.2  Stochastic Motion  42  3.3.3  Flow and Changes in Appearance  43  3.3.4  Dust Accumulation  44  3.3.5  Using Particle Systems  45  Snow on the Ground  4.2  4.3  45  3.4.1  Soil  3.4.2  Animating Sand, Mud and Snow  48  3.4.3  Eroded Terrains  48  3.4.4  Fluid Dynamics  50  3.4.5  Deformable Surfaces  50  3.4.6  Oriented Surface Particles  51  3.4.7  Multi-scale Granular Model  52  • • 45  Ad Hoc Methods  Assumptions,  4.1  11  Goals and  54 Problems  55  Assumptions and Restrictions  56  4.1.1  Scale of Snowfall  56  4.1.2  Snowfall Locality  56  4.1.3  Ground Evolution of Snow  57  4.1.4  Drawing Avalanches  57  4.1.5  Use of Commercial Software  59  4.1.6  Static World  60  4.1.7  Object Restrictions  61  Goals of the Snow Model  61  4.2.1  Snow Location  62  4.2.2  Snow Stability  67  4.2.3  Snow Surface Appearance  67  4.2.4  Snow Under the Influence of Wind  68  Overview of Problems  70  v  5  Accumulation of Snow on Occluded Surfaces  70  4.3.2  Surface Representation  72  Level of Detail Representations  72  Proportional Accumulation  72  Spatially Adaptive Representations  73  Temporally Adaptive Representations  74  Surface Smoothness  74  4.3.3  Snow Stability  74  4.3.4  Edge-Crossing Snow  77  Approaches  79  5.1  Volume-Based Paradigm  80  5.1.1  Advantages  80  5.1.2  Disadvantages  81  5.2  5.3 6  4.3.1  Surface-Based Paradigm  83  5.2.1  Advantages  83  5.2.2  Disadvantages  84  Hybrid Paradigm  86  A Computer Model of Fallen Snow  89  6.1  Overview of the Solution  90  6.2  Entities  9!  6.2.1  World  92  6.2.2  Model  92  6.2.3  Faces  92  6.2.4  Launch Site  93  6.2.5  Subdivision Area (or Launch Area)  93  6.2.6  Edge Group (or Subdivision Mesh)  95  6.2.7  Drops  95  6.2.8  Flake  96  6.2.9  Flake Group  96  6.2.10 Snow Planes  96 vi  6.3  6.4  7  6.2.11 Avalanche  97  6.2.12 Avalanche Flake  97  Overview of Phases  98  6.3.1  Storms and Storming  98  6.3.2  Subdivision Approximation Phase  98  6.3.3  Regular Storm Phase  99  6.3.4  Uncompleted Phases  100  Acquiring the Surface Model  101  6.4.1  101  Exporting the Surface Model From Alias|Wavefront  Snow Accumulation  103  7.1  Overview of Snow Accumulation  103  7.1.1  104  7.2  7.3  Computing the Snowfall Accumulation Pattern  Importance Ordering  104  7.2.1  Closeness to Camera  109  Overview of Launch Site Meshing  Ill  7.3.1  Edge Groups  Ill  7.3.2  Constrained Voronoi Diagrams as a Meshing Representation  112  7.3.3  Using Voronoi Diagrams To Determine Neighbourhoods  116  117  7.3.4  Problems With Our Meshing Representation  Increasing and Decreasing Mesh Resolution  117  7.4  Initial Particle Distribution  120  7.5  Simulating Snowflake Motion  123  7.5.1  Overview of the Algorithm  123  7.5.2  Depressions and Snow Boundary Control  125  7.5.3  Swirl Length  127  7.5.4  Choosing Parameters for Flake Motion  128  7.5.5  Efficient Intersection Bucketing  128  7.5.6  Wind and Flake Location  130  7.6  Locating Particles in the Sky  131  7.6.1  132  Allocating Mass to Snow Surfaces  vii  7.6.2 7.7  8  9  Writing in the Sky  133  Texture Generation  134  7.7.1  Flake Dusting  134  134  Computing Flake Dusting Location and Densities  7.7.2  Texture Replacement for Distant Surfaces  135  7.7.3  Problems with Texture Generation  138  Texture Interpenetration  138  Snow Transit and Flake Dusting Density  140  Snow Stability  141  8.1  Overview of the Algorithm  141  8.2  Maintaining a List of Unstable Samples  142  8.3  Angle of Repose  143  8.4  Stability Test  144  8.5  Moving Snow Over Edges  147  8.6  Stability Termination Criteria  149  8.7  Rebucketing Surfaces  150  8.8  Wind  151  Snow Final Phases  155  9.1  Surface Resampling  155  9.2  Implicit Functions  156  9.2.1  Advantages and Disadvantages of Implicit Functions .  162  9.2.2  Cornices  163  9.2.3  Vertex Reduction  163  9.3  Surface Rendering  164  9.4  Extras  165  9.4.1  Creating Snow Layers  165  9.4.2  Tracks and Surface Patterns  169  9.4.3  Water and Rain  170  9.4.4  Multiple Output Formats  171  vni  10 Results  1  10.1 Validation of the Visual Appearance of Snow  7  2  173  10.11 Validation of Visual Properties  173  10.1.2 Validation of Desirable Algorithm Properties  174  10.2 Timing and Complexity  174  10.3 Snow Images Generated Using the Complete Model  180  10.4 Oil Painting Effects  200  11 Future Work  202  11.1 Snow Properties  202  11.1.1 Precipitation 11.1.2 Density  203 :  203  11.1.3 Grain Shape  204  11.1.4 Grain Size  204  11.1.5 Liquid Water Content  204  11.1.6 Snow Temperature and Snow Metamorphism  204  11.1.7 Surface Smoothness  205  11.1.8 Stability Jamming  205  11.1.9 Stability Non-Termination  206  11.2 Issues with the Base Model  206  11.2.1 Inconsistent Normals  206  11.2.2 Polygons vs. Curved Surfaces  207  11.2.3 Weatherproofing Models  207  11.3 Rendering and Timing  208  11.4 Moving Objects  208  11.5 Nature Rubs Our Nose In It  209  12 Conclusions  211  12.1 Conclusions  211  Bibliography  213  Appendices  218 ix  A Table of Variables  218  B  219  List of Parameters  Glossary  220  The Homage  230  x  List of Tables  * 2.1  Primary Physical Characteristics of Deposited Snow  11  2.2  Typical Snow Layer Densities  12  2.3  Grain Size Classification  12  2.4  Liquid Water Content  13  2.5  Experimental Angles of Repose  9.6  Typical Mesh Reduction  . 30 164  10.7 Recap of Snow Appearance Goals  174  10.8 Recap of Desirable Algorithm Properties  175  A. 9 Description of Variables  218  B. 10 List of Parameters  219  xi  List of Figures  0.1  The author, buried in work  xx  2.1  Orographic lifting  6  2.2  Frontal lifting of air masses (incoming warm front)  6  2.3  Ice crystals grow through vapour transport  8  2.4  Location of crystal growth axis  9  2.5  ICSI classification of crystal types  10  2.6  Snow Impurities  14  2.7  Measuring snow hardness  15  2.8  The temperature profile of a snowpack  16  2.9  Molecular motion from convex to concave surfaces  17  2.10 Rounding of newly fallen snow  17  2.11 Sintering creates necks between grains  18  2.12 Recrystalization forms depth hoar crystals  19  2.13 Slush  20  2.14 Melt-freeze cluster  21  2.15 Types of mass transport due to wind  21  2.16 Turbulent wind flow can scour out snow near objects  22  2.17 Areas of wind acceleration and deacceleration  23  2.18 Sun crust  •  24  2.19 Perforated crust  25  2.20 Sun cups  25  2.21 Rime forms on the windward side of an obstacle  26  2.22 Tensile stress and tensile fractures  28  2.23 Decomposition of gravity into compressive and shear stresses  29  2.24 Slab release off a house roof  31  2.25 Reflectance as a function of wavelength  32 xii  2.26 Extinction coefficient as a function of density for fine-grained snow  33  2.27 Extinction coefficient plotted against density for a wide range of snow types  34  3.1  Initial undistorted alpine scene  37  3.2  Early spring alpine scene  37  3.3  Optical paths for the shading of snow  38  3.4  Reference pattern  39  3.5  Reference pattern  39  3.6  Examples of rendering detail  40  3.7  Snow-covered car  40  3.8  Parallel computation of falling snow  42  3.9  Falling snow  43  3.10 The failure plane  46  3.11 Considering slices as containers  46  3.12 A scoop loader is a dumping  47  3.13 A synthetic fractal terrain patch before erosion  49  3.14 A synthetic fractal terrain patch after erosion  49  3.15 One point grain and deformations  52  3.16 Correspondence between scales  54  4.1  Avalanches-in-progress are outside the scope of this research  58  4.2  Avalanche debris  58  4.3  Overview of the use of commercial software  60  4.4  Snow falls everywhere in the world  63  4.5  Areas under the bush are not directly visible to the sky yet still receive snowfall  63  4.6  A transition zone between obscured and unobscured  64  4.7  Conservation of mass  65  4.8  Snow creates bridges across gaps  66  4.9  Snow creates cornices and overhangs  66  4.10 Unstable snow moves to a stabler area  67  4.11 Snowy worlds have many smooth curves  68  4.12 Tracks in fallen snow  69  xiii  4.13 Wind shadow of an obstacle  69  4.14 Wind transport causes snow redistribution  70  4.15 An example of varying snow accumulation due to surface blockage  71  4.16 Needs of a falling snow algorithm  78  5.1  Advantages and disadvantages of a volume approach  80  5.2  Advantages and disadvantages of a surface approach  84  5.3  Advantages and disadvantages of a hybrid approach  86  6.1  The model entity  92  6.2  The face entity  93  6.3  Launch sites (white points) are the basic sampling unit  94  6.4  The subdivision area entity  94  6.5  Three geometric objects, plus the +Z axis  95  6.6  The edge group  95  6.7  The snow plane entity  96  6.8  An overview of the snow pipeline  99  6.9  Overview of required and optional phases  100  7.1  Importance ordering  108  7.2  View dependent importance ordering - from the viewpoint  109  7.3  View-dependent importance ordering - from the top  110  7.4  An edge group.  7.5  An object our meshing strategy considers "hard"  113  7.6  The scene is divided up by drops  114  7.7  A Voronoi diagram, with points shown in black  114  7.8  The Delaunay triangulation of Figure 7.7  115  7.9  The constrained Voronoi diagram clips Voronoi areas to face edges  117  •,  112  7.10 Snow on a torus  118  7.11 Top view - an obstacle causes mesh improvement in the transition zone  119  7.12 Side view - an obstacle causes mesh improvement in the transition zone  120  7.13 The final mesh with the initial launch sites shown in white  121  xiv  7.14 The final mesh with all launch sites shown in white  121  7.15 The final mesh with added points shown in red  121  7.16 A mesh with random initial launch site locations  122  7.17 A mesh with 5-time launch site optimisation  122  7.18 A mesh with near-optimal launch site optimisation  122  7.19 Side view of flake motion  124  7.20 Flake motion influenced by random walk  125  7.21 Occlusion boundaries depend on object size  126  7.22 Increased flake variability leads to more snow under an obstacle  126  7.23 As f  127  a  increases, the depression increases in size and depth  7.24 Occlusion boundaries depend on object height  127  7.25 Changing a flake's Z increment test changes the flake's direction  128  7.26 Objects are assigned to buckets within a grid  129  7.27 Locating a flake within sky buckets  132  7.28 Sky-writing input image  133  7.29 Writing on a curved object  133  7.30 An example of flake dusting on tree bark  135  7.31 Real and CGI texture comparison  137  7.32 A bare scene, before view-dependent texture replacement  138  7.33 A tree-covered scene containing a combination of full-depth and view-dependent textures. 139 7.34 All view-dependent textures (highlighted red for better visibility)  139  7.35 Texture interpenetration artifacts  140  8.1  Transition curve for fresh snow  144  8.2  Transition curve for older snow  144  8.3  Obstacle cases  145  8.4  Snow on a very steep surface accumulates behind an obstacle  146  8.5  Figure 8.4 with an invisible obstacle plane  146  8.6  Test stab checks for obstacles  147  8.7  Trajectory of avalanche  8.8  Wind on a snow globe  flakes  148 154  xv  9.1  Implicit surfaces can be used to simulate bridging  157  9.2  Volume inflation must be prevented on flat surfaces  157  9.3  Implicit functions  158  9.4  A single face, converted to implicit functions  159  9.5  The top generator is displaced from the actual snow top plane. The function is clipped to 0 outside the X Y extents of the plane  160  9.6  Parameters controlling bridge creation  160  9.7  The polygonal surface, before bridging  161  9.8  Implicit surfaces representation allows bridging  161  9.9  Implicit surfaces can form unsupported bridging and clumping  162  9.10 Snow material, showing translucence and texturing  166  9.11 A closeup of the snow material, showing translucence and texturing  166  9.12 The original image of the mailbox  167  9.13 The mailbox with a layer of snow  167  9.14 The mailbox with a second layer of snow  167  9.15 The mailbox with a third layer of snow  168  9.16 Color coding shows the layers of snow  168  9.17 Snow surfaces can be modified with additive patterns, simulating tracks and wind-ripples. 169 9.18 Snow stability algorithms can also be used to simulate water accumulation  170  9.19 Snowcovered scenes can be converted to other  171  formats  10.1 Flour compared with CGI snow  173  10.2 Timing - subdivision growth  177  10.3 Timing - intersection buckets  178  10.4 Timing - intersection buckets dropoff  178  10.5 Timing - flakes per group  179  10.6 Timing - sky buckets  179  10.7 Bare gazebo  181  10.8 Gazebo with automatically added snowfall  182  10.9 An older gazebo model without the trees  182  lO.lOAnother view (with different snowfall amounts and properties)  182  xvi  lO.HSanta's bare yard  183  10.12Frame 0  184  10. Ii!Frame 200  184  10.14 Frame 300  '•  185  10.15Frame 400  185  10.16Frame 500  186  10.17Frame 550  186  10.18Frame 650  187  10.19Frame 665  187  10.20Bare grass  188  10.21Snowcovered grass  188  10.22Snowcovered grass, different properties.  189  10.23A snowcovered haystack, inspiration by [Monet 91]  190  10.24The same haystack scene at sunset  190  10.25The bare tree. This model is not particularly detailed  191  10.26The tree with snow  191  10.27The inside of an open-air cloister  192  10.28The snowcovered cloister  192  10.29A bare sign  193  10.30The snowcovered sign  193  10.31The bare SGI logo  194  10.32The snowcovered SGI logo  194  10.33Writing using non-constant snow allocation from the sky  195  10.34A closeup of the same scene, showing "anti-aliasing" due to snowflakes  195  10.35The bare hydrant model  196  10.36The initial hydrant mesh  196  10.37Unstable, snow-covered hydrant  197  10.38The stable result  197  10.39The bare house  198  10.40The snowcovered house  198  10.41The mesh used in the house model  199 xvii  10.42A snowcovered street light, oil painting effect  200  10.43Bare grass, oil painting effect  201  10.44Snow covered grass, oil painting effect  201  11.1 Snow can cling to the underside of objects  209  11.2 Sastrugi (wind formations)  210  12.1 Before  212  12.2 After  212  B.l  231  The road ahead  xvm  Acknowledgements  "This material is based upon work supported under a National Science Foundation Graduate Research Fellowship. Any opinions, findings, conclusions or recommendations expressed in this publication are those of the author, and do not necessarily reflect the view of the National Science Foundation." The following people and groups should be acknowledged for their help in the direction and development of this thesis and the research behind it. Their contribution was greatly appreciated. •  Funding support was provided by the US National Science Foundation (NSF), National Science and Engineering Research Council ( N S E R C ) , the B C Advanced Systems Institute (ASI), the University of British Columbia ( U B C ) , and the Imager Computer Graphics Laboratory.  • Technical support was provided by my supervisor, Alain Fournier. M y supervisory committee and various assorted Imager Lab members provided encouragement and review. • The U B C Faculty of Forestry Center for Advanced Landscape Planning provided invaluable computational assistance in the form of an SGI Onyx 2. • Moral and emotional support was provided by my parents, bad hair and wet polypro were provided by the Varsity Outdoor Club. •  Erika (and Roo and Lili and Chickory and Snowball and Mandibles) reminded me of the future ahead and prevented me from going insane. Or maybe they d i d n ' t . . . N Y A H - H A - H A - H A !  •  A n d finally, spiritual support was provided by the glorious beauty of the West Coast, and the unexplored and secret places that made most of these ideas possible.  xix  Figure 0.1: The author, buried in work. Photo: Mark Grist  xx  Copyrights A n d Credits  Images from A C M publications are copyright of the Association of Computing Machinery. Images from GI are copyright of the Canadian Information Processing Society. A l l photos by Douglas Powell were kindly provided with the permission of the Geo-Images Project [Bain 00]. T h e cover photo was generated with the software described in this thesis. Unless otherwise noted, all real photographs in this thesis were taken by the author.  In the construction of my software, I used code and algorithms from several outside sources, including help from: • Geoffrey S. Heller - Provided source code for "Heller's Marching Cubes", which formed the core' of my marching cubes implementation. •  L E D A - The L E D A [Naher 96] libraries provided invaluable support for list and string handling, as well as the range tree data structure.  •  Dave Eberly - The M A G I C library [Eberly 97] provided some geometric distance routines used to generate implicit surfaces.  •  X F O R M S - This library [Zhao 00] greatly simplified the construction of the user interface.  xxi  Additional Material  The following additional resources may help the reader: •  O n l i n e list of figures - A "clickable" list of figures included in this thesis can be found at: http://www.cs.ubc.ca./nest/imager/contributions/fearing/snow/thesis/lof.html The hope is that the reader can see the original color figures undistorted by the printing process.  • T h e web site - http://www.es.ubc.ca./nest/imager/contributions/fearing/snow/snow.html contains a very brief overview of the problem, and a number of images from this thesis, including a 30-second animation. • T h e p a p e r - A summary of this thesis will be published in the S I G G R A P H 2000 conference proceedings [Fearing 00].  Chapter  1  Introduction  Photo: Douglas Powell  Computer graphics is continually striving to emulate the appearance of the real world. Some of the most interesting problems deal with "natural phenomena", which can be loosely defined as objects created by Nature. In [Fournier 89], Fournier points out some of the challenges of modelling natural phenomena, including: • very complex object descriptions that sometimes cannot be easily described using surfaces (e.g. smoke, fire) • potential for reuse in large quantities, such as leaves on a tree, or trees in a forest • motion that does not often fit easily into a system of rigid transformations, and may need to be built specifically into the model • a dependence on the natural world that makes us biologically alert to its appearance, and allows us to be especially good at picking out visual discrepancies  1  Chapter  I:  2  Introduction  Falling snow and its subsequent accumulation are a prime example of a natural phenomenon. A single snow layer is simultaneously absurdly complex on the local scale and apparently simple on the global level. A smooth, even blanketing of snow consists of a multitude of crystals, each one moving, bonding, changing, and contributing to the overall effect.  When snow moves, it can slip, fold, creep, fracture  or avalanche - travelling like a single handful of flour or an acre of solid concrete. This complexity is both intrinsically interesting, and important in the way it might give us new insights into how other natural phenomena might eventually be modelled. Snow crosses the conceptual boundary between fluids and solids, smooth and rough surfaces and predictable and unpredictable paths of motion, making it simultaneously applicable to other work, yet distinctly its own problem. To date, there has been no previous published computer graphics research on thick snow pack modelling, and only two papers on snow.pack rendering, which clearly makes it an appropriate and worthy topic of research. Beyond novelty, there are other reasons for learning about snow. In many countries, snow is a common fact of life during the winter months. For example, January snow coverage in the Northern Hemisphere has ranged between 41.7 - 49.8 million square kilometres [Robinson 99], or nearly half of the hemisphere's total land mass . A phenomenon that is so common and pervasive to 1  most of us is clearly of interest and importance. If we move beyond simulating snow for its pure challenge and ubiquity, we can see that there are many practical advantages to a working algorithm. Without an automatic model of snowfall, artists must use natural intuition to produce snow covered surfaces - an extremely time consuming task, depending on how realistic one wants the final result. A single tree might have a hundred branches, each with a complex drapery of snow, and each showering snow down onto branches below. Even if one painstakingly created every surface by hand, the image would still likely be very far from what Nature creates in a few hours. A n automatic snowfall algorithm would simultaneously save time, and open up computer generated effects to the possibilities of a whole new season. Besides the practicalities of research and application, there is another reason for investigating snowfall. Snow transforms commonplace scenes into fantastic wonderlands, greatly changing the appearance and mood of the landscape. It is our desire to be able to add snow to some of the classic images produced throughout the development of computer graphics - to see familiar scenes in a fresh, exciting way. This research concentrates on a model of snow pack accumulation and stability for computer graphics, 'High value measured in 1985. Low value measured in 1981. Land mass in the Northern Hemisphere is about 100 million square kilometres. For comparison, the area of Canada is about 10 million square kilometres.  Chapter  1:  Introduction  3  with a minor excursion into snow pack rendering. It is organised as follows: • Chapter 2 contains basic background information on snow - how it is created, how it changes, and what properties might be relevant in a simulation. • Chapter 3 provides an overview of previous computer graphics attempts at drawing snow or similar materials. • In Chapter 4, we set out our criteria for success, make some assumptions about the scope of the problem, and introduce some of the main obstacles that must be overcome. • Chapter 5 describes some competing and alternative approaches we considered and discarded en route to our final solution. • Chapter 6 introduces the algorithm and some important data entities, and gives an overview of the important phases. • Chapters 7 and 8 respectively describe the various components of our accumulation and stability algorithms, and how they relate to each other. This is the heart of our research, and contains our primary contribution to the problem. • Chapter 9 describes smooth surface generation and snow pack rendering, issues outside of the main thrust of our research but still important for creating images. In addition, there is some discussion of miscellaneous issues and side effects of our approach. • Chapter 10 helps the reader validate the model with cross-references between initial goals and the thesis sections containing validation of those goals. As well, this chapter contains a gallery of more complex examples, generated using all the components of our algorithm. • In Chapter 11 we describe some of the drawbacks associated with our approach, and indicate some areas of future research. • Chapter 12 contains the conclusions of this research. For the convenience of the reader, this thesis also includes some additional sections: • Appendix A contains list of all variables referenced in the text. • The Glossary describes the meaning of various relevant terms. • The Index provides a cross-reference for terms of interest.  Chapter 2 Snow  Fresh snow on a stream.  This chapter describes the creation and evolution of the snowpack. It contains a brief summary of the conditions required for snow, the different forms snow takes, metamorphic changes in fallen snow, and certain non-precipitated effects that are often mistaken for snow. We also discuss some material properties that are important indicators and characterisers of the snowpack. Before the reader gets too excited, we must admit up front that our model ignores, omits, and oversimplifies - almost without exception - every single one of these properties, in order to deal with computer-related issues of speed and complexity. As a result, we have included this chapter primarily to give the reader a brief understanding of the complexity of the phenomena.  1  Even the simplest and  crudest computer approximation of any natural phenomena must first start with an appreciation of the 1  See [McClung 93] and [Daffern 92] for a much more detailed discussion of snow and all its properties.  I  5  Chapter 2: Snow  ideal goal, however unapproachable it may seem. However, our approach is flexible enough to eventually allow many of these properties to play their rightful role in the simulation. To encourage the reader, we have included Section 11.1, which lists these properties and briefly discusses how they can be fully included in future implementations. This Section occurs near the end of this thesis, where the reader is aware of our solution and its implicit problems and advantages.  2.1  Snow in the A i r  This section summarises the origin and evolution of snow before it hits the ground. The atmospheric conditions that cause snowfall greatly influence the size, shape, and material properties of snow in any subsequently accumulated snowpack.  2.1.1  Precipitation  There are four major methods of large-scale air mass movement leading to snowfall. These are orographic lifting, cyclonic convergence, frontal lifting and convection. A l l of these methods cause warm air to rise, condense, and drop precipitation. Orographic lifting (lifting due to topography) is the most important moisture- transport mechanism in Western Canada, producing a large portion of the winter snowfall. Warm, moist air from the Pacific Ocean moves horizontally across the coastline, where it rides up the western side of the Coastal range. The mountains compress air into a smaller volume, causing wind acceleration and forcing the incoming air mass to lift very quickly, as shown in Figure 2.1. As the incoming air rises, the pressure decreases and the air temperature drops, eventually reaching the dew point. The air at the dew point is saturated with water vapour; it contains as much free water vapour as possible for a given temperature.  Beyond  this, water vapour begins to condense into water droplets. The warmer the air mass, the more moisture it may hold before it reaches the dew point; accordingly, maritime coastal regions usually receive more snowfall (and rainfall) than interior regions. Water droplets form around tiny (approx 1 0  - 6  m) particles called condensation nuclei, which are  generally airborne motes of dust or soil pervasive throughout the atmosphere. These droplets grow larger as free water molecules condense on the outer surface, eventually providing enough mass to leave the droplet unsupported by the atmosphere.  Chapter  2:  fi  Snow  fast rising, moist air  moutain range Figure 2.1: Orographic lifting  Orographic lifting is important because it causes such a quick upward movement of the air mass, where the faster the upwards motion, the greater the condensation and subsequent precipitation. This means that steep mountain ranges oriented perpendicular to the airflow will receive greater precipitation than more rounded mountains lying off-perpendicular. There are other methods of large scale air mass cooling and condensation. Cyclonic lifting is the circulation of air around a (potentially very large) centre of low pressure. As air spins into the low pressure centre, it forces existing air upwards.  precipitation due to lift cooling  warm air rides up over dense cold air  Figure 2.2: Frontal lifting of air masses (incoming warm front) Frontal lifting occurs when an air mass of one temperature bumps into another mass of a different temperature. The cold, dense air remains at the bottom, causing the warmer, lighter air to ride upwards.  7  Chapter 2: Snow  Precipitation is usually restricted to the boundary area of the two air masses. A moving warm front causes a shallower slope than a moving cold front, resulting in precipitation over a wider area. A n example of warm front frontal lifting is shown in Figure 2.2. Convective lifting occurs when surface air is directly warmed by the sun, causing it tp rise in place. This is usually a local effect, producing only small amounts of precipitation. Moisture transport methods operate at many different scales. Just as large air masses are pushed over coastal ranges, smaller air masses are pushed over bumps and valleys of an individual range of mountains, or even a single mountain. This often results in a wide variation in snowfall within a local area.  2.1.2  F o r m a t i o n of Ice C r y s t a l s  Snow formation starts to occur when the air temperature of the water-bearing cloud drops below 0 ° , allowing ice crystals to form around a type of foreign particle called a freezing nucleus. Like condensation nuclei, freezing nuclei also consist of dust and soil particles, although the conditions for water vapour freezing onto these surfaces are more stringent than for condensation. Only certain molecular structures will suffice, depending on temperature: McClung states that "at —10°C there are about 10 active nuclei per cubic centimeter"  [McClung 93], compared with condensation nuclei densities of 150,000 per c m  [Daffern 92]. As the temperature decreases, freezing conditions relax, increasing the number of potential freezing nuclei, eventually reaching —40° C , where water droplets can freeze without the help of a freezing nucleus. It is hypothesised that when previously frozen crystals collide and shatter, the resultant ice fragments can also form freezing nuclei, effectively speeding up the process once it has started.  2.1.3  Ice C r y s t a l G r o w t h  Once initially frozen, there are three main processes for atmospheric ice crystal growth. The first process involves vapour movement, where a pressure differential causes water molecules to migrate from water droplets towards ice crystals. Incoming water vapour molecules condense on the ice crystal, as shown in Figure 2.3. Crystal formation due to vapour motion generally causes growth along two defined axis, as shown in Figure 2.4. The basal plane forms a flat plane pierced by the c-axis, and divided into three 120° segments, each separated by an a-axis. The crystal structure of the ice grain forces all ice grain to have  3  8  Chapter 2: Snow  water droplet  Figure 2.3: Ice crystals grow through vapour transport six sides. Growth along the c-axis produces needle or column-like crystals, growth along the a-axis produces flat plates, and deposition on a combination of these two axis produces a wide variety of possible shapes. Depending on environmental conditions, the primary growth direction may switch between the a-axis and the c-axis, where the main factors determining crystal growth direction are temperature and the density of the excess water vapour on the surface of the crystal. The International Commission on Snow and Ice (ICSI) has created a classification system for crystal types, as shown in Figure 2.5. The classification system organises newly fallen snow into categories based upon general shape and creation. There are even more complex crystal classification systems [Magono 66]. The initial shape of an ice crystal can vary considerably from the final crystal shape that eventually reaches the ground. As a crystal passes through different layers of the atmosphere, it may encounter different temperatures and moisture conditions that may encourage growth of a different shape-type on the initial crystal. Turbulence in the atmosphere may repeatedly cycle a crystal through different atmospheric layers, producing an extremely complicated growth pattern.  9  Chapter 2: Snow  120°  A  c axis  a axis  a axis  a axis  Figure 2.4: Location of crystal growth axis  Crystal Growth Due to Riming  Besides vapour deposition, riming is another process for atmospheric ice crystal growth. As a crystal gains size, it increases its chance of a collision with super-cooled water droplets still in the air. These droplets freeze onto the main crystal structure, with effects ranging from a slight modification of the basic shape to a complete rounding. In very heavily rimed crystals (called graupel), the initial crystal structure is impossible to determine. Riming collisions occur both during a snowflake's gravitational descent, and during upwards motion from turbulent air - thus, as the atmosphere grows more turbulent, more riming occurs.  Crystal Growth Due to Collision  The final method of crystal growth occurs when individual ice crystals collide and matte together, producing larger super-crystals (or snowflakes). This occurs most easily with stellar and dendritic crystals (large branches with many contact points) falling in warm, moist conditions. Collisions are also thought to create small ice fragments that then act as freezing nuclei for subsequent crystal formation.  2.2  Snow on t h e Ground  As snow falls from the sky, it lands on the ground, forming a snow layer and eventually accumulating with other layers to form a snowpack. Because individual snow layers are usually deposited under varying  Chapter 2: Snow  la Columns  1,  1b Needles  1c Plates  o  1d  Hail  ih Ice pellets  turation at -3' to -9"C  or hollow  a n d betow  Weeoie-tike,  G r o w t h at high super sa-  approx  turation at  -22T  3 ' to - 5 ' C  Growth at high s u p e r s a -  Plate-like, mosrtty  turation at 0* to 3 * C  hexagonal  and  Sx-fotd  G r o w t h at h i g h super sa-  star-iike,  turation at temperatures  p l a n a r or  b e t w e e n -12*  ~8* to  -25X  to  16*C  C l u s t e r s of  Porycrystals g r o w i n g at  very small  varying environmental  crystals  conditions  A  Heavily r i m e d particles  H e a v y riming of partsctes b y a c c r e t i o n of s u p e r c o o l e d water  k  Laminar internal structure, ttansiucent or milky, g l a z e d surface  Irregular particles  ig  G r o w t h at high supe^sa-  spatial  le  11  Short prismatic crystal, softd  cylindrical  Stellar Crystals  Graupel  )  A  G r o w t h by accretion of s u p e r c o o l e d water  Transparent, mostly smaii spheroids  Figure 2.5: ICSI classification of crystal types, with permission from [McClung weather conditions, each snow layer often represents a different type of snow, with its own physical properties. Although a single snow layer is often fairly homogeneous, there are several possible counterinstances where snow properties can vary widely within a layer, such as wind deposit or avalanching from above. In between snowfalls, all layers are affected by a wide number of environmental factors that can accentuate or remove the differences between layers, and greatly change the physical properties of the entire snowpack. Snow layers continually change until they melt, evaporate, or become ice.  11  Chapter 2: Snow  2.2.1  I n i t i a l C o n d i t i o n s of Fallen Snow  The ICSI has attempted to standardise the description of snow properties. This is of extreme importance in avalanche work, where evaluating overall snow stability requires a great deal of local snowpack observation. To this end, the ICSI system characterises individual snow layers by a set of primary physical variables, given in Table 2.1. These variables are used to make the main distinctions between different types of snow. Some of these variables (most notably snow strength), are discussed in more detail in further sections. Section 11.1 describes how these properties can be hooked into future development. Table 2.1: Primary Physical Characteristics of Deposited Snow Feature Density Grain shape Grain size Liquid water content Impurities Strength (compressive, tensile, shear) Hardness index Snow temperature  Units kg see Section 2.1.3 mm diameter % by volume % by weight Pa  Symbol  depends on instrument °C  R  P F E  e j E  T  Adapted with permission from [Colbeck 90]  Density  Snow density (p) is measured as the mass of snow per cubic meter. It is also often measured as an average porosity (1 —  ), or a specific gravity (the ratio of the weight of snow to the weight of an  equal volume of water). Typical snow densities are given in Table 2.2. As can be seen from the table, the various types of snow may have densities that differ by several orders of magnitude. At the extreme end of the scale is the so-called "wild" snow, which falls in very cold and totally windless conditions, forming the one of the most porous materials in nature. As the density of the snowpack increases, there is a corresponding increase in the the number of intra-grain bonding sites per unit volume. These bonds (described in Section 2.2.3) provide coherency to the snow, leading to a direct relationship between density and snow strength. Density also affects the radiometric properties of snow, as it affects the amount of intergrain scattering and reflection.  12  Chapter 2: Snow  Table 2.2: Typical Snow Layer Densities Type of Snow wild snow dry flakes loose dry snow wind packed snow cornice snow wet snow glacier ice hard ice water  Density ( ^ ) 3 30-50  50-90 120-300 400 500-830 800 917 1000  Adapted from [Daffern 92] and [Upadhyay 95]  Grain Shape  Grain shape was previously summarised in Figure 2.5. In addition to the shapes of newly fallen crystals, there are several forms caused by snowpack metamorphism. These new shapes and the methods that produce them are described in Section 2.2.2. Crystal shape affects the reflectance properties of the snow surface - flat, faceted surface hoar reflects light in a more anisotropic way than more rounded grains.  Grain Size  Snow crystal sizes are measured according to an ICSI size classification, given in Table 2.3. In a homogeneous snow layer, the grain size is the approximate average maximal diameter of all ice crystals, measured in millimetres. As the average grain size of the snowpack decreases, there is a corresponding increase in diffuse reflection, and in the rate of extinction of visible light.  Table 2.3: Grain Size Classification Term Very fine Fine Medium Coarse Very coarse Extreme  Size (mm) < 0.2 0.2-0.5 0 . 5 - 1.0 1.0 - 2.0 2.0-5.0 > 5.0  Adapted with permission from [Colbeck 90]  13  Chapter 2: Snow  2.2.1 A  Liquid Water Content  Water content is expressed as a percentage of free water by volume. Table 2.4 provides a general idea of the moisture ranges of snow.  Table 2.4: Liquid Water Content Term Dry  Moist  Wet  Very Wet  Slush  Remarks Usually T is below 0°C, but dry snow can occur at any temperature up to 0 ° C Disaggregated snow grains have little tendency to adhere to each other when pressed together, as in making a snowball. T = 0°C. The water is not visible even at 10 X magnification. When lightly crushed, the snow has a distinct tendency to stick together. T = 0 ° C The water can be recognised at 10 x magnification by its meniscus between adjacent snow grains, but water cannot be pressed out by moderately squeezing the snow in the hands. T = 0°C. The water can be pressed out by moderately squeezing the snow in the hands, but there is an appreciable amount of air confined within the pores. T = 0°C. The snow is flooded with water, and contains a relatively small amount of air.  Approximate Range of 9 0%  < 3%  3-8%  8 - 15%  > 15%  Adapted with permission from [Colbeck 90]  The amount of water content in the snowpack has an effect on both the snowpack strength and appearance. When the snow is wet, snowpack evolution is influenced by free water moving through the snow layers, promoting rapid changes in the snowpack, reducing the cohesion between snow grains and generally reducing the snow strength. The amount of water in a particular snow layer also contributes to the total weight pressing down on the underlying snow layer, influencing the settling of lower layers, as well as the speed at which the entire snowpack moves over the ground on a slope. Wet snow also generally reflects less light than dry snow, mainly because wet snow contains so many large, rounded ice crystals (see Grain Size, above).  ( 'hapter 2: Snow  Snow Strength  The strength of snow helps determine when (or if) a particular area of snow will remain in the same place or maintain the same shape. Under the right set of strains, snow will deform, move slowly, or move quickly (avalanche). Because this property is so important to any simulation of snowpack, we discuss it on its own in Section 2.4.  Impurities  Impurities are rarely present in sufficient quantities to affect anything but the visual properties of snow, but may provide useful visual clues during rendering. Some common impurities include pine needles, dirt, and various algae and fungi that can turn snow green, yellow, or even red (Figure 2.6).  Figure 2.6: Red snow, caused by dynoflagelletes.  Snow Hardness  Snow hardness measures the resistance of the snowpack to an intruding object, where the harder the snow, the greater the snow strength. Hardness measurements are fairly subjective, and depend on the measurement instrument and measuring technique. Common instruments include fingers, fists, pencils,  Chapter 2: Snow  and the ram penetrometer (a vertical pole driven into the snow by a calibrated hammer).  Figure 2.7: The author preparing to measure snow hardness with the "face penetrometer"  Snow Temperature  Snow temperature is measured in "C, and varies with the the depth of the snowpack. The change in temperature over the change in snowpack depth is defined as the temperature gradient, and is measured in ° . The temperature gradient runs from the top of the snowpack, where the temperature is influenced by the atmosphere, to the bottom of the snowpack, where the temperature is influenced by the ground. The ground interface temperature is usually near 0° due to stored summer heat, while the air interface temperature varies dramatically, depending on the weather. Temperature gradient is usually measured in the direction of increasing temperature (down), as shown in Figure 2.8. Rocky outcroppings complicate the prediction and modelling of the temperature gradient, as they cause lateral variations in snow temperature and compress the vertical profile. The temperature gradient is extremely important in the metamorphism of the various snow layers, as it strongly influences the ratio of strengthening vs. weakening processes occurring simultaneously. In coastal climates, the ambient air temperature is usually fairly warm, and the snowpack is quite deep, implying a gentle temperature gradient.  16  Chapter 2: Snow  -20° I  |  -10° | [  0° |  snow surface  -20° -10° I I I I  outcropping causes a different gradient  0° I  rock outcropping  Figure 2.8: The temperature profile of a snowpack  2.2.2  Initial Rounding and Settling  As we have previously described in Section 2.1.3 and shown in Figure 2.5, the shape of new individual ice crystals is related to temperature, vapour density, wind speed and turbulence at the various atmospheric levels above the ground. As soon as snow hits the ground, it immediately begins to change. The primary agent of crystal change is temperature, with much lesser effects due to the natural minimisation of ice crystal surface area to volume. We are interested in these changes because they have implications for the strength, stability, crystal shape, density, and moisture content of the snowpack, all of which have been experimentally shown to influence the visual appearance of the top layer.  Metamorphism due to Implicit Crystal Shape  Ice grains start an initial process of rounding to minimise surface area, where crystal types with a large ratio of surface area to volume are the most prone to change. Instability is related to the amount and type of curvatures in a particular crystal, with a sphere being the most stable form. Concave surfaces have a lower vapour pressure than convex surfaces, encouraging water molecules belonging to complex branch-structures to move into the atmosphere before molecules belonging to gently curved structures (Figure 2.9). Thus, dendrites and stellar-type crystals metamorphize the fastest, while rounder crystals remain stable for longer. Rounding causes a reduction in the particle size, as dendritic branches are destroyed.  Chapter 2:  17  Snow  ice  crystal  ice  crystal  Figure 2.9: Molecular motion from convex to concave surfaces The destruction of interlocking branches also leads to an immediate loss of strength, as well as settling and densification. The initial rounding process does not last for very long, as curvature effects are fairly small compared to the longer-term effect of the temperature gradient. Figure 2.10 shows the rounding process.  Figure 2.10: Rounding of newly fallen snow, from [Colbeck 80]. The leftmost image shows the initial crystal, the rightmost image shows the rounded crystal. The next step is a general intra-particle transfer of water vapour. The warmest layer of the snowpack (adjacent to the ground) has the ability to contain more water vapour than the colder layers above. As this layer fills with water vapour, there is an upward pressure and upwards vapour movement. When  IS  Chapter 2: Snow  the temperature gradient is small (below about 10^-), the water vapour has a tendency to condense on larger ice particles and escape from the smaller ones, causing large particle growth at the cost of small particle shrinkage and increasing the average grain size over time.  Rounding (also called destructive  metamorphism) is the primary metamorphic process when the temperature gradient is small.  2.2.3  Sintering  As the snowpack consolidates and settles, the rounded ice grains undergo further settling via a process called sintering. Water vapour has a tendency to move from convex surfaces to concave surfaces, causing condensation and buildup at locations where ice grains touch. The touching areas develop into bridges, called "necks". The stronger and thicker the connections between ice grains, the stronger and more consolidated the entire snowpack.  Figure 2.11: Sintering creates necks between grains. Photo: E. Akitaya, from [Colbeck 90] The speed of sintering increases as the snow temperature increases (temperature gradient decreases).  2.2.4  Crystal Growth and Recrystalization  When temperature gradients are low, individual grains undergo rounding, or destructive metamorphism. There is an opposing process of recrystalization (constructive metamorphism) that occurs when the temperature gradients are high (air temperatures are low). Vapour originates from the bottom of the snowpack (where the heat is greatest), and moves towards the surface. This vapour deposits itself near the upper snowpack layers, increasing average crystal size.  19  Chapter 2: Snow  Locally, water vapour molecules move from the bottom of rounded ice grains to the top, making the bottom surface more faceted, or angular. At the extreme, faceted crystals can form depth hoar, which are often large, cup-like, hexagonal crystals, as shown in Figure 2.12. The increased angularity, pore spaces, and crystal sizes all contribute to a reduction in the internal cohesion of the snow layer. Shear failure between snowpack layers often occurs at the interface between a depth hoar layer and the layer above. Given the correct conditions (usually shallow snow depth and cold temperatures), the entire snowpack can eventually transform into depth hoar.  Figure 2.12: Recrystalization forms depth hoar crystals. Photo: K. Izumi, from [Colbeck 90] Recrystalization occurs in the direction of the temperature gradient, which is not always vertical with respect to the snowpack. Buried rocky outcroppings cause lateral temperature gradients and thus lateral recrystalization movement.  Depending on the weather, the primary snowpack metamorphism process  can alternate between rounding and recrystalization. The two processes often occur simultaneously at different snowpack layers, or at different spatial locations on the same snow-slope. The process of faceting vs. rounding depends on the growth rate of the crystals, which is in turn mainly determined by the temperature gradient, which is in turn controlled by local variations of the terrain and the air and ground temperature. A much more detailed description of the recrystalization process can be found in [Colbeck 82].  m  Chapter  2:  Snow  2.2.5  Wet Snow Metamorphism  Wet snow is even  more complex than dry snow because it contains appreciable amounts of water in  all three states. Change in very wet snow is usually driven by differences in the individual grain sizes. Larger rounded grains freeze slightly on the surface, giving off a small amount of heat, and causing smaller grains to melt. The temperature differences between the larger and smaller grains is minute, but is amplified by the heat-conductive properties of the surrounding water. The reduction of the smaller particles and the increase of the larger particles quickly inflates the average grain size of the snowpack. The weight of the snowpack above applies pressure at individual grain boundary points, reducing the melting temperature, and causing grains to melt at contact points, destroying previously formed bonds. Thus, melting snow undergoes both densification and a reduction in internal snow cohesion. When water content exceeds 15%, the snow layer is usually composed of rounded grains entirely separated and surrounded by water. This snow form is called slush (Figure 2.13), and is weak enough to avalanche at very shallow slope angles.  Figure 2.13: Slush. Photo: S. Colbeck, from [Colbeck 90] Water has a tendency to drain away from the slush layer (unless prevented by a impermeable lower layer), so extreme water contents usually do not last long. If the temperature drops during the night, the water between the ice grains freezes into a large globular mass called a melt-freeze cluster, shown in Figure 2.14. This produces a melt-water crust, or a thin layer of ice grains joined together by frozen ice bonds. During the day, the sun will warm and melt the clusters. As the cycle progresses, these clusters get larger and larger.  21  Chapter 2: Snow  Figure 2.14: Melt-freeze cluster. Photo: S. Colbeck, from [Colbeck 90]  2.2.6  W i n d and Snow  Wind is a very important agent in the large-scale transport of snow. There are three methods of snow transport by wind (Figure 2.15), generally classified by the height at which the majority of mass is moved, which is in turn influenced by the wind speed.  suspension (many meters)  saltation 10 cm vertical displacement • >.  rolling 1 mm vertical displacement  Figure 2.15: Types of mass transport due to wind When snow is dry and loose, small particles will roll or creep over the surface of the snowpack, with a typical crystal vertical displacement of 1 mm. When winds remain between 5 and 10 m/s (1836 km/h), particles are swept up and bounced along within 1-10 cm of the snow surface, in a process called saltation. Saltation and rolling generally occur when the winds are in a laminar flow. When the windspeed rises, wind motion becomes turbulent, capturing snow in eddies and carrying it a considerable  Chapter 2: Snow  22  distance above the snow surface. This method of snow transport is called suspension, and is responsible for wider scale snow movement, and loss of visibility. The windspeed threshold for suspension is around l )in/s (51 km/li). Turbulent eddies can be quite transitory, moving over the surface of the snow adding r  and removing snow. Some examples of eddy-caused wind effects are the "scour holes" around obstacles, as shown in Figure 2.16. Wind hits an obstacle, swirls around picking up snow, and then moves around the object, possibly depositing some of the snow on the lee side.  Figure 2.16: Turbulent wind flow can scour out snow near objects. Photo: Douglas Powell It should be noted that the speed at which snow is picked up from surfaces is not well understood - thus, the above given threshold values are quite approximate. In general, the more cohesive the snow, the greater the wind speed required to pick up crystals, with greater snow transport if crystals are fine and fresh. McClung says "for loose unbonded snow, the typical threshold wind speed (at a 10 m height) is 5 m/s (18 km/h). For a dense, bonded snow cover, winds greater than 25m/s (90 km/h) are necessary to produce blowing snow." [McClung 93] Although snow transport due to wind can occur anywhere, it can be especially obvious in places where there is fairly sharp ridgecrest, where deposited snow can form cornices. As air travels across the ridgecrest, it is compressed, forcing it into a smaller volume and increasing the wind speed. During this period of acceleration, wind picks up snow, scouring the windward side of the ridge. When the wind crosses the ridgeline, it is now free to expand, reducing wind speed and depositing carried snow on the lee side of the obstacle, as shown in Figure 2.17.  23  Chapter 2: Snow  s c o u r e d area  dropped snow  cornice  F i g u r e 2.17: Areas of wind acceleration and deacceleration  As wind moves snow across the surface, individual particles collide and shatter, producing very fine and very dense deposits, with a dull white color.  2.2.7  Additional Snow Forms  Certain types of environmental conditions produce snowpack with very distinctive conditions. We describe a few of these.  Sun Crust  When sun shines on fresh snow, it causes melting in a thin layer just under the surface. The free water is later refrozen to produce a strong, icy layer, or crust. Over time, this crust grows, although it rarely gets thick enough to support any weight. Figure 2.18 shows a sun crust surface.  Perforated Crust  Depending of the angle of the sun, and intervening obstacles, some ice grains receive more solar energy than others. These grains melt first, producing tiny depressions, which continue to enlarge, changing shape as they are influenced by sun direction. Figure 2.19 shows a perforated crust.  24  Chapter 2: Snow  Figure 2.18: Sun crust. Photo: E. Wengi, from [Colbeck 90]  Sun Cups  As winter progresses, natural depressions have a tendency to enlarge and increase, much like a perforated crust. At the end of the season, many of these depressions have grown quite large, forming a distinctive pattern of sun cups, as shown in Figure 2.20.  2.3  Snow-like Effects  There are a several other atmospheric processes that produce snow-like effects, but are not caused in the same manner. The two most important forms of non-precipitated snow are rime and hoarfrost.  2.3.1  Rime  Rime occurs when the low atmosphere is filled with supercooled water droplets. The wind moves these droplets laterally over the ground, where they collide with obstacles, freeze to the surface of the obstruction, and create buildups pointing into the wind. The stronger and longer the wind blows, the greater the length of the rime feathers. Figure 2.21 shows rime effects.  25  Chapter 2: Snow  Figure 2.19: Perforated crust  Figure 2.20: Sun cups. Photo: Douglas Powell  2.3.2  Surface Hoar  Surface hoar (or hoarfrost) is a layer of large ice crystals that forms on the snowpack surface, much like dew. It usually occurs in the presence of a warm layer of air lying just above the snowpack, which might occur when sun warms the snow, with a subsequent reflection of heat from the snow surface. The warmth allows the air to maintain a relatively high density of water vapour, which condenses into hoar crystals as the surface temperature drops during the night. Individual crystals are angular and faceted, and range in size between 1 mm and 1 cm. The extreme faceting of surface hoar causes a distinctive anisotropic scattering of reflected light. Surface hoar also occurs on surfaces other than snow; the same process that deposits crystals on  Chapter 2: Snow  26  Figure 2.21: Rime forms on the windward side of an obstacle. Wind blows from left to right. Photo: Adam R. Jones snowpack also results in crystals on plant surfaces, car windows, or the sidewalk.  2.4  Snow Strength  The strength of the snowpack is its ability to resist applied stress. We are concerned with three main types of applied stress: compressive stress, shear stress and tensile stress.  2.4.1  Compressive Stresses  The primary compressive stress on snow is the force of gravity. Gravitational force pushes all the grains together, rearranging them in a closer, tighter arrangement and leading to a denser, harder and stronger snowpack. This densification process is called settlement, and occurs at a wide range of speeds: "Rates of settlement vary widely in alpine snow, ranging from 10 cm per day or more in low-density newly fallen snow to perhaps 1/100 mm per day in dense snow at the end of the winter under a deep snow cover." [McClung 93] Local compressive stresses may also be exerted by an object on top of the snow (such as a car tire  Chapter 2:  27  Snow  or a foot).  2.4.2  Shear Stresses  Shear stresses are applied by gravity in a direction parallel to the snow slope, causing, lateral shear deformations (strains). This shear deformation is resisted by the shear strength of the snowpack, which is a product of snowpack internal friction and cohesion. As individual particles slide past each other, the speed at which they slide is influenced by interparticle friction, which is in turn controlled by the grain texture, water content, and snowpack depth. Grain texture describes the shape of the grain (heavily branched vs. smooth plates, etc.)  and size  of the grain (more vs. less inter-grain contact area). Increased moisture content causes lubrication of particles and reduced friction during motion, while a large snowpack depth applies additional weight that compresses particles together and increases the frictional force.  2.4.3  Tensile Stresses  Tensile stresses act to pull apart adjacent ice grains, usually in a lateral direction, and often in conjunction with shear deformation. These effects can occur when a downslope area of the snowpack has less ground friction than an upslope area, causing the bottom snow area to "pull away" (sometimes suddenly). When the failure is sudden, the top (crown) area of a slab will undergo a simultaneous tension fracture, as shown in Figure 2.22.  Decreased downslope friction is often caused by poor shear strength in the  snowpack.  2.4.4  A n g l e of Repose  The angle of the surface is an extremely important contributing factor to the overall stability of the snowpack. Shear forces operate parallel to the slope of the snowpack, while settlement (compression forces) occur perpendicular to the slope of the snowpack. These two forces work to opposite ends settlement promotes densification and slope stability, while shear stress causes slope failure. Figure 2.23 shows the gravitational force broken into settlement and shear components. As can be seen, as the angle of the slope increases, the compressive forces are reduced, the snow settles less, and the snowpack becomes weaker and less stable.. The angle of repose (also called the angle of static friction) is the angle beyond which a particular snow  Chapter 2: Snow  reduced friction along the failure plane  ground  Figure 2.22: Tensile stress and tensile fractures type is highly likely to move. The angle of repose is influenced by the variables controlling snow strength, such as the irregularity of grain geometry, temperature, cohesiveness, and water content. Table 2.5 shows some experimental angle of repose measurements for dry snow, adapted from [Kuroiwa 66].  2  From the data, one can see that as the temperature increases, dry snow bonds at a faster rate, generally leading to more cohesion and greater angles of repose. The difference in surface irregularity between natural dendritic snow and pulverised snow reflects in dendritic snow's ability to cling to steeper surfaces. Other, more rounded crystal types sluff at smaller angles, with a fair difference in the range of possible angles. The angle of repose is critical to the visual appearance of the snowpack, because it controls which surfaces can support snow. In wet snow, the amount of free moisture is often the most important determining factor, where completely saturated snow (slush) can have an angle of repose as low as 15°. It is also important to remember that as snow changes over time, it will also undergo a corresponding variation in its possible angle of repose. Once local action initiates a small avalanche, it will continue to move down the slope until the slope angle is less than the angle of kinetic friction. This angle is based upon the same variables that determine the angle of static friction, but is usually about 10° less. As snow crosses into these shallower slopes, This experiment also found that the measured angle-of-repose varied with the height of the dropped particles; primarily because the impact of the falling crystals "rounded off" areas of otherwise high cohesion. 2  Chapter  29  Snow  slope angle a  ground  Figure 2.23: Decomposition of gravity into compressive and shear stresses friction works against the movement of the snow and begins to slow it down. As avalanches gain speed, other braking methods come into play, including air drag of the advancing front, and sub-surface interaction as the avalanche ploughs under layers of snow.  2.5  Avalanches  Avalanches are generally classified into loose-snow avalanches and slab avalanches, where there is a considerable difference in the causes and effects of the two types of avalanches. Avalanche prevention and detection is very important to many people, including ski resorts, road maintenance crews, and backcountry travellers. As a consequence, the causes, nature, and triggering effects have been studied quite extensively, and in fact provide a major source of information about snow. We refer the reader to one of several excellent books on the subject [Daffern 92] [McClung 93] for material beyond the scope of this thesis.  2.5.1  Loose Snow Avalanches  Loose snow avalanches usually originate on the topmost layer of the snowpack, within an area of unconsolidated surface snow. When enough snow accumulates to overcome the angle of static friction, snow slips, and an avalanche begins. According to McClung [McClung 93], the initial mass of snow in motion can be very small (less than 1 0 m ) . The actual catalyst for avalanching in slopes at or very near the _4  3  30  Chapter 2: Snow  Table 2.5: Experimental Angles of Repose Snow type pulverised snow pulverised snow pulverised snow pulverised snow pulverised snow pulverised snow pulverised snow pulverised snow pulverised snow natural dendritics natural dendritics  grain size (mm) 0.5-0.6 0.5-0.6 0.5-0.6 0.5-0.6 0.5-0.6 0.5-0.6 . 2.0-5.0 0.84-2.0 0.42-0.84 -  temperature (°C) -35 -20 -12 -5 -3.5 +1 -35 -35 -35 -35 > -12  angle of repose (°) 45 45 47 50 55 near 90 40 41 42 63 near 90  Compiled from [Kuroiwa 66]  angle of static friction can vary, with several possibilities including the loss of cohesion due to the effects of metamorphism and the environment, sloughing from above, or excess force due to the passing of a skier. Because the triggering volume is so small, loose snow avalanches typically start at a point and spread outwards in a triangular shape (assuming a fairly homogeneous slope). Unlike slab avalanches, snow does not fail at well-defined layers in the snow, or originate along fracture lines. As the initial snow mass descends, it can sweep along other unconsolidated snow in its path. The total mass of a loose snow avalanche depends greatly on the moisture content of the snow - wet, loose snow avalanches can be quite massive and can sometimes trigger slab avalanches on lower slopes. During storms, many steep surfaces self-stabilise through a continual sloughing-off of snow. These are essentially tiny snow avalanches that move snow to lower-incline slopes before it can build up into a large, consolidated snowpack. This is important to us because it produces obvious visual results, such as snow accumulating in gullies.  2.5.2  Slab Avalanches  A slab is a large chunk of cohesive snow that usually breaks off along well- defined fracture lines, and travels over a distinct weak layer in the snowpack (Figure 2.24). If the fracture layer is particularly deep in the snow, a slab can be many meters thick, with a mass large enough to cause great destruction. As such, the causes of slab formation and slab release are of primary interest to many, and have been the subject of a great deal of experimentation and observation. However, for reasons described later in Section 4.1.4, this thesis is not concerned with modelling slab release, and so we refer the interested  31  Chapter 2: Snow  reader to the references for more detail [McClung 93] [Daffern 9*2].  Figure 2.24: Slab release off a house roof. Photo: Paul Fearing  2.6  Illumination Properties  In this section, we briefly characterise the optical properties of snow, drawing upon experimental research from the area of snow sciences.  2.6.1  Reflectance  Snow reflects a large majority of the incoming light in a fairly even manner (biased towards blue) across the visible spectrum, giving the appearance of "white" or "blue-white". Because of the generally similar treatment of the spectrum, researchers often measure albedo instead of wavelength (albedo is the percentage reflectance of incoming light integrated over all wavelengths). Albedo varies for different types of snow, from dry packed snow (0.78 - 0.93), to moist snow (0.73 - 0.78) to wet snow (0.56 - 0.72) . 3  As the grain size of the snowpack decreases, there are more scattering and reflecting boundaries, and the albedo increases. Conversely, as the grain size increases (wetter snow), the albedo decreases. The age of the snowpack is also important, as it gives an indication of the amount of surface dust and impurities on the surface layer. 3  All data from [Mellor 77].  32  Chapter 2: Snow  There is still some uncertainty as to where the majority of reflected light occurs; Mellor [Mellor 77] provides an overview of the competing experiments. In short, the prevailing view is that most of the reflectance occurs very close to the surface: "Measurements of reflectance as a function of thickness for a finite snow layer, showed that [the albedo] became effectively independent of thickness above 10 to 20 m m in snow with 0.5mm grain size". Other experiments indicate that a 5 m m depth is the minimum required to conceal any sublayer, and that sub-surface scattering is responsible for only 1 % of the reflection in fresh dry, snow. It should also be pointed out that there actually is some wavelength dependence with respect to reflectance, as shown in Figure 2.25, although these measurements are for certain grain characteristics only.  ^ ^ 1  1  n  "  our  -  SNOW  •  •  "  - 1  '  ^  •  ^  ^  ^  •-^-<\t  o— —"  \  \ \  viiifcl* spectra* niiis tr*m 0. 4D ta 0.70  .50  1  0.40  0.50 WAVELENGTH  ' \  \  V  V'  0.60 Um)  Figure 2.25: Reflectance as a function of wavelength, from [Mellor 77]  The crystal makeup of the snowpack affects the uniformity of the reflected light, especially the presence of surface hoar. These form large, flat, oriented crystals, giving the overall surface an anisotropic appearance.  Chapter 2: Snow  2.6.2  ;  ^  Extinction  Extinction measures how easily transmitted light can pass through a medium - in general attenuating exponentially according to the distance into the snowpack. Extinction is measured with an extinction coefficient, where a larger coefficient implies that light is attenuated faster. This is useful to us if we wish to render scenes where light shines through the snow and reaches the viewpoint: for example, a scene where the sun sets behind some snow-dusted trees. Extinction is wavelength dependent (Figure 2.26), and depends on a number of factors, including density and grain size (Figure 2.27).  DENSITY,  Mgm  -3  o.r  Figure 2.26: Extinction coefficient as a function of density for fine-grained snow, from [Mellor 77]. The parameter is wavelength.  34  Chapter 2: Snow  Clearice: data summarized by llobbs. 1974 tKuliiin. 1935 :.Saubercr: 1950) Bubbly ice: data summarized by Hobhs. 1974 (Kaiitin. 1935: Souherer. 1950;Pisi-  (3D  akora. 1947: PoM. 1950: Lyons and Sttiiber. 1959) Deposited wow of various types: data summarized by Mellor. 1964 Fresh Jine-grauied snow: Mellor. 1966b Blowing snow: calculatedfromdata of Budd. Dingle and Radok, 1965: by Melhr. 1966a Various types of falling snow and gratipel: O'Brien. 1969. 1910 Falling snow: Meltor. 1966a  10  r 8 o  101  io"* io ip Moss Concentration or Snow Density I M j n i )  10  1  Figure 2.27: [Mellor 77]  Extinction coefficient plotted against density for a wide range of snow types, from  Chapter 3 Previous Work  A tangled thicket of real grass.  Despite snow's everyday existence in many parts of the world, there has been little or no previous research towards a comprehensive model of snow for computer graphics. The few papers that do mention snow are generally concerned with falling snow, and usually only then en route to a description of some other technique. We should take a moment to underline the important distinction between the problem of falling snow and the problem oi fallen snow; falling snow consists of snowflakes travelling through the air towards the ground, and in any given image number on the order of hundreds or thousands. Because snowflakes can usually just be drawn as points, the interest here is in simulating a natural-looking falling motion under the influence of wind. In contrast, fallen snow is the accumulation of falling snow over hours, days, or even an entire season.  Fallen snow (otherwise known as the snow pack) consists of  countless ice particles that continually change and grow. The appearance, properties and simulation of the snow pack are many times more complex and complicated than the simulation of component snowflakes. In this section, we describe the few papers that mention snow, as well as several others that have some obvious or potential connection with the problem at hand.  36  Chapter 3: Previous Work  3.1  Geospecific Rendering of Alpine Terrain  Premoze et al. [Premoze 99] generate realistic mountainous terrains that are likely the most convincing snow-covered scenes so far. Starting with a digital elevation model enhanced with an aerial photo, they use a detailed model of snow pack evolution to add zero-thickness (ie, textured) patches of seasonal snow cover. A fair amount of work in the paper goes into preconditioning the aerial photos: normalising for sun position, removing shadows from the initial images, and classifying pixels in the image into several types of terrain cover. Once corrected, the terrain is divided into a height field of snow-water equivalent values, influenced by snow accumulation and subsequent melting. The user starts by specifying the precipitation amount and density at some base elevation, and providing an interpolating function for use at other elevations.  The base elevation is also given an  average, minimum, and maximum ambient air temperature. Temperature decreases with elevation, and changes over the simulated day. When the air temperature of a particular grid point drops below zero, the grid gains snow water equivalent based on the precipitation amount and density. The authors then compute a melt rate for each grid point based on snow hydrology concepts. The melt rate is dependent on daily ambient temperature, the snow albedo, the incoming solar radiation, and the amount of light blocked by any covering vegetation.  2D vegetation positions were obtaining  during the image classification step, while temperature and solar radiation can be computed from the simulated time and sun position. The albedo decreases as the snow gets older, rising to 0.8 for (new snow), or dropping to 0.4 for rain. Melted water disappears from the simulation. With appropriate initial environmental data, the authors can simulate an alpine area as it passes through the seasons.  Figure 3.1 shows an example of the initial undistorted scene, while Figure 3.2  shows the same scene with snow appropriate for early spring. This work produces very pleasing results, and the underlying model has many nice properties. Especially applicable are the solar influence of melting, allowing snow to remain in deep shadowed gullies long after more exposed areas have disappeared. Unfortunately, this paper generates snow on a significantly larger scale than our main area of interest. Initial satellite data comes at approximate 30m sample spacing, is mapped onto a height field with zerothickness textured snow that varies considerably by elevation. We wish to create scenes with "thick" snow accumulation on models that cannot be represented as a simple grid. Our scenes will generally be  Chapter 3: Previous Work  38  too small to have much variation by elevation. However, the model of snow evolution appears useful and applicable to our work.  3.2  Rendering Snow With Metaballs  Nishita et al. [Nishita 97] introduce a method of snow pack rendering based on multiple scattering of light within the snow volume. It is primarily concerned with snow rendering rather than modelling, and it is largely based upon previous work with clouds [Nishita 96]. However, it does produce some images at the scale we are most interested in. The general idea is to model the set of light influences on snow, including directional sunlight and incoming scattered diffuse atmospheric light as shown in Figure 3.3. Central to the approach is the notion that snow's appearance is greatly affected by the multiple scattering of light within the volume of snow. direct sunlight  Figure 3.3: Optical paths for the shading of snow, from [Nishita 97] The authors start by creating a snow surface consisting of various sized metaballs, each representing a given density function. Metaballs are placed by hand such that the centre of each metaball is flush with the underlying surface, in order to simulate the effect of increasing density with increasing depth. As light passes through grains, it exhibits strong forward scattering (i.e, is reflected within a certain set of angles in the previous direction of travel). The authors claim that the range of forward scattering angles is quite narrow. Ideally as a ray-traced view ray enters the surface, one wants to integrate the scattered  Chapter 3: Previous Work  39  light along the viewing ray, including in reflected light from the underlying object (see Figure 3.3). In order to do this efficiently, the authors divide the space up into voxels, and create a multiplescattering reference pattern that tries to model the effects on a single voxel (Figure 3.4)  Figure 3 . 4 : Reference pattern, from [Nishita 97]  The reference pattern solves a number of equations, determining the influence of nearby voxels based upon nearness, number of scatters, amount of angular scatter, etc. Because voxels are aligned with the viewing direction (as shown in Figure 3.5), the pre-computed reference pattern can then be used at multiple voxel locations, scaled by snow density.  The authors scatter small prisms and sub-metaballs throughout the snow volume in the form of a 3 D  C 'ltapt< r S: Pre vious Work  (0  texture. Each sub-ball that a viewing ray hits adds to the snow density of the particular voxel, while each hit prism acts like a tiny specular reflector (in an attempt to provide specular snow glints). Figure 3.6 shows a comparison of the various levels of rendering detail, with a noticeable difference between the opaque iso-surface (a) and the full case (d).  Figure 3.6: (a) opaque (b) single scattering with sub-balls (c) multiple scattering w/o sub-balls (d) multiple-scattering with sub-balls. All from [Nishita 97] The full image is shown in Figure 3.7.  Figure 3.7: Snow-covered car, from [Nishita 97] Because the focus is on rendering, there is little or no effort spent on creating a stable or reasonably shaped snow pack. Clearly adding metaballs by hand is not feasible for larger or more complex models, and is not appropriate for handling instabilities. Even if metaballs were added through some automatic  Chapter 3: Previous Work  41  snowing process, we suspect that the computational cost would quickly become extreme: in order to create a sufficiently detailed snow- covered environment, one needs many, many small metaballs. Figure 3.7 contained 1088 snow metaballs, and took 28 minutes to compute on an IRIS Indigo2; typical scenes could easily have 10 times as much snow detail. No matter how realistically the actual snow surface is rendered, computational limits on the number of metaballs causes distractions such as the snow on the car's tires and mirror. Secondly, there is some evidence in snow science (see Section 2.6) that suggests that the reflection properties of snow are influenced only within a very thin layer of the snow pack (10-20mm), implying that a voxel- oriented process should proceed at an extremely fine resolution near the surface only. In the car example, the authors used a voxellation of 110 for the entire world. 3  3.2.1  O t h e r Snow R e n d e r i n g W o r k  Additional work that touches indirectly on rendering the appearance of snow includes Hanrahan and Krueger [Hanrahan 93], who developed an algorithm that modelled the appearance of subsurface scattering. Although the authors mention its applicability to snow, they provide no further details or examples. Krueger [Krueger 88] described an algorithm to generate flashing and sparkles on natural objects. The authors mention that the algorithm is applicable to snow, but provide no details or examples.  3.3  Particle Systems  Snow is naturally composed of a multitude of ice grains, so it seems reasonable to try and model snowfall by creating and simulating a particle system [Reeves 83]. This paradigm assigns each particle a number of physical or pseudo-physical properties, and updates them over time according to pre-set physical or pseudo-physical rules. Particle systems have previously been used to model several other types of natural phenomena, such as fire [Sims 90], trees and grass [Reeves 85], and water [Sims 90].  3.3.1  P a r t i c l e A n i m a t i o n and R e n d e r i n g U s i n g P a r a l l e l C o m p u t a t i o n  Unfortunately, particles systems are notorious for the amount of computational power required to track thousands and thousands of individual elements. Previous particle models of falling snow are no exception. Sims [Sims 90] used a Connection Machine (4K to 64K parallel SIMD processors) to create and model a multi-purpose particle simulator. Among other effects, this was used to display falling snow  Chapter 3: Previous Work  42  under the influence of wind. An initial group of particles originated off camera, and travelled downwards at a certain preset terminal velocity. The particles were under the influence of several different motion constraints, such as vortexes and spirals, that were hand-controlled until the desired swirling effects were achieved.  Sims' particle model allowed for particle "bounce" and friction, where a zero bounce and  high friction represented flakes immediately sticking to the ground. Once snow landed, it was duplicated for visual effect, increasing the apparent amount of mass currently on the ground ( Figure 3.8). Sims shows snowfall falling onto a single horizontal surface, and does not show accumulation beyond a layer of zero thickness. This work shows that a brute-force direct approach may not be sufficient, especially when parallel computation is not an option.  Figure 3.8: Parallel computation of falling snow, from [Sims 90]  3.3.2  Stochastic Motion  Shinya and Fournier [Shinya 92] introduce more detailed and more realistic rules for motion due to wind. Their paradigm applies to several different types of graphical models, including a particle system they used to simulate falling snow. This model improves on the approach described in [Sims 90] by removing the need for "wind choreography", but still does not model snow pack accumulation. Figure 3.9 shows a scene involving snow.  Chapter 3: Previous Work  43  Figure 3.9: Falling snow, from [Shinya 92]  3.3.3  Flow and Changes in Appearance  Dorsey, Pedersen and Hanrahan [Dorsey 96] used particle systems to model the flow of water over objects in order to simulate weathering via washing and staining. Although the authors describe an absorption and deposition model for dirt-carrying water particles, the most relevant work is the initial exposure of the surface to the water source. In outside models, this is usually rain. Each water droplet is represented as a particle with mass, position, and velocity properties, created by a "rain source" that emits particles towards the model geometry. As particles travel, they are carried along in the direction of the prevailing wind, which is occasionally perturbed.  When particles eventually intersect with the surface, they are  placed in an "exposure map" (represented as a texture on the surface), and removed from the simulation. The exposure map's purpose is to arrive at an approximate distribution for the incoming water, which is in turn used as a probability function for the creation of new water particles on the surface. These particles trickle over the surface, carrying and depositing dirt, influenced by gravity, surface friction, inter-particle repulsion and surface diffusion. The force acting on each particle is projected onto the tangent plane of the underlying surface so that particles move only in the plane of the geometry. When the surface normal changes abruptly, a probability distribution determines if the particle continues on,  Chapter  3:  Previous  Work  44  or bounces away, where bouncing allows the authors to'model the splash effect of falling water. Particles fall off a surface if the difference between the velocity vector and the surface tangent exceeds a certain preset angle. Once a particle has left the surface, it falls vertically until it hits another surface. The authors state that "since computing such intersections can be computationally intensive, we precompute a table of positions where particles will land when they fall off a surface.", although they do not provide more details. The idea of the exposure map is useful for snowfall, in that it could give an initial probability distribution for potential snowfall accumulation. However, as it stands there are a number of problems with this method.  First, many particles are needed to compute a sufficient probability distribution.  There is some leeway when computing rain, because an improper distribution of water only results in slightly more staining (as the raindrops quickly trickle away). W i t h snow, an improper distribution leads to unusual accumulations in places with insufficient sampling, aggravated by the longevity of the snow pack (i.e. mistakes lie there for a while). Furthermore, as snow accumulates, some areas become fully or totally blocked (changing the exposure), which is not easily reflected with a one-time a priori exposure map. Secondly, even with a very dense exposure map, any snow accumulation based on probability must correctly handle conservation of snow mass. That is, the entire mass of snow that falls from the sky must land somewhere - again not so critical for rain, because any missing mass can be assumed to have trickled away, soaked into the ground, etc. There are also restrictions on the amount of mass that can originate from any particular area of the sky. Even with these restrictions, the idea of pre-computing an exposure distribution is potentially quite useful.  3.3.4 Hsu  Dust Accumulation  and Wong [Hsu 95] simulate the accumulation of a zero-thickness layer of dust on model surfaces.  The amount of dust that accumulates on a model is initially given by a dust prediction function, based upon surface normal, direction of dust source and surface stickiness. The prediction function is then multiplied by an external factor constant that attempts to model scraping, wind and coverage. The final result is rendered as a texture over the original model. The authors compute coverage for point P on the surface by firing rays randomly distributed over the upper hemisphere around P. The actual coverage is based on based on the weighted average of the distances to found intersections.  Chapter  3:  Previous  Work  45  This work has some parallels to potential snow accumulation, including firing rays to determine sky exposure. However, it is missing notions of dust mass conservation, stability and accumulation thickness. Still, it forms a useful basis for very thin layers of snow accumulation, such as the flake dusting described in Section 7.7.1.  3.3.5  U s i n g Particle Systems  With any appreciable amount of snowfall, the number of potential particles quickly becomes unmanageable. It is clearly not feasible to model the migration and metamorphism of individual snowflakes over a season, without some deviation from the straight particle paradigm, even if particles represent more than a single snowflake. Even so, previous particle research can provide us with some useful ideas. First, the erratic local motion of a wind-influenced falling snowflake implies that it can fall onto surfaces that are not directly visible to the sky. Particle systems might prove useful in determining which of these surfaces actually receive snow, and how this exposure changes during subsequent accumulation. So, even though it might prove intractable to simulate every snowflake in a particle system, it might be feasible to simulate just enough of them to determine where the rest of the snowfall can be placed. Secondly, wind contributes to large-scale movement of snow mass along the ground. Modelling wind movement and turbulence with a simple particle model might give one an idea of areas of thick and thin snow distributions. Finally, particle systems may prove useful for modelling the flow of snow moving over and around underlying supporting surfaces.  3.4  Snow on the Ground  Once snow has been appropriately allocated to surfaces in our target world, we must perform subsequent checks for stability and motion. There are several possible approaches that may be successfully applied to the motion of accumulated snow pack.  3.4.1  Soil  Li and Moshell [Li 93] produced the work that is probably most visually and physically connected to thick snowpack modelling, other than Nishita et al. [Nishita 97]. The authors were responsible for a dynamic soil model, allowing for volume conservation, soil slippage, and manipulation of soil with a simulated bulldozer.  Chapter 3:  Previous  Work  46  The authors use concepts from soil science to quantify soil stability, introducing soil shear strength and shear stresses to show that failure on a soil slope will occur along a failure plane, as pictured in Figure 3.10.  Figure 3.10: The failure plane, from [Li 93] The strength force s resists motion along the failure plane, and is calculated as s = cL+Wcos(a)tan(0), where c is a measure of the soil's cohesion, L is the length of the failure plane, W is the soil weight above the failure plane, a is the angle of the slope, and 6 is the angle of internal friction (also a measure of cohesion). The shear stress force r influences motion along the failure plane, is primarily caused by gravity, and is calculated as r = Wsin(a). The factor of safety is defined as ^ is the measure of stability: any ratio less than 1 implies that the soil is unstable. To determine if a particular 2D soil configuration is stable, the soil is divided into regular containers, where each container holds an average of  , as shown in Figure 3.11.  Slice l+i  Slice i  Slice i-i  hi  i  yi-i 00>},J>t  v  ^i'Vt*^ * **  '^.'•J***  *l* . V *" 1  Ax  Figure 3.11: Considering slices as containers, from [Li 93] The critical area is the top triangle of each slice. The authors attempt to solve for each slice's a, which indicates a failure plane exists. If such a failure plane exists, they solve to determine a stable slope configuration, which is the critical slope angle where the factor of safety is 1. The authors then use an equation from soil mechanics to approximate the force on the top triangle  Chapter 3: Previous Work  47  of each slice. After some rearrangement to include in volume conservation, the authors are left with n+l ordinary differential equations for n slices. Solving these equations (or a linearised version) produces a set of heights j/,- for a given time step. To extend to 3D, the authors tessellate the ground surface into a number of prisms, each surrounding an elevation post. The forces on each post are approximated by combining the 2D solutions in multiple plane directions. In their examples, the authors used a 90 x 90 grid. This paper makes some important contributions. It solves a global model of soil stability, using a fair approximation of the forces involved. It also produces some visually pleasing pictures, as shown in Figure 3.12.  Figure 3.12: A scoop loader is dumping, from [Li 93] However, there are some important differences that make this work difficult to directly extend to snow. Li and Moshell model the soil environment as a regular grid of elevation posts, which will not be possible on almost any real model of a snow environment, where surfaces will be irregularly shaped, sized, and located. Of special note are the boundary cases, where snow travels across surfaces of different orientation, possibly leaving the surface entirely. For a given surface, volume is not necessarily constant. The results of some previous computation may cause snow to fall off an edge and settle upon an underlying surface. Depending on the height and velocity of the initial snow slippage, the destination surfaces(s) may be some distance away. The number of elevation posts is also not necessarily constant. Assume snow slides off a roof onto a  Chapter  3:  Previous  48  Work  large open area. Unless the open area already had sufficient elevation posts, we may need to add more to model the lump of incoming snow. It is important to represent incoming snow, because target locations are strong candidates for subsequent sliding. Because we do not know a priori where snow may slide, the only way to avoid adding posts during the computation is to perform an extremely fine (and expensive) gridding.  3.4.2  A n i m a t i n g Sand, M u d and Snow  Sumner et al.  [Sumner 98] describe a general model of a deformable ground material that can be  compressed by rigid objects to form tracks, footprints, and marks.  The approach centers around a  regular height field, where each grid element supports a column of material. T h e user chooses the grid size and distributes mass to form the initial base surface. The authors then use a dynamic simulation to compute the position and orientation of rigid objects (for example, a foot) that push into the ground surface. A t each time step, each column performs a raycast to determine if, and where, it has intersected with an incoming rigid object. Intersected height fields are reduced until there is no interpenetration; part of the material is compressed using a compression constant, while the remaining mass is displaced outwards to adjacent columns, based upon the distance a given column is to an unaffected column. This generally has the effect of moving compressed material to a ring immediately outside the area of influence of the rigid object. T o reduce the steep appearance of this ring, an erosion step moves material based upon an angle of repose with nearby neighbours. This repeats until all adjacent grid neighbours are until a particular slope. The authors provide several interesting animations, and include some good comparisons with actual tracks in mud, and snow. Unfortunately, although the regular grid is refinable, it is not general enough to handle a wide range of models, leading to large numbers of grid elements (up to 2 million in one example). However, the compression and erosion aspects are quite applicable to snow underfoot.  3.4.3  E r o d e d Terrains  Musgrave et al. [Musgrave 89] describe fractal terrain generation, including an erosion and thermal weathering model that is quite applicable to snow stability. Terrain is modelled as a height field. Rain is deposited on vertices of the height field, sweeping up erodible material and moving it downhill. The  ChapterS:  Previous  Work  19  amount of material moved is based on water's sediment carrying capacity, the rate sediment leaves moving water to settle at a vertex, and soil softness constant that controls the conversion of stable soil into erodible sediment. Thermal weathering attempts to simulate the overall movement of (non-water influenced) talus-like material from steep slopes down to gentler location. At each time step, the authors compare the angle between neighbouring vertices with a talus angle of repose. If the angle of repose is too steep, a fraction of the difference is moved down onto the neighbour. With sufficient repetitions, excess material will be moved to stable locations. The combination of erosion and thermal weathering serve to etch water-filled river beds, and weather steep surfaces into more gentle inclines. Figure 3.13 and 3.14 shows an example of before and after erosion simulation.  Figure 3.13: A synthetic fractal terrain patch before erosion, from [Musgrave 89]  Figure 3.14: A synthetic fractal terrain patch after erosion, from [Musgrave 89] Although not directly concerned with snow, the authors show that a simple relaxation technique can produce reasonable results with similar granular materials. This method was developed for use  Chapter 3:  Previous  50  Work  on a height field, and does not immediately support interpenetration, obstacles, stationary water, or transit over cliffs and drops. However, relaxation is a very general and applicable technique, and is a good candidate for a stability algorithm. The terrain models of [Li 93] and [Musgrave 89] seem quite appropriate to snow, with sufficient extensions for non-uniform meshing, obstacles, and transport over edges.  3.4.4  '  Fluid Dynamics  Another possible approach would be to treat accumulated snow like a viscous fluid, where the viscosity (stiffness) depends on temperature and moistness (i.e. stiff water). Thus, it might be appropriate to use concepts from fluid dynamics to model the stability and motion of the snow on the ground. Fluid dynamics have been modelled extensively in other disciplines, but there have also been a few papers concerning computer graphics applications. Stora et al. [Stora 99] animate lava flow with a particle model influenced by forces determined from a model of temperature, viscosity, and heat. Kass and Miller [Kass 90] describe water as a regular grid of height fields, implement some simplified equations for shallow water motion, and then build an approximation of the discrete sampling of the continuous, partial- differential equation that must be solved. The algorithm claims to be quite stable, real time (for a 30 x 30 grid) and linear in the number of height field samples. Of particular note is an example where the authors show drops of rain falling, landing and trickling down surfaces to accumulate into puddles. Unfortunately, fluid approaches have not demonstrated their ability to handle outdoor models of the complexity and size we desire.  3.4.5  Deformable Surfaces  Terzopoulos and Fleischer [Terzopoulos 88] modelled inelastic deformations incorporating the idea of viscoelasticity, plasticity and fracture. An inelastic deformation will not return to the original reference shape, and will usually occur given sufficient force or temperature on an object. Viscoelasticity characterises creep, where an object deforms over time under constant force. In addition, the rate of deformation depends on the duration of the disturbing force, unlike elastic deformation where deformation is proportional only to force. The authors give silly putty as an example of  Chapter 3:  51  Previous Work  viscoelastic behaviour; bouncing it provides a large, sudden disturbing force that does not cause much deformation, while a slow gentle stretching causes elongation. Characterising plasticity allows materials to deform elastically up to a point, beyond which the rate of deformation may change. Finally, fracture describes a limiting stress level at which the object cannot deform any more, and breaks. The authors implement this by modelling the object as a grid of particles connected with by combinations of the different deformation types, and then solving over time using large systems of simultaneous ordinary differential equations.  This paradigm is continued by Terzopoulos et al.  allowing deformation properties with temperature dependence.  [Terzopoulos 89],  This allows object melting and heat  conduction. Piatt and Barr [Piatt 88] introduce a model for deformable solids, where objects are represented as a finite element model, subject to constraining functions (such as collisions) and the minimisation of goal functions. The solution requires solving large numbers of differential equations. Deformable surfaces and fluid flow models are possibly applicable to the transformation of snow on the ground, and a close simulation of the physical properties of snow may ultimately prove to produce most accurate results. Unfortunately, our particular problem domain contains many complicating factors that reduce the desirability and ease of use of these algorithms. First, we must perform calculations on a very large number of irregular surfaces,. likely without the help of a regular grid. Within each surface, snow mass continually leaves over edges and arrives from the sky and surfaces even higher up. Secondly, we expect the number of interesting areas to be very large, and we will likely need either an extremely fast and simple algorithm, or hierarchical control over stability detail, or both. Of particular concern are the general instability and efficiency of numerical equations with so many variables and the number of particles required to properly model a complex object. In [Terzopoulos 88], the most complex example was on a 180 x 127 mesh, while [Terzopoulos 89] simulated melting on a 5 x 5 x 10 grid. A n d finally, for reasons described in Section 4.3.3, the stability of snow is not as important to us as the initial distribution, and thus must be allocated a proportionally smaller amount of computational work.  3.4.6  Oriented Surface P a r t i c l e s  The previously described models allow for long range attractive forces (grouping particles into objects) and short range repulsive forces (controlling inter-particle spacing). Szeliski and Tonnesen [Szeliski 92] modify this paradigm by adding an orientation to each particle, such that particles tend to align themselves into locally spherical or planar surfaces rather than volumes.  Chapter 3:  52  Previous Work  This allows surfaces to be created from fairly sparse data, by growing particles to fill in gaps within the previously determined local plane. If surfaces, are stretched or bent, the same method can be used to provide a smooth interpolation across the deformed area. Surfaces can also be easily joined, merely by bringing the appropriate sets of particles close enough to interact. All of these properties seem quite applicable to snow. One could generate seed particles through a wind influenced snowfall-decent algorithm, and then join up the target locations to form a continuous surface. The fact that surfaces need not be directly connected to supporting surfaces might easily allow for snow bridges and overhangs. However, in order to describe snowfall accumulation, we would likely have to use several types of surfaces (i.e. top and bottom) in order to represent volume accumulation.  3.4.7  Multi-scale Granular M o d e l  Luciani et al. [Luciani 95] introduce a multi-scale physical model for granular materials designed to simulate such granular phenomena as piling, arching, and avalanching." The authors create a set of so-called spherical "intermediate-scale" grains, that can interact with the ground and with other particles to produce deformed aggregates. A particle's new shape depends upon a collective interaction with its neighbours, leading to complex deformations. Individual grain interactions are computed using non-linear equations to simulate a visco-elastic collisions, allowing squashing, clumping, collisions, and interlocking (Figure 3.15), and the formation of piles and arches. Using nonlinear methods allows for multiple collisions between grains, including contact oscillation and damping. This compares with an impulse (direct bounce) model, where a collision occurs at a single point in zero time. To model interactions with the ground, the authors cover the base surface with small, fixed intermediate-scale particles that provide surface friction, depending upon radius and spacing.  o  o  F i g u r e 3.15: One point grain and deformations, from [Luciani 95]  Chapter 3:  Previous  Work  53  In order to allow a finer scale, the authors refine each physically modelled intermediate scale particle with a cluster of much smaller particles that are constrained by the higher level solution. Special care must be taken to correctly handle intermediate particle motion so that a discrete stream of flowing large grains can be replaced with a continuous stream of tiny grains (as opposed to discrete groups of tiny grains). To model the small-scale particles, the authors create a grid of phyxels (physical elements) where each phyxel is represented as mass point. Interactions are modelled with a linear viscoelastic force equation, connected to each of the four phyxel neighbours and the ground. Each intermediate-scale particle affects a number of phyxels, creating a deformation field that may move as the intermediate-scale particle moves. To draw small-scale particles, the effect of several distortion field phyxels are blended and mapped to a color scheme. The attraction between phyxel elements and intermediate scale grains may be controlled and thresholded so as to influence the small scale distortion shape of the larger particle's motion.  One of the  important contributions of this work is the idea that physical properties can be modelled at a number of different scales, using a different solution method with a different discretization steps in time and space. Combining the results avoids having to work at the finest level of detail across all properties one wishes to model. In this case, intermediate-scale particles are tracked using non-linear interactions, with a minimum element of a single particle. They are then replaced with a linear model of smaller particles, where a minimum element is a spatial location (Figure 3.16). This research has some potential usefulness towards a solution for snow. Clearly, it seems advantageous to model interactions at various scales using an appropriate level of physical detail. However, to apply this to snow, this paper must be extended to 3D. The authors claim that "..the achievement of 3D intermediate-scale simulations is straightforward with no additional cost. The refinement model can also be extended to 3D". Yet, they provide no proof, examples, or algorithms. As well, we expect a multi-scale situation opposite from the one described. The authors simulate the dynamical motion of large chunks of particles, and then replace each with groups of fine particles. With snow, we will likely need to replace many fine particles with larger groups that then can be checked for stability. Finally, computing interparticle interaction forces a fairly significant computational cost, even if only performed at the scale of intermediate grains. T h e authors found that a 900 grain 2D simulation required 6.6 seconds per step, running on a SGI Indigo Extreme. Even at a very coarse resolution, a 3D  Chapter  3:  Previous  54  Work  («o ««o ««o  Figure 3.16: Correspondence between intermediate scale non-linear simulation and the final small-scale simulation, from [Luciani 95]  simulation will have hundreds or even thousands of times more particles.  3.5  A d Hoc Methods  Perhaps one of the most famous examples of computer generated snow are the popular commercials showing polar bears drinking Coca-Cola [Shifflett 97].  Although these show scenes from a flat snow  covered Arctic background, they were created by hand generating and texturing a deformed mesh, without following any underlying model.  Chapter 4 Assumptions, Goals and Problems  A real tire swing in the snow.  Chapter 2 introduced some of the main properties of snow, and described the central mechanisms of snow metamorphosis, hopefully giving the reader an appreciation of the complexity and inter-dependence of the many variables affecting a snow layer. We suspect that snow is far too complex a material to simulate completely, especially given the lack of previous research. Thus, in this Chapter we set out some of the problems we must solve, and then make a number of restricting assumptions and simplifications in order to make the overall task more tractable. This essentially "sets the stage" for the later introduction of our chosen methods for snow simulation. We start by making a number of assumptions about the snow pack, the environment, and the size of the problem, in an attempt to reduce the overall complexity of the task. Next, we discuss the visual goals of our snow modelling system. The algorithm described in Chapters  55  6, 7, 8 and 9 should  56  Chapter /,: Assumptions. Goals and Problems  attempt to produce these visual effects first, and other effects later, or as a side effect. These are the essential characteristics that say "snow" to the viewer (or at least the subset that we have chosen to be important). We then continue on to discuss some of the problems our chosen snow distribution, stability, and rendering algorithms must overcome, setting out some of the criteria we must address in order to eventually arrive at an approach.  4.1  Assumptions and Restrictions  This section describes some of the assumptions and restrictions we have imposed on our model, primarily out of the need for implementation and algorithmic simplicity when dealing with huge, outdoor scene. T o briefly summarise, we restrict the problem to a static scene,-with homogeneous snowfall falling in a local scale and accumulating into a single layer. We avoid slab avalanching, and both atmospheric and ground evolution of the snowpack.  4.1.1  Scale of Snowfall  We restrict our snowfall model to simulating sizes where the visual depth of the snowpack is appreciable relative to the size of the surroundings. This implies the viewer is close to the objects, within a- "real world" distance measured in meters or tens of meters. This is a rather vague definition, but it is meant to avoid large world scales  1  where snow can be  modelled solely by covering the scene with a white texture. We are primarily interested in scenes where the depth of the snow is of clear visual importance, and cannot be ignored. This does not prevent us from using a larger scale if we wish; if we choose mountains as our environment, our model will treat them on the scale of houses (i.e., on a scale too-small for weather variation).  4.1.2  Snowfall L o c a l i t y  We assume that snowfall is homogeneous across the scale of our world model. That is, we assume that snow falls on all parts of our model world with the same properties and duration, etc, no matter what the size of the our target world. We are specifically concerned with giving snow properties associated 'We use "scale" in the sense of "general size", and not in the sense commonly associated with maps (ie, 1:100 scale).  Chapter  4: Assumptions,  Goals and Problems  57  with a single, uninterrupted snow storm. We assume snow lands on previously bare surfaces and the duration of the storm is not long enough for any significant changes in the accumulated snowpack. The homogeneous snowfall restriction does not prevent an uneven distribution of settled snow due to slopes and wind.  It merely states that all snow enters the model world in the same way.  In the  real world, snowfall will vary depending on local atmospheric altitude, conditions and the terrain. For example, more snow often falls on the windward side of a mountain range, and the tops of mountains receive more snow than in the valley floor. Our model produces a constant snowfall no matter what the underlying terrain, fitting with our previous assumption that we are modelling snow on the scale of a human standing on the ground. Practically, these restrictions can easily be relaxed (see Section 11.1) to include a much wider range of parameters and sky-variant snow, but we impose this loose restriction now to highlight that we are not concerned with simulating an appropriate weather model as well.  4.1.3 For  Ground Evolution of Snow  the purposes of this research, we assume that snow does not undergo evolution on the ground,  essentially restricting us to simulating freshly fallen dry snow. This ignores melting, settling, and other types of snow metamorphosis. We impose this restriction primarily for simplicity; otherwise the scope of the problem becomes unmanageable.  4.1.4  Drawing Avalanches  When our model indicates that snow is unstable, we need to adjust the model until it reaches a state of stability. In the real world, this means that snow avalanches. This immediately gives rise to the problems: "What does an avalanche look like?", and "What properties should we assign avalanched snow deposited on top of fallen snow?"  Chapter  4: Assumptions.  Goals and Problems  58  Figure 4.1: Avalanches-in-progress are outside the scope of this research. This image shows a small powder avalanche.  Figure 4.2: Avalanche debris is outside the scope of this research.  Chapter  J,: Assumptions,  59  Goals and Problems  During stability calculations on fresh-falling snow, we make the assumption that avalanched snow has the same properties as the underlying snow. This is a fair assumption: for a single layer of newfalling snow, any avalanches that occur will consist of a single type of snow, as yet unchanged from any other fallen snow. The avalanches that occur will be dry snow sloughs, and not large slab movement. This allows us to assign avalanched snow to the new location as if it had fallen from the sky in an unusually large clump. We assume that any surface disturbance is small scale (because the moving snow is unconsolidated) and quickly covered up by the next layer of snow. These assumptions do not hold as soon as the snow properties change or if the snow contains more than one layer. In reality, as fallen snow consolidates, its properties may change across a scene. This means that avalanching snow will have different properties from snow at its final destination. Even the very process of avalanching changes the properties of snow, causing large chunks and balls that leave trails, indentations, and look much different from the smooth, underlying snow layer (Figure 4.2). Things get even more complex when slabs of snow (containing multiple snow layers) avalanche on top of a pre-existing set of layers.  This essentially requires simulating a wave made up of layers of  non-homogeneous material. For these reasons we make the critical assumption that our model only has a single layer of snow. We expect that using only a single layer will result in no apparent visual errors in locations where no avalanches have occurred.  4.1.5  Use of Commercial Software  This research is primarily concerned with modelling the appearance of snow surfaces. To that end, the natural goal is to produce a set of pictures that can be viewed and checked for errors.  In order to  save effort, we have decided that we will use commercially available software to perform as many of the modelling and rendering tasks as possible. In our case, we have chosen Alias|Wavefront as our primary support tool, although it is important to note that almost all common modelling and rendering packages are equally usable. In Figure 4.3, we give an overview of our design pipeline. We start by creating models and scenes using the construction tools available in Alias|Wavefront.  By using a commercial package, we ensure  that there exists a large database of precreated scenes to choose from, greatly minimising the amount of time it takes to create realistic looking test scenarios. If we created our own modelling tools, they would not be as flexible or easy to use, nor would there by any pre-existing models. Once we have created our scene, we export it to Wavefront ".obj" format [Murray 94]. The ".obj"  Chapter J : Assumptions, t  Coals and Problems  Commercial modelling  1,1)  Commercial modelling  Figure 4.3: Overview of the use of commercial software format is easy to parse and convert, and is supported by almost all major commercial modelling packages, allowing us to generate snow independent of the platform, vendor, and format of the model. The original scene contains a great deal of information we do not need to export, such as textures and lighting. Armed with the scene geometry, we simulate snowfall according to our snow model. This is the primary goal of our research. We then create a set of surfaces that model the added snow, and export these surfaces into Wavefront ".obj" format, along with some additional information about the rendering properties of the snow. The scene database we used to compute the snowfall is only partially complete, so we discard it. Back in the commercial software, we load the new snow surfaces into the complete scene database, creating and rendering a new scene with all the required properties. We argue that snow looks most realistic in a realistic setting, which is most easily achieved using all the lights, textures, and effects already available. Using commercial software as the primary rendering engine also makes this research immediately applicable to the general public. For example, we might be able to add snow to models sent to us over the WWW, or even create a complete plug-in that operates seamlessly with existing commercial software. 4.1.6  Static W o r l d  Our current assumption implies that the world geometry is static. As nice as it would be, we do not allow for drooping tree branches, for example. We also assume that our primary goal is to render a single, static image representing the end of the snow storm. We specifically do not plan to create animations of gradual snow accumulation, primarily because rendering avalanches-in-progress is out of our desired scope of interest. For example, consider a sequence of stability slices. During any single stability check,  Chapter  4-  Assumptions,  Goals and  Problems  61  there may be a large-scale movement of snow, including the possibility of snow travelling over any edge. Real snow movement produces tracks, snow clouds, and snow chunks, all of which are complex problems in their own right, and none of which we plan on rendering. This means that large bumps of snow will mysteriously move across snow surfaces, possibly disappearing over an edge and magically reappearing lower down. Of course, some environments will be stable enough that small amount of this can happen without too much visual disturbance. Note that although we do not plan to animate the accumulation of snow, we do wish to fly through and around the final scene to allow for the maximum use in animation and visualization.  4.1.7  Object  Restrictions  Using Alias|Wavefront does, not particularly restrict the set of available input surfaces, because this system supports many major different surface representations.  However, we will assume that for the  purposes of our algorithms, all surfaces are triangles. We make this assumption primarily for simplicity, and because it is easy to implement.  During the model export phase, Alias|Wavefront allows several  different levels of triangulation, so that non-polygonal models can be easily transformed to meet our assumptions.  4.2  Goals of the Snow Model  It is not our purpose to transport an entire model of snow theory from hydrology into computer graphics. First, snow is one of the most complex materials around, at least for the purposes of characterising engineering properties. A given type of snow can change over time, weather, and temperature, leading to subsequent changes in behaviour. There can often be large differences between individual snow layers, and even within a single layer. Secondly, to the author's knowledge, no comprehensive model of snow actually exists to be copied or adapted. And finally, we're not in the business of environmental snow simulation - our goal is instead to produce pictures (admittedly, we hope that by simulating what really happens, we can do a better job). So instead, we will list a number of the visual characteristics of snow that we want our  results  to have. These are the cues that say "snow". It should be mentioned here that  snow has hundreds of visual characteristics, and the set chosen here is based largely on the preference and beliefs of the author. T h e reader may judge if the set of visual goals chosen is sufficient to express  62  Chapter J : Assumptions, Goals and Problems t  most snow-covered scenes. For each of these goal, we make reference to a later section that describes the algorithm chosen as a solution. As well, Chapter 10 includes a table that cross references our goal against result images containing that property. In a later section, we talk about the goals of our implementation and algorithms (ie, speed, adaptability, etc).  4.2.1  Snow Location  There are a number of properties associated with an initial snowfall, without stability calculations (implying that all slopes are too gentle to cause snow movement). These restrictions are fairly general, and allow us many simulation options. • Initial snow location. Under windless conditions, snow falls on all surfaces that have an unobstructed view of the sky directly overhead. This seem obvious, but this property has been ignored before (Figure 3.7 - especially on the window and back of the car). Of importance to us is the idea that A L L eligible surfaces get snow - not just the largest or the luckiest. We cannot tolerate bare patches in the snow cover due to insufficient time or processing power. Our approach is discussed in Section 7.5. •  Partial coverage. Under windless conditions, snow also falls on surfaces that do not have a direct, unobstructed view of the sky directly overhead. Snow flutters, spins and twists on the way down, causing small deposits to accumulate in sheltered areas.  Note that as snow accumulates, the  partial coverage of a surface may change, as it becomes increasingly blocked by larger and larger snow deposits. Figure 4.5 shows an example of an area not visible to the sky (under the bush) that still receives snowfall. Our approach is discussed in Section 7.5. •  Partial snow amount. Under windless conditions, the amount of snow an obscured surface receives is related to the proximity of surfaces with unobscured snowfall. This essentially says that there needs to be some kind of snow depth transition between obscured and unobscured areas.  The  transition may be obvious, as shown in Figure 4.6, or confused, blurred and hard to detect. Our approach is discussed in Section 7.5. •  Conservation of snow mass. Under windless conditions, the amount of snow that falls on obscured surfaces is removed from visible surfaces that might have otherwise received that snow. Again, this  Figure 4.5: Areas under the bush are not directly visible to the sky yet still receive snowfall.  64  Figure 4.6: A bush leans out over a rock wall. Note the transition in depth between the unobscured right-hand part of the picture, and the obscured area deep under the bush. See also how this particular drop-off is not linear.  Chapter 4'- Assumptions.  Coals and Problems  65  Figure 4.7: Conservation of mass. Note how the depth of snow at the edge of the obstacle is less than elsewhere. The missing depth ends up under the obstacle. is an obvious natural restriction, but must be reiterated for the non-physical realm of computer graphics. If the sky produces 10 cm of snow, that depth must be allocated somehow between visible and occluded surfaces: if all visible surfaces receive 10 cm, then there is none left over for partially occluded surfaces. This prevents approaches where we just give all upwards facing surface 10 cm of snow, as the total apparent mass is then greater than that produced from the sky. Figure 4.7 shows an example of mass conservation. Our approach is discussed in Section 7.6.1. • Surface bridging. As snow falls, it can create small bridges across closely positioned surfaces. These bridges close up variations in the surface, and present new incoming snow with a smoother surface to build on. As an example, consider snow falling on a sewer grate. The spaces between the grate slowly close up as snow falls, eventually disappearing altogether. (Figure 4.8) Our approach is discussed in Section 9.2. • Cornices and overhangs. As snow falls, it may lean out over the supporting surface or form cornices (Figure 4.9). Our approach is discussed in Section 9.2.2.  Chapter  \:  Assumptions,  Coals and  Problems  66  Chapter 4  :  Assumptions, Coals and Problems  67  Figure 4.10: Unstable snow moves to a stabler area. In this photo, snow has slid down the slide to pile up at the bottom.  4.2.2  Snow Stability  There are also properties associated with initial snow stability. Our overall approach is discussed in Section 6.1. • Instability. Depending on characteristics of the surface and the snow, a particular accumulation of falling snow can be classified as stable or unstable. Unstable snow has a vector of instability, that generally indicates the direction moving snow will travel. • Snow restabilization. Locally unstable snow will move in the direction of instability until global stability is reached. The final destination of snow may be one that previously had little or no snow (Figure 4.10).  4.2.3  Snow Surface Appearance  We can give some visual properties of new fallen snow. The appearance of snow should be a major consideration, but is given short shrift here because it is outside our research concentration. • Snow surfaces are predominantly composed of smooth curves. While there are places for sharp surfaces, snow covered worlds generally have few discontinuities (Figure 4.11). Our approach is  Chapter J : t  Assumptions,  Goals and  68  Problems  Figure 4.11: Snowy worlds have many smooth curves, discussed in Section 9.2. • Uneven local snow distribution. The final snow surface has the appearance of bumps and depressions larger than the scale of snowflakes, where this is independent of the surface distribution. Our approach is discussed in Section 9.3. • Tracks. The shape of the final snow surface may be modified by user-added primitives, that simulate local indentations or tracks (Figure 4.12). Our approach is discussed in Section 9.4.2.  4.2.4  Snow Under the Influence of W i n d  We can provide some properties of snow under the influence of wind. • Initial distribution. The initial distribution of snow is affected by wind, such that the the amount of snow a particular surface receives depends upon wind shadowing by intervening objects (Figure 4.13). Our approach is discussed in Section 7.5.6. • Surface redistribution. Snow redistributes itself across surfaces based upon purely local snow transport model based upon a homogeneous wind field. Moving snow may redistribute across edgos.  and to locations  that previous!)  d i d not  have  snow.  N o t e that  we specifically do not plan  Chapter  £  Assumptions,  Goals and  Problems  Chapter J : Assumptions, Goals and Problems f  70  Figure 4.14: Wind transport causes snow redistribution. on modelling non-homogeneous wind flow, or turbulence (computational wind dynamics is a large body of research in itself). Figure 4.14 shows an example where the windward side of a ridge has been swept bare of snow. Our approach is discussed in Section 8.8. Starting with Chapter 6 , we introduce a snow model that chooses a particular approach and attempts to address these specific problems.  4.3  Overview of Problems  This section describes some of the major obstacles and problems that must be overcome before arriving at a final model. This section differs from the previous in that we are specifically interested in the goals and requirements of the snow algorithms we implement, and not in the visual results they produce. Chapter 10 contains table that cross-references goals of our algorithm with our actual implementation. 4.3.1  A c c u m u l a t i o n of S n o w on O c c l u d e d Surfaces  Snowflakes do not fall straight down. Even under nearly windless conditions, they spiral, twist, and spin as they descend, allowing individual snowflakes to land upon surfaces that are not directly visible to the sky. Consider the side view of a bookshelf, shown in Figure 4.15. The very top of the bookshelf is  Chapter  J,: Assumptions,  71  Goals and Problems  exposed to the sky, and collects the majority of the snow. At the extreme front edge of the bookshelf, most of the snow mass will be caught on the top surface, but some will swirl around underneath this surface. Some of this remaining snow mass will be caught on the next shelf, and so on. However, at the extreme back edge, no surfaces collect snow because they are blocked by the back of the bookshelf. .• sky  surface sees sky directly snow may bounce or stick to surface  surface occluded from the sky still receives snow  but  surface occluded from the sky and does not receive snow surface multiply occluded from the sky still receives snow side view of object  ground  Figure 4.15: An example of varying snow accumulation due to surface blockage Thus, in order to determine if an obscured surface will receive snow, we have to determine if there exists a 3D path to the sky. We can constrain the possible path volumes somewhat by imposing rules of maximal snowflake motion, but this still a very complex task. A snow-receiving surface can be buried under hundreds, or even thousands of higher surfaces, making it much harder to deposit snow based upon object height or sky visibility. We might be able to get away without simulating coverage of occluded . surfaces, except that 2  B y "occluded", we really mean that there is no vertical straight-line path between a surface and the sky. See Figure 4.15. 2  Chapter  72  Assumptions. Goals and Problems  we expect this property to be very important when snowing on trees.  Trees consists of many small-  diameter needle, twig and branch surfaces that partially block tree surfaces lower down. If we simulated straight-line fall only, then only the very upper part of the tree would receive any snow. Our approach to handling and simulating snowflake motion is described in Section 7.5.  4.3.2  Surface Representation  In this section, we describe some of the needs of our surface representation algorithm.  Given some  accumulation of snow, and an arbitrary planar surface representing our model, we want to be able to describe the snow surface in such a way that we can provide detail sufficient enough to describe any snow shadowing due to occlusion. Our representation should provide us with several important capabilities, described below.  Level of Detail Representations  In order to conserve memory, reduce the size of our final surface, and improve simplicity, we desire a surface representation that allows us to provide various levels of detail at various places within our world model. This allows us to use the same surface primitives to describe a wide range of snow accumulation masses. For example, we wish to be able to represent large, flat areas with a low level of detail, while complex, curved surfaces might use finer and finer levels of subdivision. T h e finest level of detail should be user-definable. We would also like the ability to vary the initial surface subdivision based upon the properties of the surface, such as simulating near-vertical surfaces with only a few samples. The idea here is that snowfall on this surface is unstable, and will immediately fall off as soon as we check for stability, reducing the work we need do to at the cost of undersampling. See Section 7.3.2 for our approach to multi-scale meshing.  Proportional Accumulation  A multi-scale representation of a snowfall surface (as above) allows us to represent the results using a non-constant level of detail. It should also reduce the amount of work required to achieve these results, by an amount proportional to the size of the subdivision. This allows us to build an accumulation on a large flat roof faster than one on a complex aggregate of tree surfaces, and to specifically link algorithm  Chapter /,: Assumptions. Goals and Problems  73  speed to surface representation coarseness. See Section 7.2 for our approach.  Spatially Adaptive Representations  Given that we do work proportional to our a priori subdivision size, and that subdivision size affects the level of detail of our final surface, we also want the ability to change the subdivision size on the fly. This allows us to increase the sampling rate in places where we determine the snow accumulation function is particularly complex. For example, we may initially represent a roof with only a few surface primitives;.however as soon as the snowfall starts, we may need to redefine our representation to account for a complex pattern of snow shadowing produced by a nearby tree. As well, stability calculations may result in snow falling off surfaces and landing on objects below, possibly requiring more subdivisions to represent the deposit pattern. To prevent the proliferation of subdivisions (and thus work), we would also like to be able to remove unneeded subdivisions in areas where the surface detail and work is sufficiently similar. One obvious time for subdivision collapse is immediately after the first stability calculation: at this point we have included all subdivisions required for occlusion, and all of those added to handle the first pass of snow redistribution. We may need more subdivisions as more and more snow redistributes, but most of the initial subdivisions have been placed. If no occlusions exist, and surfaces have not been showered with redistributed snow from above, we can potentially combine similar subdivisions into a large scale representing a bigger area and less work. If we can add new subdivisions, then we have the flexibility  to increase the resolution of the snow field in a very local way. For example, we could place a  set of new samples shaped like a footprint, allowing us to make tracks in the snow. If we placed those samples before snow accumulated, we could tag them as giving off heat (for example), and modify the snow mass accumulation accordingly. This could simulate snow falling on a steaming manhole cover, or on the hood of a car warmed by the engine. We would also like the ability to provide spatial detail based upon view-specific parameters, such as reducing the sampling rate as a function of distance from the eye position. See Section 7.2 for our approach.  Chapter  J: t  Assumptions,  Goals and Problems  74  Temporally Adaptive Representations  We expect that building a snow accumulation model will be extremely expensive, no matter which surface representation we choose. Thus, we wish a great detail of control over levels of detail; such that we can both a priori and on-the-fly change the amount of work we need on each pass. To give us an even greater level of control, we would like a surface representation (and associated algorithm) that allows individual locations to be allocated an amount of work proportional to the confidence of the accumulation. Consider a complex surface with many fine subdivisions that continually receives the same amount of snow on allocation pass after allocation pass.  As time progresses, we want to be able to reduce  work spent on this object by increasing the distance between allocation steps, so that a good deal of snow is added automatically. This property basically says that as we gain confidence in our initial snow distribution, we can sample less in time (at the cost of all of the drawbacks usually associated with potential undersampling). See Section 7.2 for our approach.  Surface Smoothness  Finally, we want a surface sampling method that gives us some flexibility in connecting individual areas together to form asmooth surfaces. Our surface representation should prevent inter-surface cracks, allow general rounded areas, and bridging over open spaces. See Section 9.2 for our approach.  4.3.3  Snow Stability  Once we have allocated an initial amount of snow to the environment, we compute snow stability. These two major phases alternate until all snow has been added.  Our stability calculation has two main  requirements: stability determination, and unstable snow redistribution. T h e first part is responsible for locating all subdivisions (surface or otherwise) that have more snow mass than they can physically support. Unstable snow moves from these locations and flows downhill until a stable situation is reached, possibly crossing subdivision and surface boundaries, and undergoing any number of drops. If the surfaceto-surface drop is sufficiently large, the falling snow may be affected by wind. Because new snow is continually sloughing off from surfaces, the more times we run the stability algorithm, the more accurate our created surfaces will be. For example, if we wait until the very end  Chapter  J,: Assumptions.  Goals and Problems  75  to compute stability, our computations may produce a large volume of snow on the move. This may shock-load lower surfaces, cleaning them bare, when more realistically all surfaces would have remained largely stable due to a small, slow, movement of mass. All other things being equal, we prioritise correct snow location over correct snow stability (and subsequent movement). Sky-originated snow must always land someplace, and its final destination is often immediately obvious, while secondary effects due to unstable snow are often harder to see.  We claim that it is visually better to place shifting snow in  the wrong place than it is to leave a surface improperly bare. In order to determine the correct place for unstable snow, the viewer has to infer back the series of higher surfaces that might have produced that initial snow. Depending on the environment, accumulated snow can be very stable (or completely stable). Even unstable areas will release only periodically. Thus, a snow movement model may be needed only infrequently or not at all. We describe the properties that our snow stability determination algorithms should have. •  Partial stability. For a given unstable area, we must compute the mass of unstable snow and the mass of stable snow belonging to a particular subdivision. This essentially means that during gentle accumulation, the top surface of the snow sloughs off, but some of the underlying snow remains. A little dusting of snow usually remains on even the steepest slopes.  •  Non-determinism. We wish to avoid a completely deterministic model of stability, where accumulation beyond a critical point always results in snow movement. In reality, there are often small scale local effects (such as variable solar exposure or proximity to an out-cropping) that affect snow properties and thus snow movement. Our model should include some small random factor both over time and over area.  The actual stability calculations involve a summation of the forces on the snow mass.  Different  forces may be included (or not), or given greater (or less) physical importance than actually exist. At the minimum, stability calculations should include some notion of: •  Level-of-detail. We wish to be able to interactively control the level of detail we use to perform our stability calculations, such that stability calculations may be at least as coarse as our surface representation, and possibly coarser.  •  Angle of repose. As the angle of repose increases, the shear force increases and the compressive force decreases.  Chapter  •  J: t  Assumptions,  Goals and Problems  76  Snow cohesion properties. Snow with high cohesion (intra-particle bonding) resists being pulled apart more than snow with low cohesion.  •  Surface friction. The rougher the surface, the greater the friction at the snow interface.  •  Mass of snow present. The snow mass influences friction and the compressive and shear forces.  •  Shock loading. The arrival of moving snow mass should influence an area's stability, proportional to the momentum vector of the incoming snow.  •  Obstacle force. The motion of snow should be influenced by stationary obstacles.  Once snow has been classified as unstable, it starts to move and must be redistributed. Any redistribution algorithm should have these properties: •  Movement termination. Moving snow should be affected by a model of gravity and surface friction. This allows streams of snow to accelerate or deaccelerate as the surface angle of the underlying snow changes. Deacceleration due to friction allows moving snow to eventually come to a stop.  •  Multiple neighbour distribution. Unstable snow may leave a subdivision to several adjacent neighbours, depending on their closeness to the gradient of steepest decent and the velocity of moving snow. This redistribution may include falling over an edge. The amount of snow that travels to each neighbour may vary.  •  Primary velocity. We assume that unstable snow leaving a subdivision travels in a primary direction, called the vector of unstable snow motion s„. This direction is determined by combining the masses and velocities of all sources of incoming snow with the mass of the static snow, as well as some idea of friction, and the angle of repose of the travel direction.  •  Available neighbours. We assume that moving snow can always travel on to at least one other neighbour, even if that neighbour is higher than the original surface, or has been previously checked for stability. This implies that snow can travel up and over obstacles if it is moving fast enough. Consider a subdivision located next to a wall at the base of a large slope. A large amount of snow sloughs off the slope and reaches the low point at the base of the wall. By our available neighbour property, this snow must go somewhere, so it surges up the side of the wall. This may result in snow travelling over the wall, or falling back to its original location.  Chapter  •  J,: Assumptions,  Goals and Problems  77  Adaptive subdivision. Moving snow may end up in an area with very large subdivisions (such as snow falling off a roof onto a flat surface).  Without knowing all stability calculations a priori,  we have to assume that an arbitrarily complex pattern of snow may fall from above at any point in erly assign sufficient detail to just-arrived snow mass, we may be required to resubdivide the surface. Even more onerous, we may have to subdivide in the midumber and area of a subdivision • affects its ability to propagate snow. This last property is a very tricky requirement, because it means that the set of all possible neighbours and snow locations may change on the fly. It is extremely likely that this requirement will rule out global methods of stability determination. The available neighbour property will also cause us trouble, because it violates the idea that snow always travels downwards.  If this principle held, we could check for stability in descending order of  surface height. Now, because moving snow can ride up over a surface we've already checked, we may have to process surfaces out of order or multiple times. One can think of fast moving snow sloshing up and down alternate sides of a valley like liquid in a bowl, obviously quickly damping. The adaptive subdivision property is also quite demanding. This implies the set of stability sites is not constant within a single stability pass, and is especially problematic if the available neighbour property requires multiple passes to a single location. See Chapter 8 for our approach.  4.3.4  Edge-Crossing Snow-  Snow mass will occasionally fly off a surface (roof, tree branch, etc.)  and be deposited lower down.  Falling snow follows its previous velocity vector, affected by gravity, air friction, wind, and any surfaces it may hit on the way down.  Any falling snow redistribution algorithm should have the following  properties: •  Area distribution. We assume that a single mass vector of falling snow will spread out into an area following some distribution (such as a Gaussian). For some set of given snow properties and falling distances, there is a maximal spread of this area distribution.  •  Surface coverage. Every surface within the radius of this distribution should receive some part of the original snow. This problem is somewhat similar to the potential snow mass accumulation  Chapter  J: t  Assumptions,  Goals and Problems  78  problem, except that we can assume that only surfaces visible to the area distribution can receive snow. • Surface resolution, [f the area distribution is sufficiently small compared to the area of the target surface, we want to create a new surface sampling such that we can capture detail about the fallen "bump" (Figure 4.16).  new bump  possible new subdivisions Figure 4.16: Needs of a falling snow algorithm It should be noted that we specifically do not draw falling snow.  Chapter 5 Approaches  Real snow-covered yard.  Now that we have defined our problem, we must decide upon an appropriate set of computer techniques to produce a good model of snowfall. Obviously, there are many different ways to proceed, each with its own set of strengths and weaknesses. Because this research is one of the first to attempt this problem, it is useful to briefly discuss and describe some competing approaches, even though we do not implement all of them in our final solution. Comparing and contrasting these methods will highlight the strengths of our chosen algorithm, as well as providing ideas for future researchers trying alternative methods. We believe there are three major paradigms of note, differentiated primarily by the method used to compute snow distribution and the way snow on the ground is maintained. The volume-based paradigm locates and computes snow stability as a large number of volume elements, not necessarily connected to a surface. The surface-based paradigm represents snow as mass directly allocated to subdivided world surfaces without the use of particles to compute initial exposure. The hybrid-based model combine  79  Chapter  5:  80  Approaches  elements of the first two approaches, alternately representing snow as volume elements and as surfaceallocated masses.  5.1  Volume-Based Paradigm  Within the volume-based approach, we represent snow as a set of volume elements (such as particles [Reeves 83], voxels, or blobbies [Blinn 82]) that interact with each other and the environment. Throughout the snow location and stability phases, the snow accumulation remains represented as a volume, although the discretization of that volume may change. We consider previous work by Nishita [Nishita 97] to be a prime example of this type of approach. Figure 5.1 summarises some of the advantages and disadvantages of the volume-based method. -. volumes are natural : sampling elements o  stability requires interaction with many neignbours  qOO  forms cornices layering is easy A \ large number of elements  accumulation blocks snow discrete chunks are visible  potential holes  forms bridges  Figure 5.1: Advantages and disadvantages of a volume approach  5.1.1  Advantages  Volumetric models do not have to be associated or linked with individual surfaces, in the same way that a surface-based or hybrid paradigm might. Snow can accumulate anywhere, allowing shapes that do not require direct support by a particular surface. Snow can lean far out over an edge, create snow bridges,  Chapter  5:  Approaches  81  cornices, bulges, and generally pile up in a fashion more directly related to reality. This contrasts with underlying-surface methods, where bridges and unsupported surfaces must be treated as exceptions, and specifically created in a special way. In order to solve the partial occlusion problem, we need to determine some notion of (not necessarily straight-line) sky exposure. As snow accumulates, a given surface may be blocked from the sky by a height of snow on another surface, gradually changing the exposure of the first surface and its subsequent rate of accumulation. At the limit, a pile of snow elsewhere in the model may eventually completely block a surface from any further accumulation. Models in the volume-based paradigm will generally account for this problem, mainly because the position of a new element depends on the position of previously accumulated elements. Of course, the down side of this is that as more and more volume elements arrive in the environment, the cost of correctly locating the next one also increases. This location cost can be quite large, depending on the number of volume elements in the scene, although simplification and bounding box solutions exist. For example, we could represent a multi-particle aggregate surface by its convex hull (or an even tighter wrapping), computing a local destination for the final particle at the last possible minute.  Volumetric  models can better account for layering and intra-layering effects, because we do not specifically have to keep track of layer mixing. Assuming we have a sufficiently discretized model, individual volume elements can be given general layer properties as they arrive. If a surface sloughs off, falls, and mixes, we can determine the properties of the result by looking at the relationship of individual particles or volume elements. We obtain the final renderable snow surface by polygonizing surfaces belonging to the volume aggregate (possibly with [Wyvill 86]) with smoothing or crack filling as appropriate. Volumes may also be ray-traced directly [Kajiya 84].  5.1.2  Disadvantages  Probably the biggest disadvantage of the volumetric model is the potential for very large data sets. In order to get a sufficient surface resolution, the volume discretization may have to be very small. For example, assume we have a 10m by 10m scene, containing a house and several trees. To obtain a bruteforce snow deposit of 1 m, to a surface resolution of 1cm, we need 100 million volume elements. These sizes are not unreasonable, and are certainly not at the outer limits of the areas and depths we are interested in.  Chapter  5:  82  Approaches  There are two main problems with numbers of this size. If we deposit each volume element by tracing it from the sky, we must perform on the order of 100 million particle "locates". To compute a single locate, we may have to (worst-case), test the majority of surfaces at multiple locations along a particle's spiral path, including a large number of previously fallen particles. Secondly, the final position of each particle must be maintained, along with up to 100 million others.  Even at a few bytes per particle,  this immediately exceeds available memory and computational resources. It is likely that the only way a volume-only paradigm will be feasible is if we can perform some data elision by dropping particles with different masses, volumes, or some other form of weighting.  At the very minimum, we need to  prevent wide, flat, inactive surfaces from receiving the same discretization of volume particles as much more interesting areas in the model. If our world is completely flat without occlusion, we may drop 100 million particles to absolutely no end. This clearly calls for some scheme where surfaces receive an adaptive number of particles, based upon some criteria of "surface interest". In a surface-only model, we have the option of accumulating the entire snowpack at once (or in very large steps), which is much harder to do with a volume-only model. We could drop very large particles on large, wide open surfaces. But what happens if it turns out that our particle lands on a very small nearby leaf? The discretization of the volume accumulation must be matched to the size of the target surface. Like most sampling schemes, increasing the discretization increases the chance for error. If we are willing to trade view-independence (under our assumption, we are not), we could send different sized particles according to the distance to the viewer; however there is no guarantee that snow instability will not cause these large distant particles to return to the viewer and dwarf nearby smaller volume elements. If we assume that particles are spherical, then any accumulation of snow will contain many gaps and holes left by locally uneven snowfall, inflating the perceived snow depth, unless somehow "settled". Once volume-elements have fallen, they must be checked regularly for stability. In the volume-only paradigm, this implies we must perform some inter-particle or inter-volume force computation, similar to the n-body computation commonly performed elsewhere. The n-body problem is worst- case  0(N ) 2  operation, although there are several practical algorithms to perform it faster [Appel 86] [Barnes 86] and some obvious bounding box speedups to restrict interaction to a fairly local neighbourhood. Unfortunately, because individual particles may move (potentially a long ways) the "local neighbourhood" of interaction and intersection may quite large. This is expensive, although simplification algorithms exist [Milenkovic 96]. Even worse, we need to compute the stability algorithm many times (although with less than N particles all times but the last).  Chapter  5:  A pproach.es  83  To provide a size comparison, the largest gravitational N-bocly simulation ever clone [Warren 97] contained 17.2 million bodies, required the help of parallel supercomputers, and is considered a "Grand Challenge" sized problem. Clearly, a brute force stability algorithm with 100 million particles is totally out of the question. Obviously, we can do an immense amount of culling, reducing inter-particle interaction to a few nearby neighbours over the entire path of a particle's travel. Even so, we suspect that we cannot perform stability calculations on anything more than a few tens of thousands of particles. We have stated elsewhere that the snow location problem is more important to us than the snow stability problem, thus it does not seem practical to expend less computational resources on the former than the latter. Thus, we can greatly simplify volume-based stability calculations without worrying too much about accuracy beyond some minimum level. The most obvious way to gain speedups is to bucket nearby particles into a single volume, possibly associated with a surface, and then compute forces only on grouped volumes. This method starts to resemble the surface-based or hybrid paradigms. To summarise, volume-based paradigms are useful for solving the partial occlusion problem, creating bridges and unsupported accumulations, and modelling inter-layer interaction. However, once volume elements have landed, the computational expense of tracking and updating them is likely more than we can afford.  5.2  Surface-Based Paradigm  The surface-based paradigm represents snow based upon some subdivision of triangles of the underlying model, where all.snow in the environment must "belong" to a particular surface. Partial occlusion is computed without the use of volume-type primitives like particles. Figure 5.2 summarises some of the advantages and disadvantages of the surface-based method.  5.2.1  Advantages  Snow height (and the distribution of future snow accumulation) is stored implicitly, and only computed when needed.  Thus, as snow accumulates, we must recompute surface visibility or else potentially  receive improper snow accumulations. This is one of the big tradeoffs of a surface-based representation: we shorthand many objects into one to gain simplicity and speed, at the loss of accuracy. The resolution of the volume-based model depends on the number of volume elements used to represent snow, while the surface-based paradigm gains resolution through the number of subdivisions on an individual surface.  Chapter  5:  84  Approaches  A A A  rays sample sky coverage simple neighbour interaction no curved paths  n  o  cornices  no wind effects  "  does not handle accumulation blocking mass depends on sky coverage  subdiv.s.on resolution depends on need  bridges need special handling  F i g u r e 5.2: Advantages and disadvantages of a surface approach  We expect that the number of surfaces in a model will be much smaller than the number of particles a volume model might need; thus we should be able to represent snow on the ground in a much more compact way.  5.2.2  Disadvantages  Probably the biggest drawback with surface-only models is computing the initial exposure to snow, without the use of particles or volume elements as a sampling tool. One possible approach would be to compute some form of snow "form factor", parallelling structures used in radiosity computations [Cohen 93]. The general idea here is to precompute all inter-surface interactions, and then solve a single set of equations to simultaneously determine the amount of snow accumulation. A form factor here is a measure of the exposure a particular surface has to snow originating from the sky (or another surface). The big advantage is that we can add very large amounts of snow for much less than the cost of tracing individual particles. The main problem here is that snow does not travel in a straight line, (due to small and large scale wind) implying that completely sky-occluded surfaces can still receive snow. Furthermore, the amount  Chapter  5:  Approaches  85  of snow an occluded surface gets is dependent upon the exposure of the (potentially many) occluding surfaces, mainly because each previously encountered surface will remove some of the swirling snow. The amount removed depends partially upon the distance from the previous blocking surface. The gradual buildup of snow changes the occlusion relationship between surfaces, forcing the sky form factors to be recomputed every so often. The form factor idea is most useful if we sent out particles, instead of straight-line rays (suggesting a hybrid surface/volume approach). Each surface generates a (possibly adaptive number) set of particles, which are traced towards the sky. The snow a surface receives is dependent upon the number of particles that reach the sky, and the available snow mass in a particular sky patch. Using particles places this method within the hybrid paradigm. The form factor idea can be extended to the snow stability calculation, where the approach would be to precompute the probability and amount of snow a given surface a would contribute to another surface b and then solve multiple simultaneous equations to produce the final result. Unfortunately, this is extremely hard to do in practice. Depending on the configuration of snow-dumping surfaces above surface a, outgoing snow may travel in any direction, at any speed. Thus, the probability of a certain amount of a's snow reaching surface b is non-constant, and depends upon the inputs to a. These, in turn, depend upon a higher set of input snow sources, and are always changing as neighbours continuously vary in snow depth. Additionally, there is not necessarily a single path of travel between surface a and surface 6. As snow travels over a set of intermediate surfaces, there is a probability it will stop, continue, or be deflected in a new direction with a new speed. The number of intermediate surfaces may be very large, may increase as new neighbours are added on the fly, and may include drops. All of these factors make precomputing a snow form factor extremely difficult. Furthermore, this approach is also very inefficient.  Any given surface has the potential to receive  snow from any other given surface. For most distant surfaces, wind is the primary transport mechanism, although two surfaces can exchange snow during a huge avalanche that rushes up the sides of large obstacles. The problem here is that even though two surfaces may potentially receive snow, the actual probability that they do so may be very, very small, leading to large numbers of inter-surface form factor computations that may be used only once, or not at all. For example, most of the surfaces belonging to two different trees will rarely interact. In this case, it is easier to compute the interaction only when it happens. Stability form factors make more sense in locations where there is continual flow of snow between  Chapter  5:  SO  Approaches  surfaces. Realistically, we expect this to be fairly small set of surfaces because even at the most unstable locations, there will be a fair number of snow accumulation increments for every step causing snow movement. The infrequent nature of the stability calculations make it easier to compute stability with a sequential, as-it-happens algorithm. Snow layering can also be represented in a surface-only model. However, if there is appreciable mixing during snow movement, the number of layers may increase very significantly (and the depth of each layer decrease, until the layers represent volume elements at the limit). This research does not deal with snow layering, but it is useful to think about as an eventual goal.  5.3  Hybrid Paradigm  As we have seen, both the surface-based and the volume-based paradigms have some advantages and disadvantages.  We would like to combine the best of both to produce a hybrid model.  Figure 5.3  summarises some of the advantages and disadvantages of the hybrid-based method. volumes are natural : sampling elements simple neighbour interaction  ©  no cornices  handles curved path becomes part of surface fewer subdivisions  subdivision resolution depends on need Figure 5.3:  A ,i  accumulation blocks snow coverage is uniform  bridges need special handling  Advantages and disadvantages of a hybrid approach  The partial occlusion problem indicates that volume-like particles are best suited for computing snow  Chapter- 5:  Approaches  87  exposure, because they can easily handle wind fields, swirled paths, and inherently perform sampling. Secondly, when snow flies off an edge, it must be moved to its final destination based upon gravity, wind and air friction (also likely a particle system). It would be nice to have the same algorithm for all in-air particle movement, either due to initial snowfall or unstable movement. This is more consistent, and easier to implement. In the volume model, we argued that the number of interactions required for a volume-only stability calculation far exceeded the accuracy we needed. This strongly argues for a stability model based upon surfaces, where snow is represented abstractly. Secondly, it is to our advantage to only remember particles while they are in transit, meaning that particle storage is constant, rather than increasing proportional to the number of dropped particles. We propose to use a surface model to represent snow, where individual volume elements are added to something like a surface height-field. As particles arrive, the appropriate locations grow in depth, possibly subdividing to account for differences in snow distribution. In order to handle changes in snow occlusion due to increasing snow depth, we propose to replace each model surface with a non-zero snow surface. This snow surface will slowly rise away from the original surface, blocking incoming particles as it accumulates. Particle intersections are computed against the snow surfaces, and not the original surfaces (assuming no upward motion due to wind), meaning that the intersection cost does not increase proportional to the number of particles dropped. Because particles are subsumed into the surface representation, there is no need to settle particles so that they are densely packed. Stability calculations are performed solely by looking at subdivision-to-subdivision interactions. This means that stability calculations are performed only when required, and moving snow only has a series of local interactions.  Snow that flies over an edge is likely best handled using the same distribution  method we use for the initial exposure. One disadvantage of using a surface representation for snow is that it is not quite so obvious how to perform bridging, overhanging, or other forms of snow buildup that aren't directly supported by an underlying surface. We can use our previously introduced scenario for a very general work comparison with the volumeonly approach. We assume that our 10m by 10m area is divided into 1,000,000 subdivision surfaces to achieve a 1cm surface resolution. If we compare with the 100 million particles sent in the volume-only example, we are left with a budget of 100 particles per surface subdivision. This particle budget must be spread out over all partial occlusion computations: i.e, we might send out 20 particles every 20cm of  Chapter  i>:  Approaches  88  snow depth. We expect our model paradigm to be cheaper, if only because we can skimp on the number of times we recompute our snow occlusion (at a cost of accuracy), the hope that areas of low interest will not need all of their allocated particles. As well, we can obviously reduce the amount of work required by using fewer surfaces; using 10,000 subdivision surfaces gives us an equivalent work budget of 10,000 particles per surface. The more large, flat polygons present in the scene, the fewer surfaces we need.  Chapter 6 A Computer Model of Fallen Snow  Our first tracks.  Beginning with this chapter, we describe the core algorithms of our research - a model of snow that can be used to create realistic-looking images. By necessity, this model will be a gross simplification of what actually happens in Nature, partly because snow is extremely complex, and partly because we can claim success if we can provide a framework that allows later refinement through improved experimental data. To soothe the reader, we have separated our approach into several distinct chapters. This chapter contains a global overview of the algorithm, as well as information on obtaining input worlds, and a description of some important data entities.  Additionally, we describe the relationship between the  phases the algorithm passes through before completing.  89  Chapter  6: A Computer  Model of Fallen  Snow  90  In Chapter 7, we tackle one of these major phases (snow accumulation), explaining our approach to determining how much snow falls where. Chapter 8 contains the other major phase, snow stability. Finally, Chapter 9 contains some miscellaneous steps required for a final solution but outside our main areas of interest and exploration.  6.1  Overview of the Solution  As described in Section 4.1.5, our snow-adding algorithms are part of a larger pipeline involving commercial modelling and rendering programs. Our particular area of concern is determining how much sky-origin snow should be placed where, and to a lesser extent, determining the final destination of unstable snow. To determine how much snow a particular surface should receive, sites on each surface launch a series of particles aimed upwards toward the sky. As particles flutter upwards, they are checked for intersection with intervening surfaces, where a "hit" indicates that a particle is somehow blocked, and cannot contribute snow to its source surface. As particles reach or are blocked from the sky they slowly build a picture of a surface's sky occlusion. We extrapolate the sky occlusion statistics for a particular surface to arrive at the snow mass an area should receive. Surfaces are subdivided whenever the particle tracing indicates that there is an interesting occlusion boundary, and likewise combined whenever we are fairly confident that a surface is either completely exposed or completely occluded. At first glance, our general approach towards generating snow is counterintuitive: we concentrate on surfaces instead of snowflakes, and shoot flakes upwards rather than let them fall downwards. It is our belief that our chosen method gives us a greater degree of control than other approaches, specifically in being able to precisely limit the amount of work we need to generate a particular look. As soon as we have generated a mass accumulation picture that meets some resource criteria (compute time, number of samples, size of sample or some other importance-driven function), we can add an appropriate (and arbitrary) amount of snow. This generates a complete set of 3D snow surfaces that rise off the base model. All newly added snow is then subjected to a stability test. Snow-covered sample points are sorted by approximate height, and compared with connected neighbours using a simple model of snow material  Chapter  6:  A Computer  Model of Fallen  Snow  91  properties. Snow from a locally unstable sample is sluffecl off to nearby downhill neighbours, until the initial sample is either stable, is free of snow, or is blocked by an obstacle or a growing pile of snow. Occasionally, a sample will have no downhill neighbours, causing snow to avalanche over an edge or drop-off.  Falling snow is simulated as a small number of particles, and is traced down through the  model until it eventually comes to rest. During stability, launch site resolution may increase in areas where snow falling from above generates a pattern that is not sufficiently represented by the existing site density. The movement of snow down a slope may affect the stability of both uphill and downhill neighbours, so the stability algorithm performs performs multiple passes over the sorted list of launch sites until the overall snow motion is small or the program runs out of allocated time. Since the addition of a layer of blocking and obscuring snow changes the previously computed mass accumulation pattern, we can repeat the accumulation and stability steps as desired, increasing accuracy at the cost of computation time. Once the last pass is complete, we generate a set of polygons representing the completed 3 D snow surface. Using the mass pattern computed from the accumulation and stability steps, we also generate a set of textured areas representing surfaces that have some snow coverage or flake dusting, but do not retain enough mass to be represented as a solid layer. As a final, optional step, snow volumes are replaced with several different types of implicit functions, based on snow properties, type and location. This closes up cracks, forms bridges across gaps, creates cornices, and improves the smoothness of the overall result. We polygonize the resulting isosurface and perform mesh reduction to reduce the overall number of polygons. The final model can then loaded into commercial software for rendering, manipulation, or as the starting point for a new layer of snow.  6.2  Entities  This section contains an overview of some of the objects and entities that exist in our solution, essentially a very high level object-oriented data dictionary. Later sections will describe how these objects interact.  Chapter  6:  A Computer  Model of Fallen  Snow  92  Figure 6.1: The model is a collection of connected and non-connected polygons that eventually accumulate snow.  6.2.1  World  The "world" is the volume containing the entire snow-generating system, including the sky, the ground, the wind, the original input model and any allocated snow. The world consists of a rectangular, axisaligned box completely containing and bounding everything of interest. The "sky" is the topmost Z plane of the box. The bottom Z plane is the "ground", and is used to detect and catch unstable snow (such as might fall off a roof) that would otherwise disappear into the void.  6.2.2  Model  The "model" is the set of input polygons that collectively form the base snow surfaces for later snow collection, and may consist of many different connected and non-connected components. Figure 6.1 shows an example input model, shown in Alias|Wavefront. Models are (ideally) obtained from animators and artists.  6.2.3  Faces  The "face" is the primary structure used to represent model geometry. Each face is a planar triangle, with one or more bounding vertices potentially shared with neighbouring faces. The initial connectivity between faces originates from the structure of the mode, as decided by the creating artist (Figure 6.2).  Chapter 6: A Computer Model of fallen Snow  93  k  Figure 6.2: The model is composed of a number of (potentially connected) faces. Note how the pine needles in this model are not actually connected to the branches. Each face has a normal /normal,  a  computed area  f a, are  and a slope angle  fi  s ope  (where f [ s  ope  = 0 is a  flat surface). The position, orientation and number of the initial faces remains constant throughout the snow simulation. Faces with a downward normal are considered too steep to retain snow.  6.2.4  Launch Site  A "launch site" is the basic unit of sampling in the snow algorithm, where each launch site represents the mass accumulation pattern of some amount of surrounding surface area. Every snow-eligible face is given at least one launch site, with additional launch sites added based upon importance properties of the surface (Figure 6.3).  6.2.5  Subdivision Area (or Launch Area)  Each launch site is responsible for an immediately surrounding area, termed the "subdivision area", or "launch area". As launch sites are added or removed, a subdivision area may change in size and shape. The amount of snow a launch site receives is spread out on its associated area. Section 7.3.2 describes the representation used to divide a face into areas, while Figure 6.4 shows area boundaries.  Figure 6.4: Subdivision areas (blue boundaries) surround launch sites.  Chapter  6:  A Computer  Model of Fallen  Snow  6.2.6  Edge Group (or Subdivision Mesh)  95  An edge group is a set of launch sites and their associated subdivision areas, bounded on all sides by drops. Edge groups are used to isolate connected components - within a mesh, launch sites handle avalanches and generate snow planes differently from launch sites that are not in the same mesh. Different meshes (usually, but not always) indicate a geometric discontinuity, such as that shown in Figures 6.5 and 6.6. All launch sites on a single face belong to the same mesh. A mesh may span multiple faces.  I  El  1  6.2.7  Figure 6.5: Three geometric objects, plus  Figure 6.6: Each object forms a different  the +Z axis.  edge group.  Drops  A "drop" is the edge of a single face, and is used to denote the outermost boundaries of an edge group. A sequence of drops isolates one edge group from any adjacent edge groups, and are generally located where the change in adjacent surface normal is large. When travelling snow reaches a drop, it is converted to airborne particles influenced by gravity. Although drops are conceptually edges where there are no adjacent faces (i.e., snow drops off the edge) drops may actually split connected, adjacent, coplanar faces where there is no change of elevation. Certain degenerate models may have unconnected drop edges.  Chapter  6: A Computer  Model of Fallen  Snow  Drops are discussed further in Section 7.3.1.  6.2.8  Flake  The flake is the basic sampling unit of sky occlusion. Flakes wiggle and swirl upwards from launch sites towards the sky plane.  6.2.9  Flake Group  In order to determine how much snow falls on a launch site, each site repeatedly fires groups of flakes at the sky. The running percentage of successful flakes over all groups helps determine the amount of snow a given launch site should be allocated.  6.2.10  Snow Planes  The 3 D model of snow is composed of a number of snow planes, generated by interpolating adjacent launch sites in an edge group. There are two general types of snow planes - top snow planes and edge snow planes. Top snow planes are triangular, and connect nearby launch sites. Edge snow planes are quadrilateral, connect two launch sites, and are always vertical. For a given subdivision mesh, the union of the faces, the top snow planes and the edge snow planes produces a closed volume bounded by polygons (Figure 6.7).  Figure 6.7: A single triangular face, with subdivided top and side snow planes that form a closed volume. Note: all surfaces are triangles, despite the appearances of the drawing.  Chapter  6.2.11  6:  A Computer  Model of Fallen  Snow  97  Avalanche  An avalanche is a small delta of unstable snow that travels from one launch site to another, or from a launch site to some intersection point on a drop. When an avalanche intersects a drop, it is converted into several avalanche flakes. A single launch site may generate several simultaneous avalanches.  6.2.12  Avalanche Flake  When an avalanche hits a drop, it is converted into a number of particles called "avalanche flakes". These particles are traced downwards, under the influence of gravity and wind, until they reach a stable destination or disappear into the void. Avalanche flakes simulate snow showering from above.  Chapter  6.3  6:  A Computer  Model of Fallen  Snow  9S  Overview of Phases  Figure 6.8 shows an overview of the snow pipeline. T h e input consists of a 3 D model, complete with lighting and material rendering properties.  The  polygonized 3 D model is passed on to the snow algorithm, along with user-specified snow parameters. The computed snow surfaces may be passed through an optional smoothing phase (not shown in Figure 6.8), that produces bridging. The polygonalized base model is discarded, and the final snow surfaces are loaded back into to the original, unmodified scene, along with appropriate snow material rendering properties. The complete scene is then rendered, using the combined world and snow surfaces and material rendering properties. The snow accumulation and stability algorithms can be run in a series of phases, increasing accuracy at the cost of additional time. Figure 6.9 gives an overview of the optional and required phases available to the user.  6.3.1  Storms and S t o r m i n g  A "storm" consists of a number of (possibly non-uniform) discrete time steps, with a corresponding continuous function modelling rate of snow deposit over time. On each time step, the snow accumulation algorithm determines how much snow accumulates on the model surface, while the snow stability algorithm determines the destination of unstable snow. It is possible to turn snow stability off for illustration purposes or for models where stability is not important.  6.3.2  S u b d i v i s i o n A p p r o x i m a t i o n Phase  This phase happens once, at the start of the storm.  The subdivision approximation phase tries to  compute most of the useful and interesting surface mesh areas before snow is added. The hope is that all interesting occlusion areas can be determined, sampled, and represented with a sufficiently dense mesh while the number of surfaces is low. As soon as a layer of snow is added to the world, the additional snow top and snow side surfaces can slow down the intersection algorithm considerably, and degrade the sampling effect received in return for each particle. Without snow surfaces present in the world, we test for intersection and occlusion using using entire face planes. The subdivision upon the faces does not affect the number of intersections we must perform, thus making it cheaper to increase mesh density on the fly. Thus, if the timing budget is small, this  Chapter  6:  A Computer  Model of Fallen  Snow  99  ALIAS world lighting (INPUT)  world geometry (INPUT)  world materials (INPUT)  USER snow parameters complete world scene (OUTPUT)  polygonized world geometry  (INPUT)  SNOW ALGORITHMS  resolve stability  snow surface properties  generate smoothed snow surfaces  discard polygonized world geometry  ' ALIAS  snow material  snow surfaces  y world & snow geometry  world & snow materials  world lighting  / complete world scene with snow (OUTPUT)  Figure 6.8: Ari overview of the snow pipeline phase is usually given the majority of available duration.  6.3.3  Regular Storm Phase  After the initial subdivision approximation phase, there can be a large number of subsequent regular storm phases. Each regular storm phase is allocated a certain, variable amount of clock time to compute snow allocation for an amount of simulated time. All regular storm phases keep track of many more surfaces than the initial subdivision approximation phase, since now the model maintains all resampled snow surfaces created by joining launch sites, plus vertical faces surrounding edge groups.  Chapter  (>: A Computer  Model of Fallen  Snow  100  phase is optionally repeated n times  Subdivision Approximation Phase Snow Accumulation  Snow Stability  Regular Storm Phase Snow Accumulation  (no snow planes included)  Snow  Finished  Stability  (snow planes included)  phase is optionally skipped  Figure 6.9: Overview of required and optional phases. Although the hope is that most of the mesh variation occurs in the cheaper first phase, regular storm phases can still undergo any required changes in mesh resolution. Most of the mesh resolution variation is due to unfinished computation from the first phase, or a changing accumulation pattern due to blocking by the rising level of snow. Given a limited timing budget, regular storm phases usually are given a short duration, and in many cases, skipped altogether. Changes in the mass accumulation pattern because of rising levels of snow are subtler to detect, and thus may be less visually important. In this case, the subdivision approximation phase is usually sufficient if the amount of snow on the ground is fairly small compared to the total size of the model.  6.3.4  Uncompleted Phases  In all accumulation phases, each launch site must have a chance to fire at least once before any mesh refinement occurs. In the subdivision approximation phase, if the user has not allocated enough time to the do this, the algorithm continues to work until every launch site has at least some basic information. Otherwise, there would be areas of the initial mesh with absolutely no knowledge of snow coverage. In the regular storm phase, we also test all existing launch sites before we allow any additional mesh increases. The rationale here is that the subdivision approximation phase already determined a good pattern of launch sites, and the additional particles are better spent increasing the accuracy of already-present occlusion boundaries over adding new boundaries. If there is not sufficient time to test all existing launch sites once, the phase terminates anyway, with some sites untested. This is not a problem, since each site already has an idea of snow accumulation from one or more previous phases. Once all existing launch sites have been tested at least once, mesh resolution variation may occur.  Chapter  6.4  (i: A Computer  Model of Fallen  Snow  101  Acquiring the Surface Model  At least as important as the snow is the construction and appearance of the underlying bare surface. If the initial scene looks unrealistic, it is especially challenging to make the resulting snowy version look any more natural. Unfortunately, it is quite hard to find complete models of natural-looking outdoor scenes. We ended up cobbling together a number of independent models and adding lighting, texturing and animation sometimes with less-than-perfect artistic results. Most models were obtained from various sites offering public archives of free 3D meshes [3DCafe 00] [Avalon 00]. Models were not selected for any particular criteria, and were generally combined ad-hoc to produce a sufficiently complex scene. One pleasant side effect of snow: a fresh layer can often hide, and obscure problems with the initial base scenes. For example, the poorly parameterised grass texture in the foreground of Figure 10.11 has been hidden in Figure 10.12.  6.4.1  E x p o r t i n g the Surface M o d e l F r o m Alias|Wavefront  We obtain our initial underlying world geometry from ".obj" [Murray 94] files exported from Alias|Wavefront. The .obj file provides a polygonal description of the world, with minimal information about lights, shaders and materials. This representation was chosen because it is simple, and contains a minimum of extra data that must be parsed and discarded. Objects are described as a list of common vertices, followed by faces composed of a list of enumerated vertex references. Practically, using this format does not limit us with respect to the original model source - Alias| Wavefront contains import/export features that perform conversions from a wide variety of other programs and formats. During export, Alias| Wavefront provides a number of features that allow us to modify our geometric world descriptions, specifically the polygonal approximation of curved surfaces, and the triangulation of polygons. This provides us with a model description where every face is a planar triangle; an uncomplicated description that will greatly simplify later steps. Curved surfaces may be polygonized with several methods to meet various curvature tolerances, with a corresponding variation in the number of initial polygons. The user is required to select an appropriate curvature representation, keeping in mind the available memory resources and the uses of the model. It is important to remember that all original surfaces are used in the final rendering, and the polygonal  Chapter  6: A Computer  Model of Fallen  Snow  replacements are only needed for the purposes of snow allocation (see Section 11.2.2)  102  Chapter 7 Snow Accumulation  Snow accumulates underneath a real bench.  7.1  Overview of Snow Accumulation  Peculiar to snow is the idea of "flake flutter", where the path of falling ice crystals are affected by crystal shape and atmospheric micro-turbulence. These local disturbances can prevent falling snow from descending in a straight line, instead allowing flakes to sidestep blocking obstacles and land underneath on surfaces that have no direct exposure to the sky. Thus, simulating and modelling an accumulation pattern is akin to raytracing for light, except that we are interested in path (instead of straight-line) visibility. Where an obstacle, such as a porch or a bush, blocks the ground underneath, the flake flutter effect eventually produces an occlusion boundary between completely blocked and unblocked areas. An example of this can be seen in Figure 7.1(a), where snow accumulates well underneath the overhang of the bush. Over billions of flakes, these occlusion boundaries exhibit a smooth drop-off, where the shape 10:',  Chapter  7: Snow  Accumulation  104  of the curve and amount of snow under an object depends on the size, shape, and number of blocking occlusions, the closeness of the occlusion to the ground, and the magnitude of the fluttering effect. In objects with many occluding components (such as a pine tree) the occlusion boundaries are still present, but are much less pronounced. Most falling snow accumulates on the uppermost layer of branches, but some accumulates on the next layer, and most lower branches and the ground get at least a small dusting. This contributes to the visual impression that snow is everywhere in a scene, and not just sitting on the uppermost surfaces exposed to the sky.  7.1.1 Our  C o m p u t i n g the Snowfall A c c u m u l a t i o n P a t t e r n  goal is to generate an accumulation pattern for every surface in the model, where the amount of  snow each surface receives is proportional to the occlusion factors described above.  Our approach is to  allow launch sites on each surface to emit a series of particles aimed upwards towards a "sky" bounding plane. As particles flutter upwards, they are checked for intersection with intervening surfaces, where a "hit" indicates that a particle is somehow blocked, and cannot contribute snow to its source surface. A "miss" means that the particle made it through or around all blocking obstacles and reached the sky. As particles reach or are blocked from the sky they slowly build a picture of a given launch site's sky occlusion. New launch sites are added whenever the particle tracing indicates that there is an interesting difference with nearby neighbours. Likewise, launch sites can be removed whenever nearby neighbours are consistently confident that they are either completely exposed or completely occluded. As soon as we have generated a mass accumulation picture that meets some resource criteria (compute time, number of samples, size of sample or some other importance-driven function) we can add an appropriate (and arbitrary) amount of snow. This generates a complete set of 3 D snow surfaces that rise off the base model. •Since the addition of a layer of blocking and obscuring snow changes the previously computed mass accumulation pattern, we can repeat the accumulation step as often as desired, improving the approximation at the cost of computation time.  7.2  Importance Ordering  The rationale for shooting upwards generally arises from the need for control: the idea that each individual surface can locally influence its resolution by deciding how many launch sites it needs, and how  Chapter  7: Snow  Accumulation  105  many particles each site should shoot. Since our sampling rate is orders of magnitude less complete than Nature's, prioritising the few samples we do have allows us to make better use of them. This ensures that even the tiniest surface is guaranteed at least a rough estimate of snow accumulation. This is a major advantage over potential approaches that drop blobby particles, since small surfaces are often missed at the expense of covering large ones. Figure 10.23 shows how our approach covers individual blades of grass, despite being in the middle of a very large snowy field. Each launch site is given an importance ordering used to determine order of site testing, to determine the number of particles to shoot per site, and to decide if more sites are needed nearby to improve the resolution. As long as the allocated time has not expired, the most important launch site shoots a small batch of particles, gets a new importance based on the results, and is placed back in sorted order. The importance ordering is a heuristic weighting based on the following factors (weighted in Algorithm 1). •  Completeness. Launch sites with no previous chances to shoot are more important than sites that have had at least one chance, ensuring a crude global approximation exists before any further refinement begins.  •  Area. As the area of a launch site increases, particles from a single site will pass through less of the volume immediately overhead. To prevent missing occlusions, large areas may need more particles per launch site and more initial sites.  Occlusion boundaries in large areas are more  visually obvious, and so gain preferential allocation of new refinement sites. This is measured with a rough metric involving "flutter area". This an area based upon a parameter used to control the amount of lateral motion of a snowflake; a large flutter area implies that snowflakes will reach more areas of the sky. See Section 7.5.1 for more on this parameter. •  Neighbourhoods. If the particle hit percentage of two neighbouring sites is sufficiently different, it implies that there is a nearby obstacle causing some kind of occlusion boundary. Both sites gain importance, asking for more particles to improve knowledge of the shape, orientation and magnitude of the boundary.  If the neighbours are sufficiently different and important, a new  refinement site may be added to the perturbed midpoint. Likewise, launch sites that are the same as all nearby neighbours become less important, and may be candidates for removal. •  Effort. If all other factors are equal, launch sites should use approximately the same number of particles, aiming for consistency of confidence.  •  Limits. The user can set several parameters that limit the approximate scale of the finest allowable  Chapter 7: Snow Accumulation  106  increase in resolution. This prevents launch sites from increasing indefinitely along very complex occlusion boundaries. If all sites have been resolved to this limit, the phase can terminate early. • Steepness. Launch sites that are too steep to support much snow are less important than other sites. If a surface is quite steep, all of the accumulated snow will usually (unless blocked) slide off to eventually come to rest on more stable surfaces. If flake dusting (see Section 7.7.1) on steep surfaces is important to the look of the scene, this ranking can be turned off. • Camera. When optionally enabled, sites closer to the camera receive more particles, greater refinement, and improved accuracy at the cost of imposed view dependence. See Section 7.2.1 for more discussion on this. • User. Importance ordering allows users to arbitrarily tag surfaces as being "boring" - useful for ignoring areas that will eventually be occluded or matted out. Likewise, "interesting" areas can be highlighted for additional processing and better accuracy. For example, a user might indicate that areas near the center of the model are more important than boundary areas. Algorithm 1 gives a (simplified) pseudo-code description of the importance ordering used for most of the images in this thesis. This particular weighting performs a stratified ranking, where launch sites with no information get the highest priority. After the first pass, the importance of all sites drops considerably, leaving sites with at least one neighbour requiring subdivision at the top of the priority queue. Since there will likely be a large number of these sites, they are additionally distinguished by coverage area, repose, and optionally camera distance, where larger areas gain much more attention than small ones. Once a new subdivision is created, it leaps to the head of the queue (since it has shot no flakes), is processed, and then drops down based upon area. Any time a subdivision occurs, it shatters an area into much smaller pieces, sending it down the list considerably. The actual weighting is much less important than the idea that some launch sites get freer access to a limited sampling budget than others, based on criteria important to the user for a particular scene. Figure 7.1(a) shows the occlusion boundary under a real snow-covered bush, illustrating the type of visual effects we want sampling to determine. After 10 seconds, the importance ordering has found the boundary, and generated an initial approximation. Spending an additional 90 seconds results in more subtle improvement, refining launch sites of less visual interest. Note that the meshing shows launch areas, but not launch sites. There is no stability in this example.  Chapter 7: Snow Accumulation  Determination of Launch Site Importance Procedure Computelmportance(launchSite) / / a launch sites with no shots and mass from an avalanche mass-importance = 1 if has mass, but no flakes sent; 0 otherwise / / do all untested launch sites before creating any new ones first-pass-importance = 1 if no shots sent; 0 otherwise / / rank sites that need subdivision neigh-importance = 1 if a neighbour needs subdivision; 0 otherwise / / rank larger launch areas as more important coverage-importance = ratio of area/flutter area / / i f everything else is equal, try and make sure / / all launch sites do the same number of passes pass-importance = M A X - P A S S E S - number-passes / / surfaces that are too steep go last repose-importance = 0 if too steep to support snow; 1 otherwise / / optional, based on distance to camera camera-importance = distance to camera / maximum distance / / these weights are arbitrary, but should be >> other weights, / / t o force these to be the most important importance = first-pass-importance * 50000 + mass-importance * 10000 + / / these will be the same for a large number of sites neigh-importance * 500 + pass-importance + / / these actually do most of the distinguishing coverage-importance -fcoverage-importance * repose-importance + coverage-importance * camera-importance return importance end Procedure  Algorithm 1: Determination of Launch Site Importance  107  108  Chapter 7: Snow Accumulation  (e) Snow after 100 seconds  (f) Mesh after 100 seconds  Figure 7.1: (a) A bush leaning out over a wall provides a real example of the flake-flutter phenomena, (b) Initial meshing of a crude bush model. No measurement of the real bush was done, (c) After 10 seconds, importance ordering has found the general shape of the boundary - surfaces are created by resampling launch sites within each launch area, (d) The denser mesh reflects the more interesting launch areas. A significant amount of refinement occurs behind the bush and is not visible from this viewpoint, (e) After 100 seconds, the boundary shape is essentially the same, due to the importance ordering of launch areas, (f) The denser mesh after 100 seconds.  109  Figure 7.2: When optional camera importance is enabled, surfaces closer to the camera are more likely to subdivide. This image shows the viewpoint from the camera; the next Figure shows an alternative viewpoint. For clarity, snow has been removed, and obstacles have been made transparent. Flake dusting is still visible around the edges.  7.2.1  Closeness to Camera  Although our algorithm is view independent, it offers the opportunity for increased visual accuracy when the position of the camera is fairly fixed with respect to the scene objects. When camera importance is enabled, surfaces closer to an approximate camera position are given a higher importance than surfaces more distant from the camera. This means that there will be more subdivision, increased resolution and better accuracy in areas that are more likely to be viewed (Figure 7.2 and Figure 7.3). More distant surfaces receive fewer increases in resolution, leading to larger errors and coarser approximation, but faster results. The more distant a surface, the more likely the approximations will be missed or the entire surface occluded. The decision to enable or disable camera importance generally depends on how freely the camera is expected to roam throughout the model.  Chapter  7: Snow  Accumulation  110  Figure 7.3: The previous scene viewed from the side. The ball at the left edge of the scene is the camera position. Note that although areas closest to the camera received the most mesh subdivision, the distant side of the model also received some subdivision. Flake dusting is still visible around the edges of the scene.  Chapter 7: Snow Accumulation  7.3  O v e r v i e w of L a u n c h Site  111  Meshing  Launch site surfaces are represented as triangles, generated from the original (potentially non-polygonal) base models.All up wards-facing triangles in the approximation of the underlying model are initially allocated at least one launch site. Additional launch sites are allocated based on the importance ordering of the surface, user-set resolution parameters, and the magnitude of the flake-flutter. In order to properly allocate snow, each launch site must be responsible for some non-overlapping portion of the surface, ideally the area immediately surrounding the sample point. We have chosen a strategy based upon Voronoi diagrams [Mulmuley 94] [O'Rourke 94], although there are numerous other valid meshing possibilities. Launch sites are connected in a constrained Delaunay triangulation, where each launch site is responsible for its own immediately surrounding Voronoi area.  Voronoi areas are  computed within a single face, and clipped to the edge of the triangle for maximal surface independence. Advantages of this approach include fast point-in-area tests and neighbour location, and the ability to quickly generate triangulations for intersection testing.  Figure 7.1 (b) shows an example of a sparse  initial mesh undergoing the addition of more and more launch sites, shown in Figures (f).  7.1 (d) and  Note how neighbouring constrained-Voronoi areas vary in size at the transition zones, and mostly  minimise extreme angles. In practice, many surfaces are small and isolated (such as the brush and pine needles in Figure 10.12), and meshes are reduced to the trivial case of one or two samples in a triangle. Significant meshing occurs on large, connected surfaces, such as the ground.  7.3.1  Edge Groups  Launch sites and their associated meshes are additionally divided into edge groups, which are isolated world objects, projected into the X Y plane, bordered by the X Y silhouette edges. Launch sites in the same edge group can be neighbours, share information, and be resampled, subdivided and coalesced with nearby edge group members. Launch sites in different edge groups cannot directly communicate. This simulates a significant change in angle-of-repose between eligible and non-eligible snow supporting surfaces. Projecting objects into X Y implies that launch sites can only be placed on surfaces with an angle of repose of [0..90)°. Since edge group silhouettes are not necessarily convex, we must do some additional processing to "break" constrained Delaunay neighbour links that cross a silhouette boundary, as described in Section  Chapter 7: Snow Accumulation  112  Edge groups are used primarily for avalanche resolution, denoting sharp boundaries where snow may slide off from one edge group to another. It is convenient to think of edge groups as places where a physical drop occurs, but this is not always so. Edge groups may be coplanar and adjacent, depending on the initial geometry of the model. They may also interpenetrate. A single edge group may also be arbitrarily broken into smaller edge groups, although this is inefficient since moving snow across group boundaries is more expensive than moving snow within the same edge group. Figure 7.4 shows how a sphere is converted into an edge group.  Figure 7.4: An isolated object (a), bordered by X Y silhouette edges (in red) forms an edge group after projection on the X Y plane - top view (b), side view (c). Our particular meshing strategy means that we have trouble with certain types of connected models that overlap in Z, such as a helix. However, this can be fixed by either splitting the model's natural object hierarchy, or increasing the number of edge groups, ultimately reaching the level of a group per polygon, if needed. Figure 7.5 shows a model that our meshing algorithm considers hard. Note that although the knot was split into 200 evenly sized and space edge groups, boundaries between the groups are not visible in the final result. Figure 7.6 shows an additional example of a scene divided into edge groups, where separating drops are shown with red lines. 7.3.2  C o n s t r a i n e d V o r o n o i D i a g r a m s as a M e s h i n g R e p r e s e n t a t i o n  Voronoi diagrams have been used extensively in a number of different areas of research, so we will refer the reader to [Mulmuley 94] [O'Rourke 94] for more discussion. Basically, given a set of sample points, the Voronoi diagram divides the plane into regions where the area of a given Voronoi subdivision i is  Figure 7.5: An object our meshing strategy considers "hard". Knot model courtesy of [Scharein 98]  Chapter 7: Snow Accumulation  114  Figure 7.6: The scene is divided up by drops closer to sample point i than any other point. Voronoi edges express a location where the distance to two points are the same. Figure 7.7 shows a Voronoi diagram.  Figure 7.7: A Voronoi diagram, with points shown in black. The Voronoi diagram captures a number of different proximity relationships between points. These include the closest distance to a neighbour, the set of neighbours, point location, and responsibility area. All of these are useful. As an added bonus, the Voronoi diagram can be converted into the Delaunay triangulation (Figure 7.8) [Mulmuley 94] [O'Rourke 94] which has additional desirable properties. Assume we have a face /,, and a set of j arbitrarily chosen launch sites (see Section 7.4 to see how  Chapter  7: Snow  115  Accumulation  Figure 7.8: The Delaunay triangulation of Figure 7.7. the sites are actually chosen) points. We locate each sample point on the face, and create a Voronoi diagram V. Each region Vij has the following properties: • The region Vij contains the sampling point fij. We want each launch site to be connected to an area, so we can compute "  a r  e  " = height (given constant density)  • All of the area in Vij is closest to point fij. This allows us to throw down some points, and get back a launch site subdivision where there is a strong connection between the location of the launch site and its geometric area. • All areas of V are convex. This allows us to compute the area of Vij, necessary to compute snow height for the point. • For j partitions, there is a O(logj) [O'Rourke 94] algorithm for finding the nearest neighbour of a (new) point. This is useful for simulating snow travelling over an edge group boundary and flying into space. When snow eventually intersects a destination face, we can quickly compute the closest neighbours to the intersection point, in order to reassign the incoming mass. • The nearest neighbours of any sample point j associated with region Vij are those points that have edges on Vij. Once we've computed the initial V, we can compute this for free. This property will be useful during snow stability. If we determine that a point is unstable, and it needs to shift mass, the obvious first destinations for the new snow will be the immediate neighbours. • Voronoi diagrams can be computed incrementally. This allows us to add (and remove) additional sample points on the fly. Each added point requires a recomputation of some of the existing sub-areas, and a subsequent proportional reallocation of snow mass to the new point.  Chapter  •  7:  Snow  Accumulation  116  Finally, one of the most useful properties - incremental launch sites can be added anywhere (outside of a few Delaunay degenerate cases). This allows us to cluster new launch sites in areas of real interest, contrasting with more uniform mesh representations.  The dual of the Voronoi diagram is the Delaunay triangulation, which also has some nice properties. •  We can compute a triangulation of the surface directly from the Voronoi diagram. This gives us a set of snow faces where all launch sites /,-.,• form the subface vertices. If we wish to draw the resulting snow surface as a connection of planar polygons, this gives us a fast way of building a surface from a connected set of launch sites.  •  The Delaunay triangulation maximises the "angle fatness" over the triangulation. Basically, this means that an angle belonging to any given triangle is as large as possible. If we triangulated a different way, we might get much narrower angles, leading to very long, skinny triangles. Long, narrow triangles are harder to put a bounding box around, can cause numerical trouble, and do not represent a launch site's "neighbourhood" very well.  We use Voronoi/Delaunay support described in [Naher 96], which implements incremental computations.  7.3.3  U s i n g V o r o n o i D i a g r a m s To D e t e r m i n e N e i g h b o u r h o o d s  Our meshing representation actually uses a constrained-Voronoi algorithm that takes into account boundaries and edge conditions.  A l l launch sites in an edge group are initially placed into a single Voronoi  diagram. Because edge groups are not necessarily convex in shape, all Voronoi neighbour connections are clipped to the outline of the edge group. These clipped neighbours are not actually removed, but are just tagged as being invalid. We would ideally like each representative Voronoi area to lie entirely on a single face, providing us with nice planarity, site importance based on properties of only a single face, and local cohesion (in case the face moves). However, we also need neighbourhood information across faces in order to produce smoother resampling. So, we actually keep two Voronoi diagrams: one over all elements in an edge group, and one over all elements in a face. The face-level Voronoi diagram is used to compute launch site areas, which are clipped to the boundaries of each individual face. It is not actually stored, but created on-the-fly anytime new launch sites are added to the face. The face clipping step can be seen in Figure 7.1, or in Figure 7.9, where adjacent Voronoi areas do not match neighbours clue to the presence of a clipping edge.  Chapter  7: Snow A cc u. in ula tion  117  Figure 7.9: The constrained Voronoi diagram clips Voronoi areas to face edges. The model has 8 initial faces that tile the square.  Problems With Our Meshing Representation  Neighbour clipping can lead to problems in certain geometries, particularly in edge groups containing holes. Although holes are considered an edge group boundary and are handled properly by the clipping procedure, they can lead to unusual neighbour configurations (i.e. across the hole), which are generally removed. Figure 7.10 gives an example of one of the worst cases - a torus (invisible in the picture). This object is considered a single edge group with a large hole in the middle, however the edge clipping properly removes the interior from consideration.  7.3.4  Increasing and Decreasing Mesh Resolution  One of the nicest advantages of our meshing scheme is that we can generate points wherever we need them, leading to fairly smooth transitions between sparse and dense points. Figure 7.11 shows the top view of a very sparse initial mesh that has received additional samples around the transition zone. The centre circle is an obstacle, shown in the side view Figure 7.12. Any time a new launch site is added, the meshing representation is adjusted, potentially affecting the areas of nearby neighbours. The new site gains snow mass from any nearby neighbours that have lost area, where the snow transferred is proportional to the area loss. Thus, the new launch site starts  Chapter  7:  Snow  118  Accumulation  Figure 7.10: Snow on a torus. The results were not passed through the implicit function step (see Section 9.2) life with snow strongly related to the area it now represents. For practical concerns, we prevent launch sites from being placed some distance A , ,  t e  from any  existing sites. When this user parameter is small, it prevents round-off error from creating sites with zero representative area.  When this user parameter is large, it implicitly limits the resolution, by  preventing any change in the mesh that is below some user-selected distance. The criteria for increasing mesh density is given below. Remember that launch sites are considered for resolution change in order of importance. Adding a new point reduces a site's area, forcing it down the list. A launch site cannot increase if: • the launch area is less than a user-specified absolute minimum. • the ratio of launch area to flutter area is less than a user-specified value. • some launch sites have not yet been tested during this pass. A launch site increases resolution in the direction of neighbour n if: • all launch sites have been tested during this pass. • neighbour n is sufficiently far away such that the addition of a new launch site is at least A n s  e  from both existing sites. This enforces a minimum resolution and prevents runaway site creation.  Chapter  7:  Snow Accumulation  119  F i g u r e 7.11: Top view - an obstacle causes mesh improvement in the transition zone. The initial mesh started with only 8 launch sites.  •  the accumulation coverage difference is greater than a threshold.  The mesh resolution can be decreased incrementally. Snow mass is allocated to proportionally nearby neighbours based on new area received from the remove of the old site. A launch site cannot decrease if: •  some launch sites have not yet been tested this pass.  •  at least one neighbour n is has a sufficiently different accumulation coverage and height, expressed as a percentage. The implication here is that neighbouring sites must be both the same height and the same accumulation coverage.  A launch site is removed if: •  all launch sites have been tested this pass.  •  the user enables mesh reduction.  •  all neighbours n have sufficiently similar accumulation coverage and height.  Chapter  7: Snow  Accumulation  120  Figure 7.12: Side view - an obstacle causes mesh improvement in the transition zone. The initial mesh started with only 8 launch sites. Figure 7.13 shows a plane, with initial sites shown in white (the final mesh is shown in blue). Figure 7.14 shows the final result (all sites in white, final mesh in blue). Figure 7.15 shows all sites that were added. This particular example is trying to increase resolution under four obstacles.  7.4  Initial Particle Distribution  Each launch site must shoot a minimum number of particles to ensure some basic occlusion information depending on the projected area of the surface, where smaller areas require fewer particles. As surface size decreases, a snowflake with a constant flutter area (see Section 7.5.1) is likely to pass through more of the volume above the underlying surface. As well, errors associated with a smaller surface are less likely to be noticed. With a very large surface, a single snow-flake's flutter area does not pass through as much of the volume, and so more flakes are required to cover test the entire surface. The number of initial launch sites is user-specified as a constant times a ratio of launch site area to the flutter area. The area of each face then determines how many initial launch sites it receives, between a minimum of 1 and a user-set maximum to limit complexity. As the flutter area decreases, more initial samples are needed to properly cover the entire volume over a face.  We place initial launch sites on a face  ('hapter  7: Snow  Accumulation  Figure 7.15: The final mesh with added points shown in red.  121  Chapter  7: Snow  Accumulation  122  Figure 7.16: A mesh with random initial  Figure 7.17: A mesh with 5-time launch  launch site locations.  site optimisation. The launch sites are more evenly spaced.  gure 7.18: A mesh with near-optimal launch site optimisation. Note that the launch sites are very evenly spaced, but too many launch sites are on face edges, leading to degenerate mesh elements.  Chapter  7: Snow  Accumulation  with a pseudo-random location technique that tries to  123  partially  maximise the distance between all of a  launch site's neighbours. The algorithm randomly selects a point on the face to be the candidate for the first launch site. We use our edge group information to determine the minimum distance to an existing neighbour, potentially on a different face. We then repeat the random location step, usually from 0 to 5 times, and select the final location that maximises the minimum distance to the closest neighbour. Note that it is possible, given existing neighbour information stored in our mesh to compute the best location, i.e. the new sequentially added launch site is optimally distant from all existing neighbours. This is a desirable property, except that launch sites tended to be forced to face edges, causing thin, narrow triangles. We compromise by not picking the optimal, but instead picking the best from a small number of choices. This generally results in fairly even spacing without undesirable mesh elements (it would be possible to improve this by negatively weighting the edges). Figures 7.16, 7.17 and 7.18 show the results of changing the number of mesh location steps. Note that the random process results in variance in launch site area, even though the average is equal to the ratio of flutter area.  7.5  Simulating Snowflake Motion  Snowflake motion forms the heart of our snow allocation phase, providing the initial accumulation data about the world. Our algorithm has to balance computational feasibility with the ability to mimic or produce observable snow effects. We have to emphasise up front that we have no experimental data on how flakes of various sizes and shapes move when dropped from a significant height. In the absence of anything better we have made some assumptions about snowflake motion (normal distribution etc), and provided some parameters that animators can change. However, these distributions are built so as to be easily replaceable with real data if and when it arrives.  7.5.1  O v e r v i e w of the A l g o r i t h m  We simulate snowflake motion with a series of straight-line vectors approximating a curved path, where vector length and end position are determined with a random walk process based upon a circle of radius f , and step resolution is influenced by the importance ordering. Each straight-line vector represents r  a single "swirl" of the flake, and must be checked for surface intersection along its entire length, as  Chapter  7: Snow  Accumulation  conceptually shown in Figure  124  7.19.  F i g u r e 7.19: Side view offtake motion.  The value of f  r  approximates the magnitude of the lateral "wiggle" of a flake as it travels upwards.  At each step, the value of f  r  (i.e., the X Y length of the vector) is randomly chosen from a normal  distribution curve with a mean / f  r  M  and a standard deviation f . a  The direction of the vector of X Y length  is also chosen randomly. T o save space without too much loss of accuracy, we truncate the normal  curve to 3 standard deviations on either side of the / . M  We have no real reason to believe that flake  motion is based upon a normal curve, but it serves us as a placeholder for real experimental data. As f  r  increases, the "area of effect" of a flake widens, generally blurring occlusion boundaries and  making it less obvious where bumps and depressions came from. Figure 7.20 shows the X Y projection of several flakes originating from the same point. As fr approaches zero, flakes wiggle less and less, eventually travelling in a straight vertical line at the limit of f^ = 0 and f  a  = 0. At this point our flake motion algorithm essentially duplicates ray-casting  towards the sky, resulting in sharp occlusion boundaries and no partial occlusion. Every sample with a unoccluded path to the sky directly overhead will receive the full snow depth, while samples directly blocked from the sky will receive nothing. The only intermediate-depth areas are those resulting from a surface triangulation between adjacent launch sites. Since the "area of effect" of a flake is only directly above the launch site, and there can potentially be only a few samples per polygon, it is very easy to  Chapter  7: Snow  Accumulation  125  'random!.data' -•— 'random2.data' • 'random3.data' 'randomd.dala' - K — 'randoms, da la' 'random6.dala'  5 8  -30  X - absolute position (om)  Figure 7.20: Flake motion influenced by the random walk, projected into X Y plane. Flake origin is at [45, -31] miss overhead occluders. Although unable to produce partial occlusion effects, these raycast-style rays actually are quite useful for providing a first-cut approximation of overall snow locations. Each launch site need only fire a single particle, and move it a single time. Since the majority of computational effort during this phase is spent in repeatedly determining intersections, this produces results much faster than the full algorithm.  7.5.2  Depressions and Snow Boundary Control  As fn increases, the area of effect of flake motion also increases, leading to two different visible results: depressions in areas not blocked from the sky, and snow accumulation in occluded areas. Depressions occur in unoccluded areas when particles wander laterally, hitting distant objects underneath or on the side. As  increases, the depression area widens, as flakes from further and further  away can collide with a local obstacle. The size and depth of the depression also depends on the surface  Chapter 7: Snow Accumulation  126  area of the object potentially available to block flake motion. A shorter object thus causes less of a depression than a taller one, as depicted in Figure 7.21. The thinnest object has the least chance of blocking lateral snow motion, leading to depressions that most closely match the obstacle sky occlusion boundary. The widest object blocks more flakes, leading to depressions that extend much further from the obstacle. Mass conservation is broken in these images (and subsequent), because for illustration purposes obstacles cannot accumulate snow.  Figure 7.21: As an obstacle gets smaller, the depression area around an obstacle is reduced. Fp, = 4 cm, f = 1cm a  In areas directly occluded by an obstacle, greater flake reach implies that one or more flakes may miss the obstacle and contribute snow. Figure 7.22 shows how increased flake variability produces more and more snow directly under an obstacle.  Figure 7.22: Increased flake variability leads to more snow under an obstacle. to 4 cm to 7 cm from left to right. F = lcm  increases from 1 cm  a  The variability of the wiggle is controlled with f„, with larger values increasing the range of flake's effect. Figure 7.23 shows the effect of increasing f  a  results similar to increasing  for a constant / . Increasing variability produces M  except that effects are not so evenly distributed. This parameter can  127  potentially be thought of as an "air turbulence" factor.  Figure 7.23: As f  a  increases, the depression increases in size and depth. Here  = 4cm, j  a  increases  from 1 cm to 2 cm to 4 cm from left to right.  The vertical proximity of an obstacle to a launch site also has an influence on the snow boundary. As obstacles are positioned closer to the snow surface, flakes have a reduced chance to spread in all directions. More will be blocked, leading to a more pronounced boundary, as shown in Figure 7.24.  fe.  Figure 7.24: Flakes have less chance to disperse when obstacles are close, leading to sharper occlusion boundaries.  = 4 cm, j  a  = 1cm  Note that for illustration purposes, all of these examples use a fixed number of launch sites (i.e., no resolution increase). In most scenes, boundary areas that both hit and miss obstacles are tagged as interesting, and scheduled for a clarifying increase in mesh resolution.  7.5.3  Swirl L e n g t h  Because we must repeatedly check for intersection after each update, the distance of the swirl vectors will influence our computational speed. Swirl distance is initially set by the user to reflect the general scale of the results required. It is increased heuristically as launch sites loose importance, requiring fewer  Chapter  7: Snow  Accumulation  128  and fewer intersections. An intersection computation is valid no matter what this increment is; however the z-height will change the motion of the snowflake so that it wiggles less and may travel in a slightly different direction (as shown in Figure 7.25).  Figure 7.25: Changing a flake's Z increment test changes the flake's direction  7.5.4  C h o o s i n g Parameters for F l a k e M o t i o n  We were unable to find any experimental data describing the motion of a snowflake falling under windless conditions. This is not particularly surprising, since there are so many potential experimental variables to control (atmospheric and otherwise) adding to the challenge. As a result, we were essentially forced to guess at some reasonable parameters, helped by comparing with some of the parameterised images shown above, and real images. As a general rule, more surfaces receive snow when flake variability is increased, leading to a "denser" looking scene, where snow is literally everywhere.  7.5.5  Efficient Intersection B u c k e t i n g  After obtaining a list of flakes from a launch site, we want to move them towards the sky, continually checking to see if their movement path intersects with anything. This is quite similar to a ray-tracing computation. Unfortunately, because flakes do not move in a straight line, we must repeatedly perform  Chapter  7: Snow  Accumulation  129  the intersection computation for discrete parts of the snowflake movement.  In order to determine if  a flake is blocked, we must check it for intersection with all faces (the model), plus the top and side polygons belonging to all active snow surfaces. In order to reduce the number of potential intersection tests, we try to immediately throw away large numbers of polygons that are not located near the flake's path of motion. We start by dividing the world X Y plane into a regular grid of buckets, where each bucket contains a list of all polygons that have some portion within the bucket boundaries. This requires some small storage duplication in the case of objects that overlap multiple buckets (Figure 7.26). It is not immediately clear what the optimal bucket size should be. If buckets are too small, than it is likely that a single flake must be checked against many buckets. If buckets are too large, then a single bucket will contain more polygons to be checked. In most large models, it is usually best to increase the number of buckets while memory is still available (see Figure 10.3 for timing). object b u c k e t s  Figure 7.26: Objects are assigned to buckets within a grid  Within each bucket, we compute the minimum and maximum Z values of the polygon as it passes through the bucket bounding box, giving an active Z range p . We then insert all Z ranges in a range tree r  [Willand 85]. For a tree containing n ranges, it takes 0(log n) 2  per insert and delete, and 0(log n 2  + k)  to return a list of the k elements that overlap the query range. Once we have all intersection polygons bucketed, we can begin testing for snowflake intersection. We assume that snowflakes travel in many, small linear increments. For each linear increment, we determine  Chapter  7: Snow  130  Accumulation  the set of buckets that must be tested for potential intersection, and query a range tree against the flake's Z-range. Each query on bucket i returns a list of a list of ki overlapping ranges in time 0(log n{ + ki). 2  Because polygons may overlap different buckets', we need to prevent a single polygon from being tested multiple times. In order to do this, we concatenate all lists into a master list, tagging polygons as they are checked for intersection, and tossing out polygons that have already been previously checked in the list. This intersection method will handle our expected situation, where: • the intersection cost of totally occluded surfaces initially starts high, and quickly lowers as the distribution shoots fewer and fewer particles less and less often • most of the surfaces that receive snow are located at the top of the model, are likely to be totally unoccluded. Flakes originating from these snow piles must be tested against immediately neighbouring piles, and then generally have fast and free motion upwards through the empty sky. One problem with this scheme is that when all surfaces have shot their flakes, we must potentially rebucket all polygons in preparation for the next pass. During the very first accumulation phase, rebucketing is only needed upon completion, when snow mass is added. During later accumulation phases, object polygons carrying snow may have subdivided or coalesced, leading to an intra-phase change of height. Although it is possible to perform a selective update, we choose to rebuild the buckets after each time step. Unfortunately, it leads to some "lag" where the actual height of the snow is different than the surface bucketed for intersection testing. Section 8.7 discusses this a bit more. It is important to point out that there are a number of other alternatives to our particular intersection approach. See [Cleary 88] and [Wyvill 90] for alternative possibilities borrowed from ray-tracing.  7.5.6  W i n d and Flake Location  Wind effects on snow produce some of the most compelling and interesting scenes. At each step of the flake's travel, we can include in a wind influence based upon the flake's location and distance travelled. The "wind influence" is essentially a velocity vector for every point x,y,z in space. Although we can accept a wind velocity field of nearly arbitrary complexity, all of our examples are computed using a simple everywhere-constant field. This provides an ad-hoc first approximation. Effects such as wind shadowing, eddying, wind acceleration around obstacles and other complex effects  Chapter  7: Snow  131  Accumulation  are dependent upon the wind model, and the method used to compute it - likely a computational fluid flow package. Figure 10.23 shows the influence of the slight breeze to the right. Other wind examples can be seen in the section on stability and wind (Section 8.8).  7.6  Locating Particles in the Sky  When a launch site reaches the head of the importance queue, it shoots a batch of particles towards the sky.  Batch size is user definable, but generally within the order of 10-15 flakes. Particles originate from the launch site's snow surface, potentially reaching the sky plane unimpeded  and contributing to the growth of the parent. We use a simple bucketing and filtering scheme to allocate the successful flakes to the total mass of the sky's available snow, while ensuring that small local areas of sky do not over-contribute. This is important, since the number of particles hitting any particular area of the sky may vary dramatically depending on the complexity of the underlying surfaces. We must ensure that a large concentration offtakes (say, directly above a tree), draws the same total snow as would the sky above a sparse flat surface. Furthermore, importance ordering implies that not all launch sites shoot the same number of particles. We start by dividing the sky into a grid of constant size buckets, each of area sky . area  When a flake reaches the sky, we locate it within the appropriate bucket.  If a flake's represented area is greater than the area sky , area  we place copies of the flake in surrounding  buckets, such that each bucket contains an appropriate amount of the flake area - essentially equivalent to spreading a flake's area over the sky, as shown in Figure 7.27. Without such a filtering scheme, a single flake (representing a very large area) might land in adjacent buckets in adjacent turns; because allocated snow is based on the sum of projected flake area, a change in bucket location would greatly affect the snow allocated to nearby flakes. Sky buckets only conceptually store individual flakes. For efficiency reasons, each sky bucket actually keeps a list, sorted by launch site, that maintains the total flake area sent from that site. The flake data structure used to find and compute intersection is lost (actually, reused for the next launch site) as soon as all flakes in a group have been resolved.  Chapter  7: Snow  Accumulation  132  Figure 7.27: Locating a flake within sky buckets  7.6.1  Allocating Mass to Snow Surfaces  When the snow accumulation phase finishes, all sky buckets are allocated some mass based upon the arbitrary depth of snow desired. Each bucket 6 computes a mass per area value, based on available mass of b and the summation of all representative flake areas extending into b. An individual launch site / then receives new mass proportional to the summation of the representative area of all flakes belonging to / that hit b. The height computation is affected by snow density (at the moment this is constant for all launch sites in the model, but could potentially vary). A single launch site may receive snow from multiple buckets. Flake area filtering is done at the end of the accumulation phase, when a given launch site cannot change in area due to added or removed refinement sites. After this time step's contribution of snow has been allocated, each still-active flake now contains a summed accumulation of snow mass. This mass is added to the launch site (along with the mass from other flakes sent by this launch site this turn), resulting in an increase in height. The sky is cleared of all flakes, and we move onto the next time step. Since a launch site's accumulation pattern may change with the addition of blocking snow, it is sometimes useful to split the desired snow depth up and run the accumulation phase more than once. Depending on the time allocated for each phase, lower-importance launch sites may not get a chance to shoot particles every pass. To allow fair mass allocation to those launch sites, we keep flake information in the sky until replaced by a "fresher" set of shot particles from a new pass.  Chapter  7: Snow  Accumulation  133  Unreplaced flakes remain in their sky buckets, and on the next turn will act exactly as if they had been sent upwards from the original launch site and had successfully hit the sky. As long as these flakes are kept within their sky buckets, they will contribute as usual to the snow mass allocation of a parent launch site.  7.6.2  Writing in the Sky  In the proposal phase of this research, we made an initial assumption about the homogeneity of the incoming snow, i.e., that the sky generates snow evenly. As it turns out, this assumption can be removed quite easily, and can even be exploited to generate some interesting and useful effects. By multiplying each constant-mass cell in the sky grid by a distribution factor, we can modify, in a complex way, the underlying accumulation of snow.  2001 Figure 7.28: Cells in the sky grid are factored by this input image (black = 0.0, white = 1.0)  Figure 7.29: Non-constant allocation of snow mass to sky bucketing can be used to "write" the SIGG R A P H 2000 logo with snow Figure 7.29 shows a scene where the sky's constant distribution was multiplied by the input image  Chapter  7: Snow  Accumulation  .134  shown in Figure 7.28. The obstacle occlusion of all flakes in the scene remains unchanged, but now the sky generates uneven amounts of snow. This is the ultimate in anti-aliased writing (ie, with snowflakes!) One obvious use of this property is for artistic effect: a scene may require a thick layer of snow on the background, but a lesser layer in the foreground. Other potential uses include modifying the sky distribution to simulate melting due to sun exposure or ground heat.  7.7  Texture Generation  One of the whole reasons for this research is to produce models of snow at a scale where the thickness is visually important. We can then handle scenes that cannot be done via texturing and "just turning distant scenes white". However, this does not mean that we cannot use texturing to augment our current method. The idea behind texture generation is to automatically replace 3D surface models with a textures in instances where a 3D surface is either not noticeable or not practical. The following sections describe several ways textures are used in our current implementation.  7.7.1  Flake Dusting  In many instances, accumulated snow is not thick enough to completely obscure the underlying surface, appearing instead as a light "dusting"of flakes. This phenomena often occurs in areas of low snowfall, high instability, or on surfaces with microtexture bumps, such as tree bark (Figure 7.30) Since it is not practical to model dusting as thick 3D objects, we use already-computed snow occlusion percentages to generate procedural noise textures of the appropriate averaged dusting density. Dusting textures are semi-transparent, textured polygons oriented to float slightly in front of the original model. Figure 7.31 compares the texture dusting of a real and a computer generated sign. Another example is shown in Figure 7.5.  Computing Flake Dusting Location and Densities  The surface resampling phase (see Section 9.1) connects launch sites to produce triangulated snow surfaces, along with vertical enclosing edge surfaces. Texture processing occurs on a surface-by-surface basis, and generally results in one of the following conditions: • No texturing, surface retained. If all snow vertices are of a sufficient depth, no texture is added  Chapter  7: Snow  Accumulation  135  Figure 7.30: An example of flake dusting on tree bark. and the original resampled surface is kept. • Texturing added, surface removed. The snow surface is everywhere less than the replacement delta /, and so can be removed and replaced by a texture. • Texturing added, surface kept. The snow surface is partly less and partly more than the replacement delta I, A new texture is added to handle the thin areas, while the original surface is retained to reflect the thick areas. The resulting two surfaces may overlap or occlude in different areas, but will be properly drawn in the rendering phase. The texture replacement algorithm is expressed in algorithm 2, given some texture replacement limit I.  Textured surfaces are not sent on to the implicit surfaces phase. 7.7.2  T e x t u r e R e p l a c e m e n t for D i s t a n t Surfaces  Our current texturing method replaces very thin snow surfaces with textured equivalents. The replacement metric is measured in absolute terms, i.e., mm of snow depth, and occurs throughout the model. However, it is easy to imagine a relative replacement measure, where surfaces are replaced by textures (coverage = 100 %), in areas where the surface depth is apparently small relative to some fixed viewpoint.  Chapter  7: Snow  Accumulation  136  A p r o c e d u r e t o replace r e s a m p l e d surfaces w i t h t e x t u r e s based o n c o r n e r s k y coverage a n d snow d e p t h Procedure ReplaceTexture(surface, density) compute the 3-vertex max coverage; compute the 3-vertex average coverage; compute the 3-vertex max depth; compute the 3-vertex min depth; compute the 3-vertex average depth; / / all 3 corners are sufficiently deep if (min depth > /) density = N/A; / / not a texture return Case 1; / / The next three cases handle surfaces that are totally occluded from the sky. / / Thus, it's not immediately obvious what the density should be. / / a surface can't see the sky, and has no snow else if (max coverage < 0) and (max depth < 0) density = N/A; / / not a texture return Case 1; //a  surface can't see the sky, but has snow, all less than the limit else if (max coverage < 0) and (max depth < /) density = (average depth/1); / / remove the surface, fudge the coverage density return Case 2;  //a  surface can't see the sky, has snow, but only some corners less than the limit else if (max coverage < 0) and (min depth < /) and (max depth > /) density = MIN(1, (average depth/1)); / / keep the surface, fudge the coverage density return Case 3;  / / a surface can see the sky, and all corners are less than the limit else if (max coverage > 0) and (max depth < /) density = average coverage; / / remove the surface return Case 2; //a  surface can see the sky, and some corners are less than the limit else if (max coverage > 0) and (min depth < /) and (max depth > /) density = average coverage; / / keep the surface return Case 3; end Procedure  Algorithm 2: Determining Texture Density and Replacements  Chapter  7: Snow  Accumulation  137  Figure 7.31: (a) A real sign covered with real snow, (b) A computer generated sign covered with computer generated snow. Note how dusting density increases near the top and edges in both models. Since textured surfaces are not sent on to the implicit surfaces phase, there is a very real reduction in total polygons for each surface we can replace. Additionally, replacing distant surfaces with texturing improves the results of a given marching cubes pass (See Section 9.2) by reducing the overall volume that must be marched. The downside of texture replacement is the introduction of view dependency: textures are replaced based on apparent size, and as soon as the viewer moves that apparent size changes but the textures will not. This technique is very powerful for situations where you are willing to give up camera independence in return for a good deal of simplification in areas that will always remain too distant for the error to be observable. The replacement algorithm is standard: • for each corner of a resampled polygon, we determine the distance from the camera. • given some information about the viewpoint, camera, and eventual rendering size, we can determine the apparent height of each corner's snow depth in terms of image pixels. • if the apparent depth of all corners of a resampled surfaces are less than a user-defined cutoff (for example, 0.5 of a screen pixel), we can replace the entire surface with a texture of 100% density. An obvious possibility for this type of replacement is a dense snow-covered forest. The camera can fly around the close trees and get the full benefit of the thickness, and as long as the camera remains sufficiently far from from the distant trees, everything appears reasonable. Figure 7.32 shows a bare surface with some trees located at varying distances from the camera.  Chapter  7: Snow  Accumulation  138  Figure 7.32: The many distant branches in the scene make it a good candidate for view-dependent texture replacement. Snow was added to the surface with view-dependence texture replacement on, replacing 16,888 thick surfaces out of a total of 44,314 surfaces. The pixel replacement limit was 2.0 pixels. Figure 7.33 shows the results. Figure 7.34 shows the only the replacement textures, highlighted red for better visibility.  7.7.3  Problems with Texture Generation  Our current texture generation algorithm has some drawbacks, as described in the following sections.  T e x t u r e Interpenetration  Before the snow generation phase, we export our initial model as a set of polygonal surfaces. Snow occurs on the polygonal surfaces (not the original curved surfaces) so there will be some discrepancy between the underlying and snow surfaces. This discrepancy is not usually very visible for snow, since the thickness of the snow is usually greater than the difference between the real and approximated curves. However, since the textures have no depth, differences between real and approximated curves are more  Figure 7.34: All view-dependent textures (highlighted red for better visibility)  Chapter  7: Snow  Accumulation  140  Figure 7.35: Texture interpenetration artifacts. visible. Figure 7.35 shows the original sphere with a polygonal texturing approximation that completely encompasses the sphere in a band. Note how the textured surfaces interpenetrate the original surface, producing an unpleasing tile-like artifact on the left. There are two obvious solutions: increase the resolution of exported surfaces (at the cost of a higher polygon count), and locate textured surfaces further away from the base surface (leading to a greater chance the viewer will see "floating" artifacts).  Snow Transit and Flake Dusting Density  Ideally, flake dusting should be influenced by snow that passes over surface, and not just the skydependent snow accumulation pattern. For example, when an avalanche sweeps across an area, the flake dusting density should increase, due to snow transit residue, and settling of airborne snow accompanying motion. If the avalanching snow arrived from some height above, then flake dusting should be increased at the target site, to simulate a plastering or "splatting" effect. Neither of these possibilities are implemented in our current approach.  Chapter 8 Snow Stability-  several low angle slides. The right slide exposed the ground.  8.1  Overview of the Algorithm  The stability phase of the system is responsible for re-distributing recently accumulated snow mass into a configuration that is stable, according to surface and snow properties. It can be run at intermittent times as computational power and desired accuracy permit, usually immediately after snow accumulation. The stability phase sorts snow-covered launch sites by accumulation height, and then compares each site with all adjacent connected neighbours using a simple model of snow material properties.  Snow  from a locally unstable site is sluffed off to nearby downhill neighbours, until the initial sample is either stable, is free of snow, or cannot transfer snow due to a blocking base surface or snow obstacle. Occasionally a launch site near an edge-group boundary (qv) will have no downhill neighbours, causing  141  Chapter  S: Snow  Stability  142  snow to avalanche over an edge or drop-off. Falling snow is simulated as a small number of particles, and is intersection-traced down through the model to eventually come to rest. During stability, launch site resolution may increase in areas where snow falling from above generates a pattern that is not sufficiently represented by the existing sample density. The movement of snow down the slope may affect the stability of both uphill and downhill neighbours, so the stability algorithm performs multiple passes over the sorted list of launch sites until the overall snow motion is small or the program runs out of time.  8.2  Maintaining a List of Unstable Samples  All launch sites are initially sorted by absolute Z height plus accumulation, and placed in a list of unresolved sites u\. The list is examined in decreasing Z order, immediately resolving unstable launch sites as they are discovered. The resolution of a single launch site s may affect a number of nearby neighbours: lower sites may receive new snow from s, while the loss of snow from s may create unstable angles with previously stable higher neighbours. Affected samples also include sites receiving edge-transit snow from s, or sites newly created to improve mesh density. If not there already, all launch sites affected by s, including s, are placed in a new sorted list U2- At the completion of an entire pass through u i , the list is destroyed and replaced with 112, and the entire pass is repeated until termination. The length of u\ is not guaranteed to decrease on each pass, and in fact may increase, or undergo large fluctuations: consider a large amount of very unstable snow on a wide flat surface. On the first pass, the vast majority of interior samples are considered stable, since they are at the same height as their neighbours. The band of instability exists only at the edges, where unsupported snow avalanches off into the void. As edge sites lose mass, adjacent interior neighbours are affected, and the area of instability widens. After multiple passes, the area of instability may reach deep into the centre of the surface. Fortunately, the erosion of snow from the edges towards the centre is very physically plausible. Repeatedly passing through the list can be very slow, as a nearly stable initial situation degenerates completely before regaining stability. Furthermore, it is not immediately clear how many passes are required before the stable situation is regained. However, the obvious alternative is much worse: consider sorting the unstable list by Z, and only moving downwards in Z when all higher launch sites are stable. This implies that stable areas are  Chapter  8:  Snow  Stability  l/,3  guaranteed, known, and will not affect unstable lower areas. The problem is that this type of algorithm drives a large wave of unstable snow downwards from the top, potentially overwhelming launch sites at the bottom of the model. Moving very large volumes of snow between lower launch sites generally leads to poor results, as the unstable snow can tower ordersof-magnitude over the initial accumulation. Additionally, if the stability phase is terminated early, the top parts of the model are stable, while the lower areas are extremely chaotic. The top-to-bottom multipass method works much better: by having all parts of the model move snow simultaneously, the snow moves more naturally and there is less chance that lower sites are overwhelmed. Early termination means that as a long as at least one pass was completed, the model often looks reasonable. It may not be completely stable, but at least all launch sites have moved at least one big chunk (usually the majority) of unstable snow.  8.3  Angle of Repose  Despite the wide range of physical factors influencing real snow, for simplicity we base our stability test mainly on the angle of repose (AOR) of a particular snow type. The AOR measures the static friction of a pile of granular material, and is one of the major parameters influencing our scene. It can range [Kuroiwa 66] [McClung 93] from near 90° in fresh dendritic snow to 15° in extreme slush conditions. For a given type of snow, we use a transition curve that models the probability of stability over a range of angles around the AOR." Increasing the width of the transition curve gives a stability solution with bumpier surfaces and increased variation at snow boundary edges near the critical angle. A narrow curve generates smoother surfaces with less variation. There is no reason the curve has to be linear, or decreasing. Figure 8.1 and Figure 8.2 show two different AOR curves used to compute other figures in this thesis. The selection of curve is one of the major influences that an animator has over the overall snow appearance. The AOR is based on the relative heights of accumulated snow, and not on the fixed angle of launch sites on the underlying surface. As snow drains from one launch site to another, the A O R changes continually. This means that launch sites on very steep surfaces may still support snow if the A O R of neighbouring sites is low enough, possibly because snow is blocked from moving away. Figure 9.18 shows an example of this using water (AOR = 0°) filling a fountain basin. The basin sides are too steep to support water, so mass avalanches towards the basin bottom. As the basin fills, this downward  Chapter  8: Snow  Stability  144  -i  0  10  20  1 1 r 'repcurve.data'  30 40 50 60 X - Angle of Repose  70  80  90  0  10  20  30 40 50 60 X - Angle of Repose  70  80  90  Figure 8.1: Transition curve for angle of repose,  F i g u r e 8.2: Transition curve for angle of repose,  simulating (fairly) fresh snow.  simulating older snow.  This curve was  This curve was used to  produce Figure 10.38  used to produce Figure 10.8  movement is blocked by the rising water level. Eventually, the basin fills to the brim, leaving a stable flat surface supported by the steep sides of the bowl.  8.4  Stability Test  The actual stability test for a single iteration on launch site s can be described as follows: 1. compute A O R between s and all neighbours  lower than s  2. for each i with large enough A O R , perform an obstacle test between s and n,i. if there is a non-snow obstacle in the way, the avalanche is blocked, and the neighbour n,- is ignored (Figure 8.3  (i))  ii. if there is a vertical snow surface (an edge group boundary) in the way, there is an interpenetrating surface carrying snow between s and n;, so the avalanche is also blocked (Figure 8.3 (ii))  Chapter  8:  Snow  L{5  Stability  iii. if there is a non-vertical snow surface in the way, there is an interpenetrating surface B between s and n  u  where the interpenetrating surface could potentially receive the snow  destined forre,-.Replace n; with the closest launch site on B (Figure 8.3 (iii)) 3. for all neighbours ni still in contention, shift a delta of snow from s 4. repeat steps 1 to 3 until all there are no unstable neighbours left, or s is bare of snow Figure 8.3 illustrates some of the obstacle cases. S  S  S  B  B  (ii)  (i)  ( i i i )  Figure 8.3: An illustration of stability test obstacle cases i , ii, and iii. The obstacle test (step 2) checks to make sure that avalanche motion is not blocked by interpenetrating surfaces or snow belonging to other objects. Interpenetrating surfaces usually are the result of nearby edge groups that are under, near, or inside the edge group of the site s under test. For example, the pine tree models in Figure 10.8 have many interpenetrating needles. The more snow that is added at once, the more likely interpenetration will be a problem: note how the leftmost foreground tree has two or three long unstable spikes on the lower branches, the result of the stability algorithm "jamming". The best solution is to either reduce the amount of snow added per pass, or rework the models so they are more physically plausible (ie, in this case the underlying needle surfaces are not actually connected to the branch surfaces). If an obstacle is found, snow is blocked and forced to pile up unless there is an alternative escape direction or snow rises above the intervening obstacle. Figure 9.18 shows how blocked water rises above the level of the basin sides, transferring to the top of the basin edge, eventually overflowing into the next basin.  Chapter  8:  Snow  146  Stability  Step 2 is expensive. Practically, we achieve large speedups by allowing the user to reduce the frequency of this step - from every test, to every pass, to once per stability phase, with corresponding decreases in accuracy. The most infrequent testing is usually sufficient for models where there is little inter-object penetration, although some blocking due to rising snow will be missed. Figure 10.12 was computed using the fastest method. Figure 10.23, containing thousands of interpenetrating and closely spaced grass blades, was computed using the slowest method. Any time there is a non-snow obstacle between two adjacent neighbours, we can optionally improve the way snow builds up against the obstacle by adding a new launch sites just before the intersection point. The number of "obstacle launch sites" added per pass, and the proximity to an obstacle before a new site is added are both under user control. Figure 8.4 shows how stable snow accumulates behind a blocking obstacle, even when the underlying face is too steep to support any snow. Figure 8.5 shows the same scene, with the obstacle faces removed for illustration. Note how the sides of the snow build-up are consistent with the low angle of repose of the snow (35°). Obstacle points were automatically added along the boundary of the blocking obstacle to improve the sharpness of the transition.  Figure 8.4: Snow on a very steep surface accu-  Figure 8.5: Figure 8.4 with an invisible obstacle  mulates behind an obstacle. Angle of repose= 35°  plane.  Chapter  8.5  S: Snow  Stability  147  Moving Snow Over Edges  If an unstable launch site has no downhill neighbours, it is next to an edge. Before snow cascades over an edge into the air, we perform an intersection test with a very short vector oriented in the direction of avalanche motion. The purpose of this tiny "test stab", shown in Figure 8.6, is to ensure that there is no obstacle immediately blocking edge transit snow. If an intersection is found, then some surface or nearby snow is sufficiently close to the avalanche origin to block movement. Blocked avalanches continue to accumulate until the original launch site has enough snow to pass over the obstacle.  Figure 8.6: Test stab checks for obstacle intersections that might block edge transit snow. The right image shows an avalanche that is. prevented from moving due to an obstacle. If no intersection is found, the avalanche heads over the edge and is approximated as a few (usually < 5) avalanche particles moving on a simple projectile trajectory, based on the motion of the avalanche from the origin to the edge (Figure 8.7). The initial velocity is determined by assuming the avalanche started at the launch site at rest, and accelerates under gravity and a snow frictional parameter p^ until it reaches the edge. This is a very large simplifying assumption in that an avalanche could have been accelerating for quite a while before reaching a launch site adjacent to an edge, but is instead assumed to have started from rest, quite near the edge. Each particle is assigned a proportional amount of the avalanche mass, and is perturbed slightly in both velocity and direction to encourage spreading.  Chapter  8:  Snow  148  Stability  avalanche at rest  collisions stop horizontal motion  Figure 8.7: Downwards trajectory of avalanche flakes. Avalanche particles are tracked downwards using simple projectile motion, with the potential inclusion a wind vector. We determine collisions using the same intersection-test bucketing scheme as used for upwards-travelling flakes. If we encounter a surface that cannot support snow, we set all X Y velocity and acceleration to zero, i.e., an avalanche hits a wall, and falls straight down to the base as shown in Figure 8.7. Particles are tracked downward until they intersect a surface that supports launch sites or the top snow polygon of a surface that supports launch sites. If a particle reaches the ground plane without encountering an obstacle, it has travelled into the void and all mass associated with it is removed permanently from consideration. Once a particle reaches an intersection, we then use the properties of our Voronoi meshing to quickly determine the nearest eligible launch site on the destination surface, and contribute incoming snow to the height of the target site. Depending on user-set distance threshold parameters, new launch sites may be created if existing launch sites are not dense enough to  Chapter  8: Snow  Stability  J/,9  capture the pattern of falling snow. This particular approach has some problems: in order to properly express the pattern of falling snow from above, we wish to decrease the distance threshold to increase the number of launch sites created. However, this directly affects the launch site count, reflecting on the overall speed and complexity. We ideally wish to strike a practical balance between adding new sites and site inflation. In most of our examples, we have generally lowered the tendency for incoming particles to create new launch sites. This leads to the creation of stalagmites, where particles hit at different spots, but the large scale of the meshing lumps them together at the same spot. The problem is discussed as future work in Section 11.1.7.  8.6  Stability Termination Criteria  A single pass of the stability algorithm reaches completion when it runs out of time, when the unresolved list W i becomes empty, or when all avalanches in the last pass moved only a very small amount of snow. In most scenes, the first few passes through ui resolve a majority of the unstable snow, with subsequent passes handling smaller and smaller avalanches. Forced early termination may leave unstable areas, but all launch sites will usually have avalanched at least once. Our multi-pass approach avoids driving a large wave of snow downwards in a single pass, which leads to chaotic results on early termination.  Because of lagging rebucketing (see Section 8.7), the final "stable" solution may actually contain a  number of confused, improperly stable/unstable launch sites. Additionally, any speed-accuracy tradeoffs (such as missed obstacle testing) may have also generated some improper sites. If the stability phase completes before the alloted time expires, we do a complete rebucket, and re-run the entire stability phase with refreshed obstacle testing. A l l sites are retested for stability, in case they somehow escaped from the active list. The extra phase usually fixes a few missed sites and completes in a fraction of the first pass. Termination depends on: •  Time.  Since the stability algorithm is time-driven, all stability computation stops dead when  the allocated time expires. This mean that areas of the model have unstable snow, although the magnitude of unstable snow is usually small if all sites have had at least one pass. •  Minimal snow movement. The algorithm keeps track of the maximum height of snow moved by any launch site over a complete pass through the list. If this amount ever dips below a user-set  Chapter  8:  Snow  Stability  150  cutoff, the stability pass terminates, before continuing with a complete rebucket and check-pass, described above. The cutoff limit is usually set at 0.5% of the accumulated snow height of the previous accumulation phase. • Empty List. If the active stability list is empty, there are no more sites to resolve. The algorithm continues with a complete rebucket and check-pass before terminating.  8.7  Rebucketing Surfaces  Our flake/surface intersection method depends on allocating possible blocking surfaces to a bucketing grid and a Z range tree, in order to reduce the number of candidates tested. After the first subdivision approximation phase, some of the surfaces in the intersection bucketing will be snow surfaces, implying that they will vary in position based on changing snow height. Any time a surface moves in height, it must be updated in all appropriate Z range trees. All phases begins with a complete rebucket, allowing surfaces to start off in the correct Z range. The subsequent effect of moving snow surfaces in Z then varies, depending on the phase. In subdivision accumulation phase, there are no snow surfaces, so any points added or removed can be treated as part of the base model, and ignored. After each subsequent accumulation phase, the only major change in snow height occurs at the end, when all launch sites are given allocated sky-mass. The simultaneous rise of all sites means that it is most efficient to clean out the entire intersection bucketing scheme and rebucket every surface in all the range trees. Additional launch sites may be added during this phase, but they essentially subdivide existing snow polygons and do not significantly change in position. Thus, rebucketing changing snow polygon is not an issue in these phases. However, during the stability phase, snow surfaces are changing continuously, (and possibly dramatically) with the motion of avalanches. Each avalanche results in the raising or lowing of a number of nearby neighbours, each which must be rebucketed in the Z range trees to properly allow correct occlusion testing. We have decided that it is too costly to replace and rebucket each snow polygon exactly as it changes. Instead, we have made the implementation decision to rebucket only every so often, and put up with a lag between the real position of a snow surface and the perceived position of a snow surface as maintained in the intersection bucketing. The lag parameter is set by the user to a number that represents "the number of stability checks before a complete rebucket". The smaller the lag parameter,  Chapter  S: Snow  Stability  151  the more accurately matched real and perceived surfaces will be, to a limit of I. where the problem goes away but it takes a long time to execute. Experience has shown that the lag parameter can be kept reasonably high and still obtain good results. The fire hydrant example in Figure 10.38 had a lag value of 40,000 checks for approximately 2,400 launch sites. As long as a launch site is dealing with neighbours, the lag is not a problem. Neighbourhood angle of repose is computed on the fly based upon snow and area, without use of the intersection-bucketing. Lag begins to cause problems in several specific cases: • Obstacles that interpenetrate, as shown in Figure 8.3. A given stability test may show that there is a surface blocking neighbour-to-neighbour movement, when in reality that snow has moved away and progress is unimpeded (Figure 8.3, case ii). The lag results in a "stable" result instead of an "unstable" result. The alternative case means that some snow can leak by an interpenetrating surface that has not yet been rebucketed properly to show up as an obstacle. • Avalancheflakesfall from above and hit a surface that is no longer there. This results in particles generally landing closer to the source than expected, as they hit a phantom surface early. The phantom intersection prevents the flake from completing the proper amount of horizontal motion. • Avalanche flakes "stab" and hit or miss phantom surfaces as part of edge drop jumping (see Section 8.5). This can result in either a stable result (a stab hits a surface which is not really there, when instead it should have been allowed to move) or an unstable result (a stab misses all surfaces, when it should have hit an up-to-date surface and been prevented from moving). When rebucketing occurs, lag problems generally go away or resolve themselves as long as the launch site has not been completely removed from stability consideration. This is one of the prime reasons why we perform multiple complete snow passes - the complete rebucketing before each one usually catches the few sites that cause problems.  8.8  Wind  Wind is a major factor in the large-scale transport of snow, producing some very compelling and interesting effects. Our model handles wind in both snow accumulation and snow stability phases.  Chapter  8:  Snow  152  Stability  During snow accumulation, wind influence is easily included by modifying a flake's direction and distance by a velocity vector.  Wind velocity vectors can be approximated with a constant direction,  or much more accurately computed offline. The appearance of the wind effects is strongly related to the resolution and accuracy of the wind field.  The foreground haystack in Figure 10.23 shows the  asymmetrical accumulation effects of a very slight breeze to the left, where the wind influence is globally constant. T h e appearance of the wind effects is strongly related to the resolution and accuracy of the wind field. During stability, we widen our single-site stability test to include neighbours that are within 90° of the downwind direction. Snow transport is then dependent on the neighbour's angle with respect to the local wind vector, the duration of the wind influence, and the carrying capacity of a given wind velocity, based on [Mellor 77]. The instability vector is moved according to the rules of Section  8.4, including  obstacle testing. We use a simple heuristic to compensate for the different number of times each launch site may be stability tested. Each neighbour is given a flux maximum in grams per square cm/second, which is converted (using the simulated time of the phase) to a total mass that can potentially be moved for the launch site's given surface area. This maximum limits the amount of snow each launch site can both send and receive during a phase, ideally leaving sites with the same amount of snow, except in situations where outgoing snow is blocked (ie, against an obstacle). Each stability trial results in a small reduction in the flux maximum, where the reduction is scaled by an arbitrary "resolution" N. The idea is that if all neighbours move their wind capability N times, the flux maximum will be reached. Of course, some neighbours will reach their ./V limit faster than others, thus not contributing wind movement snow on later passes. Unfortunately, if the stability phase is terminated early, some neighbours that have not reached their N limit are prevented from moving all allowable wind transport snow. Furthermore, the evenness of wind transport is directly related to the number of nearby neighbours and their configuration with respect to the wind direction. If we happen to have an unlucky configuration of neighbours such that all neighbours are nearly perpendicular to wind direction, then a given launch site may move less snow than other nearby sites with a better neighbour configuration. Thus, snow movement works best with a very dense mesh. All of these flaws result in a wind stability algorithm that works best for obstacles, where moving snow quickly piles up against an obstacle, and stays there. It does not work well for even transport on a widely varying mesh (a field, for example).  Chapter  S: Snow  Stability  Figure 8.8 shows an example of wind and stability effects using a simple, globally constant wind vector. This particular example has some artifacts not related to wind, such as insufficiently detailed texturing most obviously seen in (c).  Figure 8.8: (a) The initial scene without snow, (b) The scene with snow, but no wind, (c) A globally constant wind blows snow against the wall. Textures have been erroneously set to float too high off the base surface, (d) Even stronger wind. Much of the snow has blown completely away.  Chapter 9 Snow Final Phases  A rime encrusted tree in a large windscoop.  This section describes some of the extra steps (outside of snow accumulation and stability) required to get finished, rendered images. In particular, it describes a smoothing step required to produce bridging. It additionally describes some miscellaneous advantages to our approach.  9.1  Surface Resampling  After snow allocation, each launch site is elevated by recently accumulated snow mass divided by the current launch site area. The polygonal top snow bounding surface is then the constrained-Delaunay triangulation of elevated launch sites, with corner vertices set to the minimum of adjacent neighbours. The resulting polygonal mesh has vertices that are either launch sites (with snow height) or vertices  155  Chapter  9:  Snow Final  Phases  156  (with average snow height). Since the constrained-Delaunay triangulation is constructed over all launch sites in a face group, adjacent neighbours on different faces are still resampled properly. Additional vertical planes are included around edge group boundaries to close the surface down to the base plane. This resampling occurs any time snow is added, and results in the snow planes placed in the intersection bucketing. In order to maintain boundary discontinuities, the resampled surface must have at least the same edge boundaries as the original surface. This is important to give the appearance that the underlying face is completely covered; otherwise it produces "shrinkage", where the edges of the underlying face are not completely obscured by covering snow. It is not clear that maintaining edge discontinuities is always a good thing: our strategy is to reduce "surface mismatch" during the implicit surface phase.  If unconnected surfaces are sufficiently close,  we can make them blend via implicit functions, instead of through an explicit resampling step across different edge groups.  9.2  Implicit Functions  Importance-ordering accumulation algorithms are surface-based, implying that snow can only accumulate on supporting objects.  To allow for unsupported snow, such as gap bridging, edge bulges and wind  cornices, we perform an additional (optional) conversion step using implicit functions.  As shown in  [Blinn 82] and others, implicit functions can be used to smoothly blend surfaces together, based purely upon distance and a blobbiness parameter. For example, Figure 9.1 shows an example of two implicit surfaces that form a bridge structure.  By manipulating our implicit function equations, we allow inter-  surface blending, allowing both bridges and a general smoothing and rounding of previously sharp edges (partially due to the re-conversion of blended surfaces back to polygons). Implicit functions have the important disadvantage of volume inflation, i.e., the apparent increase in mass due to function blending. Volume inflation is clear in Figure 9.1; as the main bodies of the sphere come closer, a bridge is formed, increasing the interior volume of the system considerably from just the sum of the two initial spheres. Volume inflation is a well known problem in implicit surfaces [Desbrun 95] [Bloomenthal 97] but most solutions are too complex to be practical for models with hundreds of thousands of surfaces. Additionally, we have the advantage that most of our surfaces are so small that completely accurate volume conservation is not a necessity.  Chapter  9:  Snow Final  Phases  157  Figure 9.2: Volume inflation must be prevented on flat surfaces; the line-like bulge shows the edge of two adjacent surfaces. Our approach is to use a different type of implicit functions for the top surface and the edge surfaces, with the goal to match the top surface as closely as possible (no inflation) but to allow mass inflation at the edges (bridge creation). We want to minimise volume inflation on the top of the snow surface to allow us to correctly generate "flat" surfaces, since otherwise connected polygons form bumps where neighbouring edges touch, and peaks where multiple vertices touch. (Figure 9.2). In order to reduce these discontinuities and apparent volume inflation at function boundaries, we use known adjacency information to shrink and clip implicit functions so that the isosurface is coincident with the polygonal top surface. A small variable-radius line generator function blends cracks between adjacent functions, and smoothes over sharp creases in the snow (Figure 9.3). This allows a very limited neighbour-to-neighbour blending, producing a small smoothing effect. It is not quite satisfactory, since in some cases the depressions are still visible. Figure 9.4 shows this sort of artifact on a single face subdivided into four mesh areas by a single, central launch site. However, gaps are often minimised sufficiently to be destroyed during mesh reduction [Garland 97] after polygonisation.  Chapter  9: Snow Final  158  Phases  b l e n d i n g t o p  s u r f a c e s  l i n e s  i s o s u r f a c e  (b)  ( a )  Figure 9.3: (a) Side view of adjacent snow volumes, (b) Side view of adjacent top and edge generator functions, with crack-filling blending lines. The "top surface" implicit function is based on a distance function from Blinn [Blinn 82]: F(x,y,z) =  T-exp(-S d -B) 2  s  where d is the distance from (x,y, z) to the closest point on the generator plane, T is a threshold scale factor (set to 1), and B is a fixed "blobbiness" constant. This function is truncated where effect is small (generally, a threshold of 0.05), giving it a fixed distance of effect. Other, more efficient functions could easily be used [Wyvill 94] (and probably should be). R can be thought of as the radius of the function, where any point (x, y, z) that is distance R from the generator surface will result in F(x,y,z) = 1. If we choose 1 as our threshold for being inside or outside the function (a common practice), this means that R defines the boundary of the surface. We want this boundary to be as closely matched to the actual polygonal boundary as possible, so we move the virtual plane representing the generator function downwards by a factor ^ R C0S  O  (dependent on the  AOR), as shown in Figure 9.5. R is variable, in that multiple values produce the same isosurface. However, the choice of R affects the rate of function falloff beyond the isosurface (larger R results in a wider influence), so we need to choose some reasonable value. We allow the user to select an influence value R  max  that gives the position  where the F(x,y, z) influence goes to zero (actually truncated to a small positive amount). Beyond this distance, a generator function has no effect. Given a maximum influence R  max  and fixed B, we can then  solve the inverse of F(x,y, z) to find the appropriate R that drops off to reach zero influence at  R xma  In order to prevent wide scale volume inflation where surfaces join, the effect of each top generator  Chapter  9: Snow Final  159  Phases  Figure 9.4: A single face, converted to implicit functions. There are four mesh areas - the lower right corner show some crack-filling artifacts. function is restricted to the volume directly above the snow volume. This means the function extends upwards from the underlying snow polygon, and so may blend with other functions (not part of the same edge group) located directly above it. This produces bridging effects when, for example, a tree branch comes very close to the ground or when nearby surfaces interpenetrate. However, there is no blending with adjacent neighbours - not exactly what we want. In addition, the uphill edge of a generator surface has a small gap that varies depending on the neighbour. We fill this gap with a line generator (using the equation above) that is displaced at the average, interpolated R of the two neighbours. The "edge surface" function is one-sided, extending away from the edge surface. The strength and shape of this function directly influences the creation of bridges, including determining how close surfaces must be before bridges are formed. As a default, all edge surfaces create a small bulge centred at c = | of the height of the edge, with a radius of b  max  bi  m n  — 0.3 of the height at the maximum and a radius of  — 0.05 height at the top and base of the edge (Figure 9.6). This implies (very approximately) that  2 adjacent street grates separated with by 1 inch space will just start to form a connecting bridge under about l i inches of snow, given our default parameters. Several examples show the edge function in operation. Figure 9.4 shows the edge-rounding of a very simple surface. Figure 9.7 and Figure 9.8 show a comparison of a pure polygonal and an implicit function version of the same snow surfaces. Note that in the implicit function version, surfaces are slightly smoother, and an unsupported bridge has formed.  Chapter  9: Snow Final  160  Phases  Figure 9.5: The top generator is displaced from the actual snow top plane. The function is clipped to 0 outside the X Y extents of the plane.  edge surface top surface h = height  snow ^ c = bulge center  t  Figure 9.6: Parameters controlling bridge creation  Figure 9.7: The polygonal surface, before bridg-  Figure 9.8: Implicit surfaces representation al-  ing.  lows bridging  Chapter  9:  Snow Final  162  Phases  Figure 9.9: Implicit surfaces can form unsupported bridging and clumping. Figure 9.9 shows an example where snow on many closely-spaced pine needles has formed unsupported bridges and clumps. We have to admit our bridge parameters are very ad-hoc. However, given data for actual snow (perhaps a side picture of a bridge) we believe this approach is flexible enough to allow us to match the bridge closest touch point, and thickness of bridge. Note this method is geometric and scale-independent, and does not use any properties of the snow (such as the mass above it). One final step is required before we can render our snow scenes in Alias|Wavefront. We must somehow convert our implicit surfaces into a form importable by a commercial package. We have chosen to sample our volumes using the marching cubes algorithm [Wyvill 86] [Lorensen 87], in particular a version that polygonalises the isosurface in 0(n ) space [Watt 92], This generates a (potentially huge) 2  set of triangular facets, subsequently reduced with a decimation algorithm (Section9.2.3). The separate implicit step allows repeated re-execution at various resolutions.  9.2.1  Advantages and Disadvantages of Implicit Functions  The implicit surface generation step has some advantages, but many problems. We consider it more as post-processing step to get apparent visual results, rather than than a central step in our algorithm.  Chapter  9: Snow Final  163  Phases  It provides the ability, with the proper manipulation of parameters, to obtain a fairly wide range of results from a single set of snow surfaces. In particular, the blobby factor can be increased beyond optimal (where optimal is the value selected to match the flat surface) causing large, rounded, very smooth surfaces. Manipulation of this parameter can produce results that are quite different from those generated in the previous accumulation and stability phases - essentially invalidating physical plausibility of those phases. Implicit functions potentially allow us to add animal tracks, wind ripples, and other patterns to snow surfaces by "stamping" the snow surface with appropriately scaled negative functions, although in practice these effects are better done in the main modelling software. By interrupting the pipeline before the implicit function step in Section 9.2, we obtain polygonal results with no bridging or smoothing effects and a much lower polygon count. These compact intermediate results are more appropriate for scene setup and real-time viewing, and may actually be sufficient for the final image. Figure  7.1,  9.18,  10.23 and  10.42 were computed without the smoothing step.  As well, intermediate polygonal results can be used as the underlying model for a completely new snow accumulation run, producing the effect of true snow layers (Section 9.4.1).  9.2.2 The  Cornices  edge function already allows overhangs. We would like to extend this to allow cornices, which  are overhangs that are related in size and shape to the wind parameters. This requires mapping wind parameters to an appropriate (c, b  , 6  max  m j n  ) edge parameter, where stronger wind implies larger  b  ,  max  and c is likely placed at the topmost part of the edge surface. Unfortunately, we did not implement this, although we do not foresee too many difficulties. One complication: some cornices have sharp "curls" that may be hard to match in shape exactly, forcing us to use approximated shape.  9.2.3  Vertex Reduction  Since the work involved in generating snow surfaces is proportional to the number of input surfaces in the base model, it makes sense to reduce the size of the input through mesh reduction. We reduced the number of faces in the input model with Qslim [Garland 97], a popular (and implemented) mesh reduction algorithm. This program requires a maximum allowable error for vertex  Chapter  9:  Snow Final  16 J,  Phases  replacement, which we conservatively chose as 1% of the length of the smallest edge in the model. Limiting maximum vertex error provided reasonable initial compression with little or no apparent visual difference compared to an uncompressed version. Mesh reduction is even more useful as a post-process to the marching cubes step. Function polygonisation can generate easily generate millions of faces, especially with an 0(n ) (n = march step) 2  algorithm generally limited only by disk space. Execution time to march most models was generally under 10 minutes, on an SGI Onyx. We found that it was most efficient to run the marching cubes algorithm at a very high resolution (N > 800), and then use mesh reduction to provide a model at the detail we needed. We used the percentage of the march cube size as an error metric. Unfortunately, the mesh reduction implementation has some memory problems that enforce a hard limit on the number of surfaces that can be simultaneously reduced. We got around this problem by cutting the polygonisation into smaller strips, passing each strip to the reducer, and then reassembling strips. This worked well in almost all cases, with a few exceptions - see Figure 10.22, where join artifacts can be seen laterally in the foreground. Post process mesh reduction resulted in very significant shrinkage of the number of surfaces exported to the commercial rendering program. Sample reduction ratios are shown in Table 9.6, although for admittedly some of our largest models. Reduction depends on the number of large flat surfaces in the scene, and the number of strips (see above - fewer strips allows better reduction). Table 9.6: Typical Mesh Reduction Figure Figure 10.40 Figure 10.8 Figure 10.28  9.3  Marched Surfaces 6,232,544 7,409,067 5,695,876  Reduced Surfaces 372,225 775,625 752,457  Error Metric (% of march edge) 30% 5% 2%  Surface Rendering  It is outside of the scope of this thesis to perform snowpack rendering per se, although this is certainly an interesting and underexplored problem. However, in order to generate images we need something that looks like snow; we settled on the rendering abilities of Alias| Wavefront,shader, although other programs  Chapter  9:  Snow Final  165  Phases  will have different shaders and abilities. We attempted a (very rough) mapping of snow parameters to shader parameters, simulating fairly fresh snow. The  primary coloration comes from a solid texture of type "Rock" (an Alias predefined texture)  that simulates a mixing of two grains of different color. We used a grain size of 1 m m (medium-large grain size) where one grain is a slightly wetter snow than the other. We very approximately mapped computer R G B values R to midrange wavelength (R to 700pm, G to 535pm and B to 500pm),  1  and  roughly interpolated values from the reflectance chart in Figure 2.25 to arrive an approximate "fresh snow" R G B = [0.86 , 0.92, 0.93]. This is slightly lower than the true curve, which tended to produce whiteouts unless the lighting was just right.  The wetter grain was mapped to the "dry on top, wet  below" curve, and given an R G B = [0.79, 0.85, 0.89]. The shader was given a diffuse reflectance of 0.85, a translucence limit of 0.5 cm (completely opaque beyond this limit), and a fractal bump mapping with a maximal surface displacement of 0.5 cm. Figure 9.10 shows the snow material, and how the translucence properties allow partial visibility of objects buried in the snow, generally to a depth of 0.5 cm. The scale of the entire snowman is quite small - less than 30cm. Figure 9.11 is a closeup of the material. Snow is very sensitive to changes in the lighting, and the scene appearance depends greatly on the illumination. For example, Figure 10.23 shows sharp polygonal edge artifacts in the foreground. However, in Figure 10.24, the identical model under different lighting conditions, the polygonal edges are still present but not visible.  9.4 9.4.1  Extras Creating Snow Layers  Although snow layers are beyond the scope of this thesis, it is possible to create them in ad hoc fashion by using the results of one snow creation pass as input to the next snow creation pass. At  the moment, it is somewhat time consuming, since it requires the user to combine both the  existing model and the new layer, before re-exporting it to the start of the snow process. However, there 'Much better ways of doing this include mapping the wavelengths to the peaks of R G B phosphors in the output monitor, or even better, integrating the snow reflectance curve into a color space dependent on observer sensitivity (CIE X Y Z space), and then mapping this into the RGB phosphor space. This coidd be interesting, since it might be that snow's colours cannot actually be fully represented on a computer monitor. However, for our purposes we take some large shortcuts - justified somewhat by the flatness of the wavelength-dependence curve shown in Figure 2.25  Chapter  9:  Snow Final  Phases  166  Figure 9.10: Snow material, showing translucence  Figure 9.11: A closeup of the snow material, show-  and texturing.  ing translucence and texturing.  is nothing preventing this step from being done automatically, allowing many thin layers to pile up on top of each other. This sort of layering implies that once a layer has been put down it cannot be changed in any way, since it becomes part of the base scene. Therefore, there are no inter-layer interactions. Figure 9.12, 9.13, 9.14 and 9.15 shows such a layering sequence. Note how the resulting model overhangs from the base surface. Figure 9.16 highlights the various layers (blue = first, green = second, red = third) in the final image. The first and second layers have a very high angle of repose (88°), while the last layer has a lower repose (75°). We are not particularly trying to match the real image, except in approximate final height of snow.  Chapter  9:  Snow Final  Phases  167  snow.  Chapter  9:  Snow Final  Phases  168  Figure 9.15: The mailbox with a third layer of  Figure 9.16: Color coding shows the layers of snow.  snow.  Blue is the base layer, green is the second layer, and red is the top layer.  Chapter  9: Snow Final  Phases  9.4.2  Tracks and Surface P a t t e r n s  169  We initially wanted to add macro-roughness effects to snow surfaces, effectively acting as a displacement map to the already computed implicit surfaces. These roughness effects are simulated by "stamping" negative implicit surfaces on existing implicit surfaces, potentially based on the underlying properties of the snow. The blending nature of implicit surfaces prevents some of the edge discontinuities common with texture mapping. Figure 9.17 how a snow surface can be stamped with multiple patterns. The original (flat) implicit surface was modified with a noise function on top of a ring.  Figure 9.17: Snow surfaces can be modified with additive patterns, simulating tracks and wind-ripples. Although this has the potential to generate tracks (one of our stated goals), in practice it is not very efficient. First, we actually need a displacement map of depth, which is hard to obtain for complex real tracks (such as Figure 4.12) because dark areas do not necessarily equate with deep areas. With a 3D model of the track-producing object, we can generate an accurate displacement map by rendering for depth. However, in this situation it is easier to add the tracks as a post-processing step in Alias]Wavefront, using either displacement maps or the "warp" feature [Alias—Wavefront 96]. Both  no  Chapter 9: Snow Final Phases  features do something similar, providing better control and nicer user interfaces. Additionally, neither technique requires the very high marching cubes resolution needs to capture small track details. The main advantage of performing track stamping within our program is the ability to scale pattern amplitude based on some parameter (ie, wind speed, amount of avalanched snow, etc). We did not implement this ability, although we believe it is fairly easy to do using our framework.  9.4.3  Water and R a i n  By setting AOR = 0° and flake-flutter f = 0 the basic snow algorithm can also simulate the accumulation r  of water, from the sky as rain or elsewhere.  Figure 9.18 shows an example of an empty fountain  slowly filling up with water. Only the patch of sky shown as a red square has any mass to contribute, approximating how water appears at a spout, fills the first basin, and overflows to lower basins. As might be expected, the time to converge to a stable solution is dependent on the angle of repose, with AOR = 0° the worst possible case. We gave Figure 9.18 40 hours to complete for approximately 32,000 launch sites - an order of magnitude more time than our most complex models.  Figure 9.18: Snow stability algorithms can also be used to simulate water accumulation. Water from the red patch fills the first basin before overflowing into subsequent basins.  Chapter  9:  Snow Final. Phases  By changing a few simple parameters, the basic snow algorithm can also simulate the accumulation of rainfall, including puddles and damp spots. A "damp spot" is anywhere rain initially fell before draining off. This corresponds to choosing a single flake dusting texture of appropriate coloration. Since the depth of fallen rain is much less than the depth of fallen snow, we can't rely on errors being "buried" as much, nor can we tolerate stability jamming. The rain algorithm allows us to, for example, produce a small raindrop on every (cupped) leaf on a tree without having to drop millions of particles from the sky. This is a fairly interesting benefit and is worth further study and exploration.  9.4.4  M u l t i p l e O u t p u t Formats  Since our final output is a generic set of polygons, we aren't restricted to Alias|Wavefront as the only possible rendering engine. In fact, the "obj" format is either directly supported or easily convertible to almost every other popular graphics 3D format. As an example, Figure 9.19 shows a screen capture of an SGI "ivview" window displaying a snowcovered knot in Inventor [Manual 94] format. The knot can be rotated and viewed from any direction.  Figure 9.19: An Inventor window that allows a snow-covered torus to be rotated in real-time.  Chapter  10  Results  A real tree. Compare this with the bare tree in Figure 10.12  This chapter serves several distinct purposes. In the first section of this chapter, we recap the visual goals outlined in Chapter 4 and stated as our criteria for visual success. For each previously stated goal, we point to one or more images somewhere in the thesis that provide an example of how these goals were met. Next, we recap some of the desired goals of the implementation and algorithm, and again cross-reference with figures that illustrate how these goals were met. In the second section, we provide some simple timing numbers to give the reader an idea of the bottlenecks and timing dependencies of the approach. Finally, we provide a number of complete examples using our entire snow pipeline.  172  Chapter 10: Results  10.1  173  Validation of the Visual Appearance of Snow  Validation of snow-covered scenes is hard, in that snow observed outdoors is the result of uncontrollable and unknown environmental factors. Creating artificial snow is beyond our capabilities as a graphics lab, so instead we restrict validation to observation, asking the question: "does our algorithm produce phenomena and/or effects that are observable in nature?" However, we were able to perform a few simple experiments to show that our snow stability algorithms are at least plausible. We substituted sifted flour for snow, to improve controllability and show that our algorithms work for materials other than snow. Figure 10.1 shows a side-by-side comparison of real and computer generated scenes. Figure 7.31 shows an additional side-by-side validation image of flake dusting.  Figure 10.1: A real flour-covered scene (a) and a computer generated scene (b) compared to show that our stability algorithms are at least plausible. Our experimental setup was fairly ad-hoc: despite our best efforts, flour was distributed unevenly around the base of the real sphere.  10.1.1  Validation of Visual Properties  The following Table 10.7 restates desirable properties of the visual results of any snow algorithm, and cross-references with a select ion of figures and notes that demonstrate that goal. Unfortunately, we are hard pressed to validate our results much more than this. Real snow-covered scenes are very hard to experimentally quantify, and there are few available metrics that can be used  174  Chapter 10: Results  to compare real and computer generated scenes.  There is little previous computer graphic work on  snow that can be used as a baseline. At a minimum, we claim that our algorithm is able to at least produce some of the phenomena observable in real snow, but we cannot claim that we can mimic physical properties. At this time, we must rely upon the viewer to make a decision on the success or failure of our algorithms.  Properties  Table 10.7: Recap of Snow Appearance Goals (Section 4.2) D e s c r i b e d I n Some R e p r e s e n t a t i v e F i g u r e s  Snow location. Initial snow location. Partial coverage. Partial snow amount. Conservation of snow mass. Surface bridging. Cornices and overhangs.  4.2.1 4.2.1 4.2.1 4.2.1 4.2.1 4.2.1 4.2.1  Snow stability. Instability. Snow restabilization. Snow Surface Appearance Curved surfaces. Uneven local snow distribution. Tracks. Snow Under the Influence of Wind Initial distribution. Surface redistribution.  4.2.2 4.2.2 4.2.2 4.2.3 4.2.3 4.2.3 4.2.3 4.2.4 4.2.4 4.2.4  10.1.2  10.14 (under bench), most others 7.1 7.21, 7.23 7.5(less snow on bottom of knot) 10.8 (foreground tree) 9.15 (overhang) cornice not directly supported 8.4 10.1 (snow under ball) 7.5 most, due to sampling noise 9.17, but otherwise not directly supported 10.23 8.8  V a l i d a t i o n of Desirable A l g o r i t h m P r o p e r t i e s  The following Table 10.8 recaps some of the desirable properties of the implementation and operation of the snow algorithm, and cross-references them with illustrative figures or notes.  10.2  Timing and Complexity  Figure 10.2 gives a feel for the progression of the subdivision phase in a situation where the goal is to create as many subdivisions as possible. Four models (the bush shown in Figure 7.1, the hydrant shown in Figure 10.38, the gazebo shown in Figure 10.8, and cylinder shown in Figure 7.22) were each given 15  Chapter  10:  175  Results  Table 10.8: Recap of Desirable Algorithm Properties (Section 4.3) Properties Accumulation of snow on occluded surfaces. Surface Representation Multiple level of detail Proportional accumulation Spatially adaptive meshing Temporally adaptive meshing Surface smoothing Snow Stability Partial stability Non-determinism Stability Forces Level of detail  Described In 4.3.1 4.3.2 9.2 4.3.3 4.3.3 4.3.3 4.3.3 4.3.3  Angle of repose Snow cohesion Snow frictional properties Surface frictional properties  4.3.3 4.3.3 4.3.3 4.3.3  Shock loading Obstacle forces  4.3.3 4.3.3  Some Representative Figures  10.23 (haystack vs field) 10.23 (field needs less work) 10.40 (near tree vsflatfield) 7.1 (less sampling under obstacles) 9.8 10.26(not all snow falls off) due to width angle of repose curve multi-resolution mesh gives multi-resolution stability 9.18 not used in the model not used in the model 10.40 (low friction on roof causes distant snowfall, in this model) not achieved by our model 8.4  minutes to execute on a 250MHz SGI Onyx. Each model used 30x30 intersection bucketing, and 60x60 sky bucketing, shooting 10 flakes per flake group. Of additional consideration is the number of launch sites per launch-site eligible faces (eligible faces: cylinder = 8, bush = 128, hydrant = 2224, gazebo = 28050). In models with a very small number of faces (i.e. the cylinder with 8) there are a very large number of sites per face (> 1000), and the bottleneck is the time it takes to recompute subdivision areas within a single face. In the larger models, new launch sites are much more spread out among faces, so that the number of sites per face remains low. As as result, the larger models actually perform better - the very complex gazebo model created more subdivisions in the same time than the simple cylinder model. We suspect that subdivision creation is approximately an 0(log n) algorithm, which is the time 2  it takes to interactively insert and delete from the Voronoi diagram. Figure 10.3 shows that increasing the number of intersection buckets improves the number of subdivisions created in afixed15 minutes period (other conditions as above). However, Figure 10.4 shows that the rate of improvement slows as buckets reduce in size, likely because of the greater duplicated  Chapter  10:  Results  176  faces covering more than one bucket. Figure 10.5 shows the effect of increasing the number of flakes in a flake group (intersection buckets = 60x60, sky buckets = 60x60). Clearly adding more flakes reduces the total number of launch sites produced, but the reduction is not linear, i.e., a ten times increase in flakes produces only a 4-5 times reduction in sites. Figure 10.6 shows the number of sites created for a varying number of sky buckets (intersection buckets = 60x60, flakes = 10 per group). From this, it appears that the number of sky buckets does not influence computational speed in as long as memory is not an issue.  Chapter  10:  177  Results  'bush.data'.f^'hydrant.data'' 'gazebo.aata' 'cylinQer.data' 15*loa(xi'15*IOQ(x) — •  .a'' 35000  £  3  20000  300  400 500 X - Time (seconds) 15 minute max  600  F i g u r e 10.2: The total number of created launch sites depends on the number of launch sites per face.  Chapter  10:  Results  Figure 10.3: More total launch sites can be created with more intersection buckets.  Figure 10.4: As more intersection buckets are added, the benefit slowly drops off.  178  Chapter  10:  Results  179  40000  10000  10  20  30  90  100  Figure 10.5: As the number of flakes per group increases, fewer total sites can be created. 'skybuck.data' -•— I  3 1  3»0»  20000  Figure 10.6: The number of sky buckets is not a primary influence on launch site creation.  Chapter  10.3  10:  180  Results  Snow Images Generated Using the Complete Model  This section contains a compilation of snow images. It shows the addition to snow to a wide selection of different models. Figure 10.7 shows an initial gazebo, with a snow-covered version shown in Figure 10.8. Figure 10.9 shows an older version of the gazebo model, without the trees. Another view (with different snowfall amounts and properties) is shown in Figure 10.10. Figures 10.12 to 10.19 are from the animation "Santa's Suitcase". They show how the snowcovered model can be viewed from many directions and easily used for animation.  The bare model contains  hundreds of thousands of polygons, and is a good example of how our approach handles very complex scenes. Figure 10.11 shows the bare yard. Figures 10.21 and  10.22 demonstrate snow with several different angles of repose. The original  image is shown in Figure 10.20. Note how adding snow makes the original dark and bland bare grass image much more interesting and visually appealing. The faint dark lines on the foreground snow of Figure 10.22 are artifacts produced by the surface simplification algorithm. Figures 10.23 and  10.24 show examples of the multi-scale ability of the algorithm. Our approach  handles both individual blades of hay, and the very large field surrounding the haystack. This particular model has a very large amount of base surface interpenetration, where an individual polygon may pass through many other polygons. Figure 10.25 and Figure 10.26 show a gnarled tree model covered in snow.  This model is not  particularly detailed - many of the branches are prism-shaped, instead of smoothly rounded, and as result snow is unable to cling to certain areas. Figures 10.27 and 10.28 show the before and after of the inside of an roofless, open air cloister. Figures 10.30 and background.  10.29 show the before and after of computer generated sign in front of a real  The initial model has many problems, including foreground grass that floats about the  ground plane. Figures 10.31 and  10.32 show the before and after snow-covering of the Silicon Graphics logo.  Figures 10.33 and  10.34 show an example of how non-constant snow allocation from the sky (see  Section 7.6.2) can be used to write on surfaces. Figure 10.33 shows the writing from a distance. The closeup in Figure 10.34 shows how natural the patches of snow look. Figure 10.35 shows a bare hydrant. Figure 10.36 illustrates initial mesh. The model before stability  181  Chapter 10: Results  is shown in Figure 10.37, while the final result in shown in Figure 10.38. Figure 10.39 shows an initial house model. This scene is quite poorly conditioned (see Section 11.2.1) where many normals are pointing in improper directions. This leads to artifacts in the result (Figure 10.40), especially around the leftmost, lower floor railings. Figure 10.41 shows the mesh used to construct the results. Note the change in level-of-detail, especially under the trees.  Figure 10.7: Bare gazebo.  Figure 10.9: An older gazebo model without the  Figure 10.10: Another view (with different snowfall  t  amounts and properties).  r e e s  Chapter  10:  183  Results  Figure  10.11:  Santa's bare yard.  Figure 10.13: Frame 200  Figure 10.15: Frame 400  Chapter  10:  Results  186  Figure 10.17: Frame 550  187  Chapter 10: Results  Figure 10.19: Frame 665  Figure 10.20: Bare grass.  Figure 10.21: Snowcovered grass.  Chapter  189  10: Results  Figure 10.22: Snowcovered grass, different properties.  Chapter  10:  190  Results  Figure 10.24: The same haystack scene at sunset.  Chapter 10: Results  Chapter  10:  Results  Figure 10.28: The snowcovered cloister.  Chapter  10:  Results  193  Figure 10.30: The snowcovered sign  Figure 10.32: The snowcovered SGI logo.  Figure 10.34: snowflakes.  A closeup of the same scene, showing "anti-aliasing" due to  Figure 10.35: The bare hydrant model.  Figure 10.36: The initial hydrant mesh.  Chapter  10:  197  Results  Figure 10.37: Unstable, snow-covered hydrant.  Figure 10.38: The stable result.  Chapter 10: Results  F i g u r e 10.40: The snowcovered house. This scene is a good example of a model with poor surface normal organization.  198  Figure 10.41: The house mesh. Note how the mesh density changed under the trees.  Chapter 10: Results  10.4  200  Oil Painting Effects  Several popular image viewing programs have the ability to post-process images to add an "oil-painting effect" [Adobe 95] that essentially smoothes and clumps areas of similar color. This does not really have much to do with snow, except that we rather like the appearance of some of these images, since they tend to smooth out some small discontinuities and errors. Figure 10.43 shows Figure 10.20 with an oil painting effect.  Figure 10.44 gives the result (same  initial image as Figure 10.22). Figure 10.42 shows a very simple black and white street lamp scene. The lights accumulate snow (a flaw).  Figure 10.42: A snowcovered street light, oil painting effect.  Figure 10.44: Snow covered grass, oil painting effect.  Chapter 11 Future Work  More research on snow is needed. Much, much more. O h yeah.  This section describes areas of future work that augment and build upon our solutions.  11.1  Snow Properties  Snow is one of the most complex naturally occurring materials on the planet, often existing and interacting in solid, liquid and gaseous phases simultaneously. Its behaviour is very hard to predict, even on a large scale where hydrology and avalanche science models exists. When we first proposed to work on snow, we had plans to construct a very detailed model of the snow pack, encompassing as many physical variables as possible. Unfortunately, we ran right into the biggest roadblock of the thesis: snow scenes are generally extremely complicated, consisting of hundreds of thousands of snow surfaces on top of an already very large and complicated bare scene. In order to compute reasonable results in a reasonable time frame, we first needed to come up with some way of addressing the speed and complexity issues. The result was our importance ordering particle scheme, which we consider the prime contribution of this  202  Chapter 1.1: Future Work  203  thesis. However, as a side effect of changing thesis direction, many snow properties fell by the wayside. We'd like to add more of these back in to our model, following the original intent of the thesis. Some of the extremely important physical properties and effects we have ignored include snow compression and packing, layers, slab avalanches, snow creep, snow pack metamorphosis, melting, snowpack reflectance properties, and solar effects. In Chapter 2 we discussed many of the properties important to the formation and evolution of real snow. We briefly summarise how future versions could include these properties into our overall approach. We are confident that the work done to create an efficient, multiresolution, importance-based framework will allow us to improve the physical model without too much effort.  11.1.1  Precipitation  The various types of precipitation that result in snowfall generally operate on a very large scale. They can be ignored in many cases, since the scale of most models is too small to receive much variation in snow depth due to the environmental effects of mountains, weather fronts, etc. However, varying precipitation effects can actually be simulated through use of the sky-bucketing (Section 7.6.2). For example, assume we had a model of a mountain range and a map of current wide-scale precipitation (possibly obtained from a weather broadcast). By modifying sky mass by the precipitation map, we can produce snowfall patterns that reflect the current weather conditions (independent of altitude effects). This can be done without any modification to the current program.  11.1.2  Density  Our current implementation has hooks for snow density computation, in that each launch site actually computes height based on a density value and a mass, rather than just a mass. However, we do not currently vary density across the mesh, instead setting it to a constant value for a given snow type. In future implementations, we would like to increase density in areas where snow arrives from above, thus compressing depth and reducing the stalagmite problem. Our avalanche model could also instability in areas of densely packed snow. We believe density (and subsequent settling and compaction) to be one of the most important missing properties of our model, primarily because different snow types can have an order-of-magnitude variation in value, greatly affecting the apparent height and strength.  Chapter  11:  Future  Work  204  Density also affects some of the reflection properties of snow, which could be simulated through careful measurement and mapping to a shader.  11.1.3  G r a i n Shape  The variation of external temperature plays an important role in generating ice crystals of differing shapes. Although we have no data on this, we expect that grain shape would affect the fluttering motion of individual flakes - the greater surface area of plate-shaped flakes would likely be affected by microturbulence more than spherical or needle-shaped flakes. With appropriate experimental data, these effects can be included in our program by mapping flake type to the flake parameters described in Section 7.5.  11.1.4  G r a i n Size  Grain size affects the rendering properties of snow, which could be simulated through a careful shader mapping. We expect (although again we have no data) that grain size affects also affects flutter parameters, since larger particles have both more surface area and more mass.  11.1.5  Liquid Water  Content  Liquid water content plays a very important role in snow stability, in that it partially determines the angle of repose of various snow types. In a sense, we have encompassed this by providing a user-specifiable angle of repose curve. We could potentially improve things by having a given snow type automatically generate a predefined AOR curve, although there are also advantages to giving the animator direct control.  11.1.6  Snow T e m p e r a t u r e and Snow M e t a m o r p h i s m  We believe that our method is flexible enough to capture individual snow types, given appropriate experimental data. Snow metamorphism is the transition between those types, influenced in a major way by the temperature gradient. At minimum, it is possible to recompute the entire snow simulation with interpolated snow parameters, resulting in the apparent transition of snow between types. It is also possible to modify snow parameters during the simulation, but this is not quite sufficient, since it results in changing the entire snowpack, instead of just the most recent top layer.  Chapter  11.1.7  11: Future  Work  205  Surface Smoothness  Avalanching real snow distributes snow in a much wider and more complex cloud that we currently model with our few particles. Our model tends to cluster edge-fallen snow in a much smaller area than in real life, leading to snow stalagmites, such as those near the foreground wall in Figure 10.12. We would like to use a more detailed snow-cloud dispersal model to smooth out these artifacts. Related to this is the idea that we might want multiple angles of repose for multiple areas of our scene. At the moment, we use a single global angle of repose, since all snow is considered homogeneous. But by introducing even a few small changes (ie, snow that powdered out over an edge is less cohesive than slightly settled ground snow) we introduce snow property variation that both adds modelling complexity, and allows us to improve smoothness allowing spatially uneven stability. As a general observation, our snow surfaces are not as smooth as those that occur in nature. This is due in part to insufficient allocated time, discontinuities in the snow representation, poor base model creation and structuring, limitations of the implicit step, and the avalanching problem described above. Although it is possible to artificially enforce smoothness, through surface simplification, high-tolerance mesh reduction, or some other averaging, we have not done so in this thesis. We feel that by improving some parts of our algorithm (especially edge-fallen snow), we can improve the overall "gentle curve" look of our scenes.  11.1.8  Stability Jamming  Adding arbitrary amounts of snow to the scene can lead to a problem we term "stability jamming". If the amount of added snow is large, new snow surfaces rise a considerable distance above their base surfaces, potentially interpenetrating other snow surfaces located on objects above. These cases are usually handled by our obstacle-finding stability tests, but occasionally an interpenetrating surface will be termed as "stable", usually due to a particular configuration of surfaces and snow. An example of this can be seen in Figure 10.8, where the foreground left pine tree has several tall snow spikes hidden in the branches. This illustrates the amount of snow that was added in a single pass, and shows how these particular needles have "jammed" and are unable to move, even though they are clearly unstable. Stability jamming goes away as the user increases the number of accumulation and stability phases, and decreases the snow allocated to each. As a snowpile slowly rises up, any obstacle above begins to shadow it more and more, leading to a slowing of accumulation that eventually reaches zero at the base of the obstacle above.  Chapter  11.1.9  11:  Future  Work  206  Stability Non-Termination  Our stability algorithm does not always terminate, given unlimited time. Occasionally, it reaches a state where all sites are stable except for a very small number (usually < 10) that cannot be resolved properly, even with up-to-date rebucketing and no speed shortcuts for obstacle testing. Unresolvable sites are generally caused through round-off error or degenerate base surfaces. These sites are often related to stability jamming, above.  11.2  Issues with the Base Model  Many of our biggest problems originate from the underlying base models themselves. These models have usually been taken from free repositories, and were creating by a wide variety of authors using a wide variety of programs and techniques. As a result, models have all sorts of interpenetrations, shortcuts, and errors that can confuse our algorithm. For example Figure 6.2 shows a closeup of the tree model used many times in the thesis. Many of the "needles" are not connected to the branch, or each other, making it hard for snow to clump where they join. Model errors like this are surprisingly common, and caused us much grief.  11.2.1  Inconsistent N o r m a l s  Our algorithm assumes that objects consist of "thick" materials with the normals correctly pointing outwards from each surface. Only surfaces with normals pointing upwards can receive snow: downwards facing normals can't retain snow, although they block it and reflect it. "Thin" surfaces (represented by only one layer of polygons) are usable, as long as normals are consistent. For example, a sphere should have all surface normals pointing away from the centre. We found that one of our biggest problems was obtaining models that had normals that met our assumptions. Almost all of the the models taken from the W W W had no rhyme nor reason to the orientation of normals, and in some cases were aggressively disorganised. Since many of the models were very large, correcting normals by hand was impractical. We ended up building a helper raytracing tool that attempted to orient normals in a consistent direction based on the number of surface intersections. The tool gets it right most of the time, but there are many cases where the normal is flipped, leading some major visual artifacts. Figure 10.40 was one  Chapter  11: Future  Work  207  of the most poorly structured models. Even after normal flipping, there are many, many surfaces with inconsistent normals - note gaps in the snow cover on the lower leftmost porch. We don't really consider this a flaw of the snowing algorithm, although it is of immense practical importance. Surface normals are used for a wide variety of other algorithms, and any shortcuts in model generation will cause problems for many other techniques as well.  11.2.2  P o l y g o n s vs. C u r v e d Surfaces  We snow on a polygonal representation of the underlying model surface, which may be constructed out of non-polygonal primitives. When snow surfaces are loaded back into the original scene, there may be some mismatching, depending on how close the polygonal representation of the model matches the true surface. In most cases, mismatches are obscured by a layer of snow, and are thus rendered nearly invisible given sufficient snow. However, these mismatches are much more visible with very small amounts of snow - in particular the textures used to simulate flake dusting. Textures generated via the dusting process were meant to lay flat against the polygonal snow surface. If they are too close to the surface, they cut through the true curved representation, as illustrated in Figure 7.35. If they are moved further away from the surface to prevent interpenetration, they start to interfere with the object's silhouette, as shown on the blades of grass in Figure 10.22. The immediate solution is to increase the accuracy of the polygonization, and vary the distance textures float above the surface. The problem here is that the a priori accuracy is not known beforehand, so in the worst case a scene might require several iterations before a sufficient approximation is determined.  11.2.3  Weatherproofing Models  It is quite common for digital artists to take many time-saving shortcuts during the creation of 3D models. One of the most prevalent ways to save time and energy is to build 3D models that are valid only from a few viewing directions. The idea here is that there is no point in modelling surfaces that will never be seen by the camera. Figure 10.40 shows such a model. Unfortunately, without back blocking surfaces to prevent it, snow will happily accumulate inside the model, potentially filling it. Depending on the model, this extra snow accumulation ranges from completely hidden to visible (ie, through one of the windows in the front) In all cases, it results in  Chapter  11:  Future  Work  208  additional, unnecessary computation time. Even complete models must be fairly tightly constructed. Since every surface shoots particles, it is likely that surfaces near cracks or poor joins will eventually let particles escape, allowing snow to creep in where it is unwanted. We don't consider this a flaw in the algorithm at all - instead, it is a reflection of what happens in the real world. If you leave the door open, the snow and rain are going to come inside!  11.3  Rendering and Timing  T i m i n g results are not fully applicable to our importance ordering scheme, as models are usually allocated a running time convenient to the user. However, the timing bottleneck of snow as a useful effect is the rendering phase, which is outside the scope of our current work. and  Large models such as Figure 10.8  10.12 were given overnight for snow accumulation, yet required weeks to raytrace animations of  several hundred frames. Rendering is aggravated by aliasing in moving scenes - such as the distant, tiny, white snowpatches resting on distant, tiny, dark needles shown in Figure 10.12. We are interested in physically realistic, multi-resolution snow shaders or rendering models that are fast and accurate. Improving rendering times would provide our model with wider applicability.  11.4 Our  Moving Objects  initial problem statement restricted us to static images, however most of the tools are in place to  both create time-lapse animations and handle moving objects. Unfortunately, the easiest way to do this is to generate a number of snow models, one for each timestep, which immediately leads to large bottlenecks in the rendering phases (where human intervention is required). We could improve things by allowing snow models to pass into the rendering pipeline automatically, thus reducing the size of data created. Our stability algorithm does not allow time-lapse avalanches. However, in time-lapse photography, it is most likely any pictures taken over a period of hours would miss any very short duration avalanches occurring during that time.  Chapter 11: Future Work  11.5  209  Nature Rubs Our Nose In It  It is important to realize that even after categorising and understanding snow's physical properties, Nature still has the ability to make us go: "Wow! How did that happen?". Figure 11.1 shows snow clinging to the underside of a ski pole pointing straight upwards, violating all of the thesis' angle of repose assumptions. Figure 11.2 shows an example of complex wind formations called sastrugi, that are likely inexpressible with any current simulation method.  Figure 11.1: Snow can cling to the underside of objects. It is the author's personal belief that the unpredictability and complexity of snow is what makes it so beautiful. A simple walk to the park in a snowstorm turns into a voyage of discovery, as snow interacts with the everyday city to form fantastic shapes and strange scenes. We expect that even after much additional research, we (or other authors) will never be able to completely model all the complexity and beauty of snow using a computer. And in some ways - that makes us glad.  Chapter 11: Future Work  210  F i g u r e 11,2:  Sastrugi (wind formations).  Chapter 12 Conclusions  12.1  Conclusions  In this thesis, we have described some basic properties of snow, and hopefully given the reader an understanding and appreciation of the complexities of this material. Snow has been sparsely researched in computer graphics despite its commonplace occurrence, and this alone makes it an interesting and worthwhile topic. Our particular contribution to the field is a new algorithm for the creation of snow covered images, using a novel particle location scheme that allows surfaces to independently control the sampling effort needed to determine accumulation. Separability of surface accumulation produces many useful side effects, including importance ordering, adaptive refinement, reasonable degradation upon early termination, and greater control over the final result. Our placement algorithm allows us simulate effects such as accumulation under obstacles, flake dusting, wind, rain, and "snow-writing".  211  Chapter  12:  Conclusions  212  Efficiently computing multi-resolution accumulation patterns across very large outdoor scenes turned out to be the key problem of the thesis, overwhelming many other snow-related concerns. As a result, we consider our importance ordering particle scheme to be the main contribution to both the thesis and the existing body of research. We have also presented a simple but effective model of snow stability that handles avalanches, edgetransit snow, obstacles supporting and blocking snow, materials other than snow, and mass transport due to wind. Our algorithm handles nearly arbitrary surface configurations, with support for snowbridges and cornices. Integration with commercial software allows us to snow upon existing models in a variety of formats using various levels of detail, providing greater flexibility, power, and ease of use. Finally, we have shown that our approach is able to handle large, complex outdoor scenes consisting of hundreds of thousands of surfaces. It is our hope that this work will open up an entire new season to computer graphics, and will stimulate other researchers to explore the natural, glorious beauties of winter.  Figure 12.1: Before.  Figure 12.2: After.  Bibliography  [3DCafe 00]  3DCafe. 3D CAFE, W W W address: http://www.3dcafe.com, Pictures Multimedia, March 20 2000.  [Adobe 95]  Adobe. Adobe Photoshop 3.0 User Guide. Mountain View, C A , 1995.  Platinum  Adobe Systems Incorporated,  [Alias—Wavefront 96] Alias—Wavefront. Alias 8.5 Overview. Silicon Graphics Limited, Mountain View, C A , 1996. [Appel 86]  Appel, A . " A n efficient program for many body simulation". SI AM nal on Scientific  [Avalon 00]  Avalon.  and Statistical  Avalon:  Computing,  The Original  Public  Jour-  6(1):85-103, 1986. 3D Archive,  W W W address:  http://avalon.viewpoint.com, Viewpoint, March 20 2000. [Bain 00]  Bain, D . G . The Geo-Images Project, W W W address: h t t p : / / W W W Geography.Berkeley.EDU/GeoImages.html, Department of Geography, University of California at Berkely, March 20 2000.  [Barnes 86]  Barnes, J . and Hut, P. " A hierarchical O(nlgn) force calculation algorithm". Nature, 324:446-449, 1986.  [Blinn 82]  Blinn, J . " A Generalization of Algebraic Surface Drawing". A CM actions  on Graphics,  Trans-  l(3):235-256, 1982.  [Bloomenthal 97]  Bloomenthal, J . Introduction to Implicit Reading, Massachusetts, 1997.  Surfaces.  [Cleary 88]  Cleary, J . and Wyvill, G . "Analysis of an Algorithm for Fast Ray Tracing using Uniform Space Subdivision". Visual Computer, 4(2):65-83, 1988.  [Cohen 93]  Cohen, M . and Wallace, J . Radiosity  and Realistic  Morgan Kaufmann,  Image Synthesis.  Aca-  demic Press Professional, Toronto, Ontario, 1993. [Colbeck 80]  Colbeck, S. Dynamics 1980.  [Colbeck 82]  Colbeck, S.  Growth  of Snow and Ice Masses. Academic Press, Toronto, of Faceted  Crystals  in a Snow  Cover.  U . S . Cold  Regions Research and Engineering Laboratory, Hanover, N.H., 1982. [Colbeck 90]  Colbeck, S., Akitaya, E . , Armstrong, R., Gubler, H . , Lafeuille, J . , Lied, K., McClung, D . , and Morris, E . International  Classification  for  Seasonal  Snow on the Ground. International Commission on Snow and Ice (IAHS), World Data Center A for Glaciology, U . of Colorado, Boulder, C O , U S A , 1990. 213  Bibliography  [Daffern 92]  21/,  Daffern, T . Avalanche  Safety for Skiers and Climbers.  Rocky Mountain  Books, Calgary, Alberta, 1992. [Desbrun 95]  Desbrun, M . and Gascuel, M . - P . "Animating Soft Substances with Implicit Surfaces".  Computer  Graphics  (Proc.  SIGGRAPH),  29:287-290,  August 1995. [Dorsey 96]  Dorsey, J . , Pederson, H . , and Hanrahan, P. "Flow and Changes in A p pearance". Computer  Graphics  (Proc. SIGGRAPH),  30:411-420, August  1996. [Eberly 97]  Eberly, D .  MAGIC:  An Object-Oriented  W W W address: http://www.cs.unc.edu/ of North Carolina, Aug 15 1997.  Library  for Image  Analysis,  eberly/index.html, University  [Fearing 00]  Fearing, P. "Computer Modelling of Fallen Snow". Computer (Proc. SIGGRAPH), 34:to appear, July 2000.  [Fournier 89]  Fournier, A . "The Modelling of Natural Phenomena". In Proceedings of Graphics Interface, pages 191-202. Canadian Information Processing Society, 1989.  [Garland 97]  Garland, M . and Heckbert, P. "Surface Simplification Using Quadric Error Metrics".  Computer  Graphics  (Proc.  SIGGRAPH),  Graphics  pages 209-216,  August 1997. [Hanrahan 93]  Hanrahan, P. and Krueger, W . "Reflection from Layered Surfaces due to Subsurface Scattering". Computer  Graphics  (Proc. SIGGRAPH),  27:165-  174, August 1993. [Hsu 95]  Hsu, S.-c. and Wong, T.-t. "Simulating Dust Accumulation". IEEE puter  [Kajiya 84]  Graphics  (Proc. SIGGRAPH),  Com-  15(l):18-22, January 1995.  Kajiya, J . and Herzen, B . "Ray Tracing Volume Densities". Graphics  [Kass 90]  and Applications,  18(3):165-174, July  Computer  1984.  Kass, M . and Miller, G . "Rapid, Stable Fluid Dynamics for Computer Graphics". Computer  Graphics  (Proc. SIGGRAPH),  24(4):49-55, August  1990. [Krueger 88]  Krueger, W . "Intensity Fluctuations and Natural Texturing". Graphics  [Kuroiwa 66]  (Proc. SIGGRAPH),  Computer  22(4):213-220, August 1988.  Kuroiwa, D . , Mizuno, Y . , and Takeuchi, M . "Micrometrical Properties of Snow". In International  Conference  on Low Temperature  Science  (Physics  of Snow and Ice), volume 1, Part II, pages 722-751. Institute for Low Temperature Science, Aug 1966.  [Li 93]  L i , X . and Moshell, M . "Modeling Soil: Realtime Dynamic Models for Soil Slippage and Manipulation".  27:361-368, August 1993.  Computer  Graphics  (Proc.  SIGGRAPH),  Bibliography  [Lorensen 87]  215  Lorensen, W . and Cline, H . "Marching Cubes: A High Resolution 3D Surface Construction Algorithm". Computer  Graphics  (Proc.  SIGGRAPH),  21(4):163-169, July 1987. [Luciani 95]  Luciani, A . , Habibi, A . , and Manzotti, E . " A Multi-Scale Physical Model of Granular Materials". In Proceedings of Graphics Interface, pages 136137. Canadian Information Processing Society, 1995.  [Magono 66]  Magono, C. and Lee, C . "Meteorological Classification of Natural Snow Crystals". Journal  of the Faculty  of Science (Hokkaido  poro, Japan), 2(4):321-335, November 1966. [Manual 94]  Manual, O. I. R. Open Inventor  C++ Reference Manual.  University,  Sap-  Addison-Wesley  Publishing Co, Don Mills, Ont, 1994. [McClung 93]  McClung, D . and Schaerer, P. The Avalanche taineers, Seattle, Washington, 1993.  [Mellor 77]  Mellor, M . "Engineering Properties of Snow". 19(81):15-66, 1977.  [Milenkovic 96]  Milenkovic, V . "Position-Based Physics: Simulating the Motion of Many Highly Interacting Spheres and Polyhedra". Computer Graphics (Proc. SIGGRAPH), 30:129-136, August 1996.  [Monet 91]  Monet, C .  Wheatstacks,  Snow Effect,  Morning.  Handbook. Journal  The Moun-  of Glaciology,  Painting: oil on canvas,  J . Paul Getty Museum, Los Angeles, 1891. [Mulmuley 94] [Murray 94]  Mulmuley, K . Computational  domized Algorithms.  Geometry:  An Introduction  Through  Ran-  Prentice Hall, Englewood Cliffs, N J , 1994.  Murray, J . and VanRyper, W .  Encylopedia  of Graphics  File  Formats.  O'Reilly and Associates, Sebastopol, C A , 1994. [Musgrave 89]  Musgrave, F . , Kolb, C , and Mace, R. of Eroded Fractal Terrains".  Computer  "The Synthesis and Rendering Graphics  (Proc.  SIGGRAPH),  23(3):41-50, July 1989. [Naher 96]  Naher, S. and Uhrig, C . The LEDA  User Manual,  Version  R 3.3. Fach-  bereich Mathematik und Informatik, Martin-Luther Universitat HalleWittenberg, Halle, Germany, 1996. [Nishita 96]  Nishita, T., Dobashi, Y . , and Nakamae, E . "Display of Clouds Taking into Account Multiple Anisotropic Scattering and Sky Light". Computer Graphics  [Nishita 97]  (Proc. SIGGRAPH),  30:379-386, August 1996.  Nishita, T., Iwasaki, H . , Dobashi, Y . , and Nakamei, E . " A Modeling and Rendering Method for Snow by Using Metaballs". In Proc. EUROGRAPHICS, volume 16. European Association for Computer Graphics, 1997.  Bibliography  216  [O'Rourke 94]  O'Rourke, J . Computational Geometry New York, New York, 1994.  [Piatt 88]  Piatt, J . and Barr, A . "Constraint Methods for Flexible Models". puter Graphics  [Premoze 99]  [Reeves 83]  (Proc. SIGGRAPH),  in C. Cambridge University Press, Com-  22(4):279-288, August 1988.  Premoze, S., Thompson, W . , and Shirley, P. "Geospecific Rendering of Alpine Terrain". In Eurographics Rendering Workshop. European Association for Computer Graphics, June 1999. Reeves, W .  "Particle Systems - A Technique for Modelling a Class of  Fuzzy Objects". Computer  Graphics  (Proc. SIGGRAPH),  17(3):359-376,  1983. [Reeves 85]  Reeves, W . and Blau, R. "Approximate and Probabilistic Algorithms for Shading and Rendering Structured Particle Systems". Computer Graphics (Proc. SIGGRAPH), 19(3):313-322, July 1985.  [Robinson 99]  Robinson, D . Northern  Hemisphere  Snow Cover  Charts.  National Snow  and Ice Data Center, http://wwwnsidc.colorado.edu/NSIDC/EDUCATION/SNOW /snow_Robinson.html, As of Dec 1, 1999. [Scharein 98]  Scharein, R. G . Interactive Topological Drawing. P h D Thesis, Department of Computer Science, The University of British Columbia, 1998.  [Shifflett 97]  Shifflett, T . "Rhythm and Hues Studios". Personal communication, January 1997.  [Shinya 92]  Shinya, M . and Fournier, A . "Stochastic Motion - Motion Under the Influence of W i n d " . In Proc. EUROGRAPHICS, pages 119-128. European Association for Computer Graphics, 1992,  [Sims 90]  Sims, K . "Particle Animation and Rendering Using Data Parallel Computation". Computer  Graphics  (Proc. SIGGRAPH),  24(4):405-413, August  1990. [Stora 99]  Stora, D . , Agliati, P.-O., Cani, M.-P., and Gascuel, J.-D. "Animating Lava Flows". In Proceedings of Graphics Interface, pages 203-210. Canadian Information Processing Society, 1999.  [Sumner 98]  Sumner, R., O'Brien, J . , and Hodgins, J . "Animating Sand, M u d and Snow". In Proceedings of Graphics Interface, pages 125-132. Canadian Information Processing Society, 1998.  [Szeliski 92]  Szeliski, R. and Tonnesen, D . "Surface Modeling with Oriented Particle Systems". Computer  Graphics  (Proc.  SIGGRAPH),  26(2):185-194, July  1992. [Terzopoulos 88]  Terzopoulos, D . and Fleischer, K .  "Modeling Inelastic Deformation:  Viscoelasticity, Plasticity, Fracture".  GRAPH),  22(4):269-278, August 1988.  Computer  Graphics  (Proc.  SIG-  Bibliography  217  [Terzopoulos 89]  Terzopoulos, D . , Piatt, J . , and Fleischer, K . "Heating and Melting Deformable Models". In Proceedings of Graphics Interface, pages 219-226. Canadian Information Processing Society, 1989.  [Upadhyay 95]  Upadhyay, D . Cold Climate Toronto, 1995.  [Warren 97]  Warren, M .  Grand  Hydrometeorology.  Challenge-scale  John Wiley and Sons,  Applications:  Cosmology,  WWW  address: http://www.acl.lanl.gov/Applications/astronomy/warren.html, Los Alamos National Laboratory, April 15 1997. [Watt 92]  Watt, A . and Watt, M . Advanced  Animation  and Rendering  Techniques.  Addison-Wesley Publishing, Don Mills, Ontario, 1992. [Willand 85]  Willand, D . "New Data Structures for Orthogonal Queries". SIAM nal of Computing, 14(l):232-253, 1985.  [Wyvill 86]  Wyvill, B . , McPheeters, C , and Wyvill, G . "Data Structures for Soft Objects". The Visual Computer, 2(2):227-234, 1986.  [Wyvill 90]  Wyvill, G . and Trotman, A . "Ray-Tracing Soft Objects". In Graphics International'90, pages 469-475, 1990.  [Wyvill 94]  Wyvill, B . "Explicating Implicit Surfaces". In Proceedings of Graphics Interface, pages 165-173. Canadian Information Processing Society, 1994.  [Zhao 00]  Zhao, User  T.  and  Interface  Overmars, Toolkit  M.  XF0RMS  for X, W W W  Library:  address:  oldenburg.de/handbuch/oldoku/xforms/forms.html,  A  Jour-  Computer  Graphical  http://koko.offis.uni-  May 5 2000.  A note about references containing W W W addresses: Because of the fleeting nature of the Web, it is important to realize that the content and location of these references may change at any time. Listed references contain a "last-valid" date, which indicates the last date the page was used (this is not a last-modified date). Note that the author entry for a Web reference may actually be the page maintainer, instead of the creator of the content.  Appendix A Table of Variables  Variables in boldface are generally accepted by the ICSI. Other variables are introduced for the purpose of this thesis.  Name  Description  P  snow density, measured in symbol for grain shape (see Figure 2.5) ice grain size, measured in mm liquid water content of snow, measured as a % by volume snow impurities, measured as a % by weight compressive, tensile and shear snow strength, measured in P a hardness index snow surface temperature, measured in ° C direction of snow instability radius of an individual swirl of a flake mean of the radius of a flake swirl standard deviation of the radius of a flake swirl A Voronoi diagram on face i, with a sub-area around point j centre of implicit surface bridge bulge maximum radius of implicit surface bridge bulge minimum radius of implicit surface bridge bulge minimum distance between an existing launch site and any new site the Z range of a plane within an intersection bucket area of a single sky bucket list of unresolved stability sites updated list of unresolved stability sites parameter of snow kinetic friction parameter controlling top implicit surface radius blobbiness parameter for top implicit surface maximum radius of effect, top implicit surface implicit function strength for a given 3D point  F E  9 J S R T fr fa Vij C  bmax bmin A site Pr area Ui U  2  Uk  R B Rmax  F{x,y,z)  Table A.9: Description of Variables  218  Appendix B List of Parameters  Table B.10 contains a selected list of the most important parameters used by the implementation to control the appearance of the snow. To help the reader, the final column indicates the general result if a parameter is increased. Name  Description  I n c r e a s i n g Affects..  buckets skybuckets starting-ratio  number of intersection buckets number of sky buckets ratio of flake area to starting area  max-coverage-diff  subdivide if coverage diff exceeds this  min-neigh-distance subdivide-time initial-depth storm-time stability-resample  subdiv stop if neighbours are this close time for initial subdivision phase subdivision snow depth time for next step of snow rebucketing rate during stability  m i n / m ax-edge-fl akes texture-replacement snow-density flake-mu flake-sigma  number of flakes used over an edge replace very thin snow with textures snow density mean of the radius of a flake swirl standard deviation of the radius of a flake swirl centre of implicit surface bridge bulge maximum radius of implicit surface bridge bulge points at a snow type with a given A . O . R angle of repose curve  runs faster, more memory runs faster, more memory denser mesh, runs slower more memory stops sooner, sparser mesh less subdivision. stops sooner, sparser mesh denser mesh, more shots deeper snow denser mesh, more shots surfaces lag from actual position better avalanche spreading more textures, surface floating less height for same mass wider depressions wider depressions  if-c if-b-max/min snow-type A . O . R curve  Table B.10: List of Parameters  219  bulge overhangs at the top more volume inflation N/A snow clings to more surfaces  Glossary  Underlined terms are directly related to physical snow, as opposed to the implementation. angle of repose - The limit of slope angle beyond which a particular type of snow becomes unstable, and highly likely to move. angle of kinetic friction - The slope angle beyond which a particular type of moving snow will be likely to continue moving. angle of static friction - Another term for angle of repose. avalanche - A n implementation entity representing a delta of unstable snow that moves from launch site to launch site. avalanche flake - A n implementation entity representing unstable snow travelling over a drop and falling through the air. average porosity - A measure of snow density, expressing the ratio of air in the snow compared to ice. a-axis - One of three axis located on the basal plane on a crystal, each separated by 120°. Growth along these axis is symmetric, and will produce flat, plate-like ice crystals. basal plane - One of the axis of crystal growth. c-axis - A n axis of ice crystal growth that pierces the basal plane in a single place. Growth along the c-axis produces long, needle-like ice crystals. compressive stress - Caused by an applied force (usually gravity, but possibly an object), acting on the snowpack perpendicular to the surface. condensation nuclei - Tiny foreign particles of dust or soil floating in the atmosphere. As air cools, these particles form a centre for free water molecules to condense around. constructive metamorphism - A method of snow metamorphism. When the temperature gradient is high, there is a subsequent rise in water vapour pressure (within the top-bottom layers of the snowpack), leading to crystal growth near the top layer. This makes crystals more angular and faceted. Constructive metamorphism occurs mainly when the temperature gradient is large (i.e. the surface layer is cold). convection - Describes a wider-scale movement of air masses. The sun warms air, causing it to rise and possibly drop precipitation. This effect is quite local, and is usually negligible compared to the other major types of air lifting. cornice - Large mass of wind-deposited snow often found on the lee side of ridges and crests. 220  221  Glossary  crown - The line of slab fracture, denoting the uppermost plane of tensile failure. The crown is usually at right angles to the ground plane. cyclonic lifting - Describes a very large-scale movement of air masses, causing an overall lifting of moisture and subsequent precipitation. The general motion of the atmosphere causes a cyclical rotation around centres of low pressure. As air enters the central low pressure zone, it pushes up existing air, causing a slow lifting. Cyclonic lifting can extend across an extremely wide area (an entire low pressure zone, as seen on a weather map). density - Snow density is measured as the mass of snow per cubic meter. depth hoar - Large (generally fragile and cup-shaped) crystals caused by recrystalization in the snow. Conditions for depth hoar general involve large temperature gradients (i.e. very cold air temperatures). destructive metamorphism - A method of snow metamorphism. Water pressure differences between different grain curvatures slowly round off crystal branches, leading to an initial decrease in individual grain size. After that, water vapour moving between layers transfers mass from small particles to large particles. Destructive metamorphism occurs mainly when the temperature gradient is small (i.e., the surface layer is warm). dew point - A t the dew point, the pressure of the air is such that free water vapour molecules condense into liquid water droplets. drop - A n implementation entity representing the boundary of an edge group, where travelling snow flies over the edge edge group - A n implementation entity representing a group of connected launch areas surrounded by drops. elastic deformation - A n object is subject to a elastic deformation if it returns to its original or reference shape after removal of the deforming force. face - A n implementation entity representing a model surface. firn snow - Very dense, old snow that has undergone many cycles of thawing and freezing. flake - A n implementation entity representing a particle sent towards the sky. flutter area - Value  used as a metric to compare rough sky area coverage of  flakes. freezing nuclei - Tiny foreign particles of dust or soil floating in the atmosphere. These particles have special molecular properties that encourage the formation of ice crystals around their surfaces. Freezing nuclei are much rarer than condensation nuclei, and their activeness depends on the temperature (a lower temperature means than more nuclei meet the right crystalization conditions) frontal lifting - Describes a large-scale movement of air masses, leading to lifting and precipitation. Warm front lifting occurs when warm air rides up over top  222  Glossary  of colder, denser air. Similarly, a cold front occurs when cold air ploughs underneath a mass of warm air. Cold fronts are generally more local than warm fronts, producing lifting and precipitation over a narrower area. graupel - When snow crystals are heavily affected by riming in the atmosphere, they can gain so much additional ice that their original structure is lost. These types of grains are called graupel. hoarfrost - Another name for surface hoar (q.v.) ice glaze - A layer of rime plastered on a snow slope, causing a weak sliding layer, ice grains - Crystals of frozen water. implicit surface - A surface that is represented with an equation F(x) = 0. In our particular application, all implicit surfaces are three-dimensional, i.e, F(x,y,z) = 0. laminar flow - Non-turbulent wind flow. launch area - A n implementation entity representing the area controlled or represented by a single launch site. launch site - A n implementation entity representing the base unit of snow accumulation sampling. lee - The side of an obstacle protected from the wind. melt-freeze cluster - Globules of ice crystals formed during periods of warm days and cold nights. Melting crystals release free water, which surrounds the grains and subsequently freezes around them during the night, metamorphism (in snow) - Generally describes the two competing processes (constructive and destructive) that work to simultaneously strengthen and weaken the snowpack. model - A n implementation entity, representing the set of polygons that collectively form the surface eligible for snow. natural phenomena - Very loosely defined as objects created by natural means. See [Fournier 89] for a more detailed descriptions of natural phenomena modelling as it relates to computer graphics. necks - The connection between two sintered grains of ice. The neck grows due to the motion of water vapour. As neck size increases, the snowpack strength increases. object - A n implementation entity representing a group of related surfaces in the original object. orographic lifting - Describes the large scale motion of an atmospheric front. Horizontally moving air hits a mountain range, and is forced to rise parallel to the slope. When frontal direction and mountain slope are perpendicular, air rises the fastest. perforated crust - Bumpy snow surface caused by uneven distribution of solar radiation, leading to preferential melting on a small scale.  223  Glossary  plasticity - The ability to retain a shape caused by deformation due to pressure, pore spaces - Pore spaces are the empty air in between ice particles. In general. the greater the density, the smaller and fewer the pore spaces, recrystalization - Another name for constructive metamorphism (q.v.) rime - a buildup of ice crystals on the windward side of obstacles. Supercooled water droplets, carried by the wind, hit and freeze onto obstructions, sometimes growing into long feathers riming - A method of atmospheric crystal growth. As a crystal moves about in the air, collides with water droplets that subsequently freeze onto the original crystal. In especially turbulent air, the crystal may gain so much additional ice that its original structure is unrecognisable, rolling - A method of snow transport due to wind, occurring very close (1 mm) to the surface. Rolling is the minute motion of individual crystals pushing and resettling against each other across the snowpack. rounding - Another name for destructive metamorphism (q.v.) saltation - A method of mass transport due to wind. Ice crystals are grabbed by the wind and lifted in a short jump, usually within 1-10 cm of the surface. Each jump disturbs new crystals, which may make their own jumps. saturation point - The point, for a given air temperature, where the atmosphere can no longer absorb free water molecules. As the temperature decreases, the saturation point decreases. Thus, cold air can contain less moisture than warm air. scour holes - Holes in the snow caused by turbulent air. These are usually found on the windward side of objects, where air hits the obstacle, picks up snow, and then swirls around behind. Some of this snow may be found deposited on the lee side of the obstacle. settlement - Vertical compression in the snowpack due to the applied stress of gravity. Settlement leads to a denser, harder snowpack. shear stress - Caused by an applied force (usually gravity, but possibly an object), acting on the snowpack parallel to the surface. Shear stress can cause the snowpack to fracture along a plane of failure. - Single Instruction Multiple Data. SIMD is a parallel processing paradigm where many different processors operate in lockstep.  SIMD  sintering - A process of snow metamorphism. When the temperature gradient is low, water vapour has a tendency to move from convex areas to concave areas, building up ice necks between adjoining ice grains. This strengthens the snow pack. snowpack - A general descriptive term for all the snow that has accumulated upon a surface. snow polygon - A n implementation entity, representing either the top or the side of an accumulation of snow upon a surface.  Glossary  snow placement - One of the three main phases in our overall snow handling algorithm. Snow placement is primarily concerned with distributing a given mass of snow about the world, depending on surface orientation, surface occlusion, conservation of mass, and wind. The snow placement phase does not determine if a particular mass of snow can remain physically stable (stationary) on a given surface. Snow placement answers the question: "How much falls where?" snow stability - Another of the three main phases in our overall snow handling algorithm. Snow stability determines if the snowpack is locally or global stable, where stable is defined as stationary. Stability is influenced by such factors as slope angle, mass of deposited snow, cohesiveness of snow, etc. Once a given mass of snow is determined as unsupportable, the snow stability phase is responsible for reallocating it to reach stable state. Snow stability answers the questions "How much snow can a surface support?" and "Where does the unsupported snow go?" snow rendering - Another of the three main phases in our overall snow handling algorithm. Snow rendering is responsible for converting static information about the deposited, reallocated snowpack into a final format visible to humans. This requires generating surfaces of the appropriate shape, colour, and reflective properties, and using these to produce images. snow density - Measures the mass of snow per cubic meter. subdivision mesh - Another name for an edge group. sun cups - Large depressions in the snow surface, generally occurring late in the season. sun crust - A thin icy crust, caused by melting and subsequent refreezing of the snowpack. surface hoar - A layer of large, delicate crystals that form on snow surfaces due to a process of water vapour condensation. Surface hoar layers are structurally very weak. suspension - A method of mass transport due to wind. When wind speeds are high enough, snow may be picked up and carried long distances, travelling well above above the snowpack. Snow is suspended in the turbulent air. temperature gradient - The temperature gradient is a measurement of the change in snow temperature over the change in measured distance. It is expressed in degrees C / m , and is conventionally measured in the direction of increasing temperature. Because snowpack temperatures can change laterally, technically this is a vector (although it is commonly used only for depth) tensile stress - Caused by an applied force acting to pull apart the snow. Usually caused when one part of the snow has a tendency to move faster or easier than another part, leading to a tensile fracture between the two. viscoelasticity - Has both viscous properties (i.e, flows) as well as elastic properties (returns to original shape after deformations).  224  Glossary  voxel - A 3D volume partition. wind crust - A layer of wind-deposited snow, generally containing small, densely packed grains. w i n d w a r d - The side of an obstacle unprotected from the wind. wild snow - A n extremely light, non-cohesive form of new snow, formed during a snowstorm with cold temperatures and a complete absence of wind. world - A n implementation entity, representing all model objects, all snow surfaces, and volume bounded by the sky and ground planes.  225  Index  entities, 91 fr,  erosion, 48, 142  124  exposure map, 43 Abstract, ii  extinction, 33  A C M , xxi albedo, 31, 36  face, 92, 111, 116  Alias, 59, 61, 101  falling snow, 35, 77, 141  angle of repose, 75, 106, 111, 141, 143, 144, 146, 170  flake dusting, 134, 140, 196  effect of water content on, 204  flake flutter, 120  effect of water content on , 13  flake group, 96, 105, 131, 175 timing, 176  experimental data, 27, 204 explanation, 27  flakes, 96, 104  previous work, 46, 49  flour, 173  animations, 188  fluid dynamics, 50  A5I, xix  flutter area, 124  assumptions, 56  flux maximum, 152  avalanche, 28, 57, 97, 111, 208  form factor, 84  debris, 57, 59  freezing nuclei, 7, 9  loose snow, 29  friction, 42, 43, 46, 52, 76, 143  previous work, 52  frontal lifting, 6  slab, 30  funding Imager, xix  avalanche flake, 97  N5ERC, xix basal plane, 7  NSF, xix  blending function, 157, 159  future work, 202  bridging, 65, 81, 87, 156, 162  . density, 203 grain shape, 204  camera dependence, 73, 82, 106, 109, 135  grain size, 204  commercial software, 59  precipitation, 203  condensation nuclei, 5  temperature, 204  convective lifting, 7  water content, 204  cornice, 22, 65, 81, 87, 156, 163 cyclonic lifting, 6  gazebo, 183 grain shape, 12  Delaunay triangulation, 111, 114, 116  role in future work, 204  density, 11, 203  grain size, 12, 31  role in future work, 203  role in future work, 204  depressions, 125  Graphics Interface, xxi  depth hoar, 19, 25  grass, 188  dew point, 5  graupel, 9  Douglas Powell, xxi  ground, 91  drops, 77, 95, 111, 147 dry snow  hardness, 14  density, 11  haystack, 196  visual effect of grain size, 13  height field, 36, 47-50, 87  dust, 44  hybrid paradigm, 86  edge group, 95, 111, 112, 116, 156  ice density, 1 1  elevation model, 36  226  227  Index  ice crystal a-axis, 8  previous research involving snow, 36 role in snowfall, 5  c-axis, 8 formation of, 7, 16 growth, 7 axial growth, 8 collision, 9  natural phenomena challenges of, 1 definition of, 1 snow as an example of, 2  riming, 9  neighbours, 105, 111, 116  tern perature effects, 7  normals, 93. 206  turbulence. 8  N S E R C , xix  vapour motion, 7  N3F, xix  shape classification, 8 structure, 8 turbulence. 22  oriented surface particles, 51 orographic lifting, 5  vapour motion, 16, 18 Imager Lab inability to create real snow, 173 thanks, xix impending doom picture of, 15 implicit surfaces, 156, 159, 162, 169 source code, xxi importance ordering, 104, 131 impurities, 14 intersection bucketing, 128, 156, 175 rebucketing, 130 timing, 130, 175  parallel computation, 42 particle systems, 41, 43, 51, 80, 82, 87, 147 advantages, 45 disadvantages, 45 path tracing, 71, 129 perforated crust, 23 phases, 98 uncompleted, 100 pipeline of snow creation, 98 polar bears, 54 polygonisation, 156, 162 porosity, 11  lag, 149, 150  portability, 60  launch area, 93  precipitation  launch site, 93, 95, 100, 105, 155  methods of, 5  lava, 50  role in future work, 203  lighting, 165  scale of, 57 scales of, 7  marching cubes, 81, 162  previous research, 35  source code, xxi mass inflation, 157  QSlim, 162  melt-freeze, 20  Qslim, 163  melting, 51, 57 previous research, 36 mesh reduction, 164 meshing, 111 metaballs, 38, 80 model, 92 acquiring, 101 dynamics, 208 errors in the underlying, 206 exporting, 101 generated by commercial software, 60 poorly constructed, 199 static, 60 triangulation, 61, 207 mountains, 56  rain, 43, 170 ram penetrometer, 15 random walk, 124 range tree, 129, 150 rebucketing, 149, 150 reflectance, 31, 38 rendering, 60, 164 avalanches, 78 extinction, 33 goals, 67 impurities, 14 lighting, 165 oil painting effects, 200 reflectance, 31, 38  228  Index  timing, 208 restrictions, 56 riming, 9, 24  assumptions, 59 environmental influences. 10 water content in, 13  rounding, 16, 18, 19  mass  sampling, 93  meshing, 105, 106, 109, 112, 116, 117, 120  reallocation, 132 sastrugi, 209 scale, 56  holes, 117 metamorphism, 16, 57  multi-scale attributes, 53, 105 of D E M data, 38 of precipitation. 7, 203  constructive, 18 previous work, 36 on the ground, 9, 45  surface representation. 72  partial coverage, 62  temporal adaptivity, 74  partial occlusion, 70, 86, 104, 106, 125  scattering, 38,41 effects of density, 11  physical properties, 11 density, 11, 203  scour-holes, 22  extinction, 33  settling, 87, 203  grain shape, 12, 204  shock-loading, 75  grain size, 12, 204  sintering, 18, 19  hardness, 14  sky, 91, 96, 104  mass, 65  snow from the, 56 sky bucketing, 131, 132 timing, 176 used for sky writing, 133 used to simulate precipitation, 133, 203 slush, 13, 20, 28, 143 snow  reflectance, 31, 36 strength, 14 temperature, 15, 204 water content, 13, 204 porosity, 11 problem advantages of the solution, 2  accumulation, 70, 103  artistic uses, 2  computing, 104  assumptions, 56, 57, 59-61  goals, 62, 174  complexity of, 2, 5, 55  air-borne, 5  goals, 61, 62, 67, 68, 172  avalanches, 59, 60, 205  justification for study, 2  coverage of land mass, 2  novelty of, 2  due to mountains, 5, 57, 203  overview, 70, 89, 90  due to precipitation, 5  previous research on, 2, 35  flake motion, 123 flake dusting, 207 flake flutter, 103, 105, 111, 123  restrictions, 56 simplification of, 5, 55, 89 rendering, 164, 165 effect of grain shape, 12  effects of, 125  effect of grain size, 12, 165  overview, 123  effect of water content, 13, 36  parameters, 128  future work, 204  role in future work, 204  previous work, 36, 38, 41  swirl length, 127  specific gravity, 11  wind, 130  stability, 67, 74, 85, 141, 142  importance ordering, 104, 106  adaptive subdivision, 76  impurities, 14  drops, 77, 141, 147  in Western Canada, 5  due to wind, 68, 151  launch area, 105, 111, 112, 116, 117, 120  goals, 67  launch site, 105, 111, 112, 116, 117, 120  jamming, 145, 205  launch sites  level-of-detail, 75  initial location, 120 layers, 9, 19, 30, 57, 81, 86, 165  list of samples, 142  Index  229  multipass, 143  role in future work, 204  neighbours, 76, 143, 144, 148  temperature gradient, 15  obstacle testing, 144, 152  test stab, 147  obstacles, 147  texture interpenetration, 138  speedups, 145  texture replacement, 135  termination, 76, 149, 206  textures, 134, 165, 207  water, 170  timing, 100  strength, 26, 46  tracks, 48, 68, 73  compressive stress, 26  translucence, 165  shear stress, 27  triangulation, 61, 92, 101, 111, 116  surface  turbulence, 8, 42  multi-scale, 72 smoothness, 205 tracks in, 48, 68, 169 validation, 173  undersampling, 72 user interface source code, xxi  writing with, 133 snow plane, 96  validation, 173  snowflake  vertex reduction, 163  dendritic, 9  viscosity, 50  needles, 8, 204  volumetric models, 38, 79, 80  plates, 8, 204  Voronoi diagram, 111, 112, 116, 175  shape classification, 8 soil, 45, 48 strength, 46 source code implicit surfaces, xxi L E D A , xxi marching cubes, xxi user interface, xxi specular reflectance, 40, 41 stability phase, 141 staining, 43 storm, 98, 99 assumptions, 57 subdivision goals, 73 subdivision approximation, 98 subdivision area, 93 subdivisions, 73 sun crust, 23 sun cups, 24 sun position, 36 surface inter penetration, 144, 145, 196 paradigm, 83 reduction, 163 resampling, 155 smoothness, 67, 74, 143, 156, 205 surface hoar, 12, 25 temperature, 15, 204 as it affects crystal growth, 7, 16, 18 previous work, 36  water, 43, 50, 145, 170 water content, 13 role in future work, 204 wavelength, 165 weathering, 43, 49, 207 wet snow, 20, 28, 30, 31 density, 11 visual effect of grain size, 13 wild snow, 11 wind, 21, 68, 130, 147, 151 eddies, 21, 152 flux maximum, 152 saltation, 21 winter, 36 January snow coverage, 2 world, 91 writing, 133, 196  The Homage  When all is said and done, it is always best to let Nature have the last word.  230  2S1  Index  Figure B . l : The road ahead.  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items