A Fast Level Set based Algorithm for TopologyIndependent
Shape Modeling *
Ravikanth Malladi,1 James A. Sethian,2 and Baba C. Vemuri1
1Department of Computer & Information Sciences
University of Florida, Gainesville, FL 32611
2Department of Mathematics
University of California, Berkeley, CA 94720
Abstract
Shape modeling is an important constituent of computer vision as well as computer graphics
research. Shape models aid the tasks of object representation and recognition. This paper
presents a new approach to shape modeling which retains the most attractive features of existing
methods, and overcomes their prominent limitations. Our technique can be applied to model
arbitrarily complex shapes, which include shapes with significant protrusions, and to situations
where no a priori assumption about the object's topology is made. A single instance of our
model, when presented with an image having more than one object of interest, has the ability
to split freely to represent each object. This method is based on the ideas developed by Osher
& Sethian to model propagating solid/liquid interfaces with curvaturedependent speeds. The
interface (front) is a closed, nonintersecting, hypersurface flowing along its gradient field with
constant speed or a speed that depends on the curvature. It is moved by solving a "Hamilton
.1 ....1. type equation written for a function in which the interface is a particular level set. A
speed term synthesized from the image is used to stop the interface in the vicinity of the object
boundaries. The resulting equation of motion is solved by employing entropysatisfying upwind
finite difference schemes. We also introduce a new algorithm for rapid advancement of the front
using what we call a narrowband updation scheme. This leads to significant improvement in the
time complexity of the shape recovery procedure. An added advantage of our modeling scheme
is that it can easily be extended to any number of space dimensions. The efficacy of the scheme
is demonstrated with numerical experiments on low contrast medical images.
*Submitted to the special issue of the Journal of Mathematical Imaging & Vision on Topology and Geometry in
Computer Vision.
1 Introduction
An important goal of computational vision is to recover the shapes of objects in 2D and 3D from
various types of visual data. One way to achieve this goal is via modelbased techniques. Broadly
speaking, these techniques, as the name suggests, involve the use of a model whose boundary
representation is matched to the image to recover the object of interest. These models can either
be rigid or nonrigid. In the former case, we have what are popularly known as correlationbased
template matching techniques, while the later involves a dynamic model fitting process to the
data. In this paper, we present a new approach to dynamic shape modeling which retains the most
attractive features of existing methods, and overcomes their prominent limitations. For the rest of
this paper we will use the word shape model for a boundary (surface) representation of an object
shape.
Shape reconstruction typically precedes the symbolic representation of surfaces and shape models
are expected to aid the recovery of detailed structure from noisy data using only the weakest
of the possible assumptions about the observed shape. To this end, several variational shape
reconstruction methods have been proposed and there is abundant literature on the same [2, 22,
28, 3, 29, 14]. Generalized spline models with continuity constraints are well suited for fulfilling
the goals of shape recovery (see [3, 4, 26]). Generalized splines are the key ingredient of the
dynamic shape modeling paradigm introduced by Terzopoulos et al., [27]. Incorporating dynamics
into shape modeling enables the creation of realistic animation in computer graphics applications
and for tracking moving objects in computer vision. Following the advent of the dynamic shape
modeling paradigm, there was a flurry of research activity in the area, with numerous application
specific modifications to the modeling primitives, and external forces derived from data constraints
[10, 30, 31, 32, 7]. However, the aforementioned schemes for shape modeling have two serious
limitations the dependence of the final surface shape on the initial guess made to start the
numerical reconstruction procedure, and an assumption that the object is of a known topology.
The first of these deficiencies stems from the fact that the nonconvex energy functionals used in the
variational formulations have multiple local minima. As a consequence of this feature, the numerical
procedures, for convergence to a satisfactory solution require an initial guess which is "reasonably"
close to the desired shape. Existing shape representation schemes have an additional shortcoming
in that they lack the ability to dynamically sense the topological changes that may occur in the
shape of interest during the shape recovery process. In this paper, we present a modeling scheme
that makes no assumption about the object's topology, and leads to a numerical algorithm whose
convergence to the desired shape is completely independent of the shape initialization. Our method
can also recover shapes whose topology changes over time, e.g., the cell boundary in cell division.
The framework of energy minimization has also been used successfully in the past for extracting
salient image contours edges and lines. Kass et al. [10] used energyminimizing "snakes" that
exhibit an affinity toward image features such as edge points and edge segments, while maintaining
piecewise smoothness. The weights of the smoothness and image terms in the energy functional
can be .1,1 1. .1 for different kinds of behavior. Snakes, also referred to as active contour models,
are restricted examples of the more general techniques of matching deformable surface models to
image data by means of energy minimization [27]. In this modeling scheme, one seeks to design
energy functionals whose local minima comprise the set of alternative solutions available to high
level processes. In the absence of a welldeveloped highlevel mechanism to make a choice among
these solutions, an interactive approach is used to explore the alternatives. By adding suitable
energy terms to the minimization, the user pushes the model out of a local minimum toward the
desired solution. In the problem area of automatic segmentation of noisy images, snakes perform
poorly unless they are placed close to the preferred shapes. In a move to make the final result
relatively insensitive to the initial conditions, Cohen [6] defines an inflation force on the active
contour. This force makes the model behave like an inflating balloon. The modified contour model
will be "11i... 1" by a strong edge and will simply pass through a spurious edge since the force
exerted by it is too weak relative to the ambient inflation force.
Although the inflation force prevents the curve from getting 1 .1.p1. '" by isolated spurious edges,
the active contour model cannot segment complex shapes like the one shown in figure l(b). More
over, despite a good initialization, the active contour model, due to its arclength and curvature
minimization properties, cannot be forced to extrude through any significant protrusions that a
shape may posses. One possible solution to this problem is to embed the snake model, which is
an instance of a 1D thinplatemembranespline, in an adaptive environment wherein the material
jj iP
(a) CT image (b) DSA image (c) Shapes with holes
Figure 1: Test bed for our topologyindependent shape modeling scheme.
parameters controlling the relative strengths of elasticity and rigidity are adapted to the data(see
[21]). The merits of such an approach are suspect since it is not always possible to derive criteria
upon which to base the adaptation algorithm. So the problem is one of accurately modeling '/'i.i
cations and protrusions in complex structures. In [27] it has been shown that multiple instances
of deformable models are required to handle shapes with several distinct parts. This can be very
cumbersome, for it involves excessive user interaction and presumes that the shape has already
been segmented into its constituent parts. Instead, we propose a dynamic shape modeling method
that will start with a single instance of the model and will automatically "sprout" branches during
the evolutionary process. Once the shape has been segmented from the image, its constituent part
structure may be inferred using the algorithm described in [12].
Most existing surface modeling techniques require that the topology of the object be known before
the shape recovery can commence. However, it is not always possible to specify the topology of
an object prior to its recovery. As a result, most shape recovery schemes make strong assumptions
about the object topology. Unknown topology is an important concern in object tracking and
motion detection applications where the positions of object boundaries are tracked in an image
sequence through time. During their evolution, these closed contours may change connectivity
and split, thereby undergoing a topological transformation. A heuristic criterion for splitting and
merging of curves in 2D which is based on monitoring deformation energies of points on the elastic
curve has been discussed in [20].
In the context of static problems, more recently, particle systems have been used to model sur
faces of arbitrary topology [25]. Smoothness and continuity constraints are imposed by subjecting
the particle system to interaction potentials which locally prefer planar or spherical arrangement.
Particles can be added and deleted dynamically to enlarge and trim the surface respectively, while
the system dynamics strive continually to organize the particles into smooth shapes. The result
is a versatile method with applications in surface fitting to sparse data and 3D medical image
segmentation.
The scheme described in this paper can be applied to static as well as dynamic situations, where
no a priori assumption about the object's topology is made. A single instance of our model, when
presented with an image having more than one shape of interest (see figure 1(c)), has the ability to
split freely to represent each shape [16, 17]. We show that by using our approach, it is also possible
to extract the bounding contours of shapes with holes in a seamless fashion (see figure 14).
1.1 Overview
In this subsection we briefly outline our scheme to model complex shapes. This method is inspired
by ideas first introduced in Osher & Sethian [19, 24] to model propagating fronts with curvature
dependent speeds. Two such examples are flame propagation and crystal growth, in which the
speed of the moving interface normal to itself depends on transport terms modified by the local
curvature. The challenge in these problems is to devise numerical schemes for the equations of the
propagating front which will accurately approximate these highly unstable physical phenomena.
Sethian [23] has shown that direct parameterization of the moving front may be unstable since it
relies on local properties of the solution. In contrast, a method which preserves the global properties
of the motion is sought. Osher and Sethian [19, 24] achieve this by viewing the propagating surface
as a specific level set of a higherdimensional function. The equation of motion for this function is
reminiscent of an initial valued HamiltonJacobi equation with a parabolic righthand side and is
closely related to a viscous hyperbolic conservation law.
In our work, we adopt these level set techniques to the problem of shape recovery. To isolate a
shape from its background, we first consider a closed, nonintersecting, initial hypersurface placed
inside (or outside) it. Following the above level set approach, this hypersurface is then made to
flow along its gradient field with a speed F(K), where K is the curvature of the hypersurface. As in
[19], we adopt a global approach and view the (N 1) dimensional moving surface as a level set of a
timedependent function y of N space dimensions. The equations of motion written for this higher
dimensional function are then amenable to stable I ,i,,i,.I. 'fying numerical schemes designed to
approximate hyperbolic conservation laws. Topological changes can be handled naturally in this
approach, since a particular level set {y = 0} of the function y need not be simply connected.
However, there are two problems that need to be surmounted before we can use this design for
shape recovery. First, it is required that we stop the hypersurface in the neighborhood of the
desired shape. We do this by synthesizing a ,, ,I,'., speed function from the image data. Secondly,
we have to construct an extension of this speed function to other level sets {y = C} in order to
give a consistent meaning to the imagebased speed term at all points in the image (see figure 2).
In the following sections we outline a possible solution to these problems.
We note that this work on interface motion and hyperbolic conservation laws as discussed in
[19, 23, 24], has been applied in the area of computer vision for shape characterization by Kimia
et al. [11, 12], who unify many diverse aspects of shape by defining a continuum of shapes (reac
tion/diffusion space), which places shapes within a neighborhood of other similar shapes. This leads
to a hierarchical description of a shape which is suitable for its recognition. The key distinguishing
feature of our work from that of Kimia et al., is that they assume the boundary of the object shape
to be known, while we recover it from noisy data. In other words, they show that by evolving a
known shape boundary, explicit clues can be derived towards the goal of developing a hierarchical
shape description. In contrast, we start with an arbitrary function y and recover complex shapes
by propagating it along its gradient field. Shape characterization, as in [12], can be achieved once
the object shape is extracted.
An important contribution of this paper over our previous work in [17] is the design and imple
mentation of a faster numerical technique for solving the governing equation. Our new algorithm
exploits the fact that the front, which is a particular level set {y = 0} of a higher dimensional
function, can be advanced by updating the function y at a small set of points. This scheme is an
alternative to updating the y function at all the points in the computational domain. The set of
points at which the update procedure is applied belong to a narrow band lying on either side of
the level set {i = 0}. Since the narrowband update strategy involves only a fraction of the total
number of points, a significant saving in time is realized, making our method a very attractive
alternative to other shape recovery schemes. A complete discussion of the narrowband techniques
for interface propagation may be found in [1].
In summary, we present a novel scheme for shape modeling which can be used in both computer
vision and computer graphics applications. Given the reconstructed shape, our approach can also
be extended to decipher the constituent part structure for highlevel processing. The remainder of
this paper is organized as follows: section 2 introduces the curvaturedependent front propagation
problem and establishes a link between HamiltonJacobi equations and a hyperbolic conservation
law. In section 3 we explain our level set algorithm for shape recovery and in section 4 we outline
the details of our narrowband approach. Finally in section 5 we present some experimental results
of applying our method to some noisy and low contrast medical images. We close with a discussion
of advantages of our approach and direction of future research in section 6.
2 Front propagation problem
In this section we present the level set technique due to Osher and Sethian [19]. For details and an
expository review, see Sethian [24]. As a starting point and motivation for the level set approach,
consider a closed curve moving in the plane, that is, let 7(0) be a smooth, closed initial curve in
Euclidean plane R2, and let 7(t) be the oneparameter family of curves generated by moving 7(0)
along its normal vector field with speed F(K), a given scalar function of the curvature K. Let
x(s, t) be the position vector which parameterizes 7(t) by s, 0 < s < S.
One numerical approach to this problem is to take the above Lagrangian description of the
problem, produce equations of motion for the position vector x(s, t), and then discretize the pa
rameterization with a set of discrete marker particles laying on the moving front. These discrete
markers are updated in time by approximating the spatial derivatives in the equations of motion,
and advancing their positions. However, there are several problems with this approach, as discussed
in Sethian [23]. First, small errors in the computed particle positions are tremendously amplified
by the curvature term, and calculations are prone to instability unless an extremely small time step
y z = (xy,t=O)
(a) (b)  A
y
C . X (t) Level 0
Set
X
Z (xyt)
Figure 2: Level set formulation of equations of motion (a) & (b) show the curve 7 and the surface
(x, y) at t = 0, and (c) & (d) show the curve 7 and the corresponding surface (x, y) at time t.
is employed. Second, in the absence of a smoothing curvature (viscous) term, singularities develop
in the propagating front, and an entropy condition must be observed to extract the correct weak
solution. Third, topological changes are difficult to manage as the evolving interface breaks and
merges. And fourth, significant bookkeeping problems occur in the extension of this technique to
three dimensions.
As an alternative, the central idea in the level set approach of Osher and Sethian [19] is to
represent the front 7(t) as the level set { = 0} of a function t. To motivate this approach, we
consider the example of an expanding circle. Suppose the initial front 7 at t = 0 is a circle in the
xyplane (figure 2(a)). We imagine that the circle is the level set {y = 0} of an initial surface
z = (x, y, t = 0) in R3 (see figure 2(b)). We can then match the oneparameter family of moving
curves 7(t) with a oneparameter family of moving surfaces in such a way that the level set { = 0}
always yields the moving front (see figures 2(c) & 2(d)).
In the general case, let 7(0) be a closed, nonintersecting, (N 1) dimensional hypersurface. Let
y(x, t), x E RN, be the scalar function such that (x, 0) = d(x), where d(x) is the signed distance
from x to the hypersurface 7(0). We use the plus sign if x is outside 7(0) and minus sign if x is
inside. Each level set of y flows along its gradient field with speed F(K). The gradient Vy(x, t)
is normal to the (N 1) dimensional level set passing through x. Now, we derive the equation of
motion fo f function y.
Consider the motion of some level set {i = C}. Following the derivation in [19], let x(t) be the
trajectory of a particle located on this level set, so
y(x(t), t) = C. (1)
The particle speed Ox/Ot in the direction n normal to 7(t) is given by the speed function F. Thus,
Ox
S.n = F, (2)
where the normal vector n is given by n = Vy/  Vi By the chain rule we get,
Ox
t+ V = 0 (3)
and substitution yields
t + F I Vq = 0, (4)
with an initial condition i(x, 0) = d(x). We refer to equation (4) as the level set "Hamilton
.1.. .. 1." formulation. Note that at any time, the moving front 7(t) is simply the level set {i(x, t) =
0}. There are several advantages to this approach. First, since the underlying coordinate system
is fixed, discrete mesh points used in the numerical update equations do not move, resulting in a
stable computation. Topological changes in the front can be handled naturally by exploiting the
property that the level surface {y = 0} need not be simply connected. y(x,t) always remains a
function, even if the level surface {i = 0} corresponding to the front 7(t) changes topology, or
forms sharp corners. The geometric and differential properties of 7(t) are captured in the function
i and can be readily extracted. As an example, if x E R2, the curvature is given by
(.' 2xyxy + .2)
K = (5)
(y + 2)3/2
This approach can also be easily extended to higher dimensions and appropriate expressions can
be obtained for the mean curvature and the Gaussian curvature [19].
By substituting F(K) = 1 EK as a typical speed function in equation (4), the equation of
motion becomes
t+ Vi = EK V (6)
Equation (6) resembles a HamiltonJacobi equation with viscosity, where .....iy" refers to the
secondorder parabolic righthand side. This equation can be solved using the stable, entropy
satisfying finite difference schemes, borrowed from the literature on hyperbolic conservation laws
(see [19]).
3 Shape recovery with front propagation
In this section, we describe how the level set formulation for the front propagation problem discussed
in the previous section can be used for shape recovery. There is a fundamental difference between the
problems of front propagation and shape recovery. In the former, the front represents a solid/liquid
interface (crystal growth) or a boundary separating burnt and unburnt regions (flame propagation).
In these cases the computation is alive as long as there remains a physical domain into which the
front can be moved. For example, the flame front can be moved as long as there is a region to be
burnt and it hasn't crossed the physical domain in which the solution is sought. On the contrary, in
shape recovery the front represents the boundary of an evolving shape. Since the idea is to extract
object shapes from a given image, the front should be forced to stop in the vicinity of the desired
object boundaries. This is analogous to the force criterion used to push the active contour model
towards desired shapes. We define the final shape to be the configuration when all the points on
the front come to a stop, thereby bringing the computation to an end.
Our goal now is to define a speed function from the image data that can be applied on the prop
agating front as a stopping criterion. In general the function F can be split into two components:
F = FA + FG. The term FA, referred to as the advection term, is independent of the moving front's
geometry. The front uniformly expands or contracts with speed FA depending on its sign and is
analogous to the inflation force defined in [6]. The second term FG, is the part which depends on
the geometry of the front, such as its local curvature. This (diffusion) term smooths out the high
curvature regions of the front and has the same regularizing effect on the front as the internal defor
mation energy term in thinplatemembrane splines [10] (see the figure (9)). We rewrite equation
(6) by splitting the influence of F as
S+ FA I V I +FG I V l = 0. (7)
First consider the case when the front moves with a constant speed, i.e. F = FA. To this if we
add a .,. l'.., speed term synthesized from the image, such that their sum tends to zero near large
image gradient locations, we will achieve our goal of bringing the front to a stop in the neighborhood
of object boundaries. To this end, we define a negative speed FI to be
FA
F(x, y) = (M A { VG, I(x, y) M2} (8)
where M1 and M2 are the maximum and minimum values of the magnitude of image gradient
 VG, I(x, y) 1, (x, y) E Q. The expression G, I denotes the image convolved with a Gaussian
smoothing filter whose characteristic width is a. Alternately, we could use a smoothed zerocrossing
image to synthesize the negative speed function. The zerocrossing image is produced by detecting
zerocrossings in the function V2G, I, which is the original image convolved with a Laplacian
ofGaussian filter whose characteristic width is a. The equation of motion with the addition of
imagebased speed is
S+ (FA + FI) V = 0. (9)
FI is called an extension of FI to points away from the boundary 7(t), i.e. at points (x, y) E
(0 7(t)), and is equal to FI on 7(t). We shall return to the issue of extension shortly. The value
of FI lies in the range [FA, 0] as the value of image gradient varies between M1 and M2. From
this argument it is clear that the front gradually attains zero speed as it gets closer to the object
boundaries and eventually comes to a stop.
In the case when the front moves with a speed that is a function of local curvature, i.e. FG f 0,
it is not possible to find an additive speed term from the image that will cause the net speed of
the front to approach zero in the neighborhood of a desired shape. Instead, we multiply the speed
function F = FA + FG with a quantity kI. The term kI, which is defined as
1
ki(x, y) = (10)
(, y) 1+ VG, I(x, y)' (10)
has values that are closer to zero in regions of high image gradient and values that are closer to
unity in regions with relatively constant intensity. More sophisticated stopping criteria can be
synthesized by using the orientation dependent steerablee" filters [9]. The modified equation of
motion is given by
t + ki(FA + FG) I V = 0. (11)
We now come to an important juncture in our discussion. The imagebased speed term, be it F1
or ki, has meaning only on the boundary 7(t), i.e. on the level set {y = 0}. This follows from
the fact that they were designed to force the propagating level set {f = 0} to a complete stop in
the neighborhood of an object boundary. However, the equation of motion (9) is written for the
function y, which is made up of infinitely many level curves. In other words, equations (9) & (11)
control the evolution of a family of level sets. Therefore, it is imperative that the net speed used in
the evolution equation has a consistent physical meaning for all the level sets, i.e. at every point
(x, y) E Q. Speed functions such as FG which are functions of geometric properties of the surface
z = i(x, y), can be readily computed at any (x, y) E Q. However, FI is not such a function. It
derives its meaning not from the geometry of y but from the configuration of the level set {y = 0}
in the image plane. Thus, our goal is to construct an imagebased speed function FI that is globally
defined. We call it an extension of FI off the level set {y = 0} because it extends the meaning of
FI to other level sets. Note that the level set {f = 0} lies in the image plane and therefore FI must
equal FI on {i = 0}. The same argument applies to the coefficient ki.
If the level curves are moving with a constant speed, i.e. FG = 0, then at any time t, a typical
level set {y = C}, C E R, is a distance C away from the level set {y = 0} (see figure 3). Observe
that the above statement is a rephrased version of Huygen's principle which, from a geometrical
standpoint, stipulates that the position of a front propagating with unit speed at a given time t
should consist of only the set of points located a distance t away from the initial front. On the
other hand, if FG # 0, the level sets will violate the property that they are a constant distance
away from each other. However, they will never collide and cross each other if the speed function
F = FA + FG is continuous (see [8]). With the above remarks in mind we state the following:
Property 1 Any external (imagebased) speed function that is used in the equation of motion
written for the function y should not cause the level sets to collide and cross each other during the
evolutionary process.
X
Figure 3: Huygen's principle construction
Recall that the function i(x,t) has been initialized to d(x), where d(x) is the signed distance
from a point x E Q to the boundary 7(0). Since we cannot attribute any geometric meaning to the
function F1 (ki) at points away from the level set {f = 0}, we look for a meaning that is consistent
with property (1). Therefore, the question to ask is: what is the value of FP (or ki) at a point (x, y)
lying on a level set {y = C}? We answer this question in the following construction (see figure 4).
Construction 1 The value of FI (ki) at a point P lying on a level set {f = C} is exactly the
value of F1 (ki) at a point Q, such that point Q is a distance C away from P and lies on the level
set { = 0}.
It is easy to see that FI reduces to F1 on {i = 0}. We use the same construction to determine
the value of ki at a point (x, y) lying on some level set {y = C}. Note that if the definition of
a speed function adheres to construction 1, then it will also be consistent with the property 1.
Thus, having ascribed the intent of pseudodifferential equations (9) & (11) in the context of shape
modeling, we can use finite difference schemes to solve them numerically. Since y can develop
corners and sharp gradients, numerical schemes borrowed from hyperbolic conservation laws are
used to produce stable upwind schemes. Moreover, the equations of motion can be solved on a
uniform mesh and the level sets can be moved without their explicit construction.
4 Numerical solution
In this section, almost without a change, we present the arguments discussed in Sethian [24]. For
complete details of the following scheme, we refer the reader to Osher & Sethian [19, 24].
I x
Figure 4: Extension of imagebased speed terms to other level sets
The equation (6) poses an initial valued problem. It is rewritten here as
Y + (y + 2)1/2 = V (12)
with y(x,y,t = 0) = + distance from (x,y) to 7(0). As shown in Sethian [23], for E > 0, the
parabolic righthand side diffuses sharp gradients and forces y to stay smooth at all values of t.
For E = 0, the boundary moves with unit speed, and a corner must develop from smooth initial
data. Once a corner develops, it is not clear how to propagate the front in the normal direction,
since the derivative is not defined at the corner. A variety of .i.." solutions which propagate
the curve beyond the occurrence of a singularity are possible. Of all such weak solutions, one is
interested in the one that is the limit of smooth solutions as E 0. This particular weak solution
can be selected with the help of a socalled i. iI ropy condition", see [23]: If the front is viewed as
a burning flame, then once a particle is burnt it stays burnt. Thus, approximations to the spatial
derivative are sought that do not artificially smooth sharp corners and which pick out the correct
entropy solution when singularities develop. The schemes given in [19, 24] are motivated by the fact
that the entropy condition for propagating fronts is identical to the one for hyperbolic conservation
laws, where stable, consistent, entropysatisfying algorithms have a rich history.
First, consider the onedimensional version of the level set equation, with constant normal velocity
FA = 1. Then the equation (6) becomes a standard HamiltonJacobi equation
t + H(,x) = 0, (13)
where H(P ) = ( )1/2, and with a given initial value of p. Let u = Ox. Taking the derivative
with respect to x, equation (13) becomes
ut + [H(u)] = 0, (14)
where H(u) = (u2)1/2. Equation (14) is a scalar hyperbolic conservation law in one space variable.
Solutions can develop discontinuous jumps, known as shocks, even with smooth initial data. In order
to make sense of the solution after shocks form, an integral version of the conservation law which
admits discontinuous solutions is studied. Both sides of equation (14) are integrated in an arbitrary
interval [a, b] to produce
d
d I u(x, t)dt = H[u(a, t)] H[u(b, t)]. (15)
at c1 a
u is known as a weak solution of the conservation law if it satisfies the above integral equation.
Note that u need not be differentiable to satisfy the integral form of the conservation law.
When will a numerical algorithm approximate the correct, entropysatisfying solution to equation
(15)? The answer is found in the following definition:
Definition 1 Let u7 be the value of u at a mesh point iAx at time nAt. A threepoint difference
scheme is said to be in conservation form if there exists a function g(ul, u2) such that the scheme
can be written in the form
Un+1 _ g(u7, u +1) g(U 1, u)
At z Ax (16)
where g(u,u) = H(u).
This definition is natural; the scheme must approximate the hyperbolic conservation law, subject
to the consistency requirement g(u, u) = H(u). In order to guarantee that the scheme picks out the
correct entropysatisfying weak solution, monotonicity is required, i.e., that u +1 be an increasing
function of the arguments 0u_ u and un 1. The main fact is: A conservative, monotone scheme
produces a solution that satisfies the entropy condition. Equation (16) is a scheme for the slope
u, which must be converted into a scheme for p itself. First write equation (13) with a forward
difference in time as,
+1 =. AtH(u). (17)
Since the numerical flux function g approximates H, the solution to equation (17) may be approx
imated by
+1 = tg(ui1/2, u i+1/2)
Atg(D; , Dt, ), (18)
where g is an appropriate numerical flux function and the standard definitions of the forward and
the backward difference operators have been used, namely,
1
D  1= ,X 1
Ax
D+t (19)
Ax
Finally, an appropriate numerical flux function g is required. In the special case where H(u)
may be written as a function of u2, i.e., H(u) = f(u2) for some function f, one can use the
HamiltonJacobi flux function given in [19]:
g(i1/2, i+1/2) = g(D0;, D+0i)
= f((max(D O, 0))2 + (min(D+O, 0))2). (20)
This conservative monotone scheme is upwind method in that it differences in the direction of
propagating characteristics. This is important, since it imposes boundary conditions on the walls
of a finitesized computational box. An upwind scheme automatically differences in the outward
flowing direction at the box walls if the boundary is expanding, thus information flows out. In the
case when FA = 1, so that f(u2) = (2)1/2, equation (18) reduces to
+1 _=. At{(max(D q0, 0))2 + (min(D+, 0))2}1/2. (21)
This algorithm produces the correct entropysatisfying weak solution to the propagating front
problem.
The above discussion can be extended to problems in more than one space dimension (see Osher
& Sethian [19]). Recall that the original intent was to solve equations (9) and (11). In two
dimensions, the scheme given in equation (21) is extended by differencing in each direction to
produce the following numerical scheme for equation (9):
+1 = A t[FA + (Fi);,j{(max(D,' ,0))2 + (min(D+.' ,,0))2
+(max(D i,0))2 + (min(D .' i,0))2}1/2. (22)
Similarly, the numerical scheme for equation (11) is,
+1 = AtFA(k~)i,j{(max(D.' i, 0))2 + (min(D .' i,0))2
+(max(D .' 0))2 + (min(D.' i, 0))2}1/2 AtFGk, I V I (23)
The second term FGk I Vi  is not approximated in the above equation; one may use a straight
forward central difference approximation to this term.
4.1 Narrowband updation scheme
In this section we first consider a straightforward approach to solve the equations (22) & (23) and
subsequently show that the time complexity of our algorithm can be improved by applying the
update procedure to a small set of points in the domain. In other words, we propagate the front
by solving the governing equation at a small number of points. These points lie in a narrow band
around the level set {y = 0}.
It should be observed that by updating the level set function on a grid, we are moving the level
sets without constructing them explicitly. A particular level set {i = 0}, which models the shape of
interest, is forced to a stop in the vicinity of desired object boundaries. This task is accomplished
by applying an imagebased speed term on the zero set (equation (11)). An extension function
is constructed to attribute a consistent meaning to the imagebased speed term at points away
from the zero set (see discussion in section 3). Therefore, the straightforward algorithm consists of
advancing from one time step to the next as follows:
Algorithm 1
1. At each grid point (iAx,jAy), where Ax and Ay are step sizes in either coordinate directions,
the extension of imagebased speed term is computed. This is done in accordance with the
construction described in previous section; i.e., by searching for a point q which lies on the
level set {i = 0}, and is closest to the point (iAx,jAy). The value of imagebased speed
term at the current point is simply its value at the point q.
2. With the value of extended speed term (k')i,j, and calculate .'+ using the upwind,
finite difference scheme given in equation (23).
3. Construct an approximation for the level set { = 0} from +1. This is required to visualize
the current position of the front in the image plane. A piecewise linear approximation for the
front 7(t) is constructed as follows. Given a cell C(i,j), if max(.' ;, i+,j,' ,+1, i+l,,+,j+1) <
0 or min(.' ;, ii+1,j, i+1, i+1,j+1) > 0, then C(i,j) 7(t) and is ignored, else, the entry
X
Figure 5: A narrowband of width 6 around the level set {y = 0}.
and exit points where y = 0 are found by linear interpolation. This provides two nodes on
7(t) and thus, one of the line segments which form the approximation to 7(t). The collection
of all such line segments constitutes the approximation to the level set { = 0}, which is used
for future evaluation of the imagebased speed term in the update equation (23).
4. Replace n by n + 1 and return to step 1.
It is easy to see from the above algorithm that the most expensive step is the computation of the
extension for imagebased speed term. This is because at each grid point, we must search for the
closest point lying on the level set {i = 0}. Moreover, if FG = 0, then the stability requirement for
the explicit method for solving our level set equation is At = O(Ax). For the full equation (23), the
stability requirement is At = O(Ax2). This could potentially force a very small time step for fine
grids. These two effects, individually and compounded, make the computation exceedingly slow.
A marginal improvement in performance can be achieved by evaluating the extension only once
every k iterations. In other words, we take k steps in time without recomputing the imagebased
speed (ki)ij. Alternately, we could downsample the image and perform our calculations at a lower
resolution. On the other hand, we run the risk of losing accuracy while computing on a coarse
grid. One approach is to embed our level set algorithm in a multiresolution framework. Instead,
we improve the time complexity of our algorithm by updating the level set function y only at a
small set of points on the grid.
The basic idea behind this new scheme stems from the observation that the front can be moved by
updating the level set function at a small set of points in the neighborhood of the zero set instead
of updating it at all the points on the grid. In figure (5) the bold curve depicts the level set { = 0}
and the shaded region around it is the narrowband. The narrowband is bounded on either side
by two curves which are a distance 6 apart, i.e., the two curves are the level sets {i = 6/2}. The
value of 6 determines the number of grid points that fall within the narrowband. Since, during a
given time step the value of .' is not updated at points lying outside the narrowband, the level
sets {y >1 6/2 } remain stationary. The zero set which lies inside moves until it collides with the
boundary of the narrowband. Which boundary the front collides with depends on whether it is
propagating inward or outward; either which way, it cannot move past the narrowband's boundary.
Therefore, as a consequence of our update strategy, the front can be moved through a maximum
distance of 6/2, either inward or outward, at which point we must rebuild an appropriate (a new)
narrow band. We reinitialize the y function by treating the current zero set configuration, i.e.,
{i = 0}, as the initial curve 7(0) (see section 2). Note that the reinitialization procedure must
account for the case when {y = 0} changes topology. This procedure will restore the meaning of
y function by correcting the inaccuracies introduced as a result of our update algorithm. Once a
new y function is defined on the grid, we can create a new narrowband around the zero set, and
go through another set of, say 1, iterations in time to move the front ahead by a distance equal to
6/2. The value of I is set to the number of time steps required to move the front by a distance
roughly equal to 6/2. This choice depends on some experimentation. Thus, a faster algorithm for
shape recovery consists of the following steps:
Algorithm 2
1. Set the iteration number m = 0 and go to step 2.
2. At each grid point (i,j) lying inside the narrowband, compute the extension kI of image
based speed term.
3. With the above value of extended speed term (ky );j, and calculate using the
upwind, finite difference scheme given in equation (23).
4. Construct a polygonal approximation for the level set {f = 0} from + A contour tracing
procedure is used to obtain a polygonal approximation. Given a cell (i,j) which contains
7(t), this procedure traces the contour by scanning the neighboring cells in order to find the
next cell which contains 7(t). Once such a cell is found, the process is repeated until the
contour closes on itself. The set of nodes visited during this tracing process constitutes the
polygonal approximation to 7(t). In general, to collect all the closed contours, the above
tracing procedure is started at a new, as yet unvisited, cell which contains the level set
{i = 0}. A polygonal approximation is required in step 2 for the evaluation of imagebased
speed term and more importantly, in step 6 for reinitializing the y function.
5. Increment m by one. If the value of m equals 1, go to step 6, else, go to step 2.
6. Compute the value of signed distance function y by treating the polygonal approximation
of {i = 0} as the initial contour 7(0). As mentioned earlier, a more general method of
reinitialization is required when {y = 0} changes topology. Go to step 1.
The issue of boundary conditions needs to be addressed at this time. In our original method, the
level set function y was initialized and updated on a square grid. To update .' at a boundary
point (i,j), it is necessary to specify the boundary conditions in order to numerically approximate
the derivatives. We have chosen freeend boundary conditions, i.e., if a derivative computation at
a given boundary point (i,j) involves a point that falls outside the computational domain, then
that derivative is set to zero. In our new approach, since we only update y at points lying in the
narrowband, the issue of specifying boundary conditions for points lying on the edge of the band
becomes pertinent. With our relatively simple speed motion, the freeend boundary condition is
adequate, however, in more complex applications such as crystal growth, and flame propagation,
accurate specification of boundary conditions is imperative.
We now show that this new faster approach provides a correct approximation to the propagating
front problem. In figure (6), we show the result of applying narrowband updation algorithm to a
star shaped front propagating with speed F = K, where K is the curvature as in equation (5). The
calculation was done on a unit box with 64 points in either direction, and a time step of At = 0.0003
was employed. The width of the narrow band has been set to 6 = 0.075, and the y function was
recomputed once every (1 =) 40 time steps. In figure 6(a), we show the initial curve along with the
level sets {y = i0.04}, where i E [1..5]. After 40 narrowband updations (figure 6(b)), only the
level sets {i <1 0.08 } move and the rest remain stationary. We note the inconsistency between
the level sets lying on either side of the narrowband, making the reinitialization step necessary
in order to restore the meaning of the y function. Following the reinitialization step, another 40
update steps are applied (figure 6(c)), which "dIT,. " the high curvature regions of the front even
further. In subsequent figures, the results of repeatedly applying the same strategy are shown.
Finally, in figure 6(f), the peaks and troughs on the front get completely diffused, and it attains a
smooth circular configuration after 4 reinitialization steps and a total of 200 time steps.
5 Experimental results
In this section we present several shape recovery results that were obtained by applying the narrow
band level set algorithm to image data. Given an image, our method requires the user to provide
an initial contour 7(0). As we shall see, there is absolutely no restriction on where the initial
contour can be placed in the image plane as long as it is inside a desired shape or encloses all the
constituent shapes. This feature is of paramount importance in the context of automatic shape
recovery. Our front seeks the object boundaries by either propagating outward in the normal
direction or propagating inward in the negative normal direction. This choice is made at the time
of initialization. Note that after the specification of initial shape of 7(0), our algorithm does not
require any further user interaction.
The initial value of the function y i.e., i(x, 0) is computed from 7(0). We first discretize the level
set function Q on the image plane and denote .' i as the value of Q at a grid point (iAx,jAy),
where Ax and Ay are step sizes in either coordinate directions. In our implementation, since we
usually work with 2kX2k images, the computational domain is a square with Ax = Ay = h. We
define the distance from a point (i,j) to the initial curve to be the shortest distance from (i,j) to
7(0). The magnitude of .' i is set to this value. We use the plus sign if (i,j) is outside 7(0) and
minus sign if (i,j) is inside. Once the value of .' i is computed at time t = 0 by following the above
procedure, we use the update equations from the previous section to move the front.
We now present our shape recovery results in 2D. First we consider a 256 X 256 CT (computed
tomography) image of an abdominal section shown in figure 7(a). Our intention is to recover the
shape of the stomach in this particular slice. The function y has been discretized on a 12_'8 X 12_'8
mesh, i.e., calculations are performed at every second pixel. In figure 8(a), we show the closed
(b) After 40 iterations
(c) After 80 iterations
(e) After 160 iterations
(d) After 120 iterations
(f) After 200 iterations
Figure 6: Narrowband updation algorithm applied to a starshaped front propagating with speed
F = K. Calculations were done on a 64 X 64 grid with a time step At = 0.0003. t was
recomputed after every 40 time iterations.
(a) Initialization
(a) Original image
(b) Imagebased speed term
Figure 7: Imagebased speed term ki(x y) = iii xy, with a = 3.25, synthesized from the
CT image.
contour that the user places inside the desired shape at time t = 0. The function y is then made to
propagate in the normal direction with speed F = ki(1.0 0.025K). We employed the narrow
band updation algorithm to move the front with a time step size set to At = 0.0005 and the y
function recomputed after every 50 time steps. Figure 7(b) shows the imagebased speed term
which is synthesized according to equation (10). This term when applied on the propagating front,
acts as a stopping criterion, thereby causing the front to come to rest in the neighborhood of the
desired shape. Note that in figure 7(b), ki(x, y) values lying in the interval [0..1] have been mapped
into the interval [0..255]. In figures 8(b) through 8(e) we depict the configuration of the level set
{y = 0} at four intermediate time instants. The final result is achieved after 575 time iterations
and is shown in figure 8(f). It should be noted that our method does not require that the initial
contour be placed close to the object boundary. In addition, observe how the front overshoots all the
isolated spurious edges present inside the shape (see figure 7(b)) and settles in the neighborhood of
edges which correspond to the true shape. This feature is a consequence of eK component in the
speed which diffuses regions of 1 '"ii curvature on the front and forces it to attain a smooth shape.
As mentioned in section 3, smoothness of the front can be controlled by choosing an appropriate
curvature component in the speed function F = 1 eK. The objective of our next experiment is
(b) After 100 iterations
(c) After 175 iterations
(e) After 450 iterations
(d) After 300 iterations
(f) After 575 iterations
Figure 8: Recovery of the stomach shape from a CT image of an abdominal section. Narrow
band computation was done on a 1'2' X 1'' grid the front was made to propagate with speed
F = ki(1.0 0.025K) and the time step At was set to 0.0005. y was recomputed once every 50
time steps.
(a) Initialization
(a) e = 0.05 (b) e = 0.25 (c) = 0.75
Figure 9: Smoothness control in shape recovery can be achieved by varying the curvature component
in the speed F = ki(1.0 EK).
to demonstrate smoothness control in the context of shape recovery. In figures 9(a) through 9(c),
we show the results of applying our narrowband shape recovery algorithm to an image consisting
of three synthetic shapes. Initialization was performed by drawing a curve enclosing each one of
the three shapes. We compute the signed distance function y(x, y) from these curves. The level
sets are then made to propagate with speed F = k(1.0 EK). First, as shown in figure 9(a), we
perform shape recovery with the value of E = 0.05. The process is repeated with different values
of E; 0.25 in figure 9(b) and 0.75 in figure 9(c). Clearly, with every increment in the value of E,
the level set {y = 0} attains a configuration that is relatively smoother. This is analogous to the
controlledcontinuity provided by thinplatemembrane splines [26, 10].
In our third experiment we recover the complicated structure of an arterial tree. The real image
has been obtained by clipping a portion of a digital subtraction angiogram. This is an example
of a shape with extended branches or significant protrusions. In this experiment we compare the
performance of our scheme with the active contour model to bring the limitations of the later into
focus. We first attempt to reconstruct the complex arterial structure using a snake model with
",/1./',o, forces [6]. In figures 10(a) through 10(i), we show a sequence of pictures depicting the
snake configuration in the image. We present the final equilibrium state of the snake in figures
10(c), 10(f), & 10(i) corresponding to three distinct initializations, one better than the preceding.
In all three cases the active contour model, even after 1000 time iterations, barely recovers the main
stem of the artery and completely fails to account for the branches. Two prominent limitations of
the snake model immediately come into focus. The first is the dependence of final result on the
initial configuration. This is due to the existence of multiple local minima in the (nonconvex) energy
functional which the numerical procedure explicitly minimizes. The second feature is the inability
of snake model to attain a stable shape with protrusions. Observe how in the third case, despite a
good initialization (figure 10(g)), the snake .pi,." back into a relatively II idpl. " configuration
in figure 10(h). This inadequacy stems from snake's arclength (elasticity) and curvature (rigidity)
minimizing property. Snake prefers regular shapes because shapes with protrusions have very high
deformation energies. Now, we apply our level set algorithm to reconstruct the same shape. After
initialization in figure 11(a), the front is made to propagate in the normal direction. We employ
the narrowband updation scheme with a band width of 6 = 0.075 to move the front. It can be
seen that in subsequent frames the front literally II. ." into the branches and finally in 11(f)
it completely reconstructs the complex tree structure. The advantages of our scheme are quite
apparent from this example. Since our front advancement process does not involve optimization of
any quantity, the shape recovery results we obtain are independent of initialization. In addition, a
 ",/. instance of our shape model pl.iii" branches and recovers all the connected components
of a given shape. All calculations were carried out on a 64 X 64 grid and a time step At = 0.001
was used.
In the next experiment, we depict a situation when the front undergoes a topological transforma
tion to reconstruct the constituent shapes in an image. The image shown in figure 12(a) consists
of three distinct shapes. Initial curve is placed in such a way that it envelops all the objects.
The front is then advanced in the direction of negative normal. Alternately, we could perform the
initialization by placing a curve in each one of the individual shapes and propagating them in the
normal direction. We choose the former option. The level set {i = 0} first wraps itself tightly
around the objects (see figures 12(c) & 12(d)) and subsequently splits into four separate closed
curves (figure 12(e)). While the first three closed segments of {y = 0} recover the three distinct
shapes, the one in the middle (see figure 12(e)), since it does not enclose any object, eventually
disappears. Figure 12(f) shows the final result. Again it should be noted that a single instance
of our shape model dynamically splits into three instances to represent each object. To show the
working of our algorithm, in figures 13(a)(i) we show the level sets { = iC}, i E [5.. + 5], with
C = 0.05. The function y was discretized on a 64 X 64 grid and At was set to 0.001.
Lastly, we show that our approach can also be used to recover shapes with holes. We do this by
applying our method to extract the shapes of handwritten characters. The characters "A" and "B"
in the figure (14) are examples of shapes with holes. Recovery of handwritten character shapes
finds application in the area of character recognition. These shapes can also serve as an input stage
to algorithms computing medial axis transform [13, 5] or skeletons [18, 15]. The attractive feature
here is the seamless way in which the outer and inner boundaries of a given shape can be recovered
without separate initializations. In figure 14(a), we show the initial contour which encloses all the
characters. This contour is then made to propagate inward with a constant speed. Figures 14(b)
shows an intermediate stage in the front's evolution and in 14(c), the front splits into four separate
contours. The calculation comes to a halt when in figure 14(d), the level set {y = 0} recovers the
outer boundaries of four separate characters. Unlike the characters "A" & "B" for which we need to
extract the inner boundary for their complete shape description, for the characters "S" & "H", the
outer boundary completely describes their shape. In the second stage of our computation, we treat
the zero set configuration in figure 14(d) as an initial state, and propagate all the fronts inward
by momentarily relaxing the imagebased speed term. This causes the zero set to move into the
character shapes as shown in figure 14(e). Finally in figure 14(f), the zero set comes to stop in the
neighborhood of holes that are present in the shapes of characters "A" & "B", thereby achieving a
complete shape recovery. The calculations for this experiment were done on a 1'2 X 1'2 grid and
the time step At was set to 0.00025.
6 Summary and discussion
In this paper we presented a novel shape modeling scheme. Our approach while retaining the
desirable features of existing methods for shape modeling, overcomes most of their deficiencies.
We adopt the level set techniques first introduced in Osher & Sethian [19] to the problem of
shape recovery. With this approach, complex shapes can be recovered from images. Unlike the
variational formulations for shape recovery which rely on energy minimization, the final result
in our method is completely independent of the initial state. This is a very desirable feature to
have, specially in applications such as automatic shape recovery from noisy images. Moreover,
our scheme makes no a priori assumption about the object's topology. Other salient features of
our shape modeling scheme include its ability to split and merge freely without any additional
bookkeeping during the evolutionary process, and its easy extendibility to higher dimensions. The
equations of motion governing our evolutionary system resemble an initial valued HamiltonJacobi
equation with a parabolic righthand side and are amenable to stable entropysatisfying numerical
solution schemes. Thus, the result is a very general shape modeling algorithm which we believe
will find numerous applications in the areas of computer vision and computer graphics.
In the context of 2D shape recovery, we force our front to come to a stop in the neighborhood
of object boundaries by synthesizing a speed function from noisy images. An extension is defined
to attribute a consistent meaning to this speed term at all points in the computational domain. A
significant improvement in time complexity is realized by adopting the narrowband update strategy.
In the narrowband scheme, the front is advanced by updating the y function in a narrow band
around the zero set. The efficacy and versatility of our method is demonstrated by a series of shape
recovery experiments on a set of synthetic and noisy medical images.
Our future research efforts will be focused on extending this method to other application domains.
It is easy to envision a similar scheme to reconstruct the surface structures of 3D objects from 3D
medical image data. However, it is not clear how to adapt the level set formulation in its present
form to reconstruct surface shapes from sparse range data. The issue of time complexity is an
important one, more so in 3D. We have seen in section 4 that the stability requirement for solving
the level set equations forces a very small time step for fine grids and that the evaluation of the
extension is computationally very expensive. The narrowband strategy has drastically reduced
the computational cost in 2D, but implementing a similar scheme in 3D could prove cumbersome.
One approach is to embed out level set algorithm in a multiresolution framework, where the loss of
accuracy at rapidlyconverging coarse grid computation is compensated by highly accurate slowly
converging fine grid calculation. Alternately, a parallel implementation can be devised. We are
currently pursuing such schemes.
References
[1] D. Adalsteinsson, "A narrowband approach to level set techniques for propagating interfaces,"
LBL Report, Oct. 1993.
[2] R. M. Bolle and B. C. Vemuri, "On threedimensional surface reconstruction methods," IEEE
Trans. on Pattern Analysis and Machine Intelligence, vol. PAMI 13, No. 1, pp. 113, 1991.
[3] T.E. Boult and J.R. Kender, "Visual surface reconstruction using sparse depth data," in Proc.
IEEE Conf. on Computer Vision and Pattern P ..*. I,'/'l,,. June 1',., pp. 6876.
[4] A. Blake and A. Zisserman, Visual Reconstruction, MIT Press, Cambridge, MA.
[5] H. Blum, "A transformation for extracting new descriptors of shape," Models for the Perception
of Speech and Visual Form (W. WathenDunn, ed.), Cambridge MA: MIT Press, 1',11".
[6] L. D. Cohen, "On Active Contour Models and Balloons," Computer Vision, Graphics, and
Imaige Pi ,,i Vol. 53, No. 2, pp. 211218, March 1991.
[7] H. Delingette, M. Hebert, and K. Ikeuchi, "Shape representation and image segmentation using
deformable models," in Proceedings of IEEE Conference on Computer Vision and Pattern
P .. ,, "./',,, pp. 467472, Maui Hawaii, June 1991.
[8] L. C. Evans and J. Spruck, :.lotion of level sets by mean curvature. I," Journal of Differential
Geometry, Vol. 33, pp. 635681, 1991.
[9] W. T. Freeman and E. H. Adelson, "Steerable filters for early vision, image analysis, and
wavelet decomposition," in Proceedings of ICCV, pp. lii. 115, Osaka, Japan, 1990.
[10] M. Kass, A. Witkin, and D. Terzopoulos, "Snakes: Active Contour Models," International
Journal of Computer Vision, pp. 321331, 1',",
[11] B. B. Kimia, A. R. Tannenbaum, and S. W. Zucker, "Toward a Computational Theory of
Shape: An Overview," in Proceedings of ECCV, Antibes, France, 1990.
[12] B. B. Kimia, A. R. Tannenbaum, and S. W. Zucker, "Shapes, shocks, and deformations I:
The components of shape and reactiondiffusion space," to appear in International Journal of
Computer Vision, 1992.
[13] D. T. Lee, .1. .1Ii.1 Axis Transformation of a Planar Shape," IEEE Trans. Patt. Anal. Machine
Intell. PAMI4, No. 4, pp. 363369, 1',_
[14] D. Lee and T. Pavlidis, "Onedimensional regularization with discontinuities," IEEE Trans.
on Pattern Analysis and Machine Intelligence, vol. PAMI 10, pp. 822829, 1'P1'
[15] F. Leymarie and M. D. Levine, "Simulating the grassfire transform using an active contour
model," IEEE Trans. Patt. Anal. Machine Intell. PAMI14 No. 1, pp. 5675, 1992.
[16] R. Malladi, J. A. Sethian, and B. C. Vemuri, "A topologyindependent shape modeling
scheme," in Proceedings of SPIE Conference on Geometric Methods in Computer Vision II,
San Diego, California, pp. 246258, July 1993.
[17] R. Malladi, J. A. Sethian, and B. C. Vemuri, "Shape modeling with front propagation: A
level set approach," Center for Pure and Applied Mathematics, Report PAM589, Univ. of
California, Berkeley, August 1993.
[18] N. Mayya and V. T. Rajan, "An efficient shape representation scheme using Voronoi skeletons,"
IBM Research Report, RC 19161 (83458), Sept. 1993.
[19] S. Osher and J. A. Sethian, I i..I . propagating with curvature dependent speed: Algorithms
based on HamiltonJacobi formulation," Journal of Computational Physics, Vol. 79, pp. 1249,
[20] R. Samadani, "Changes in connectivity in active contour models," Proceedings of the Workshop
on Visual Motion, pp. 337343, Irvine California, March 1','
[21] R. Samadani, "Adaptive snakes: control of damping and material parameters," Proceedings
of SPIE Conference on Geometric Methods in Computer Vision, Vol. 1570, pp. 202213, San
Diego California, July 1991.
[22] L.L. Schumaker, "Fitting Surfaces to Scattered data," in Approximation Theory II, G.G.
Lorentz, C.K. Chui, and L.L. Schumaker, (eds.). New York: Academic Press, 1976, pp. 203
267.
[23] J. A. Sethian, "Curvature and the evolution of fronts," Commun. in AMathematical Physics,
Vol. 101, pp. 487499, 1' "
[24] J. A. Sethian, \i.ii. I.... algorithms for propagating interfaces: HamiltonJacobi equations
and conservation laws," Journal of Differential Geometry, Vol. 31, pp. 131161, 1990.
[25] R. Szeliski and D. Tonnesen, "Surface modeling with oriented particle systems," Computer
Graphics SIGGRAPH, Vol. 26, No. 2, pp. 1~', 194, July 1992.
[26] D. Terzopoulos, "Regularization of inverse visual problems involving discontinuities," IEEE
Trans. on Pattern Analysis and Machine Intelligence, Vol. PAMI 8, No. 2, pp. 413424, 1','i.
[27] D. Terzopoulos, A. Witkin, and M. Kass, "Constraints on deformable models: Recovering 3D
shape and nonrigid motion," A, /'fi '.i Intelligence, 36, pp. 91123, 1''"
i] D. Terzopoulos, "The computation of visible surface representations,"IEEE Trans. on Pattern
Analysis and Machine Intelligence, vol. PAMI 4, Vol. 10, pp. 417 .; 1I'l'
[29] B. C. Vemuri, A. Mitiche, and J. K. Aggarwal, "Curvaturebased representation of objects
from range data," Int. Journal of Iimage and Vision Computing, 4, pp. 107114, 1','i.
[30] B. C. Vemuri and R. Malladi, "Surface griding with intrinsic parameters," Pattern P..' i "'I',
Letters, Vol. 13, No. 11, pp. ,1".812, November 1992.
[31] B. C. Vemuri and R. Malladi, "Constructing intrinsic parameters with active models for in
variant surface reconstruction," IEEE Trans. on Pattern Analysis and Machine Intelligence,
Vol. 15, No. 7, pp. 668681, July 1993.
[32] Y. F. Wang and J. F. Wang, "Surface reconstruction using deformable models with interior
and boundary constraints," in Proceedings of ICCV, pp. 300303, Osaka, Japan, 1990.
I *
S~
(a) Initialization 1
I'
SYP
S
(b) 500 iterations
(c) 1000 iterations
B'
'A
(d) Initialization 2
(e) 500 iterations
S
(g) Initialization 3
Ii
(h) 500 iterations
(f) 1000 iterations
SL
(i) 1000 iterations
Figure 10: An unsuccessful attempt to reconstruct a complex shape with igli. .1i" protrusions
using an active contour model. Three different (poor) results are shown in parts (c), (f), & (i)
corresponding to three distinct initializations in parts (a), (d), & (g) respectively.
I''
I'
D
'A
Sll
S
I'
Ii
AI
f 3'7 li^
(a) Initialization
(c) After 123 iterations
(e) After 275 iterations
(b) After 60 iterations
(d) After 200 iterations
(f) After 391 iterations
Figure 11: Reconstruction of a shape with igImi. .mi" protrusions: an arterial tree structure.
Computation was done on a 64 X 64 grid with a time step At = 0.001. A narrowband updation
scheme was used with a band width of 6 = 0.075.
Yrn
IrSI
11
~*''"Z~
a.
\sy
(a) Initialization
(b) After 30 iterations
I/
(c) After 60 iterations
Ut
(d) After 100 iterations
(e) After 125 iterations
(f) After 140 iterations
Figure 12: Topological split: a single instance of the shape model splits into three instances to
reconstruct the individual shapes. Computation was done on a 64 X 64 mesh with a time step
At = 0.001.
'I
7,
V
(a) Initialization
(b) 60 iterations
(c) 89 iterations
/
r/
(d) 100 iterations
(e) 120 iterations
(f) 130 iterations
/
a
(g) 154 iterations
(h) 170 iterations
(i) 185 iterations
Figure 13: Topological split example: level sets shown are { i
/
iCI, I E [5..+5], with C
0.05.
SB
AHI
(a) Initialization
5
A
(c) 760 iterations
s B
AH
(e) Second stage
(b) 545 iterations
S
A
(d) Outer boundaries
s B
AH
(f) Inner boundaries
Figure 14: Shapes with holes: a twostage scheme is used to arrive at a complete shape description
of both simple shapes and shapes with holes. Computation was performed on 12' X 12' grid and
the time step At was set to 0.00025.
