<< . .

( 23)

. . >>

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
out to be a good choice, since steep and ¬‚at gradients as well as rapidly
φ = ’.8
changing features are avoided as much as possible.
φ =0
Implicit functions make both simple Boolean operations and more ad-
vanced constructive solid geometry (CSG) operations easy to apply. This
φ = 1.15
is important, for example, in computer-aided design (CAD). If φ1 and φ2
are two di¬erent implicit functions, then φ(x) = min(φ1 (x), φ2 (x)) is
the implicit function representing the union of the interior regions of φ1
and φ2 . Similarly, φ(x) = max(φ1 (x), φ2 (x)) is the implicit function
representing the intersection of the interior regions of φ1 and φ2 . The
Figure 1.3. A few isocontours of our two-dimensional example φ(x) = x2 + y 2 ’ 1
complement of φ1 (x) can be de¬ned by φ(x) = ’φ1 (x). Also, φ(x) =
along with some representative normals.
max(φ1 (x), ’φ2 (x)) represents the region obtained by subtracting the
interior of φ2 from the interior of φ1 .
The gradient of the implicit function is de¬ned as including x = 1, where the interface normal is N = 1, and N points to the
left for all x < 0 including x = ’1, where the interface normal is N = ’1.
‚φ ‚φ ‚φ
φ= , , . (1.1)
The normal is unde¬ned at x = 0, since the denominator of equation (1.2)
‚x ‚y ‚z
vanishes. This can be problematic in general, but can be avoided with a
The gradient φ is perpendicular to the isocontours of φ and points in the
number of techniques. For example, at x = 0 we could simply de¬ne N as
direction of increasing φ. Therefore, if xo is a point on the zero isocontour
either N = 1 or N = ’1. Our two- and three-dimensional examples (above)
of φ, i.e., a point on the interface, then φ evaluated at xo is a vector that
show similar degenerate behavior at x = 0, where all partial derivatives
points in the same direction as the local unit (outward) normal N to the
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
φ trick of adding a small > 0 to the denominator of equation (1.2) can
N= (1.2)
| φ| be a bad idea in general, since it produces a normal with |N | = 1. In fact,
when the denominator in equation (1.2) is zero, so is the numerator, making
for points on the interface.
N = 0 when a small > 0 is used in the denominator. (While setting N = 0
Since the implicit representation of the interface embeds the interface
is not always disastrous, caution is advised.)
in a domain of one higher-dimension, it will be useful to have as much
On our Cartesian grid, the derivatives in equation (1.2) need to be ap-
information as possible representable on the higher-dimensional domain.
proximated, for example using ¬nite di¬erence techniques. We can use a
For example, instead of de¬ning the unit normal N by equation (1.2) for
¬rst-order accurate forward di¬erence
points on the interface only, we use equation (1.2) to de¬ne a function N
φi+1 ’ φi
everywhere on the domain. This embeds the normal in a function N de¬ned ‚φ
≈ (1.3)
on the entire domain that agrees with the normal for points on the interface. ‚x x
Figure 1.3 shows a few isocontours of our two-dimensional example φ(x) =
abbreviated as D+ φ, a ¬rst-order accurate backward di¬erence
x2 + y 2 ’ 1 along with some representative normals.
Consider the one-dimensional example φ(x) = x2 ’ 1, where N is de¬ned φi ’ φi’1
≈ (1.4)
by equation (1.2) as N = x/|x|. Here, N points to the right for all x > 0 ‚x x
1.4. Geometry Toolbox 11 12 1. Implicit Functions

