Kluwer - Handbook of Biomedical Image Analysis Vol

important when having to interpolate the data to determine the front speed on the boundary contour. Of the three methods, this is the only one that has this capability.

The Immersed Interface Method

The immersed interface method, introduced by LeVeque and Li [74], has also been coupled to the level set method [76, 78]. Like the X-FEM described above, the immersed interface method is designed to solve elliptic equations which arise in a variety of physical applications. The advantage of the immersed interface method is that it is second-order accurate, even near the interface where jump conditions may appear.

The immersed interface method is designed to solve equations of the form

where the coefficient functions β, κ, and f may have discontinuities across an interface . The function f may also have a delta function singularity, which often arises, for example, from surface tension in multiphase flow.

The key idea in the immersed interface method is to modify the discretization of Eq. 4.58 in such a way that the jump discontinuities and singularities are accounted for, leading to a fully second-order method. At points away from the interface, where the coefficient functions and the solution are smooth, the standard central difference approximation is used. However, for grid points which are near the interface, an additional grid point is added to the usual central difference stencil to account for a second-order Taylor approximation around a point on the interface.

To illustrate how this method works, consider the one-dimensional problem

where u− is the value of u on the interval [0, α], and u+ is the value of u on the interval [α, 1]. Suppose that the point α is located between the uniformly spaced grid points xi and xi+1. The idea is to calculate coefficients γi−1, γi, γi+1, and an

is second-order accurate, with jump conditions Eqs. 4.60 and 4.61.

To determine the γi’s, Taylor expansions are taken about the point x = α to

get

These expansions are inserted into Eq. 4.62, and the u+ terms are eliminated from the equation by using the jump conditions Eqs. 4.60 and 4.61, combined with the equation

which comes from the continuity of f in Eq. 4.59. The function f on the right side

The coefficients γi−1, γi, γi+1, and Ci are now chosen so that Eq. 4.67 holds up to second order. This leads to the following equations:

Figure 4.20: Choice of stencil for (a) points not crossed by the interface and (b) points where the interface crosses the stencil. Dashed lines indicate the points used in the stencil.

These equations are solved for γi−1, γi, γi+1, and Ci, thus determining the numerical approximation corresponding to the point xi using Eq. 4.62. A similar process is followed for the approximation centered at xi+1. This results in a specialized discretization at these two points and standard central difference approximations everywhere else.

For higher dimensional problems, a similar approach is taken. At grid points not crossed by the interface, the standard central difference stencil is used (see Fig. 4.20(a)) to approximate Eq. 4.58. At grid points where the interface crosses through the stencil, an additional grid point is chosen across the interface from the center of the stencil (see Fig. 4.20(b)).

When building the specialized discretization for the stencil at grid points as in Fig. 4.20(b), a point (x , y ) is chosen for the point around which the approximation will be computed, and around which all Taylor expansions will be taken. Usually, the point (x , y ) is the point on the interface closest to the center of the stencil (in this example, point 2). Once (x , y ) is chosen, a coordinate transformation is taken so that the interface normal maps onto the x-axis. Once this coordinate transformation is completed, the computation of the stencil is similar to the one-dimensional case described above.

As noted earlier, the advantage of this method is that it is truly second-order accurate, even in the neighborhood of the interface. However, the stencil that is produced is irregular, and it sometimes can be difficult to solve the resulting linear system. Also, the choice of the points (x , y ) is somewhat arbitrary, and

it is not clear what the best choices should be. Nonetheless, the method has been used successfully in a number of applications, e.g. see the review in [76].

The Ghost Point Method

The ghost point method [50] is another method designed to solve elliptic equations with irregular and moving boundaries represented by the level set method. The idea behind this method is similar to the use of what are often called ghost points for discretizing boundary conditions in finite difference methods. In this context, ghost points are grid points located outside the computational domain, and are used to enforce boundary conditions.

The method presented in [50] is designed to solve equations of the form

in an irregularly shaped domain , where β and f are smooth functions defined on , and g is defined on ∂ , the boundary of . This is a more restrictive class of problems than can be handled by the previous two methods described, but it is a class of problems that often arises. By focusing on this simpler class, a second-order method with a simple discretization can be employed, which uses a stencil that has properties which make it easier to solve numerically than the system created by the previous methods.

