Training Material

Introduction to Level Sets

By Siri Øyen Larsen, adapted for web by Eigil Samset

Objective

The principal objective of this presentation is to show how level set methods can be applied in segmentation of medical images.

Why level set methods?

Level set methods were first introduced by Osher and Sethian in 1988. They are numerical techniques for tracking evolving interfaces, and are conveniently used in a wide variety of applications. The level set methods handle well interfaces with sharp corners, cusps, topological changes, and three dimensional complications

When applying the level set idea to image segmentation, we are seeking to detect the boundaries of some object or area in the image that we want to segment. This is done by initializing an interface somewhere in the image, and then change it by letting appropriate forces act on it until it finds the correct boundaries in the image. For this reason we will be talking about an interface that is evolving in time. Level set methods stand out from other front-tracking methods in the way they represent the boundaries, as they use an implicit representation of the interface. We will look at how this is done below.

Representation of the interface

Level set representation

As mentioned above the level set method tracks an interface, usually a contour in two dimensions or a surface in three dimensions, although the general idea works in any dimension. But rather than describing the evolution of the interface itself, the level set approach operates on a function in one dimension higher. The interface is thus described as the isocontour of this function, hence it is an implicit representation. In order to model an evolving interface we let the level set function depend on time as well as space.

Definition 1 Let $\phi=\phi(\vec x,t)$ be the level set function. Then the interface $\Gamma$ , at a given point in time, , is given as the set of points in space that corresponds to the zero level isocontour of $\phi$ , i.e. $\Gamma(t) = \{ \vec x : \phi(\vec x,t) = 0 \}$ .

Figure 1 shows an example of a simple level set function of two variables. The level set function is a surface in 3D, while its isocontours are two-dimensional curves. The evolving interface may be described by allowing one of the level set function's isocontours to represent it. Most commonly the zero isocontour is selected, as in definition 1. ¹

$\includegraphics[width=7 cm]{figures/level_set_ex1_.eps}$ $\includegraphics[width=7 cm]{figures/level_set_ex1_contour_.eps}$ $\includegraphics[width=7 cm]{figures/level_set_ex2_.eps}$ $\includegraphics[width=7 cm]{figures/level_set_ex2_contour_.eps}$

Figure 1: Level set representation of an interface at two different points in time. The level set function $\phi$ is shown on the left and the interface corresponding to $\phi = 0$ is plotted on the right. Note that $\phi$ is positive outside the interface, and negative inside. In the above example the interface first yield a single curve, while later in time the topology of the interface has changed - it has splitted into two separate curves.

Now, it should be justified why one might wish to ``complicate'' things with an extra dimension, instead of employing an explicit parametric representation. There are several reasons for this, as we will discuss next.

Parametric representation

Traditional techniques for tracking interfaces include point based methods, where the interface is given through a discrete parameteric representation. The spatial positions of points on the curve or surface at a given point in time are used to reconstruct the evolving boundary front by connecting them with line segments in 2D and (for instance) triangles in three dimensions, [1]. We thus have a representation of the contour given by values $\vec x(s,t)$ at values $s \in D$ , where is the space parameter domain, and is the time parameter. The equation of motion of an interface $\vec x$ is then given by a Langrangian formulation ²

$\displaystyle \frac{d\vec x}{dt} = \vec V(\vec x).$

(1)

Here $\vec V$ defines the velocity of the points on the interface. One difficulty regarding this approach is that the connectivity of the points on the interface must be tracked, and is likely to change during the evolution. Another consideration is to ensure that the discretization is appropriatly coarse at each point in time - it must be fine enough to reconstruct the interface, and at the same time, if the points come too close together, they may cross each other, or cause instability unless the time step is also adjusted accordingly. A remedy for such problems is to redistribute the points every few time steps, and add or remove points where this is necessary, i.e. reparametrize. However, this task quickly becomes quite complicated, especially in three dimensions.

A more serious problem arises if the evolution of the front leads to a change in its topology. The parametric approach is not capable of handling this in a general manner, unless special procedures are implemented for detecting possible splitting and merging of interfaces. This brings us to one of the main advantages of level set methods: since they represent the interface implicitly, they allow changes in topology in an elegant way. Moreover, while a method that uses an explicit representation must be modified when changing from two to three dimensions, the level set strategy works equally well in any dimension.

Slicer

