Kluwer - Handbook of Biomedical Image Analysis Vol

34 sin(3π x) sin(3π y) is computed on

Recent Advances in the Level Set Method

225

to the accepted set, A, only the immediately adjacent grid points require the approximate value of φ to be updated. Let Fmin and Fmax be the minimum and maximum values of the speed function F. For the more general ordered upwind method, all the tentative points in a radius of xFmax/Fmin around x must be updated. If the new approximate value for φ is smaller, this new value is used. This is to account for the possible highest speed direction which could allow the point x to influence grid points within this radius before the immediately adjacent grid points. The formulation for computing the approximation for φ at these tentative points uses the same type of one-sided discretization as used in the fast marching method to follow the characteristics from x.

As an example of the use of the ordered upwind method, the geodesic distance from the origin on the manifold z =

the square [− 12 , 12 ] × [− 12 , 21 ] in the x–y plane. The resulting distance isocontours are shown in Fig. 4.12

4.3.2 Improved Velocity Extensions

The velocity extension method currently in common usage was described in Section 4.2.5, and can be attributed to [3]. However, as noted in [23], the velocity extension characteristics are not supposed to be the straight line extensions that are currently constructed. While it is true that F · φ = 0 should hold at the initial interface, it does not necessarily hold off the interface.

As an example of what can happen with the current velocity extension method, consider the example of an interface consisting of two circles, with the left circle having speed 1, and the right circle having speed 2 (see Fig. 4.13). The current velocity extension method is such that the left half-plane will have

F = 1, and the right half-plane will have F = 2, with the break indicated by the dashed line in Fig. 4.13. The evolution makes a clear error when the right circle expands to the dividing line. Once the circle crosses that line, the velocity extension incorrectly changes the speed from 2 to 1. By noting the gap between successive contours, it is clear that the right-hand circle has slowed down on the left side. The reason that the velocity extension in Fig. 4.13 failed is because the characteristics of the problem were not respected. Once the interface crossed the center line, the velocity came from the left circle, while the characteristics came from the right circle. Ultimately, this happened because the velocity extension was done independent of, and prior to, the actual evolution.

226

Chopp

250
200
150
100
50
50	100	150	200	250

Figure 4.12: A contour map of the distance from the origin on the manifold z = 34 sin(3π x) sin(3π y), computed using the ordered upwind method. Reprinted with permission from [101].

The solution is to do both the fast marching method with the velocity extension at the same time:

F φ = 1,	(4.40)
F · φ = 0.	(4.41)

The discretization of these two equations is the same as before, but the solution method requires some explanation. Again, suppose the values of φ and F are already determined at xi−1, j and xi, j+1. Then Eqs. 4.40 and 4.41 become

	Fi2, j ((Dx−φi, j )2 + (D+y φi, j )2)						= 1,		(4.42)
(D− F	)(D−φ	)	+	(D+ F	)(D+φ	)	=	0.	(4.43)
x i, j	x i, j			y i, j	y i, j

These equations correspond to Eqs. 4.22 and 4.34 respectively, where the unknowns are Fi, j and φi, j , and the remainder of the terms are known. This pair

Recent Advances in the Level Set Method

227

1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1	−0.8	−0.6	−0.4	−0.2	0	0.2	0.4	0.6	0.8	1
−1

Figure 4.13: Example of two circles expanding using the current velocity extension method.

of equations results in a quartic polynomial in Fi, j which can be solved using a Newton solver or by a direct quartic polynomial solver. Once Fi, j is computed, the value of φi, j is easily computed from Eq. 4.43.

The initialization of this method uses a similar bicubic representation as was discussed in Section 4.2.3. The initialization process is based upon the following theorem from [23], and illustrated in Fig. 4.14:

Theorem 1. Suppose = {(x, y) : ax + by = c} and	F0(x, y) = dx + ey + f
for (x, y) with F0 not identically zero on , then the equations
F φ = 1,	(4.44)
F · φ = 0,	(4.45)

228	Chopp
	lines of constant ϕ
F0	lines of constant F
	(characteristics)
initial conditions	Γ
	Γ
	solution