To illustrate this method, consider first the one-dimensional problem

with ∂ = xI , and u(xI ) = uI . Assume xI lies between the two grid points xi and xi+1. For points xj in the interior of the domain, the central difference discretization, similar to the one used in the immersed interface method, is used:

At the boundary, the discretization Eq. 4.74 is again employed, but the value of ui+1 is not defined because xi+1 is outside of . Instead, a ghost value for ui+1 is computed from the boundary condition using a linear extrapolation:

For stability reasons, if θ < x, then Eq. 4.75 is replaced with ui+1 = uI . Using Eq. 4.75 in Eq. 4.74 produces the following discretization for the point near the boundary:

In multiple dimensions, this same extrapolation technique is carried out along each coordinate direction.

The resulting discretization is only first-order accurate near the boundary, but is second-order accurate overall. This is due to the confinement of the firstorder error to the nodes adjacent to the boundary. On the other hand, the linear system that comes from this discretization can be solved using faster conjugate gradient-type algorithms. Increasing the order of the extrapolation to compute ui+1 can result in a linear system that is more difficult to solve numerically, because of the non-symmetric stencil, and hence is not preferred.

This method is used primarily for its simplicity, while still yielding secondorder convergence overall. For problems where the accuracy at the boundary is critical, this is probably not the preferred method, especially if the solution is difficult to resolve near the boundary. The method has been used in a handful of applications, for example, see [124].

Comparison of the Elliptic Equation Solvers

The algorithms presented here, for solving elliptic equations in conjunction with the level set method, vary significantly in sophistication, complexity, and capability. The X-FEM approach is by far the most difficult to construct, but is also the most general, and has the greatest potential to solve challenging problems. In particular, the X-FEM approach provides a much more accurate representation of the solution near the boundary, a property that is of critical importance when the velocity of the interface depends on this very value.

The immersed interface method and ghost point method, on the other hand, are built much more easily, and still produce accurate solutions. The immersed interface method handles a larger range of equations than does the ghost point method, which is the most restrictive in this regard. Between these two methods, the immersed interface method is more accurate at the boundary, but at the expense of a more difficult system of equations to solve numerically.

The ghost point method is probably the fastest, due to its use of faster linear solvers, but an actual direct comparison has not been done. Both the immersed interface method and ghost point method will be faster than the X-FEM approach on the same mesh. However, to obtain the same accuracy near the interface, the X-FEM will not require as fine a mesh as the others, and hence can make up the difference in time by using a coarser mesh to obtain comparable results. A direct comparison of these three methods is the subject of current research.

4.3.4 Particle Level Set Method

Another modification of the level set method, called the particle level set method, was proposed by Enright et al. in [38]. In the particle level set method, the level set function is compared with the motion of particles which move along the characteristics of the same velocity field. For an interface which is passively advected using the same velocity field, the particles, in theory, should not cross the interface. By comparing the motion of the particles with the moving interface, problems with the location of the interface can be identified and corrected.

Suppose the interface velocity is determined by a velocity field v(x, t). Given this velocity, the interface speed function, F, in Eq. 4.5 is given by

Substituting this expression for F into Eq. 4.77 gives the passive interface advection equation

At the same time, the particles themselves are moving with this same velocity, v. These two evolutions are coupled together when the particles are checked to see if any has crossed the interface, which in this case indicates that a particle has moved from a point where φ > 0 to a point where φ < 0, or vice versa. At that point, the level set function is “corrected.”

In [38], a large number of particles are randomly distributed uniformly in the neighborhood of the interface φ = 0. Each particle, p, is assigned a sign, sp, to indicate whether it is starting where φ > 0 or φ < 0, and is also assigned its distance, rp, to the interface. As the evolution of the interface and the particles

proceeds, the particle locations are periodically checked to determine whether they have strayed across the level set function interface.

When a particle is determined to have strayed sufficiently far across the level set interface, the interface is reconstructed using the particle information. To do this, each particle, p, located at the point xp, is assigned a local signed distance function