abbreviated as D’ φ, or a second-order accurate central di¬erence good choice for φ turns out to be the signed distance function discussed in
the next chapter.
φi+1 ’ φi’1
≈ (1.5) The mean curvature of the interface is de¬ned as the divergence of the
‚x 2x
normal N = (n1 , n2 , n3 ),
abbreviated as Do φ. (The j and k indices have been suppressed in the
‚n1 ‚n2 ‚n3
above formulas.) The formulas for the derivatives in the y and z directions ·N =
κ= + + , (1.6)
‚x ‚y ‚z
are obtained through symmetry. These simple formulas are by no means
exhaustive, and we will discuss more ways of approximating derivatives so that κ > 0 for convex regions, κ < 0 for concave regions, and κ = 0 for a
later in the text. plane; see Figure 1.4. While one could simply use ¬nite di¬erences to com-
When all numerically calculated ¬nite di¬erences are identically zero, pute the derivatives of the components of the normal in equation (1.6), it is
the denominator of equation (1.2) vanishes. As in the analytic case, we usually more convenient, compact, and accurate to calculate the curvature
can simply randomly choose a normal. Here, however, randomly choosing directly from the values of φ. Substituting equation (1.2) into equation (1.6)
a normal is somewhat justi¬ed, since it is equivalent to randomly perturb- gives
ing the values of φ on our Cartesian mesh by values near round-o¬ error.
These small changes in the values of φ are dominated by the local approx- ·
κ= , (1.7)
| φ|
imation errors in the ¬nite di¬erence formula for the derivatives. Consider
a discretized version of our one-dimensional example φ(x) = x2 ’ 1, and
so that we can write the curvature as
suppose that grid points exist at xi’1 = ’ x, xi = 0, and xi+1 = x with
κ = φ2 φyy ’ 2φx φy φxy + φ2 φxx + φ2 φzz ’ 2φx φz φxz + φ2 φxx
exact values of φ de¬ned as φi’1 = x2 ’ 1, φi = ’1, and φi+1 = x2 ’ 1, x y x z
respectively. The forward di¬erence formula gives Ni = 1, the backward +φ2 φzz ’ 2φy φz φyz + φ2 φyy /| φ|3 (1.8)
y z
di¬erence formula gives Ni = ’1, and the central di¬erence formula can-
not be used, since Do φ = 0 at xi = 0. However, simply perturbing φi+1 to in terms of the ¬rst and second derivatives of φ. A second-order accurate
x2 ’ 1 + for any small > 0 (even round-o¬ error) gives D o φ = 0 and ¬nite di¬erence formula for φxx , the second partial derivative of φ in the x
direction, is given by
Ni = 1. Similarly, perturbing φi’1 to x2 ’ 1 + gives Ni = ’1. Thus,
for any approach that is stable under small perturbations of the data, it is ‚2φ φi+1 ’ 2φi + φi’1
≈ , (1.9)
acceptable to randomly choose N when the denominator of equation (1.2) 2 x2
vanishes. Similarly in our two- and three-dimensional examples, N = x/|x|
abbreviated as Dx Dx φ, or equivalently, Dx Dx φ. Here D+ and D’ are
+’ ’+
everywhere except at x = 0, where equation (1.2) is not de¬ned and we
de¬ned as in equations (1.3) and (1.4), respectively, and the x subscript
are free to choose it arbitrarily. The arbitrary normal at the origin in the
indicates that the ¬nite di¬erence is evaluated in the x direction. A second-
one-dimensional case lines up with the normals to either the right, where
order accurate ¬nite di¬erence formula for φxy is given by Dx Dy φ, or
N = 1, or to the left, where N = ’1. Similarly, in two and three spatial di-
equivalently, Dy Dx φ. The other second derivatives in equation (1.8) are
mensions, an arbitrarily chosen normal at x = 0 lines up with other nearby
de¬ned in a manner similar to either φxx or φxy .
normals. This is always the case, since the normals near the origin point
In our one-dimensional example, φ(x) = x2 ’ 1, κ = 0 everywhere ex-
outward in every possible direction.
cept at the origin, where equation (1.7) is unde¬ned. Thus, the origin, is
If φ is a smooth well-behaved function, then an approximation to the
a removable singularity, and we can de¬ne κ = 0 everywhere. Interfaces in
value of the normal at the interface can be obtained from the values of N
one spatial dimension are models of planes in three dimensions (assuming
computed at the nodes of our Cartesian mesh. That is, given a point xo
that the unmodeled directions have uniform data). Therefore, using κ = 0
on the interface, one can estimate the unit outward normal at xo by in-
everywhere is a consistent model. In our two- and three-dimensional ex-
terpolating the values of N from the Cartesian mesh to the point xo . If
1 2
amples above, κ = |x| and κ = |x| (respectively) everywhere except at the
one is using forward, backward, or central di¬erences, then linear (bilinear
origin. Here the singularities are not removable, and κ ’ ∞ as we approach
or trilinear) interpolation is usually good enough. However, higher-order
accurate formulas can be used if desired. This interpolation procedure re- the origin. Moreover, κ = 1 everywhere on the one-dimensional interface
quires that φ be well behaved, implying that we should be careful in how in two spatial dimensions, and κ = 2 everywhere on the two-dimensional
we choose φ. For example, it would be unwise to choose an implicit func- interface in three spatial dimensions. The di¬erence occurs because a two-
tion φ with unnecessary oscillations or steep (or ¬‚at) gradients. Again, a dimensional circle is a cylinder in three spatial dimensions (assuming that
1.5. Calculus Toolbox 13 14 1. Implicit Functions

needed for the boundary. This is easily accomplished by including the
„¦ „¦’ measure-zero boundary set with either the interior or exterior region (as
φ >0 above). Throughout the text we usually include the boundary with the
φ <0 interior region „¦’ where φ(x) < 0 (unless otherwise speci¬ed).
outside The functions χ± are functions of a multidimensional variable x. It
is often more convenient to work with functions of the one-dimensional
variable φ. Thus we de¬ne the one-dimensional Heaviside function
κ <0 0 if φ ¤ 0,
‚„¦ H(φ) = (1.12)
1 if φ > 0,
φ =0
κ >0 where φ depends on x, although it is not important to specify this depen-
interface dence when working with H. This allows us to work with H in one spatial
dimension. Note that χ+ (x) = H(φ(x)) and χ’ (x) = 1 ’ H(φ(x)) for all x,
so all we have done is to introduce an extra function of one variable H to
Figure 1.4. Convex regions have κ > 0, and concave regions have κ < 0. be used as a tool when dealing with characteristic functions.
The volume integral (area or length integral in 2 or 1 , respectively) of
a function f over the interior region „¦’ is de¬ned as
the unmodeled direction has uniform data). It seems nonsensical to be trou-
bled by κ ’ ∞ as we approach the origin, since this is only a consequence
f (x)χ’ (x) dx, (1.13)
of the embedding. In fact, since the smallest unit of measure on the Carte- „¦
sian grid is the cell size x, it makes little sense to hope to resolve objects
where the region of integration is all of „¦, since χ’ prunes out the exterior
smaller than this. That is, it makes little sense to model circles (or spheres)
region „¦+ automatically. The one-dimensional Heaviside function can be
with a radius smaller than x. Therefore, we limit the curvature so that
used to rewrite this volume integral as
’ 1x ¤ κ ¤ 1x . If a value of κ is calculated outside this range, we merely
replace that value with either ’ 1x or 1x depending on which is closer.
f (x) (1 ’ H(φ(x))) dx (1.14)
As a ¬nal note on curvature, one has to use caution when φ is noisy. „¦
The normal N will generally have even more noise, since it is based on the
representing the integral of f over the interior region „¦’ . Similarly,
derivatives of φ. Similarly, the curvature κ will be even noisier than the
normal, since it is computed with the second derivatives of φ.
f (x)H(φ(x)) dx (1.15)

