Level Set Models

A level set model[30,29] specifies a surface as a level set (iso-surface) of a scalar volumetric function, $\phi:U \mapsto \Re$, where $U \subset \Re^{3}$ is the range of the surface model. Thus, a surface $S$ is

$$S = \left\{{\mbox{\boldmath$\displaystyle s$}} \vert \phi({\mbox{\boldmath$\displaystyle s$}}) = k \right\},$$ (9)

and $k$ is the isovalue. In other words, $S$ is the set of points ${\mbox{\boldmath $\displaystyle s$}}$ in $\Re^{3}$ that composes the $k$th iso-surface of $\phi$. The embedding $\phi$ can be specified as a regular sampling on a rectilinear grid. The surfaces may propagate with (time-varying) curvature-dependent speeds. Level set methods provide the mathematical and numerical mechanisms for computing surface deformations as iso-values of $\phi$ by solving a partial differential equation on the 3D grid ($U$). That is, the level set formulation provides a set of numerical methods that describes how to manipulate the grey-scale values in a volume, so that the iso-surfaces of $\phi$ move in a prescribed manner. See Figure 6.

There are two different approaches to defining a deformable surface from a level set of a volumetric function as described in Equation (9). Either one can think of $\phi({\mbox{\boldmath $\displaystyle s$}})$ as a static function and change the iso-value $k(t)$ or alternatively fix $k$ and let the volumetric function dynamically change in time, i.e. $\phi({\mbox{\boldmath $\displaystyle s$}},t)$. Following the second approach, we can mathematically express the dynamic model as

$$\phi({\mbox{\boldmath$\displaystyle s$}},t)=k.$$ (10)

To transform this definition into partial differential equation that can easily be solved by standard numerical techniques, we differentiate both sides of Equation (10) with respect to time $t$, and apply the chain rule:

$$\frac{\partial\phi({\mbox{\boldmath$\displaystyle s$}},t)}{\parti... ...yle s$}},t)\cdot \frac{{\rm d}{\mbox{\boldmath$\displaystyle s$}}}{{\rm d}t}=0.$$ (11)

Equation (11) is sometimes referred to as a ``Hamilton-Jacobi-type'' equation and defines an initial value problem for the time-dependent $\phi$. Let ${\rm d}{\mbox{\boldmath $\displaystyle s$}}/ {\rm d}t$ be the movement of a point on a surface as it deforms, such that it can be expressed in terms of the position of ${\mbox{\boldmath $\displaystyle s$}} \in U$ and the geometry of the surface at that point, which is, in turn, a differential expression of the implicit function, $\phi$. This gives a partial differential equation (PDE) on $\phi$: ${\mbox{\boldmath $\displaystyle s$}}\equiv{\mbox{\boldmath $\displaystyle s$}}(t)$

$$\frac{\partial \phi }{\partial t} = -\nabla \phi \cdot \frac{{\rm... ...ox{\boldmath$\displaystyle s$}}, \mathrm{D} \phi, \mathrm{D}^{2} \phi, \ldots),$$ (12)

where ${\mbox{\boldmath $\displaystyle F$}}$ is a user-defined ``speed'' term which generally depends on a set of order-$n$ derivatives of $\phi$, $D^{n}\phi$, evaluated at ${\mbox{\boldmath $\displaystyle s$}}$, as well as other functions of ${\mbox{\boldmath $\displaystyle s$}}$. Typically ${\mbox{\boldmath $\displaystyle F$}}({\mbox{\boldmath $\displaystyle x$}})$ combines a data term with a smoothing term, which prevents the solution from fitting too closely to noise-corrupted data. There are a variety of surface-motion terms that can be used in succession or simultaneously in a linear combination to form ${\mbox{\boldmath $\displaystyle F$}}({\mbox{\boldmath $\displaystyle x$}})$. For the work presented in this paper, we combine a feature attraction term and a smoothing term weighted by a factor $\beta$,[28]

$${\bf F} = {\bf F}_{attr} + \beta {\bf F}_{curv}.$$ (13)

The first term ${\bf F}_{attr}$ is due to the attraction to the edges in the volume. It attracts the surface models to certain grey scale features in the input data. For instance, the gradient magnitude indicates areas of high contrast in volumes. By following the gradient of such grey scale features, surface models are drawn to minimum or maximum values of that feature. Typically grey scale features, such as the gradient magnitude are computed with a scale operator, e.g., a derivative-of-Gaussian kernel. If models are properly initialized, they can move according to the gradient of the gradient magnitude and settle onto the edges of an object at a resolution that is finer than the original volume. For this work we used the attraction force

$${\bf F}_{attr} = \nabla \vert(\nabla (G * I({\bf x}))\vert,$$ (14)

where the volume data $I({\bf x})$ is convolved with a Gaussian kernel $G$ with $\sigma \approx 0.5$, such that a positive sign moves surfaces towards maxima and the negative sign towards minima.

There are a variety of options for the curvature smoothing terms in Equation (13), and the question of efficient, effective higher-order smoothing terms is the subject of on-going research[29]. For the work presented in this paper the smoothing term uses the mean curvature $\mathcal{K}_M$ of the level set ${\bf S}$ to form a vector in the direction of the surface normal $\bf {n}$:

$${\bf F}_{curv} = \mathcal{K}_M \bf {n} = (\nabla \cdot {\bf n}){\... ...phi}{\vert\nabla \phi\vert} \right) \frac{ \nabla \phi}{\vert\nabla \phi\vert}.$$ (15)

It is weighted by a factor $\beta$, allowing the user to control the amount of smoothing, and is tuned for each dataset. The level set propagation stops when the ${\bf F}_{attr}$ and $\beta{\bf F}_{curv}$ terms cancel each other, or when the number of computational iterations reaches a user-specified value.

Level set models have a number of practical and theoretical advantages over conventional surface models, especially in the context of deformation and segmentation. 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) of degrees of freedom, and therefore they can accommodate complex shapes. Indeed, the shapes formed by the level sets of $\phi$ are restricted only by the resolution of the sampling. Thus, there is no need to re-parameterize 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 paper to give a detailed mathematical description of such a numerical implementation. Rather we shall give a short outline below and refer to the actual source code which is publicly available¹.

Equation (12) to (15) can be solved using finite forward differences if one uses the up-wind scheme, proposed by Osher and Sethian [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 Equation (10). However, this 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 iso-values 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 [24].In previous work [25] 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 sub-voxel accuracy.

 

Figure 3:

Figure 3. Level set models represent curves and surfaces implicitly using grey scale images. For example an ellipse is represented as the level set of an image (top). To change the shape of the ellipse we modify the grey scale values of the image by solving a PDE (bottom).