We have used an existing implementation of the level set method for segmentation, developed at the MIT Artificial Intelligence Lab and the Surgical Planning Lab at Brigham & Women's Hospital, an affiliate of Harvard Medical School in Boston, USA. The level set module is incorporated in a program called Slicer, which is an open-source software environment for visualization, registration, segmentation, and quantification of medical data. Slicer is a software research tool. We have used it to load the DICOM images ⁴from the MRI scanner, perform level set segmentation, and generate 3D surface models for viewing the segmented structures. Further information about the Slicer software can be found on their website, www.slicer.org.

The level set equation

The level set function $\phi: \Omega \times [0,\infty) \to \mathbb{R}$ is a scalar valued function of both space and time variables. Since we restrict our attention to the image segmentation problem, in space $\phi$ is defined on the same rectangular domain as the image, $\Omega \subset \mathbb{R}^n$ . Usually we have (a single image) or (an image volume, i.e. a set of slice images). The level set function is initialized at time zero and evolves in time until it stops, hence the function's time domain is $[0,\infty)$ .

$\includegraphics[width=6 cm]{figures/signedDistanceInit.eps}$

$\includegraphics[width=6 cm]{figures/signedDistanceFunc.eps}$

Figure 2: Initialization of the level set function as a signed distance function. The right image illustrates the signed Euclidian distance which is computed from the boundary as initialized on the left.

First of all we must build an initial value for the level set function, i.e. $\phi(\vec x,t=0)$ . This is often done using a so-called signed distance function. If the boundary we are seeking lies in 2D/3D, we might initialize the front as a circle/sphere of small radius. We then construct the initial value

$\displaystyle \phi_0 = \phi(\vec x,t=0) = \pm d,$

where $\pm d$ is the signed Euclidian distance from each point $\vec x \in \Omega$ to the initial front - assigning a positive distance if the point lies outside the region bounded by the front, and negative if inside. For points that lie on the initial interface, the distance is of course zero, hence $\phi = 0$ here, as it should. See figure 2 for an illustration.

We can now describe the motion of the front by matching it with the zero isocontour of the level set function. As outlined in [1] we must always (as the interface evolves) make sure that the level set value of a point on the front with path $\vec x(t)$ is zero, i.e.

$\displaystyle \phi(\vec x(t),t) = 0.$

(2)

Differentiating this equation with respect to time by the chain rule we obtain

$\displaystyle \phi_t + \nabla \phi(\vec x(t), t) \cdot \vec x'(t) = 0.$

(3)

Here $\phi_t$ designates the partial derivative of $\phi$ with repsect to . Now, let denote the speed that drives the evolution, more specifically, is the speed in the outward normal direction to the level set interface. Then $F = \vec x'(t) \cdot \vec n$ , where $\vec n$ is the outward unit normal to the level sets of $\phi$ , hence $\vec n=\nabla \phi/\vert\nabla \phi\vert$ . Manipulating equation 3 we get

$\displaystyle \phi_t + \frac{\nabla \phi}{\vert\nabla \phi\vert} \cdot \vec x'(t)\vert\nabla \phi\vert = 0$