The level set function is now reconstructed in two steps. First, the functions φ+ and φ− are computed where

and where P+ and P− are the sets of points which were assigned positive and negative sp respectively. The final φ function is now recovered from φ+ and φ− by the equation

where

a, |a| < |b|

absmin(a, b) = . (4.83) b, |b| ≤ |a|

There is no guarantee that the resulting reconstructed level set function will be a signed distance function, so if this is desired, a reinitialization step will be applied to reform φ into a signed distance function.

What is novel about this approach is the use of the Lagrangian and Eulerian methods to play against each other to ensure proper interface motion. However, one must carefully determine when the particle solution is correct, versus the level set evolution. This is determined by checking the local characteristics to see if they are colliding or expanding. The level set evolution tends to be better when characteristics are colliding, whereas the particle method will be more reliable when the interface is moving tangentially or stretching. Nonetheless, this combination tries to extract the positive capabilities of both the Lagrangian and Eulerian types of approaches to interface motion, while discounting the negatives.

practitioner to consider alternative numerical methods. It may or may not be the case that the level set method is still the best choice.

For a more comprehensive discussion on the level set method, the interested reader is directed to the books by Sethian [104] (which also includes the fast marching method) and Osher and Fedkiw [84].

Questions

1.What are the main advantages of the level set method?

2.What is the importance of the connection between the level set method and hyperbolic conservation laws?

3.What is the difference between the level set method and the fast marching method?

4.Why are triple junctions a problem for the level set method?

5.What is the primary purpose of reinitialization, and why is it important to do it as accurately as possible?

6.What is the alternative to using repeated reinitializations?

7.What kinds of problems can be solved by the general ordered upwind method that could not be solved by the fast marching method?

8.What is the difference between the original velocity extension and the new velocity extension methods?

9.Can the level set method be implemented using the finite element method?

10.What is the advantage of using the X-FEM over a standard finite element formulation?

11.Is the level set method appropriate for all interface propagation problems?

Bibliography

[1]Adalsteinsson, D. and Sethian, J. A., Transport and diffusion of material quantities on propagating interfaces via level set methods, J. Comput. Phys., Vol. 185, pp. 271–288, 2003.

[2]Adalsteinsson, D. and Sethian, J. A., A fast level set method for propagating interfaces, J. Comput. Phys., Vol. 118, No. 2, pp. 269–277, 1995.

[3]Adalsteinsson, D. and Sethian, J. A., The fast construction of extension velocities in level set methods, J. Comput. Phys., Vol. 48, No. 1, pp. 2– 22, 1999.

[4]Alkhalifah, T., Traveltime computation with the linearized eikonal equation for anisotropic media, Geophys. Prospecting, Vol. 50, pp. 373–382, 2002.

[5]Allaire, G., Jouve, F., and Toader, A. M., A level-set method for shape optimization, C. R. Math., Vol. 334, No. 1125–1130, 2002.

[6]Antiga, L., Ene-Iordache, B., and Remuzzi, A., Computational geometry for patient-specific reconstruction and meshing of blood vessels from mr and ct angiography, IEEE Trans. Med. Imaging, Vol. 22, pp. 674–684, 2003.

[7]Balabel, A., Binninger, B., Herrmann, M., and Peters, N., Calculation of droplet deformationby surface tension effects using the level set method, Combust. Sci. Technol., Vol. 174, pp. 257–278, 2002.

[8]Bassano, E., Numerical simulation of thermo-solutal-capillary migration of a dissolving drop in a cavity, Int. J. Numer. Methods Fluids, Vol. 41, pp. 765–788, 2003.

[9]Bazdidi-Tehrani, F., and Zaman, S., Two-phase heat transfer on an isothermal vertical surface: a numerical simulation, Int. J. Heat Fluid Flow, Vol. 23, pp. 308–316, 2002.

[10]Bertalmio, M., Cheng, L. T., Osher, S., and Sapiro, G., Variational problems and partial differential equations on implicit surfaces, J. Comput. Phys., Vol. 174, pp. 759–780, 2001.