is the integral of f over the exterior region „¦+ .
1.5 Calculus Toolbox By de¬nition, the directional derivative of the Heaviside function H in
the normal direction N is the Dirac delta function
The characteristic function χ’ of the interior region „¦’ is de¬ned as
ˆ H(φ(x)) · N ,
δ(x) = (1.16)
if φ(x) ¤ 0,

χ (x) = (1.10)
which is a function of the multidimensional variable x. Note that this dis-
0 if φ(x) > 0
tribution is nonzero only on the interface ‚„¦ where φ = 0. We can rewrite
where we arbitrarily include the boundary with the interior region. The equation (1.16) as
characteristic function χ+ of the exterior region „¦+ is de¬ned similarly as
δ(x) = H (φ(x)) φ(x) · = H (φ(x))| φ(x)| (1.17)
if φ(x) ¤ 0,
0 | φ(x)|
χ+ (x) = (1.11)
1 if φ(x) > 0,
using the chain rule to take the gradient of H, the de¬nition of the normal
from equation (1.2), and the fact that φ(x) · φ(x) = | φ(x)|2 . In one
again including the boundary with the interior region. It is often useful to
have only interior and exterior regions so that special treatment is not spatial dimension, the delta function is de¬ned as the derivative of the
1.5. Calculus Toolbox 15 16 1. Implicit Functions