$\displaystyle \phi_t + ( \vec x'(t) \cdot \vec n ) \vert\nabla \phi\vert = 0$

which brings us to the general form of the level set equation

$\displaystyle \phi_t + F \vert \nabla \phi \vert$	$\displaystyle = 0$	(4)
given $\displaystyle \quad \phi(\vec x, 0)$	$\displaystyle = \phi_0 .$

The beginning of this section dealt with how we can specify $\phi_0$ as the initial value of the level set function. In our case we will run the evolution of the levelset in three dimensions, thus we can initialize $\phi_0$ to be zero at one or several small spheres, each sphere being centered at the a location within the object to be segmented. ⁵For the purpose of covering the rest of the function domain we compute $\phi_0$ at each grid point as the signed distance from this location to the closest point on the initial interface.

The goal of the speed function is to act on the interface and ``pull it'' towards the edges of the image. We should therefore model the speed function in such a way that when the interface reaches the desired position, the speed becomes zero, hence $\phi_t = 0$ and the interface must stop. At the same time the speed function should affect $\phi$ so as to make sure it is kept smooth. Clearly, the function must be defined in a way suitable for the given application. In a later section we give a detailed discussion of the implementation that we used for our task.

Solving the equation numerically

After initializing $\phi$ and at the grid points, the interface can be moved across the grid by evolving $\phi$ forward in time by applying numerical methods to update its values. Let $\phi^n=\phi(t^n)$ represent the values of $\phi$ at a given point in time. We must update $\phi$ by finding new values after some time increment $\Delta t$ , i.e. finding $\phi^{n+1}=\phi(t^{n+1})$ where $t^{n+1} = t^n + \Delta t$ . This can be done using a simple first-order accurate method for the time discretization, the forward Euler method, [2]. We compute $\phi^{n+1}$ by approximating $\phi_t$ at time as

$\displaystyle (\phi_t)^n = \frac{\phi^{n+1} - \phi^n}{\Delta t},$

(5)

which when substituted into equation 4 and rearranging terms gives

$\displaystyle \phi^{n+1} = \phi^n - \Delta t ( F^n \vert\nabla \phi^n\vert ),$

(6)

where the superscripts in and $\vert\nabla \phi^n\vert$ denote that the respective functions are evaluated at time . Care must be taken when deciding upon a finite difference scheme to approximate the space derivatives. Since we are using a finite computational grid we are furthermore required to implement boundary conditions. We refer to the literature for further details about the numerical solution of the level set equation. Both Osher & Fedkiw's book [2] and Sethian's book [1] give thorough explanations.

Efficiency

When the interface we are seeking, as defined in Definition 1, is an object in $\mathbb{R}^n$ , we must resolve an -dimensional set in space, which may be computationally expensive if we solve the level set scheme in the straightforward manner. Instead of updating the level set function in the entire computational domain, the task can be relieved in part by focusing on the area of the domain close to the isocontour - the rest of the space may be left unresolved, [2]. This local approach to discretizing implicit representations stores and updates only the values of the level set function $\phi$ at points in a neighbourhood of the zero level set. It is known as the Narrow Band approach (see [1]) since we update the level set function only in an area shaped as a narrow band around the zero level set. Points that are outside the narrow band are given large positive/negative values if they are outside/inside the front itself. ⁶Of course, since the front is moving, the narrow band must also be advanced together with the front. Instead of doing this every time step, we stop and rebuild the narrow band only when the zero level set has reached the edge of the band. To rebuild the band requires that we compute the signed distance to the zero level set, i.e. reinitialize the distance map. Keeping $\phi$ approximately equal to the signed distance (the accuracy of the approximation depending on how often we perform reinitialization) is also desirable for another purpose: it helps to avoid numerical instabilities and to retain a smooth function. Furthermore, maintaining $\phi$ a signed distance function may simplify certain computations of geometric formulas, such as the one for curvature, since we then have $\vert\nabla \phi\vert = 1$ , [2].

The speed function

As an example, the speed function in equation 4 can be decomposed into three terms: a smoothing term , a balloon term , and an advection term , such that

$\displaystyle F = F_s + F_b + F_a.$

(7)

The function is defined on the same domain as the level set function and the image, and changes together with $\phi$ as the interface evolves. We will see that the terms of are functions of the level set function $\phi$ and the image intensity , which represents the inner and outer forces that act on the interface, respectively. The interpretation of the different components of the speed in equation 4 will be explained in the subsections to follow.

In order to understand how the interface develops with the speed function, we should look at what happens to points on the interface as they evolve in the normal direction. Remember that denotes the speed that drives the evolution of the interface in the normal direction, i.e. $F = \vec x'(t) \cdot \vec n$ , where $\vec n = \frac{\nabla \phi}{\vert\nabla \phi\vert}$ as before. Since points on the interface move in the normal direction, we have

$\displaystyle \vec x'(t) = \beta \vec n,$

(8)

for some scalar valued function $\beta$ . If $\beta$ is positive, the points on the interface move in the outward normal direction, while if $\beta$ is negative, they move in the inward normal direction. We may then derive

$\displaystyle F = \vec x'(t) \cdot \vec n = \beta \vec n \cdot \vec n = \beta,$

since $\vec n$ is a unit vector. From this we conclude that the level set equation moves points on the interface (or generally points on any level set of $\phi$ ) in the outward normal direction by a quantity corresponding to evaluated at that point. Alternatively: $\phi_t + F \vert\nabla \phi\vert = 0$ is equivalent to $\vec x'(t) = F \vec n$ . This statement will become useful in understanding how the speed function is defined and will be exploited below.

Most implementations of the level set method allow adjustment of the contribution that each of the terms of the speed function makes to the equation. This can also be done in Slicer by adjusting a coefficient for each of the speed function terms.

Smoothing term

As the name implies, this term controls the smoothness of the contour. The term involves computing the (scalar) curvature $\kappa$ of the level set function $\phi$ , and thus represents an inner energy of the interface (it is a function of $\phi$ alone).

It can be proved that motion under curvature reduces variation in the curve, yielding the wishful smoothing effect. ⁷This is also intuitive, as illustrated by the two dimensional case in figure 3. Imagine that each piece of the curve moves in the outward normal direction with a speed porportional to negative curvature, $-\kappa$ . Parts of the curve where the curvature is positive will thus experience negative speed and move inward, while negative curvature will yield outward motion.

$\includegraphics[width=8 cm]{figures/curvature1.eps}$

$\includegraphics[width=8 cm]{figures/curvature2.eps}$

Figure 3: An example of a level set function in 2D. Left: Convex regions have positive curvature, while concave regions have negative curvature. Right: The red and blue arrows show in which direction the interface move and indicate negative and positive curvature respectively.

When the level set function evolves in three dimensions, one must make a choice of what type of curvature to use, e.g. Gaussian, mean or minimal curvature. The Slicer Level Set module gives a choice of mean or minimal curvature. The minimal curvature is meant especially for segmentation of small structures with high curvatures, or tube-like structures such as arteries.

The mean curvature $\kappa_M$ may be shown to equal the divergence of the normal to the level set function, [1]. The divergence of a vector field is defined as $div ~F = \nabla \cdot F$ , where $\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z})$ . Using subscripts to denote partial derivatives ⁸, we derive the mean curvature $\kappa_M$ as