(A, B)	(A, B)

Figure 4.14: Illustration of a sample initial condition and the corresponding solution.

with φ(x, y) = 0, F(x, y) = F0

(x, y), have a solution of the form

db

−

ea

2

2

F(x, y)

X(x, y)

+

Y(x, y)

,

(4.46)

= √

a2 + b2

√

2

b

2

Y(x, y)

a

φ(x, y) =

+

tan−1

.

(4.47)

db

ea

X(x, y)

−

If db − ea = 0, where

X(x, y) =

√

b

(x − A) −

√

a

(y − B),

(4.48)

2

+ b

2

a

2

+ b

2

a

Y(x, y) =

√

a

(x − A) +

√

b

(y − B)

(4.49)

2

+ b

2

2

+ b

2

a

a

ec+ f b

af +cd

. The solution is valid in the set

R2

\ L, where

and where A = ae−bd

, B = bd−ae

L is an arbitrary line passing through the point ( A, B).

If db − ea = 0, then F0(x, y) = F0

is constant on , and the solution be-

comes

F(x, y) = F0,

±1

(4.50)

φ(x, y)

=

(ax

+

by

−

c),

(4.51)

√

2

+ b

2

F0

a

valid on all R2.

Given an initial piece of the interface, the interface is approximated using a linear function, and also the speed, F, along the interface uses a linear

Recent Advances in the Level Set Method

229

approximation. The linear approximation allows the solutions in the theorem to apply, where it is observed that the characteristics travel in circles with varying speed F, and the linear approximation of F designates a center of rotation, ( A, B), depending on where F crosses zero. This leads to a generalized form of Eq. 4.28:

φ(y) × (x − y) =	x − y 2(k · ( F(y) × φ(y)))	.	(4.52)

2F(y)

Note that Eq. 4.28 is recovered if F = 0 is assumed. Equations 4.27 and 4.52 are solved in the same manner as described in Section 4.2.3.

Using the modified velocity extension method on the earlier two-circle example produces the correct results as shown in Fig. 4.15.

Another example that illustrates the difference between the two velocity extension methods is given by an initial circle, with F varying linearly with respect to x, and near zero on the left side. The largest difference between the two methods can be seen on the side where F is small. In the old method,

1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1	−0.8 −0.6 −0.4 −0.2	0	0.2	0.4	0.6	0.8	1
−1
Figure 4.15: Two-circle example with the modified velocity extension method.

230
1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1	−0.8	−0.6	−0.4	−0.2	0	0.2	0.4	0.6	0.8	1
−1

									Chopp
1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1	−0.8	−0.6	−0.4	−0.2	0	0.2	0.4	0.6	0.8	1
−1

Figure 4.16: Comparison of the old (left) and new (right) velocity extension methods.

the interface slows down when it approaches the left side, while with the new method the interface wraps around and merges. The two solutions are shown side by side in Fig. 4.16. The characteristics for this example, represented by the lines of constant F, are shown in Fig. 4.17, illustrating the analogous solution as computed in the theorem. Note how the lines of constant F are orthogonal to the lines of constant φ, as a result of solving Eq. 4.41.

4.3.3 Coupling to Elliptic Solvers

Very often, the speed of the interface is determined by solving an associated elliptic equation, e.g. the pressure equation for incompressible fluid flow. This leads to an elliptic equation which must be solved on an irregularly shaped domain or where there is an internal boundary with jump conditions across the boundary. There are several strategies to handle this problem. When using finite elements to solve this elliptic equation, a mesh is dynamically generated so that it conforms to this irregular boundary. When using finite differences, special delta functions can be added at nodes near the interface to enforce the jump conditions, see e.g. [88].

In the context of the level set method, there are three strategies for setting up and solving the associated elliptic equation. They vary in generality,

Recent Advances in the Level Set Method

231

1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1	−0.8 −0.6 −0.4 −0.2	0	0.2	0.4	0.6	0.8	1
−1
Figure 4.17: Plot of the characteristic curves along which F is constant.

complexity, and accuracy, and provide different advantages. All three strategies are presented here.

