## Level Set Surface Models

When considering deformable models for segmenting 3D volume data, one is faced with a choice from a variety of surface representations, including triangle meshes [19,20], superquadrics [21-23], and many others [18,24-29]. Another option is an implicit level set model, which specifies the surface as a level set of a scalar volumetric function, \$ : U ^ IR, where U c IR3 is the range of the surface model. Thus, a surface S is

with an isovalue k. In other words, S is the set of points s in IR3 that composes the kth isosurface of \$. The embedding \$ can be specified as a regular sampling on a rectilinear grid.

Our overall scheme for segmentation is largely based on the ideas of Osher et al. [30] that model propagating surfaces with (time-varying) curvature-dependent speeds. The surfaces are viewed as a specific level set of a higher dimensional function \$â€”hence the name level set methods. These methods provide the mathematical and numerical mechanisms for computing surface deformations as isovalues of \$ by solving a partial differential equation on the 3D grid. That is, the level set formulation provides a set of numerical methods that describes how to manipulate the grayscale values in a volume, so that the isosurfaces of \$ move in a prescribed manner (shown in Fig. 8.1). This chapter does not present a comprehensive review of level set methods, but merely introduces the basic concepts and demonstrates how they may be applied to

Figure 8.1: (a) Level set models represent curves and surfaces implicitly using grayscale images. For example, an ellipse is represented as the level set of an image shown here. (b) To change the shape of the ellipse we modify the grayscale values of the image by solving a PDE.

Figure 8.1: (a) Level set models represent curves and surfaces implicitly using grayscale images. For example, an ellipse is represented as the level set of an image shown here. (b) To change the shape of the ellipse we modify the grayscale values of the image by solving a PDE.

the problem of volume segmentation. For more details on level set methods see [7,31].

There are two different approaches to defining a deformable surface from a level set of a volumetric function as described in Eq. (8.1). Either one can think of p(s) as a static function and change the isovalue k(t) or alternatively fix k and let the volumetric function dynamically change in time, i.e. p(s, t). Thus, we can mathematically express the static and dynamic models respectively as p(s) = k(t), (8.2a)

To transform these definitions into partial differential equations which can be solved by standard numerical techniques, we differentiate both sides of Eq. (8.2) with respect to time t, and apply the chain rule:

The static equation (8.3a) defines a boundary value problem for the time-independent volumetric function p. This static level set approach has been solved [32,33] using "Fast Marching Methods." However, it inherently has some serious limitations following the simple definition in Eq. (8.2a). Since p is a function (i.e. single-valued), isosurfaces cannot self-intersect over time, i.e. shapes defined in the static model are strictly expanding or contracting over time. However, the dynamic level set approach of eq. (8.3b) is much more flexible and shall serve as the basis of the segmentation scheme in this chapter. Equation (8.3b) is sometimes referred to as a "Hamilton-Jacobi-type" equation and defines an initial value problem for the time-dependent p. Throughout the remainder of this chapter we shall, for simplicity, refer to this dynamical approach as the level set method, and not consider the static alternative.

Thus, to summarize the essence of the (dynamic) level set approach, let ds/dt be the movement of a point on a surface as it deforms, such that it can be expressed in terms of the position of s e U and the geometry of the surface at that point, which is, in turn, a differential expression of the implicit function, p.

This gives a partial differential equation on 0: s = s(t)

d0 ds 9

91 dt ds

dt where F() is a user-created "speed" term that defines the speed of the level set at point s in the direction of the local surface normal n at s. F() may depend on a variety of local and global measures including the order-n derivatives of 0, Dn0, evaluated at s, as well as other functions of s, n, 0, and external data. Because this relationship applies to every level set of 0, i.e. all values of k, this equation can be applied to all of U, and therefore the movements of all the level set surfaces embedded in 0 can be calculated from Eq. (8.4).

The level set representation has a number of practical and theoretical advantages over conventional surface models, especially in the context of deformation and segmentation. First, level set models are topologically flexible, they easily represent complicated surface shapes that can form holes, split to form multiple objects, or merge with other objects to form a single structure. These models can incorporate many (millions) degrees of freedom, and therefore they can accommodate complex shapes such as the dendrite in Fig. 8.7. Indeed, the shapes formed by the level sets of 0 are restricted only by the resolution of the sampling. Thus, there is no need to reparameterize the model as it undergoes significant changes in shape.

The solutions to the partial differential equations described above are computed using finite differences on a discrete grid. The use of a grid and discrete time steps raises a number of numerical and computational issues that are important to the implementation. However, it is outside of the scope of this chapter to give a detailed mathematical description of such a numerical implementation. Rather we shall provide a summary in a later section and refer to the actual source code which is publicly available5.

Equation (8.4) can be solved using finite forward differences if one uses the up-wind scheme, proposed by Osher et al. [30], to compute the spatial derivatives. This up-wind scheme produces the motion of level set models over the entire range of the embedding, i.e., for all values of k in Eq. (8.2). However, this

5 The level set software used to produce the morphing results in this chapter is available for public use in the VISPACK Libraries at http://www.cs.utah.edu/~whitaker/vispack.

method requires updating every voxel in the volume for each iteration, which means that the computation time increases as a function of the volume, rather than the surface area, of the model. Because segmentation requires only a single model, the calculation of solutions over the entire range of isovalues is an unnecessary computational burden.

This problem can be avoided by the use of narrow-band methods, which compute solutions only in a narrow band of voxels that surround the level set of interest [34,35]. In a previous work [36] we described an alternative numerical algorithm, called the sparse-field method, that computes the geometry of only a small subset of points in the range and requires a fraction of the computation time required by previous algorithms. We have shown two advantages to this method. The first is a significant improvement in computation times. The second is increased accuracy when fitting models to forcing functions that are defined to subvoxel accuracy.