Shape Modeling with Front Propagation: A Level
Set Approach *
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
Developing shape models is an important aspect of computer vision research. Geo
metric and differential properties of the surface can be computed from shape models.
They also aid the tasks of object representation and recognition. In this paper we
present an innovative new approach for shape modeling which, while retaining im
portant features of the existing methods, overcomes most of their limitations. Our
technique can be applied to model arbitrarily complex shapes, shapes with protru
sions, and to situations where no a priori assumption about the object's topology can
be 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. Our
method is based on the level set ideas developed by Osher & Sethian to follow prop
agating 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. We move the interface by solving
a "Hamilton.T. ...1.1" type equation written for a function in which the interface is a
particular level set. A speed function synthesized from the image is used to stop the
interface in the vicinity of the object boundaries. The resulting equations of motion
are solved by numerical techniques borrowed from the technology of hyperbolic con
servation laws. An added advantage of this 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 synthesized images and noisy medical images.
*Submitted for publication in the IEEE Trans. on Pattern Analysis & Machine Intelligence.
1 Introduction
An important goal of computational vision is to recover the shapes of 2D and 3D objects
from various types of visual data. To achieve this goal, shape models that satisfy constraints
imposed by sensory data must be synthesized. Shape models aid the computation of cer
tain geometric and differential properties of surfaces. They also serve the purpose of an
intermediate stage in object recognition tasks, since they provide a more stable and useful
description than the original intensity images. In this paper we present a new approach to
shape modeling which, while retaining important features of the existing methods, overcomes
most of their limitations.
Shape reconstruction typically precedes the symbolic representation of surfaces. The shape
models must recover detailed structure from noisy data using only the weakest among the
possible assumptions about the observed shape. Several variational reconstruction methods
have been proposed and there is abundant literature on the same [2, 17, 23, 3, 24, 11].
Generalized spline models with continuity constraints are well suited for fulfilling the goals of
shape reconstruction (see [3, 4, 21]). Generalized splines are the key ingredient of the dynamic
shape modeling paradigm introduced by Terzopoulos et al., [22]. 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 [9, 25, 26, 27, 28, 6, 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 a strong
assumption on the object's 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 ...... l., lly" 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 during the shape reconstruction
process. Our method, which we shall describe presently, makes no assumption about the
object's topology, and it leads to a numerical algorithm whose convergence to the desired
shape is completely independent of the shape initialization.
The framework of energy minimization has also been used successfully in the problem
domain of extracting salient image contours edges and lines. Kass et al. [9] used energy
minimizing "snalk. that are attracted to the image features such as edges points and edge
segments, whereas internal spline forces impose a smoothness constraint. The weights of
the smoothness and image force terms in the energy functional can be adjusted for different
kinds of behavior. Snakes, also referred to as active contour models, are restricted examples
of the more general techniques of matching deformable models to image data by means
of energy minimization [22]. The scheme seeks to design energy functionals whose local
minima comprise the set of alternative solutions available to highlevel 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 [5] defines an inflation force on
the active contour. This new force makes the model behave like an inflating balloon. The
contour model with the above change will be stopped by a strong edge and will simply pass
through a spurious edge which is too weak relative to the ambient inflation force.
Although the inflation force prevents the curve from getting "trapped" by isolated spurious
edges, the active contour model cannot segment complex shapes with significant protrusions
like the one shown in figure (1). Moreover, 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 thinplate
membranespline, in an adaptive environment wherein the material parameters controlling
Figure 1: Digital subtraction angiogram of an arterial structure
the relative strengths of elasticity and rigidity are adapted (see [16]). 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 bifurcations
and protrusions in complex structures. In [22] 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 deciphered into its constituent parts. Instead, we propose a method that will
start with a single instance of the model and will sprout the branches during the evolutionary
process. Once the shape has been segmented from the image, its constituent part structure
can be inferred [10].
Most existing surface modeling techniques require that the topology of the object be known
before the reconstruction can commence. However, it is not always possible to specify the
topology of an object prior to its reconstruction. As a result, most reconstruction schemes
end up making strong assumptions about object topology. Vision systems which derive
quantitative models of complex object shapes by integrating different visual modalities can
not evade the issue of unknown and unpredictable topologies. Unknown topology is also 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 which is based on monitoring
deformation energies of points on the elastic curve has been discussed in [15]. More recently,
molecular dynamics has been used to model surfaces of arbitrary topology [20]. Smoothness
and continuity constraints are imposed by subjecting a particle system to interaction poten
tials 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 situations where no a priori as
sumption about the object's topology can be 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 [12].
1.1 Overview
In this subsection we briefly outline the scheme we use to model complex shapes. Our
method is inspired by ideas first introduced in Osher & Sethian [13, 19] to follow propagating
fronts with curvaturedependent 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 [18] 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 [13, 19] achieve this by embedding the surface in a higher
dimensional 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
problems of shape reconstruction. To isolate a shape from its background, we first consider a
closed, nonintersecting, initial hypersurface placed inside it. Following the level set approach
above, this hypersurface is then made to flow along its gradient field with some speed F(K),
where K is the curvature of the hypersurface. As in [13], 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 ent, ,. il'. '.If, ;.'.,, numerical schemes designed to approximate
hyperbolic conservation laws. Topological changes can be handled naturally in this approach,
since a particular level set {0 = 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
reconstruction. First, it is required that we stop the hypersurface in the neighborhood
of the desired shape. We do this by synthesizing a negative speed function from the image.
Secondly, we have to construct an extension of this speed function to other level sets {0 = C}
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 [13, 18, 19], has been applied in the area of computer vision for shape characterization
by Kimia et al. [10], who unify many diverse aspects of shape by defining a continuum
of shapes (reaction/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 object shape to be known, while we reconstruct 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 can be readily done once the object shape is extracted.
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 used for deciphering the constituent part structure. The remainder of
this paper is organized as follows: section 2 introduces the curvaturedependent front prop
agation problem and establishes a link between HamiltonJacobi equations and a hyperbolic
conservation law. In section 3 we explain our level set algorithm for shape reconstruction
and section 4 presents some experimental results of applying our method to some synthe
sized and real noisy images. We close with a discussion of advantages of our approach and
direction of future research in section 5.
2 Front propagation problem
In this section we present the level set technique due to Osher and Sethian [13]. For details
and an expository review, see Sethian [19]. 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
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 parameterization 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 one time step. However, there are several
problems with this approach, as discussed in Sethian [18]. 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 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
(a) (b)  
Y 1X
() (d) :::::::: ::::::::
0X  7 (t) Level Y=o
(x, y) at time t.
dimensions.
As an alternative, the central idea in the level set approach of Osher and Sethian [13] is to
represent the front 7(t) as the level set { = 0} of a function 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 i ./plane (figure 2(a)). We imagine that the circle is the level set {i = 0} of an initial
surface z = 0(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 {1 = 0} always yields the moving front (figures 2(c) & 2(d)).
In the general case, let 7(0) be a closed, nonintersecting, (N 1) dimensional hypersurface.
Let O(x, t), x E R', be the scalar function such that O(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 1 flows along its gradient field with speed
F(K). The gradient VO(x, t) is normal to the (N 1) dimensional level set passing through
x. Now, we derive the equation of motion for function (.
Consider the motion of some level set { = C}. Following the derivation in [13], let x()
be the trajectory of a particle located on this level set, so
0(x(t), t)= C. (1)
The particle speed ax/at in the direction n normal to 7(t) is given by the speed function F.
Thus,
ax
n = F, (2)
where the normal vector n is given by n = Vy/ I VO . By the chain rule we get,
Ox
a+ V = 0 (3)
and substitution yields
+ F  V 1= 0, (4)
with an initial condition O(x, 0) = d(x). We refer to equation (4) as the level set "Hamilton
Jacobi" formulation. Note that at any time, the moving front 7(t) is simply the level set
{y(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 {o = 0} need not be
simply connected. O(x,t) always remains a function, even if the level surface {o = 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 y and can be readily extracted.
As an example, if x e R2, the curvature is given by
( '? Q / '/ l ) 1 >
K = 20 x ;"1(5)
(;?2 + ,'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 [13].
By substituting F(K) = 1 eK as a typical speed function in equation (4), the equation
of motion becomes
+ VO I= eK VI V (6)
Equation (6) resembles a HamiltonJacobi equation with viscosity, where  ..... ty" refers
to the secondorder parabolic righthand side. This equation can be solved using the sta
ble, entropysatisfying finite difference schemes, borrowed from the literature on hyperbolic
conservation laws (see [13]).
3 Shape reconstruction 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 reconstruction. There is a fundamen
tal difference between the problems of front propagation and shape reconstruction. In the
former, the front represents a solid/liquid interface (crystal growth) or a boundary separat
ing 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 reconstruction
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 propagating 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 [5]. 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 deformation energy term in
thinplatemembrane splines [9]. We rewrite equation (6) by splitting the influence of F as
+ FA V I+FG + V0 = 0. (7)
First consider the case when the front moves with a constant speed, i.e. F = FA. To this if
we add a negative 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)= { VG, I(x,y) M2}, (8)
(Ml M2) '
where M1 and M2 are the maximum and minimum values of the magnitude of image gradient
VG, I(x, y) , (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 LaplacianofGaussian filter whose characteristic width is a. The
equation of motion with the addition of imagebased speed is
.+ (FA+FI) VO = 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 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 kz. The term ki, which is defined
as
1
kj(x,y) = (10)
1+  VG, *I(x,y) '
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. The modified equation of motion is
given by
+ ki(FA + FG) V = 0. (11)
We now come to an important juncture in our discussion. The imagebased speed term,
be it FI or ki, has meaning only on the boundary 7(t), i.e. on the level set {0 = 0}. This
follows from the fact that they were designed to force the propagating level set {y = 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) G 0. Speed functions such as FG which
are functions of geometric properties of the surface z = ~(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 {0 = 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 {q = 0} lies in the image plane and therefore FI must
equal FI on {0 = 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 Hi,.i, ,,'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' 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:
w=C
x
Figure 3: Huygen's principle construction
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.
Recall that the function y(x,t) has been initialized to d(x), where d(x) is the signed
distance from a point x EG to the boundary 7(0). Since we cannot attribute any geometric
meaning to the function FI (ki) at points away from the level set {y = 0}, we look for a
meaning that is consistent with property (1). Therefore, the question to ask is: what is the
value of F, (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 1I,.,,, on a level set { C = C} is i .i. IIl,
the value of FI (ki) at a point Q, such that point Q is a distance C ,., l,,, from P and lies
on the level set {O = 0}.
It is easy to see that FI reduces to FI on {f = 0}. We use the same construction to
determine the value of kI at a point (x, y) lying on some level set {i = 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 and experimental results
In this section, almost without a change, we present the arguments discussed in Sethian [19].
For complete details of the following scheme, we refer the reader to Osher & Sethian [13, 19].
The equation (6) poses an initial valued problem. It is rewritten here as
Y+ y)1/2 = v.( (12)
with ((x,y,t = 0) = + distance from (x,y) to 7(0). As shown in Sethian [18], 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 so
called "entropy condition", see [18]: If the front is viewed as a burning flame, then once a
particle is burnt it sui i 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 [13, 19] 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(,) = 0, (13)
=C P _^M)
VCQ (xy)
5 X
Figure 4: Extension of imagebased speed terms to other level sets
where H(9,) = _(0f)1/2, and with a given initial value of p. Let u = p. Taking the
derivative with respect to x, equation (13) becomes
ua + [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 u(x, t)dt = H[u(a, t)] H[u(b, t)]. (15)
dt 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(ui, U2) such
that the scheme can be written in the form
+ ,+i)1ii) (16)i
At 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 ul, u', and u"i. 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 0 itself.
First write equation (13) with a forward difference in time as,
+1 = A AtH(u). (17)
Since the numerical flux function g approximates H, the solution to equation (17) may be
approximated by
+l  Atg(ui1/2,Ui+1/2)
= Atg(D Oi, D+ ), (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
Dj, (19)
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 [13]:
g(ui1/2,i+1/2) = g(D;i, Df+i)
=f((max(D 0))2 + (min(D+ 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
outwardflowing direction at the box walls if the boundary is expanding, thus information
flows out. In the case when FA = 1, so that f(U2) = (u2)1/2, equation (18) reduces to
n+1 = At{(max(D;4, 0))2 + (min(D+4, 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 [13]). 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):
= At[FA + ()i,]{(max(D; ,, ,0))2 + (min(D. ,0))2
+(max(D ,0))2 + (min(DS. ,, 0))2}1/2. (22)
Similarly, the numerical scheme for equation (11) is,
i = AtFA(kI)i,j{(max(D;.' 0))2 + (min(D +' 0))2
+(max(D .' ;,0)j+)2 + min(D} ,))21/2 AtFGki  V (23)
The second term FGki I Vy I is not approximated in the above equation; one may use a
straightforward central difference approximation to this term.
4.1 Experimental results
In this section we present several shape reconstruction results that were obtained by applying
the level set algorithm to image data. We also outline some of the implementation details of
our algorithm. 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., y(x, 0) is computed from 7(0). We first discretize
the level set function y on the image plane and denote ; ; as the value of y 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 one 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 ;' 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.
It should be observed that by updating the level set function on a grid, we are moving
the level sets without constructing them explicitly. To find the position of the front and
to compute the imagebased speed terms, at each time step we need to find the level set
{1 = 0}. We construct a piecewise linear approximation for 7(t) as follows. Given a cell
(i,j), ifmax(.' i, tj, i ;i+, +,' ,j+i) < 0 or min(.' ;, I' ;+i, ;,' +l,j+1) > 0, then that
cell cannot contain 7(t), so we ignore that cell. Otherwise, we find the entrance and exit
points where y = 0 by linear interpolation. This provides two nodes on 7(t) and thus one
of the line segments which form our approximation to 7(t). The collection of all such line
segments constitutes our approximation to the level set {1 = 0}, which is used for future
evaluation of the imagebased speed term in the update equation (23). These line segments
are also used to display the current position of the front in the image plane.
The stability requirement for the explicit method for solving our level set equation is
At = O(Ax2) for the equation (23). If FG = 0, then the stability requirement is At = O(Ax).
This could potentially force a very small time step for fine grids making the computation
excruciatingly slow. Since the algorithm to extend the imagebased speed terms to other
level sets is relatively more expensive than the front propagation algorithm, we could improve
the performance by evaluating the extension only once every k iterations. In other words,
we take k steps in time without recomputing the force field (ki)i,j. Alternately, we could
downsample the image and perform our calculations at a lower resolution. In this approach
we run the risk of losing accuracy. However, we show that the results obtained with down
sampled images are very promising. We can get the best of both worlds by embedding our
level set algorithm in a multiresolution framework. We refrain from continuing this line of
discussion any further since the details of such schemes are beyond the scope of this paper.
We now present our shape reconstruction results in 2D. We first consider a 256 X 256 image
with a single shape. The function y has been discretized on a 64 X 64 mesh, i.e. calculations
are performed at every fourth pixel. In figure 5(a), we show the closed contour that the
user places around the shape at time t = 0. The function y is then made to propagate
in the negative normal direction. Figures 5(b) & 5(c) depict the configuration of the level
set {y = 0} at two intermediate time instants. The final result is achieved after 80 time
iterations and is shown in figure 5(d). It should be noted that our method does not require
that the initial contour be placed close to the object boundary. In order to quantize the loss
in accuracy which results at lower resolutions, we perform the same reconstruction on a 256
X 256 mesh. The result is shown in 5(e). In figure 5(e), we also plot the level sets {y = C}
to verify that the imagebased speed term does not violate the property 1.
In our second 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 branches and significant protrusions. In this experiment we compare
the performance of our scheme with the active contour model and bring its limitations into
focus. We first attempt to reconstruct the complex arterial structure using a snake model
with inflation forces. In figures 6(a) through 6(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 6(c), 6(f), & 6(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 light. 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 6(g)), the
snake n.ip" back into a relatively I ,i in pl " configuration in figure 6(h). This inadequacy
stems from snake's arclength (elasticity) and curvature (rigidity) minimizing nature. 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 7(a), the front is made to propagate in the normal direction. It can be seen that in
subsequent frames the front literally "fl .. " into the branches and finally in 7(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 reconstruction results we obtain are independent of initialization.
In addition, a single instance of our shape model p .ll,,I ;" branches and recovers all the
connected components of a given shape. In figures 8(a)(i) we plot the other level sets to
elucidate the process of extending the imagebased speed function to points away from the
zero set. All calculations were carried out on a 64 X 64 grid and the time step At is set to
0.001.
Lastly, in figure (9) we depict a situation when the front undergoes a topological trans
formation to reconstruct the constituent shapes in an image. The image 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. The level set {0 = 0} first wraps
itself tightly around the objects (see figures 9(c) & 9(d)) and subsequently splits into four
separate closed curves (figure 9(e)). While the first three closed segments of {o = 0} recover
the three distinct shapes, the one in the middle (see figure 9(e)), since it does not enclose
any object, eventually disappears. Figure 9(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. Again, to show the working of our algorithm, in figures 10(a)(i) we
show the level sets {y = iC}, i G [5.. +5], with C = 0.05. The function y was discretized
on a 64 X 64 grid and At is set to 0.001.
5 Concluding remarks
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 de
ficiencies. We adopt the level set techniques first introduced in Osher & Sethian [13] to
the problem of shape recovery. With this approach, complex shapes can be reconstructed.
Unlike the variational formulations for shape reconstruction which rely on energy minimiza
tion, the final result in our method is completely independent of the initial state. This is a
very desirable feature to have, specially if the problem is to recover object shapes 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 ex
tendibility 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 reconstruction, we force our front to come to a stop in the
neighborhood of object boundaries by synthesizing a negative speed term from noisy images.
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 other issue is
one of time complexity. We have seen in the previous section that the stability requirement
for solving the level set equations forces a very small time step for fine grids. On the other
hand, shape reconstruction results on a coarse grid suffer from loss of accuracy. To salvage
this situation we propose 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 slowlyconverging fine grid calculation. We are currently working on one
such scheme.
References
[1] R. Bajcsy and S. Kovacic, "Multiresolution elastic matching," Computer Vision, Graph
ics, and Image Processing, Vol. 46, pp. 121, 1989.
[2] R. M. Bolle and B. C. Vemuri, "On threedimensional surface reconstruction methods,"
IEEE Trans. on Pattern A,,.lii. 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 Recognition, June 1'I i., pp. 6876.
[4] A. Blake and A. Zisserman, Visual Reconstruction, MIT Press, Cambridge, MA.
[5] L. D. Cohen, "On Active Contour Models and Balloons," Computer Vision, Graphics,
and Image Processing, Vol. 53, No. 2, pp. 211218, March 1991.
[6] L. D. Cohen and I. Cohen, "Deformable models for 3D medical images using finite
elements and balloons," in Proceedings of IEEE Conference on Computer Vision and
Pattern Recognition, pp. Urbana Illinois, June 1992.
[7] H. Delingette, M. Hebert, and K. Ikeuchi, "Shape representation and image segmenta
tion using deformable models," in Proceedings of IEEE Conference on Computer Vision
and Pattern Recognition, pp. 467472, Maul Hawaii, June 1991.
[8] L. C. Evans and J. Spruck, "Motion of level sets by mean curvature. I," Journal of
Differential Geometry, Vol. 33, pp. 635681, 1991.
[9] M. Kass, A. Witkin, and D. Terzopoulos, "Snakes: Active Contour Models," Interna
tional Journal of Computer Vision, pp. 321331, 1'l"'
[10] B. B. Kimia, A. R. Tannenbaum, and S. W. Zucker, "Shapes, shocks, and deformations
I: The components of shape and reactiondiffusion space," Technical Report LEMS105,
Division of Engineering, Brown University, June 1992.
[11] D. Lee and T. Pavlidis, "Onedimensional regularization with discontinuities," IEEE
Trans. on Pattern An,.,lii.' and Machine Intelligence, vol. PAMI 10, pp. 822829, 1'i1i.
[12] 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, July 1993.
[13] S. Osher and J. A. Sethian, "Fronts propagating with curvature dependent speed: Al
gorithms based on HamiltonJacobi formulation," Journal of Computational P,,.:i ;
Vol. 79, pp. 1249, 1l'"
[14] W. Press, B. Flannery, S. Teukolsky, and W. Vetterling, Numerical Recipes in C, Cam
bridge University Press, Cambridge, 1'l""
[15] R. Samadani, "Changes in connectivity in active contour models," Proceedings of the
Workshop on Visual Motion, pp. 337343, Irvine California, March 1989.
[16] R. Samadani, "Adaptive snakes: control of damping and material parameters," Pro
ceedings of SPIE Conference on Geometric Methods in Computer Vision, Vol. 1570, pp.
202213, San Diego California, July 1991.
[17] L.L. Schumaker, "Fitting Surfaces to Scattered data," in Approximation Ti .. ry II, G.G.
Lorentz, C.K. Chui, and L.L. Schumaker, (eds.). New York: Academic Press, 1976, pp.
203267.
[18] J. A. Sethian, "Curvature and the evolution of fronts," Commun. in Mathematical
Pt,;, a Vol. 101, pp. 487499, li'
[19] J. A. Sethian, "Numerical algorithms for propagating interfaces: HamiltonJacobi equa
tions and conservation laws," Journal of Differential Geometry, Vol. 31, pp. 131161,
1990.
[20] R. Szeliski and D. Tonnesen, "Surface modeling with oriented particle systems," Com
puter Graphics SIGGRAPH, Vol. 26, No. 2, pp. 1Si194, July 1992.
[21] D. Terzopoulos, "Regularization of inverse visual problems involving discontinuities,"
IEEE Trans. on Pattern A,,,lii.,.' and Machine Intelligence, Vol. PAMI 8, No. 2, pp.
413424, 1'li.
[22] D. Terzopoulos, A. Witkin, and M. Kass, "Constraints on deformable models: Recov
ering 3D shape and nonrigid motion," Artificial Intelligence, 36, pp. 91123, l'l",
[23] D. Terzopoulos, "The computation of visible surface representations,"IEEE Trans. on
Pattern Al,,,1.', and Machine Intelligence, vol. PAMI 4, Vol. 10, pp. 417438, '1"'
[24] B. C. Vemuri, A. Mitiche, and J. K. Aggarwal, "Curvaturebased representation of
objects from range data," Int. Journal of Image and Vision Computing, 4, pp. 107114,
1' 1' ,
[25] B. C. Vemuri and R. Malladi, "Deformable models: Canonical parameters for surface
representation and multiple view integration," Proc. IEEE Conference on Computer
Vision and Pattern Recognition, pp. 724725, Maul Hawaii, June 1991.
[26] B. C. Vemuri and R. Malladi, "Surface griding with intrinsic parameters," Pattern
Recognition Letters, Vol. 13, No. 11, pp. 11'812, November 1992.
[27] B. C. Vemuri and R. Malladi, "Constructing intrinsic parameters with active models
for invariant surface reconstruction," IEEE Trans. on Pattern A,#,,il,.'. and Machine
Intelligence, in press.
l"] 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.
(b) After 22 iterations
(c) After 49 iterations
(d) After 80 iterations
(e) Fine grid computation
Figure 5: Shape reconstruction results at different resolutions: mesh size for parts (a), (b),
(c), & (d) is 64 X 64 and the time step At = 0.008. Computation in part (e) was performed on
a 256 X 256 mesh with a step size At = 0.00005. The other level sets shown are {y = 0.05}.
rat
(a) Initialization
I'
(a) Initialization 1
(d) Initialization 2
I'1
(b) 500 iterations
)
(e) 500 iterations
I')
mIll
(c) 1000 iterations
m
IA
(f) 1000 iterations
I!'
(g) Initialization 3
(h) 500 iterations
* .1
(i) 1000 iterations
Figure 6: An unsuccessful attempt to reconstruct a complex shape with ,*i.,licant" pro
trusions 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.
PII~
t
PI~Psj
muM
(a) Initialization
(c) After 123 iterations
(e) After '' iterations
(b) After 60 iterations
(d) After 200 iterations
(f) After 391 iterations
Figure 7: Reconstruction of a shape with ili lcant" protrusions: an arterial tree structure.
Computation was done on a 64 X 64 grid with a time step At = 0.001.
il
; W.
(a) Initialization (b) 60 iterations
(d) 170 iterations
(g) 320 iterations
(e) 200 iterations
(h) 350 iterations
(f) '' iterations
(i) 390 iterations
Figure 8: Reconstruction of a shape with protrusions: figure shows the level sets {y = iC},
i E [5.. + 5], and C = 0.05.
(c) 123 iterations
(a) Initialization
(c) After 60 iterations
(b) After 30 iterations
,t
(d) After 100 iterations
*P
(e) After 125 iterations
(f) After 140 iterations
Figure 9: 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'
(a) Initialization
(b) 60 iterations
(c) 89 iterations
/
r
(d) 100 iterations
(e) 120 iterations
(f) 130 iterations
0
(g) 154 iterations
Figure 10: Topological split example:
C = 0.05.
(h) 170 iterations
level sets shown are {y
(i) 185 iterations
iC}, i G [5.. + 5], with
/b