one-dimensional Heaviside function: interface width equal to three grid cells when φ is normalized to a signed
distance function with | φ| = 1; see Chapter 2). Then the delta function
δ(φ) = H (φ), (1.18)
is de¬ned according to equation (1.18) as the derivative of the Heaviside
where H(φ) is de¬ned in equation (1.12) above. The delta function δ(φ)
is identically zero everywhere except at φ = 0. This allows us to rewrite 0 φ<’ ,

equations (1.16) and (1.17) as πφ
1 1
’ ¤φ¤ ,
δ(φ) = 2 + 2 cos (1.23)

ˆ 
δ(x) = δ(φ(x))| φ(x)| (1.19) 0 < φ,
using the one-dimensional delta function δ(φ).
where is determined as above. This delta function allows us to evaluate
2 1
The surface integral (line or point integral in or , respectively) of
the surface integral in equation (1.21) using a standard sampling technique
a function f over the boundary ‚„¦ is de¬ned as
such as the midpoint rule. Similarly, the smeared-out Heaviside function in
equation (1.22) allows us to evaluate the integrals in equations (1.14) and
f (x)δ(x) dx, (1.20)
The reader is cautioned that the smeared-out Heaviside and delta
where the region of integration is all of „¦, since δ prunes out everything
functions approach to the calculus of implicit functions leads to ¬rst-
except ‚„¦ automatically. The one-dimensional delta function can be used
order accurate methods. For example, when calculating the volume of the
to rewrite this surface integral as
region „¦’ using
f (x)δ(φ(x))| φ(x)| dx. (1.21)
(1 ’ H(φ(x))) dV (1.24)
Typically, volume integrals are computed by dividing up the interior
with the smeared-out Heaviside function in equation (1.22) (and f (x) =
region „¦’ , and surface integrals are computed by dividing up the bound-
1), the errors in the calculation are O( x) regardless of the accuracy of
ary ‚„¦. This requires treating a complex two-dimensional surface in three
the integration method used. If one needs more accurate results, a three-
spatial dimensions. By embedding the volume and surface integrals in
dimensional contouring algorithm such as the marching cubes algorithm can
higher dimensions, equations (1.14), (1.15) and (1.21) avoid the need for
be used to identify the region „¦’ more accurately, see Lorenson and Cline
identifying inside, outside, or boundary regions. Instead, the integrals are
[108] or the more recent Kobbelt et al. [98]. Since higher-order accurate
taken over the entire region „¦. Note that dx is a volume element in three
methods can be complex, we prefer the smeared-out Heaviside and delta
spatial dimensions, an area element in two spatial dimensions, and a length
function methods whenever appropriate.
element in one spatial dimension. On our Cartesian grid, the volume of a
three-dimensional cell is x y z, the area of a two-dimensional cell is
x y, and the length of a one-dimensional cell is x.
Consider the surface integral in equation (1.21), where the one-
dimensional delta function δ(φ) needs to be evaluated. Since δ(φ) = 0
almost everywhere, i.e., except on the lower-dimensional interface, which
has measure zero, it seems unlikely that any standard numerical approxi-
mation based on sampling will give a good approximation to this integral.
Thus, we use a ¬rst-order accurate smeared-out approximation of δ(φ).
First, we de¬ne the smeared-out Heaviside function
0 φ<’ ,

H(φ) = 1 + 2 + 2π sin πφ
φ 1
’ ¤φ¤ , (1.22)

