Table Of ContentStanley Osher Ronald Fedkiw
Level Set Methods and
Dynamic Implicit Surfaces
With 99 Figures, Including 24 in Full Color
viii Preface
ofmovinginterfacesplaysakeyroleintheproblemtobesolved.Asearchof
(cid:141)levelsetmethods(cid:142)ontheGooglewebsite(whichgaveover2,700responses
as of May 2002) will give an interested reader some idea of the scope and
utilityofthemethod.Inaddition,someexcitingadvancesinthetechnology
havebeenmadesincewebeganwritingthisbook.Wehopetocovermanyof
thesetopicsinafutureedition.Inthemeantimeyoucan(cid:147)ndsomeexciting
Preface animationsandmovingimagesaswellaslinkstomorerelevantresearchpa-
persviaourpersonalwebsites:http://graphics.stanford.edu/~fedkiw
and http://www.math.ucla.edu/~sjo/.
Acknowledgments
Many people have helped us in this e(cid:1)ort. We thank the following col-
leagues in particular: Steve Marschner, Paul Romburg, Gary Hewer, and
Steve Ruuth for proofreading parts of the manuscript, Peter Smereka and
Li-Tien Cheng for providing (cid:147)gures for the chapter on Codimension-Two
Objects, Myungjoo Kang for providing (cid:147)gures for the chapter on Motion
Involving Mean Curvature and Motion in the Normal Direction, Antonio
MarquinaandFredericGibouforhelpwiththechapteronImageRestora-
tion, Hong-Kai Zhao for help with chapter 13, Reconstruction of Surfaces
Scope, Aims, and Audiences from Unorganized Data Points, and Luminita Vese for help with the chap-
ter on Snakes, Active Contours, and Segmentation. We particularly thank
This book, Level Set Methods and Dynamic Implicit Surfaces is designed Barry Merriman for his extremely valuable collaboration on much of the
to serve two purposes: research described here. Of course we have bene(cid:147)tted immensely from col-
Parts I and II introduce the reader to implicit surfaces and level set laborationsanddiscussionswithfartoomanypeopletomention.Wehope
methods.Wehaveusedthesechapterstoteachintroductorycoursesonthe these colleagues and friends forgive us for omitting their names.
materialtostudentswithlittlemorethanafundamentalmathbackground. We would like to thank the following agencies for their support during
No prior knowledge of partial di(cid:1)erential equations or numerical analysis this period: ONR, AFOSR, NSF, ARO, and DARPA. We are particularly
is required. These (cid:147)rst eight chapters include enough detailed information grateful to Dr. Wen Masters of ONR for suggesting and believing in this
to allow students to create working level set codes from scratch. project and for all of her encouragement during some of the more di(cid:2)cult
Parts III and IV of this book are based on a series of papers published times.
by us and our colleagues. For the sake of brevity, a few details have been Finally, we thank our families and friends for putting up with us during
occasionallyomitted.Thesechaptersdoincludethoroughexplanationsand this exciting, but stressful period.
enough of the signi(cid:147)cant details along with the appropriate references to
allow the reader to get a (cid:147)rm grasp on the material. Los Angeles, California Stanley Osher
This book is an introduction to the subject. We have given examples of Stanford, California Ronald Fedkiw
theutilityofthemethodtoadiverse(butbynomeanscomplete)collection
of application areas. We have also tried to give complete numerical recipes
and a self-contained course in the appropriate numerical analysis. We be-
lievethatthisbookwillenableuserstoapplythetechniquespresentedhere
to real problems.
Thelevelsetmethodhasbeenusedinarapidlygrowingnumberofareas,
far too many to be represented here. These include epitaxial growth, opti-
maldesign,CAD,MEMS,optimalcontrol,andotherswherethesimulation
x Contents
3.2 Upwind Di(cid:1)erencing . . . . . . . . . . . . . . . . . . . . 29
3.3 Hamilton-Jacobi ENO . . . . . . . . . . . . . . . . . . . 31
3.4 Hamilton-Jacobi WENO . . . . . . . . . . . . . . . . . . 33
3.5 TVD Runge-Kutta . . . . . . . . . . . . . . . . . . . . . 37
Contents 4 Motion Involving Mean Curvature 41
4.1 Equation of Motion . . . . . . . . . . . . . . . . . . . . . 41
4.2 Numerical Discretization . . . . . . . . . . . . . . . . . . 44
4.3 Convection-Di(cid:1)usion Equations . . . . . . . . . . . . . . 45
5 Hamilton-Jacobi Equations 47
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.2 Connection with Conservation Laws . . . . . . . . . . . . 48
5.3 Numerical Discretization . . . . . . . . . . . . . . . . . . 49
5.3.1 Lax-Friedrichs Schemes . . . . . . . . . . . . . . . 50
5.3.2 The Roe-Fix Scheme . . . . . . . . . . . . . . . . 52
5.3.3 Godunov(cid:144)s Scheme . . . . . . . . . . . . . . . . . 54
6 Motion in the Normal Direction 55
6.1 The Basic Equation . . . . . . . . . . . . . . . . . . . . . 55
6.2 Numerical Discretization . . . . . . . . . . . . . . . . . . 57
6.3 Adding a Curvature-Dependent Term . . . . . . . . . . . 59
Preface vii
6.4 Adding an External Velocity Field . . . . . . . . . . . . . 59
Color Insert (facing page 146)
7 Constructing Signed Distance Functions 63
I Implicit Surfaces 1
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 63
7.2 Reinitialization . . . . . . . . . . . . . . . . . . . . . . . 64
1 Implicit Functions 3
7.3 Crossing Times . . . . . . . . . . . . . . . . . . . . . . . 65
1.1 Points . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
7.4 The Reinitialization Equation . . . . . . . . . . . . . . . 65
1.2 Curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
7.5 The Fast Marching Method . . . . . . . . . . . . . . . . 69
1.3 Surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4 Geometry Toolbox . . . . . . . . . . . . . . . . . . . . . 8 8 Extrapolation in the Normal Direction 75
1.5 Calculus Toolbox . . . . . . . . . . . . . . . . . . . . . . 13 8.1 One-Way Extrapolation. . . . . . . . . . . . . . . . . . . 75
8.2 Two-Way Extrapolation . . . . . . . . . . . . . . . . . . 76
2 Signed Distance Functions 17
8.3 Fast Marching Method . . . . . . . . . . . . . . . . . . . 76
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2 Distance Functions . . . . . . . . . . . . . . . . . . . . . 17 9 Particle Level Set Method 79
2.3 Signed Distance Functions . . . . . . . . . . . . . . . . . 18 9.1 Eulerian Versus Lagrangian Representations . . . . . . . 79
2.4 Examples. . . . . . . . . . . . . . . . . . . . . . . . . . . 19 9.2 Using Particles to Preserve Characteristics . . . . . . . . 82
2.5 Geometry and Calculus Toolboxes . . . . . . . . . . . . . 21
10 Codimension-Two Objects 87
10.1 Intersecting Two Level Set Functions . . . . . . . . . . . 87
II Level Set Methods 23 10.2 Modeling Curves in 3 . . . . . . . . . . . . . . . . . . . 87
(cid:1)
10.3 Open Curves and Surfaces . . . . . . . . . . . . . . . . . 90
3 Motion in an Externally Generated Velocity Field 25 10.4 Geometric Optics in a Phase-Space-Based Level
3.1 Convection . . . . . . . . . . . . . . . . . . . . . . . . . . 25 Set Framework. . . . . . . . . . . . . . . . . . . . . . . . 90
Contents xi xii Contents
III Image Processing and Computer Vision 95 15 Two-Phase Compressible Flow 167
15.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 167
11 Image Restoration 97 15.2 Errors at Discontinuities . . . . . . . . . . . . . . . . . . 168
11.1 Introduction to PDE-Based Image Restoration . . . . . . 97 15.3 Rankine-Hugoniot Jump Conditions . . . . . . . . . . . . 169
11.2 Total Variation-Based Image Restoration . . . . . . . . . 99 15.4 Nonconservative Numerical Methods . . . . . . . . . . . 171
11.3 Numerical Implementation of TV Restoration . . . . . . 103 15.5 Capturing Conservation . . . . . . . . . . . . . . . . . . 172
15.6 A Degree of Freedom . . . . . . . . . . . . . . . . . . . . 172
12 Snakes, Active Contours, and Segmentation 119
15.7 Isobaric Fix . . . . . . . . . . . . . . . . . . . . . . . . . 173
12.1 Introduction and Classical Active Contours. . . . . . . . 119
15.8 Ghost Fluid Method . . . . . . . . . . . . . . . . . . . . 175
12.2 Active Contours Without Edges . . . . . . . . . . . . . . 121
15.9 A Robust Alternative Interpolation . . . . . . . . . . . . 183
12.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
12.4 Extensions . . . . . . . . . . . . . . . . . . . . . . . . . . 124
16 Shocks, Detonations, and De(cid:148)agrations 189
16.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 189
13 Reconstruction of Surfaces from Unorganized
16.2 Computing the Velocity of the Discontinuity . . . . . . . 190
Data Points 139
16.3 Limitations of the Level Set Representation . . . . . . . 191
13.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 139
16.4 Shock Waves . . . . . . . . . . . . . . . . . . . . . . . . . 191
13.2 The Basic Model . . . . . . . . . . . . . . . . . . . . . . 140
16.5 Detonation Waves . . . . . . . . . . . . . . . . . . . . . . 193
13.3 The Convection Model . . . . . . . . . . . . . . . . . . . 142
16.6 De(cid:148)agration Waves . . . . . . . . . . . . . . . . . . . . . 195
13.4 Numerical Implementation . . . . . . . . . . . . . . . . . 142
16.7 Multiple Spatial Dimensions . . . . . . . . . . . . . . . . 196
IV Computational Physics 147 17 Solid-Fluid Coupling 201
17.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 201
14 Hyperbolic Conservation Laws and 17.2 Lagrange Equations . . . . . . . . . . . . . . . . . . . . . 203
Compressible Flow 149 17.3 Treating the Interface . . . . . . . . . . . . . . . . . . . . 204
14.1 Hyperbolic Conservation Laws . . . . . . . . . . . . . . . 149
14.1.1 Bulk Convection and Waves . . . . . . . . . . . . 150 18 Incompressible Flow 209
14.1.2 Contact Discontinuities . . . . . . . . . . . . . . . 151 18.1 Equations . . . . . . . . . . . . . . . . . . . . . . . . . . 209
14.1.3 Shock Waves. . . . . . . . . . . . . . . . . . . . . 152 18.2 MAC Grid . . . . . . . . . . . . . . . . . . . . . . . . . . 210
14.1.4 Rarefaction Waves . . . . . . . . . . . . . . . . . 153 18.3 Projection Method . . . . . . . . . . . . . . . . . . . . . 212
14.2 Discrete Conservation Form . . . . . . . . . . . . . . . . 154 18.4 Poisson Equation . . . . . . . . . . . . . . . . . . . . . . 213
14.3 ENO for Conservation Laws . . . . . . . . . . . . . . . . 155 18.5 Simulating Smoke for Computer Graphics . . . . . . . . 214
14.3.1 Motivation . . . . . . . . . . . . . . . . . . . . . . 155
14.3.2 Constructing the Numerical Flux Function . . . . 157 19 Free Surfaces 217
14.3.3 ENO-Roe Discretization 19.1 Description of the Model . . . . . . . . . . . . . . . . . . 217
(Third-Order Accurate) . . . . . . . . . . . . . . 158 19.2 Simulating Water for Computer Graphics . . . . . . . . . 218
14.3.4 ENO-LLF Discretization
(and the Entropy Fix) . . . . . . . . . . . . . . . 159 20 Liquid-Gas Interactions 223
14.4 Multiple Spatial Dimensions . . . . . . . . . . . . . . . . 160 20.1 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . 223
14.5 Systems of Conservation Laws . . . . . . . . . . . . . . . 160 20.2 Treating the Interface . . . . . . . . . . . . . . . . . . . . 224
14.5.1 The Eigensystem . . . . . . . . . . . . . . . . . . 161
21 Two-Phase Incompressible Flow 227
14.5.2 Discretization . . . . . . . . . . . . . . . . . . . . 162
14.6 Compressible Flow Equations . . . . . . . . . . . . . . . 163 21.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 227
14.6.1 Ideal Gas Equation of State . . . . . . . . . . . . 164 21.2 Jump Conditions . . . . . . . . . . . . . . . . . . . . . . 230
14.6.2 Eigensystem . . . . . . . . . . . . . . . . . . . . . 164 21.3 Viscous Terms . . . . . . . . . . . . . . . . . . . . . . . . 232
14.6.3 Numerical Approach . . . . . . . . . . . . . . . . 165 21.4 Poisson Equation . . . . . . . . . . . . . . . . . . . . . . 235
Contents xiii
22 Low-Speed Flames 239
22.1 Reacting Interfaces . . . . . . . . . . . . . . . . . . . . . 239
22.2 Governing Equations . . . . . . . . . . . . . . . . . . . . 240
Part I
22.3 Treating the Jump Conditions . . . . . . . . . . . . . . . 241
23 Heat Flow 249 Implicit Surfaces
23.1 Heat Equation . . . . . . . . . . . . . . . . . . . . . . . . 249
23.2 Irregular Domains . . . . . . . . . . . . . . . . . . . . . . 250
23.3 Poisson Equation . . . . . . . . . . . . . . . . . . . . . . 251
23.4 Stefan Problems . . . . . . . . . . . . . . . . . . . . . . . 254
References 259
In the next two chapters we introduce implicit surfaces and illustrate a
Index 271
number of useful properties, focusing on those that will be of use to us
later in the text. A good general review can be found in [16]. In the (cid:147)rst
chapter we discuss those properties that are true for a general implicit
representation. In the second chapter we introduce the notion of a signed
distance function with a Euclidean distance metric and a (cid:141) (cid:142) sign used to
–
indicate the inside and outside of the surface.
4 1. Implicit Functions
(cid:1)= x2 (cid:138)1
1
Implicit Functions
(cid:2)+ (cid:2)(cid:138) (cid:2)+
(cid:1)>0 (cid:1)<0 (cid:1)>0
x
outside inside outside
(cid:3)(cid:2) (cid:3)(cid:2)
(cid:1)=0 (cid:1)=0
interface interface
Figure 1.1. Implicit function (cid:1)(x) = x2 1 de(cid:147)ning the regions (cid:1)(cid:138) and (cid:1)+ as
(cid:138)
well asthe boundary (cid:2)(cid:1)
1.1 Points wasteful, since the implicit function (cid:2)((cid:3)x) is de(cid:147)ned on all of n, while the
(cid:1)
interface has only dimension n 1. However, we will see that a number of
(cid:138)
Inonespatialdimension,supposewedividethereallineintothreedistinct very powerful tools are readily available when we use this representation.
pieces using the points x = 1 and x = 1. That is, we de(cid:147)ne ( , 1), Above, we chose the (cid:2)(x) = 0 isocontour to represent the lower-
( 1,1), and (1, ) as three s(cid:138)eparate subdomains of interest, alth(cid:138)o(cid:2)ugh(cid:138)we dimensional interface, but there is nothing special about the zero
r(cid:138)egardthe(cid:147)rsta(cid:2)ndthirdastwodisjointpiecesofthesameregion.Werefer isocontour. For example, the (cid:2)(cid:26)(x) = 1 isocontour of (cid:2)(cid:26)(x) = x2, de(cid:147)nes
to(cid:3) =( 1,1)astheinsideportionofthedomainand(cid:3)+ =( , 1) the same interface, (cid:1)(cid:3) = 1,1 . In general, for any function (cid:2)(cid:26)((cid:3)x) and
(1, (cid:138))ast(cid:138)heoutsideportionofthedomain.Theborderbetween(cid:138)t(cid:2)hei(cid:138)nsid(cid:3)e an arbitrary isocontour (cid:2)(cid:26)({(cid:3)x(cid:138)) = }a for some scalar a , we can de(cid:147)ne
(cid:2) (cid:4) (cid:1)
and the outside consists of the two points (cid:1)(cid:3) = 1,1 and is called (cid:2)((cid:3)x) = (cid:2)(cid:26)((cid:3)x) a, so that the (cid:2)((cid:3)x) = 0 isocontour of (cid:2) is identical to the
{(cid:138) } (cid:138)
the interface. In one spatial dimension, the inside and outside regions are (cid:2)(cid:26)((cid:3)x)=a isocontour of (cid:2)(cid:26). In addition, the functions (cid:2) and (cid:2)(cid:26)have identical
one-dimensional objects, while the interface is less than one-dimensional. properties up to a scalar translation a. Moreover, the partial derivatives
In fact, the points making up the interface are zero-dimensional. More of (cid:2) are the same as the partial derivatives of (cid:2)(cid:26), since the scalar vanishes
generally, in n, subdomains are n-dimensional, while the interface has upondi(cid:1)erentiation.Thus,throughoutthetextallofourimplicitfunctions
(cid:1)
dimension n 1. We say that the interface has codimension one. (cid:2)((cid:3)x)willbede(cid:147)nedsothatthe(cid:2)((cid:3)x)=0isocontourrepresentstheinterface
(cid:138)
In an explicit interface representation one explicitly writes down the (unless otherwise speci(cid:147)ed).
points that belong to the interface as we did above when de(cid:147)ning (cid:1)(cid:3) =
1,1 .Alternatively,animplicitinterfacerepresentationde(cid:147)nestheinter-
{(cid:138) }
face as the isocontour of some function. For example, the zero isocontour
1.2 Curves
of (cid:2)(x) = x2 1 is the set of all points where (cid:2)(x) = 0; i.e., it is ex-
(cid:138)
actly (cid:1)(cid:3) = 1,1 . This is shown in Figure 1.1. Note that the implicit
{(cid:138) } In two spatial dimensions, our lower-dimensional interface is a curve that
function(cid:2)(x)isde(cid:147)nedthroughouttheone-dimensionaldomain,whilethe
separates 2 into separate subdomains with nonzero areas. Here we are
isocontour de(cid:147)ning the interface is one dimension lower. More generally, (cid:1)
limiting our interface curves to those that are closed, so that they have
in n, the implicit function (cid:2)((cid:3)x) is de(cid:147)ned on all (cid:3)x n, and its isocon-
(cid:1) (cid:4) (cid:1) clearlyde(cid:147)nedinteriorandexteriorregions.Asanexample,consider(cid:2)((cid:3)x)=
tour has dimension n 1. Initially, the implicit representation might seem
(cid:138)
1.2. Curves 5 6 1. Implicit Functions
y interval [s ,s ], one needs to resolve a two-dimensional region D. More
o f
generally, in n, a discretization of an explicit representation needs to re-
(cid:1)
solve only an (n 1)-dimensional set, while a discretization of an implicit
(cid:138)
representation needs to resolve an n-dimensional set. This can be avoided,
in part, by placing all the points (cid:3)x very close to the interface, leaving the
(cid:2)(cid:138) (cid:2)+ restofD unresolved.Sinceonlythe(cid:2)((cid:3)x)=0isocontourisimportant,only
the points (cid:3)x near this isocontour are actually needed to accurately repre-
(cid:1)<0 (cid:1)>0
senttheinterface.TherestofDisunimportant.Clusteringpointsnearthe
inside outside
x interface is a local approach to discretizing implicit representations. (We
will give more details about local approaches later.) Once we have chosen
thesetofpointsthatmakeupourdiscretization,westorethevaluesofthe
implicit function (cid:2)((cid:3)x) at each of these points.
Neither the explicit nor the implicit discretization tells us where the in-
(cid:3)(cid:2)
terface is located. Instead, they both give information at sample locations.
(cid:1)= x2 + y2 (cid:138)1=0
Intheexplicitrepresentation,weknowthelocationofa(cid:147)nitesetofpoints
interface on the curve, but do not know the location of the remaining in(cid:147)nite set
of points (on the curve). Usually, interpolation is used to approximate the
Figure 1.2. Implicit representation ofthe curve x2+y2=1. locationofpointsnotrepresentedinthediscretization.Forexample,piece-
wise polynomial interpolation can be used to determine the shape of the
interface between the data points. Splines are usually appropriate for this.
x2+y2 1, where the interface de(cid:147)ned by the (cid:2)((cid:3)x)=0 isocontour is the Similarly, in the implicit representation we know the values of the implicit
(cid:138)
unit circle de(cid:147)ned by (cid:1)(cid:3) = (cid:3)x (cid:3)x = 1 . The interior region is the unit function (cid:2) at only a (cid:147)nite number of points and need to use interpolation
{ | | | }
opendisk(cid:3) = (cid:3)x (cid:3)x <1 ,andtheexteriorregionis(cid:3)+ = (cid:3)x (cid:3)x >1 . to (cid:147)nd the values of (cid:2) elsewhere. Even worse, here we may not know the
(cid:138)
{ || | } { || | }
TheseregionsaredepictedinFigure1.2.Theexplicitrepresentationofthis location of any of the points on the interface, unless we have luckily cho-
interface is simply the unit circle de(cid:147)ned by (cid:1)(cid:3)= (cid:3)x (cid:3)x =1 . sen data points (cid:3)x where (cid:2)((cid:3)x) is exactly equal to zero. In order to locate
{ || | }
Intwospatial dimensions, theexplicit interfacede(cid:147)nition needstospec- the interface, the (cid:2)((cid:3)x) = 0 isocontour needs to be interpolated from the
ify all the points on a curve. While in this case it is easy to do, it can be known values of (cid:2) at the data points. This is a rather standard procedure
somewhatmoredi(cid:2)cultforgeneralcurves.Ingeneral,oneneedstoparam- accomplished by a variety of contour plotting routines.
eterize the curve with a vector function (cid:3)x(s), where the parameter s is in The set of data points where the implicit function (cid:2) is de(cid:147)ned is called
[s ,s ]. The condition that the curve be closed implies that (cid:3)x(s )=(cid:3)x(s ). a grid. There are many ways of choosing the points in a grid, and these
o f o f
While it is convenient to use analytical descriptions as we have done lead to a number of di(cid:1)erent types of grids, e.g., unstructured, adaptive,
so far, complicated two-dimensional curves do not generally have such curvilinear. By far, the most popular grids, are Cartesian grids de(cid:147)ned as
simple representations. A convenient way of approximating an explicit (x ,y ) 1 i m,1 j n . The natural orderings of the x and y
i j i j
{ | (cid:6) (cid:6) (cid:6) (cid:6) }
representation is to discretize the parameter s into a (cid:147)nite set of points are usually used for convenience. That is, x < < x < x < x <
1 i 1 i i+1
••• (cid:138)
s < < s < s < s < < s , where the subintervals [s ,s ] < x and y < < y < y < y < < y . In a uniform
o i 1 i i+1 f i i+1 m 1 j 1 j j+1 n
••• (cid:138) ••• ••• ••• (cid:138) •••
are not necessarily of equal size. For each point s in parameter space, Cartesian grid, all the subintervals [x ,x ] are equal in size, and we set
i i i+1
we then store the corresponding two-dimensional location of the curve de- x = x x . Likewise, all the subintervals [y ,y ] are equal in size,
i+1 i j j+1
(cid:7) (cid:138)
notedby(cid:3)x(s ).Asthenumberofpointsinthediscretizedparameterspace andweset y =y y .Furthermore,itisusuallyconvenienttochoose
i j+1 j
(cid:7) (cid:138)
is increased, so is the resolution (detail) of the two-dimensional curve. x= y sothattheapproximation errorsarethesameinthe x-direction
(cid:7) (cid:7)
The implicit representation can be stored with a discretization as well, as they are in the y-direction. By de(cid:147)nition, Cartesian grids imply the use
exceptnowoneneedstodiscretizeallof 2,whichisimpractical,sinceitis ofarectangulardomainD=[x ,x ] [y ,y ].Again,since(cid:2)isimportant
1 m 1 n
(cid:1) (cid:215)
unbounded. Instead, we discretize a bounded subdomain D 2. Within only near the interface, a local approach would indicate that many of the
(cid:5) (cid:1)
thisdomain,wechoosea(cid:147)nitesetofpoints(x ,y )fori=1,...,N todis- gridpointsarenotneeded,andtheimplicitrepresentationcanbeoptimized
i i
cretely approximate the implicit function (cid:2). This illustrates a drawback of by storing only a subset of a uniform Cartesian grid. The Cartesian grid
the implicit surface representation. Instead of resolving a one-dimensional points that are not su(cid:2)ciently near the interface can be discarded.
1.3. Surfaces 7 8 1. Implicit Functions
We pause for a moment to consider the discretization of the one- surgery(cid:142) needed for merging and pinching is much more complex, leading
dimensional problem. There, since the explicit representation is merely a to a number of di(cid:2)culties including, for example, holes in the surface.
set of points, it is trivial to record the exact interface position, and no Oneofthenicestpropertiesofimplicitsurfacesisthatconnectivitydoes
discretization or parameterization is needed. However, the implicit repre- notneedtobedeterminedfor thediscretization. AuniformCartesian grid
sentationmustbediscretizedif(cid:2)isnotaknownanalyticfunction.Atypical (x ,y ,z ) 1 i m,1 j n,1 k p can be used along
i j k
{ | (cid:6) (cid:6) (cid:6) (cid:6) (cid:6) (cid:6) }
discretization consists of a set of points x < < x < x < x < with straightforward generalizations of the technology from two spatial
1 i 1 i i+1
••• (cid:138)
<x on a subdomain D =[x ,x ] of . Again, it is usually useful to dimensions. Possibly the most powerful aspect of implicit surfaces is that
m 1 m
••• (cid:1)
use a uniform grid, and only the grid points near the interface need to be it is straightforward to go from two spatial dimensions to three spatial
stored. dimensions (or even more).
1.3 Surfaces 1.4 Geometry Toolbox
In three spatial dimensions the lower-dimensional interface is a surface Implicit interface representations include some very powerful geometric
that separates 3 into separate subdomains with nonzero volumes. Again, tools. For example, since we have designated the (cid:2)((cid:3)x) = 0 isocontour as
(cid:1)
we consider only closed surfaces with clearly de(cid:147)ned interior and exterior the interface, we can determine which side of the interface a point is on
regions.Asanexample,consider(cid:2)((cid:3)x)=x2+y2+z2 1,wheretheinterface simply by looking at the local sign of (cid:2). That is, (cid:3)x is inside the interface
(cid:138) o
is de(cid:147)ned by the (cid:2)((cid:3)x) = 0 isocontour, which is the boundary of the unit when(cid:2)((cid:3)x )<0,outsidetheinterfacewhen(cid:2)((cid:3)x )>0,andontheinterface
o o
sphere de(cid:147)ned as (cid:1)(cid:3) = (cid:3)x (cid:3)x = 1 . The interior region is the open unit when (cid:2)((cid:3)x )=0. With an explicit representation of the interface it can be
{ | | | } o
sphere (cid:3) = (cid:3)x x(cid:3) < 1 , and the exterior region is (cid:3)+ = (cid:3)x (cid:3)x > 1 . di(cid:2)cult to determine whether a point is inside or outside the interface. A
(cid:138)
{ | | | } { | | | }
The explicit representation of the interface is (cid:1)(cid:3)= (cid:3)x (cid:3)x =1 . standardprocedurefordoingthisistocastarayfromthepointinquestion
{ || | }
For complicated surfaces with no analytic representation, we again need to some far-o(cid:1) place that is known to be outside the interface. Then if the
to use a discretization. In three spatial dimensions the explicit represen- ray intersects the interface an even number of times, the point is outside
tation can be quite di(cid:2)cult to discretize. One needs to choose a number the interface. Otherwise, the ray intersects the interface an odd number of
of points on the two-dimensional surface and record their connectivity. In times,andthepointisinsidetheinterface.Obviously,itismoreconvenient
twospatialdimensions,connectivitywasdeterminedbasedontheordering, simply to evaluate (cid:2) at the point (cid:3)x . In the discrete case, i.e., when the
o
i.e., (cid:3)x(s ) is connected to (cid:3)x(s ) and (cid:3)x(s ). In three spatial dimensions implicit function is given by its values at a (cid:147)nite number of data points,
i i 1 i+1
(cid:138)
connectivityislessstraightforward.Iftheexactsurfaceanditsconnectivity interpolation can be used to estimate (cid:2)((cid:3)x ) using the values of (cid:2) at the
o
are known, it is simple to tile the surface with triangles whose vertices lie known sample points. For example, on our Cartesian grid, linear, bilin-
ontheinterfaceandwhoseedgesindicateconnectivity.Ontheotherhand, ear, and trilinear interpolation can be used in one, two, and three spatial
ifconnectivityisnotknown,itcanbequitedi(cid:2)culttodetermine,andeven dimensions, respectively.
some of the most popular algorithms can produce surprisingly inaccurate Numerical interpolation produces errors in the estimate of (cid:2). This can
surface representations, e.g., surfaces with holes. lead to erroneously designating inside points as outside points and vice
Connectivitycanchangefordynamicimplicitsurfaces,i.e.,surfacesthat versa. At (cid:147)rst glance these errors might seem disastrous, but in reality
are moving around. As an example, consider the splashing water surface they amount to perturbing (or moving) the interface away from its exact
in a swimming pool full of children. Here, connectivity is not a (cid:141)one-time(cid:142) position.Iftheseinterfaceperturbationsaresmall,theire(cid:1)ectsmaybemi-
issue dealt with in constructing an explicit representation of the surface. nor,andaperturbedinterfacemightbeacceptable.Infact,mostnumerical
Instead, it must be resolved over and over again every time pieces of the methods depend on the fact that the results are stable in the presence of
surface merge together or pinch apart. In two spatial dimensions the task small perturbations. If this is not true, then the problem under consider-
is more manageable, since merging can be accomplished by taking two ation is probably ill-posed, and numerical methods should be used only
one-dimensional parameterizations, s and s(cid:26), and combining them into a with extreme caution (and suspicion). These interface perturbation errors
i i
singleone-dimensionalparameterization.Pinchingapartisaccomplishedby decreaseasthenumberofsamplepointsincreases,implyingthattheexact
splitting a single one-dimensional parameterization into two separate one- answer could hypothetically be computed as the number of sample points
dimensional parameterizations. In three spatial dimensions the (cid:141)interface isincreasedtoin(cid:147)nity.Again,thisisthebasisformostnumericalmethods.
1.4. Geometry Toolbox 9 10 1. Implicit Functions
While one cannot increase the number of grid points to in(cid:147)nity, desirable y
solutions can be obtained for many problems with a practical number of
grid points. Throughout the text we will make a number of numerical ap-
proximations with errors proportional to the size of a Cartesian mesh cell, N
i.e., x (or ( x)r). If the implicit function is smooth enough and well N
resolv(cid:7)ed by th(cid:7)e grid, these estimates will be appropriate. Otherwise, these N N N
errors might be rather large. Obviously, this means that we would like our
implicit function to be as smooth as possible. In the next chapter we dis-
cuss using a signed distance function to represent the surface. This turns N
x
out to be a good choice, since steep and (cid:148)at gradients as well as rapidly
changing features are avoided as much as possible. (cid:1)=(cid:138).8
Implicit functions make both simple Boolean operations and more ad- (cid:1)=0
vanced constructive solid geometry (CSG) operations easy to apply. This
is important, for example, in computer-aided design (CAD). If (cid:2) and (cid:2)
1 2 (cid:1)=1.15
are two di(cid:1)erent implicit functions, then (cid:2)((cid:3)x) = min((cid:2) ((cid:3)x),(cid:2) ((cid:3)x)) is
1 2
the implicit function representing the union of the interior regions of (cid:2)
1
and (cid:2) . Similarly, (cid:2)((cid:3)x) = max((cid:2) ((cid:3)x),(cid:2) ((cid:3)x)) is the implicit function
2 1 2
representing the intersection of the interior regions of (cid:2) and (cid:2) . The
1 2 Figure1.3.Afewisocontoursofourtwo-dimensionalexample(cid:1)((cid:3)x)=x2+y2 1
complement of (cid:2)1((cid:3)x) can be de(cid:147)ned by (cid:2)((cid:3)x) = (cid:138)(cid:2)1((cid:3)x). Also, (cid:2)((cid:3)x) = along with some representative normals. (cid:138)
max((cid:2) ((cid:3)x), (cid:2) ((cid:3)x)) represents the region obtained by subtracting the
1 2
(cid:138)
interior of (cid:2) from the interior of (cid:2) .
2 1
The gradient of the implicit function is de(cid:147)ned as including x=1, where the interface normal is N(cid:3) =1, and N(cid:3) points to the
(cid:1)(cid:2) (cid:1)(cid:2) (cid:1)(cid:2) left for all x<0 including x= 1, where the interface normal is N(cid:3) = 1.
(cid:2)= , , . (1.1) (cid:138) (cid:138)
(cid:8) (cid:1)x (cid:1)y (cid:1)z The normal is unde(cid:147)ned at x=0, since the denominator of equation (1.2)
(cid:1) (cid:2)
vanishes. This can be problematic in general, but can be avoided with a
Thegradient (cid:2)isperpendiculartotheisocontoursof(cid:2)andpointsinthe
(cid:8) number of techniques. For example, at x=0 we could simply de(cid:147)ne N(cid:3) as
direction of increasing (cid:2). Therefore, if (cid:3)x is a point on the zero isocontour
o eitherN(cid:3) =1orN(cid:3) = 1.Ourtwo-andthree-dimensionalexamples(above)
opfoi(cid:2)n,tsi.ei.n,athpeosinatmoendtihreecitniotenrfaascet,htehleonca(cid:8)l(cid:2)uneivtal(uoauttewdaardt)(cid:3)xnooirsmaavleNc(cid:3)totroththaet show similar degener(cid:138)ate behavior at (cid:3)x = (cid:3)0, where all partial derivatives
vanish. Again, a simple technique for evaluating (1.2) at these points is
interface. Thus, the unit (outward) normal is
just to pick an arbitrary direction for the normal. Note that the standard
N(cid:3) = (cid:8)(cid:2) (1.2) trick of adding a small (cid:4) > 0 to the denominator of equation (1.2) can
(cid:2) be a bad idea in general, since it produces a normal with N(cid:3) =1. In fact,
|(cid:8) | | |(cid:9)
for points on the interface. whenthedenominatorinequation(1.2)iszero,soisthenumerator,making
Since the implicit representation of the interface embeds the interface N(cid:3) =(cid:3)0whenasmall(cid:4)>0isusedinthedenominator.(WhilesettingN(cid:3) =(cid:3)0
in a domain of one higher-dimension, it will be useful to have as much is not always disastrous, caution is advised.)
information as possible representable on the higher-dimensional domain. On our Cartesian grid, the derivatives in equation (1.2) need to be ap-
For example, instead of de(cid:147)ning the unit normal N(cid:3) by equation (1.2) for proximated, for example using (cid:147)nite di(cid:1)erence techniques. We can use a
points on the interface only, we use equation (1.2) to de(cid:147)ne a function N(cid:3) (cid:147)rst-order accurate forward di(cid:1)erence
everywhereonthedomain.ThisembedsthenormalinafunctionN(cid:3) de(cid:147)ned (cid:1)(cid:2) (cid:2) (cid:2)
i+1(cid:138) i (1.3)
ontheentiredomainthatagreeswiththenormalforpointsontheinterface. (cid:1)x (cid:10) x
(cid:7)
Figure1.3showsafewisocontoursofourtwo-dimensionalexample(cid:2)((cid:3)x)=
x2+y2 1 along with some representative normals. abbreviated as D+(cid:2), a (cid:147)rst-order accurate backward di(cid:1)erence
Consi(cid:138)dertheone-dimensionalexample(cid:2)(x)=x2 1,whereN(cid:3) isde(cid:147)ned (cid:1)(cid:2) (cid:2) (cid:2)
by equation (1.2) as N(cid:3) = x/x. Here, N(cid:3) points to(cid:138)the right for all x > 0 (cid:1)x (cid:10) i(cid:138)xi(cid:138)1 (1.4)
| | (cid:7)
