LINEBURG

implicit function to be as smooth as possible. In the next chapter we dis-

N

cuss using a signed distance function to represent the surface. This turns

x

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

‚x

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

oo

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-

oo

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

inside

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,

1

’

χ (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)

ˆ

δ(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

function

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)

(1.15).

„¦

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)

2

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

distance,

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

x

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

interface

interface

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-

Copyright Design by: Sunlight webdesign