$\displaystyle \kappa_M = \nabla \cdot \frac{\nabla \phi}{\vert\nabla \phi\vert}$ where $\displaystyle \qquad \vert\nabla \phi\vert = \sqrt{\phi_x^2 + \phi_y^2 + \phi_z^2},$

yielding

$\begin{multline} \kappa_M = \frac{\phi_{xx}\vert\nabla\phi\vert -\phi_x\frac{2... ...i_z\phi_{zz}} {2\vert\nabla\phi\vert}}{\vert\nabla\phi\vert^2} , \end{multline}$

which can be simplified to

$\begin{multline} \kappa_M = \frac{\phi_x^2(\phi_{yy}+\phi_{zz}) + \phi_y^2(\ph... ...i_z\phi_{xz} + 2\phi_y\phi_z\phi_{yz} }{\vert\nabla\phi\vert^3}. \end{multline}$

We sum this up in the following observation.

Observation 1 (Smoothing term $F_s = - \omega \kappa$ ) Motion by curvature is known to have satisfactory geometric smoothing properties. When the coefficient $\omega$ is a positive constant, points on the interface will move in the direction of concavity, i.e. in the direction of negative curvature, and the total variation of the interface will be reduced.

Balloon term

The balloon term, also called the expansion term, can informally be said to control ``how hard the image is pulling on the isocontour''. It has its name from the idea that it ``inflates'' the interface like a balloon. The balloon term is a function of the image, and it often uses information about the image intensity within the object whose boundaries we are about to detect. This function should be high in areas where we are pretty sure there are no boundaries to be found, since the contour can move quickly here. When the zero level set is close to a boundary, on the other hand, the speed caused by the balloon term should be small, so that the evolution of the interface slows down and the other two speed function terms (the smoothing term and the advection term ) should take over the decision of when to stop the deformation. Note that the balloon term is constant in time - it depends exclusively on the image intensity, and not on the level set function, as the two other speed function terms do.

In Slicer, the balloon term is defined as

$\displaystyle c_{I(\vec x)} = \begin{cases}e^{-\frac{(I(\vec x)-\mu)^2}{\sigma^... ...}; \\ 1 - \tau, & \text{if $p_h \neq 0$\ and $I(\vec x) \geq p_h$}. \end{cases}$

(9)

In the above equation $I=I(\vec x)$ is the image as an intensity function, $\mu$ is the mean intensity and $\sigma$ is the standard deviation of the intensity of the object structure to segment. The parameter $\tau$ (called the probability threshold) is a threshold to allow a negative balloon term. We have a choice to set the probability high threshold, , or leave it equal to zero. All the above terms will be explained below.

$\includegraphics[width=6 cm]{figures/balloonforce1_.eps}$

$\includegraphics[width=6 cm]{figures/balloonforce2_.eps}$

Figure 4: Examples of the balloon force as function of the image intensity. Left: 0. Right: = mean intensity.

From 9 we see that the balloon force will be strong in areas of the image with intensity close to the mean, while it will slow down the speed when moving away from the mean. If is defined (user optional), the balloon term is equal to the balloon term at the mean - i.e. the highest possible - for all intensity values higher than . One choice is therefore to set equal to $\mu$ and set $\mu$ equal to the mean intensity of the part of the object that is close to its boundary, meaning that the intensity inside the object may be arbitrarily high (far above the mean) - the force will still be strong in such areas. The motivation for such a specialized implementation is a wish to segment very bright objects (white matter in medical images). Figure 4 shows examples of the balloon term as a function of image intensity.

When using the level set tool in Slicer to segment other types of objects, one can ignore setting the probability high threshold. However, in cases where we want to segment the dark areas of an image, we can exploit this feature by first negating all the intensity values in our images. Figure 5 shows an image, the negated image, and the output from the level set segmentation run on the negated imaged, but overlayed on the original.

$\includegraphics[width=0.9\linewidth]{figures/balloonforce-ill1_.eps}$

$\includegraphics[width=0.9\linewidth]{figures/balloonforce-ill2__.eps}$

Figure 5: MR images from a cryo experiement in an exvivo bovine liver. The level set algorithm is run on the left image, which is the preprocessed and negated magnitude image. On the right, we see the segmentation result on top of the original image.

The last term in equation 9 is $\tau$ , which represents the probability threshold. The exponential term (the Gaussian function) in the equation can never become negative, nor can it become greater than one, which is the maximum value reached at the mean intensity, $I=\mu$ . As mentioned $\tau$ is a threshold to allow negative balloon term, i.e. shrinking of the surface when the intensity is far from the intensity of the structure to segment. This means that the higher $\tau$ is, the more important it is to us that the speed function becomes negative, such that the forces acting on the interface are reversed when it has moved into an area with intensity values far from the mean. Here we assume that the interface is initialized inside the object. A common choice of $\tau$ is in the range .

The observation below yields a summary of the properties of the balloon term.

Observation 2 (Balloon term $F_b = c_{I(\vec x)}$ ) The balloon term is a function of the image intensity. Its purpose is to speed up the interface in areas where we know for sure that there is no edge in the image, by inducing contraction or expansion based on comparing image intensity with the intensity of the object to segment. In areas with high intensity variance, we suspect we are near the boundary, and therefore slow down the speed.

Advection term

Advection in general refers to the transport of some quantity through a vector field. The advection term in the level set segmentation problem uses information about the image to construct a vector field - it represents the outer energy working to deform the interface. This term is therefore said to control the stopping criteria of the interface evolution, i.e. detect when the zero levelset of $\phi$ has reached the boundary of the object of interest. The general form of this term is given as the (negative) inner product

$\displaystyle <tex2html_comment_mark>86 - \nabla g \cdot \vec n$

where $\vec n$ = $\nabla \phi/\vert\nabla \phi\vert$ is the normal of the level set interface and is a function of the image . Remember that $\phi_t + F \vert\nabla \phi\vert = 0$ is equivalent to $\vec x_t = F \vec n$ , thus the advection term contributes to moving points on the level set by $\vec x'(t) = \beta \vec n$ where

$\displaystyle \beta = - \nabla g \cdot \vec n.$

(10)

The function is constructed in such a way that it attracts the interface to the boundaries of the object, thus it is sometimes called an edge detector. Letting be a strictly decreasing function of the gradient of the image, it will take on small values where the image gradient is high, i.e. at the object boundaries. Suitable choices of include

$\displaystyle g(I) = \frac{1}{1+\vert\nabla \hat{I} \vert^p},$

(11)

where $\hat{I}$ is a smoothed version of the image ⁹and or . When the interface is inside the object, and the gradient field of creates a negative inner product with the normal direction of $\phi$ , then $\beta > 0$ , and the interface will move in its outward normal direction towards the object boundaries. If the level set on the other hand crosses the boundary we are seeking, the gradient field of will coincide with the direction of $\nabla \phi$ , yielding a positive inner product and $\beta < 0$ , pulling points on the level set back towards the boundary. Thus the level set interface will eventually be forced to stay on the edges of the image.

$\includegraphics[width=0.9\linewidth]{figures/advection1_.eps}$

$\includegraphics[width=0.9\linewidth]{figures/advection2_.eps}$

Figure 6: The figure shows an image of a simple disk shaped object and a level set interface (in red) initialized inside it. The blue arrows along the curve represent outward normals to the level set, i.e. $\vec n=\nabla \phi/\vert\nabla \phi\vert$ . The second figure shows a close-up look of the edge in the image. The gradient field $\nabla g$ of the edge detector function is overlayed on the image.

Figure 6 shows a simple image together with the vector field of the edge detector function as given by equation 11 with and where $\hat{I}$ is convolved with a Gaussian of size $3 \times 3$ and standard deviation 0.5. An interface $\phi = 0$ and some gradient vectors of $\phi$ along the interface are also plotted. This example serves to illustrate what is pointed out in this section.

The purpose of the advection term may be summed up as follows.

Observation 3 (Advection term $F_a = - \nabla g \cdot \vec n$ ) The advection term's task is to act as an edge detector. It adjusts the position of the interface when close to the real object boundaries.

Termination of the interface evolution

The level set interface is evolved according to the speed function until reduces to (near) zero and the interface stops. Since is composed of three terms, it is basically a problem of balance between the different speed terms. Note that the advection term and the balloon term tend to shrink the contour (by becoming negative) if it crosses the object boundaries.

The evolution is driven by repeatedly updating the level set function $\phi$ for each elapsed time step $\Delta t$ . Hence it is of course possible to stop the evolution after a predefined number of time steps, i.e. without the speed function vanishing.

References

James Albert Sethian. Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science. Cambridge University Press, second edition, 1999.
Stanley Osher and Ronald Fedkiw. Level Set Methods and Dynamis Implicit Surfaces. Springer-Verlag New York, Inc., 2003.
Steven Lobregt and Max A. Viergever. A discrete dynamic contour model. IEEE Transactions on Medical Imaging, 14(1):1224, Mar 1995.
Michael Kass, AndrewWitkin, and Dmitri Terzopoulos. Snakes: Active contour models. Int. J. Comput. Vision, 1(4):321331, 1987.
Chenyang Xu and Jerry L. Prince. Snakes, shapes, and gradient vector flow. IEEE Transactions on Image Processing, pages 359369, Mar 1998.
Karl Krissian and Carl-FredrikWestin. Fast sub-voxel reinitialization of the distance map for level set methods. Technical Report 0001, Department of Radiology, Brigham and Women's Hospital, HarvardMedical School, Laboratory ofMathematics in Imaging, January 2004. ISSN.
James Albert Sethian. Curvature and the evolution of fronts. Communication of Mathematical Physics, 101(4):487499, 1985.

¹There is no specific reason for choosing the zero level isocontour, other than that this is merely a convention. One is free to choose the level that suits the problem. As an example, note that the zero level set of the function .FNv(x, y) = x2 + y2 $B!](B 1 and the Nv = 1 isocontour of both define the same contour, namely the unit circle.

²In a Langrangian formulation the mesh is fixed to the contour, i.e. the discretization follows the contour. In an Eularian formulation on the other hand, the discretization is fixed in space, and the contour must flow through the mesh. The level set formulation is an example of the latter.

⁴ DICOM stands for Digital Imaging and Communications in Medicine. It is a file format for storing medical images.

⁵ The level set module in Slicer (cf. section 'Slicer') allows the user to pick seed points (called fiducial points in this context) interactively by pointing the mouse in the image, and also to choose the radius of the spheres centered at the fiducials.

⁶ There are several variants of the Narrow Band method. Sethian proposes in his book [1] to use Fast Marching Methods for reinitializing, while the Slicer implementation uses a method described by Krissian and Westin in [6].

⁷ For a proof of how motion under curvature decreases the total curvature, see for example [7].

⁸e.g. $\phi_x = \frac{\partial \phi}{\partial x}$ and $\phi_{xy} = \frac{\partial^2\phi}{\partial x \partial y}$

⁹ $\hat{I}$ can be computed applying Gaussian filtering or some other image smoothing filter to .