The Extended Finite Element Method

The extended finite element method (X-FEM) [29,81,121] is a numerical method to model internal (or external) boundaries without the need for the mesh to conform to these boundaries. The X-FEM is based on a standard Galerkin procedure and uses the concept of partition of unity [80] to accommodate the internal boundaries in the discrete model. The partition of unity method [80] generalized finite element approximations by presenting a means to embed local solutions of boundary-value problems into the finite element approximation.

For a standard finite element approximation, consider a point x of Rd that lies inside a finite element e. Denote the nodal set N = {n1, n2, . . . , nm}, where m is the number of nodes of element e. The approximation for a vector-valued

232			Chopp
function u(x) : Rd → Rd assumes the form
uh(x) =	I	φI (x)uI , (uI Rd),	(4.53)
	nI
	N

where the functions φI (x) are the finite element basis functions and uI are the weights.

The extended finite element method uses enrichment functions, extra basis functions which are sensitive to prescribed boundaries, to capture the boundary conditions and improve the solution in the neighborhood of regions which would otherwise require greater spatial resolution. Consider again a point x that lies inside a finite element e. The enriched approximation for the function u(x) becomes

uh(x) =

I

φI (x)uI +

J

φJ (x)ψ (x)aJ ,

(4.54)

nI N

nJ Ng

classical

enriched

where the nodal set Ng consists of nodes which are on elements cut by the boundary, for example, see Fig. 4.18. In general, the choice of the enrichment function ψ (x) that appears in Eq. 4.54 depends on the geometry, the boundary condition, and the elliptic equation being solved.

To illustrate the effectiveness of this approach, consider the following simple example. Suppose we wish to solve the radial heat equation on an annulus given

Figure 4.18: Example of choosing enriched nodes. Enriched nodes are indi-

cated by gray dots.

Recent Advances in the Level Set Method			233
by
1		ur = 0, 0 < ≤ r < L,	(4.55)
urr +
	r
ur ( ) = −10, u(L) = 0.			(4.56)
The exact solution is given by
u(r) = −10 ln(r) + 10 ln(L).			(4.57)

If we solve this equation for = 0.01, L = 9 using a standard finite element method with linear elements and with nodes at r = 0, . . . , 9, the solution for

≤ r < 1 is very unsatisfactory, as shown in Fig. 4.19. However, by using a simple enrichment function ψ1(r) = ln(r), and using this enrichment function on the first two nodes (located at r = 0, 1), dramatically better results are achieved (Fig. 4.19). Of course, refining the finite element mesh would also improve the results, but this requires remeshing as the interface (in this example the left boundary) moves. The X-FEM achieves this accuracy without remeshing.

The merits of coupling level sets to the extended finite element method were first explored in [118], and subsequently its advantages further realized in [53, 61, 82, 115, 117, 120]. The two methods make a natural pair of methods where:

1.Level sets provide greater ease and simplification in the representation of geometric interfaces.

2.The X-FEM, given the right enrichment functions, can accurately compute solutions of elliptic equations which are often required for computing the interface velocity.

3.Geometric computations required for evaluating the enrichment functions (such as the normal or the distance to the interface) are readily computed from the level set function [120].

4.The nodes to be enriched are easily identified using the signed distance construction of the level set function [115, 117, 118, 120].

Compared to the other methods to follow, this algorithm is more complex, but it is also much more general. Through the use of enrichment functions, this method provides a much better solution near the interface, providing subgrid resolution in that region without requiring additional mesh refinement. This is

234

Chopp

2.5
						Finite Element Method
						Exact Solution
						Extended Finite Element Method
2
1.5
u
1
0.5
0	1	2	3	4	5	6	7	8	9
0

							r
							(a)
	1
							Finite Element Method
	0.9						Exact Solution
							Extended Finite Element Method
	0.8
	0.7
	0.6
u	0.5
	0.4
	0.3
	0.2
	0.1
	0	ε = 0.01
		0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
	0

r

(b)

Figure 4.19: Solutions of the radial heat equation: (a) whole domain ≤ r ≤ L

and (b) across first three nodes.