1 < φ,
where is a tunable parameter that determines the size of the bandwidth
of numerical smearing. A typically good value is = 1.5 x (making the
18 2. Signed Distance Functions

2 y
x xc
Signed Distance Functions yc

Figure 2.1. xC is the closest interface point to x and y.

xC is the point on the interface closest to y as well. To see this, consider
Figure 2.1, where x, xC , and an example of a y are shown. Since xC is the
closest interface point to x, no other interface points can be inside the large
circle drawn about x passing through xC . Points closer to y than xC must
reside inside the small circle drawn about y passing through xC . Since the
small circle lies inside the larger circle, no interface points can be inside
the smaller circle, and thus xC is the interface point closest to y. The line
2.1 Introduction segment from x to xC is the shortest path from x to the interface. Any local
deviation from this line segment increases the distance from the interface.
In other words, the path from x to xC is the path of steepest descent for the
In the last chapter we de¬ned implicit functions with φ(x) ¤ 0 in the
interior region „¦’ , φ(x) > 0 in the exterior region „¦+ , and φ(x) = 0 on the function d. Evaluating ’ d at any point on the line segment from x to xC
gives a vector that points from x to xC . Furthermore, since d is Euclidean
boundary ‚„¦. Little was said about φ otherwise, except that smoothness is
a desirable property especially in sampling the function or using numerical
approximations. In this chapter we discuss signed distance functions, which
| d| = 1, (2.2)
are a subset of the implicit functions de¬ned in the last chapter. We de¬ne
which is intuitive in the sense that moving twice as close to the interface
signed distance functions to be positive on the exterior, negative on the
gives a value of d that is half as big.
interior, and zero on the boundary. An extra condition of | φ(x)| = 1 is
The above argument leading to equation (2.2) is true for any x as long
imposed on a signed distance function.
as there is a unique closest point xC . That is, equation (2.2) is true ex-
cept at points that are equidistant from (at least) two distinct points on
the interface. Unfortunately, these equidistant points can exist, making
2.2 Distance Functions equation (2.2) only generally true. It is also important to point out that
equation (2.2) is generally only approximately satis¬ed in estimating the
A distance function d(x) is de¬ned as gradient numerically. One of the triumphs of the level set method involves
the ease with which these degenerate points are treated numerically.
d(x) = min(|x ’ xI |) xI ∈ ‚„¦,
for all (2.1)
implying that d(x) = 0 on the boundary where x ∈ ‚„¦. Geometrically, d
may be constructed as follows. If x ∈ ‚„¦, then d(x) = 0. Otherwise, for
2.3 Signed Distance Functions
a given point x, ¬nd the point on the boundary set ‚„¦ closest to x, and
label this point xC . Then d(x) = |x ’ xC |.
A signed distance function is an implicit function φ with |φ(x)| = d(x) for
For a given point x, suppose that xC is the point on the interface closest
all x. Thus, φ(x) = d(x) = 0 for all x ∈ ‚„¦, φ(x) = ’d(x) for all x ∈ „¦’ ,
to x. Then for every point y on the line segment connecting x and xC ,
2.4. Examples 19 20 2. Signed Distance Functions

φ =| x | ’1
and φ(x) = d(x) for all x ∈ „¦+ . Signed distance functions share all the
properties of implicit functions discussed in the last chapter. In addition,
there are a number of new properties that only signed distance functions
possess. For example,
| φ| = 1 (2.3)
as in equation (2.2). Once again, equation (2.3) is true only in a general
„¦+ „¦’ „¦+
sense. It is not true for points that are equidistant from at least two points
φ >0 φ <0 φ >0
on the interface. Distance functions have a kink at the interface where
d = 0 is a local minimum, causing problems in approximating derivatives inside
outside outside
on or near the interface. On the other hand, signed distance functions
are monotonic across the interface and can be di¬erentiated there with
‚„¦ ‚„¦
signi¬cantly higher con¬dence.
φ =0 φ =0
Given a point x, and using the fact that φ(x) is the signed distance to
the closest point on the interface, we can write
xC = x ’ φ(x)N (2.4)

to calculate the closet point on the interface, where N is the local unit
Figure 2.2. Signed distance function φ(x) = |x| ’ 1 de¬ning the regions „¦’ and
normal at x. Again, this is true only in a general sense, since equidistant „¦+ as well as the boundary ‚„¦.
points x have more than one closest point xC . Also, on our Cartesian grid,
equation (2.4) will be only an approximation of the closest point on the
interface xC . Nevertheless, we will ¬nd formulas of this sort very useful.
φ(x) = |x|’1, gives the same boundary ‚„¦, interior region „¦’ , and exterior
Equations that are true in a general sense can be used in numerical ap-
region „¦+ , that the implicit function φ(x) = x2 ’1 did. However, the signed
proximations as long as they fail in a graceful way that does not cause an
distance function φ(x) = |x| ’ 1 has | φ| = 1 for all x = 0. At x = 0 there
overall deterioration of the numerical method. This is a general and pow-
is a kink in our function, and the derivative is not de¬ned. While this may
erful guideline for any numerical approach. So while the user should be
seem problematic, for example for determining the normal, our Cartesian
cautiously knowledgeable of the possible failure of equations that are only
grid contains only sample points and therefore cannot resolve this kink.
generally true, one need not worry too much if the equation fails in a grace-

<< . .

( 23)

. . >>

Copyright Design by: Sunlight webdesign