Citation
A topology-independent shape modeling scheme

Material Information

Title:
A topology-independent shape modeling scheme
Creator:
Malladi, Ravikanth, 1966
Publication Date:
Language:
English
Physical Description:
viii, 68 leaves : ill. ; 29 cm.

Subjects

Subjects / Keywords:
Approximation ( jstor )
Computer vision ( jstor )
Curvature ( jstor )
Equations of motion ( jstor )
Geometric shapes ( jstor )
Modeling ( jstor )
Snakes ( jstor )
Surface contours ( jstor )
Three dimensional modeling ( jstor )
Topology ( jstor )
Computer and Information Sciences thesis Ph. D
Computer vision -- Mathematical models ( lcsh )
Dissertations, Academic -- Computer and Information Sciences -- UF
Image processing -- Mathematical models ( lcsh )
Genre:
bibliography ( marcgt )
non-fiction ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1993.
Bibliography:
Includes bibliographical references (leaves 64-67).
General Note:
Typescript.
General Note:
Vita.
Statement of Responsibility:
by Ravikanth Malladi.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. §107) for non-profit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact the RDS coordinator (ufdissertations@uflib.ufl.edu) with any additional information they can provide.
Resource Identifier:
030435953 ( ALEPH )
31255145 ( OCLC )

Downloads

This item has the following downloads:


Full Text







A TOPOLOGY-INDEPENDENT SHAPE MODELING SCHEME


By

RAVIKANTH MALLADI













A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA


1993
























Copyright 1993 by

Ravikanth Malladi













To my parents, Pardhasaradhi and Vasanthavalli Malladi, my sister Lakshmi, and my uncle, Chandrasekhar Nori, who have always encouraged my academic aspirations.













ACKNOWLEDGEMENTS


First, I would like to express my gratitude to my advisor, Dr. Baba C. Vemuri, for his continuous encouragement, for providing financial support, for numerous discussions on research topics, and guidance throughout the course of this work. I am grateful to Dr. James A. Sethian for consenting to serve on my supervisory committee; it was his involvement and insightful suggestions that made progress possible. I thank him for being generous with his time and providing me with more ideas than I could implement.
I would also like to thank the other members of my committee, Drs. Gerhard X. Ritter, Andrew F. Laine, Paul A. Fishwick, and Murali Rao, for their help and constructive criticism. Special thanks are due to Randy Fischer for his help in acquainting me with the video equipment and the KSR machine.

I owe a great deal to my friend Niranjan Mayya for his constant support and understanding throughout the past five years. An early acquaintance with Mark Burch has grown into a strong friendship which I shall always cherish. I thank my other "old" friends Amitabh Mungae, Pamela Rolfe, and Vidhya Gangar for making my life in Gainesville a very enjoyable one. More recently, in Roxana Dastur, Eman Anwar, Pankaj Agarwal, and Chitrank Sheshadri, I found excellent friends with whom I have enjoyed numerous lively conversations.

Financial support for this work was provided in part by the Department of Radiology and the Department of Computer and Information Sciences at the University of Florida.
















TABLE OF CONTENTS


ACKNOWLEDGEMENTS ............................ iv

ABSTRACT ........................................... vii

CHAPTERS

1 INTRODUCTION ............................... 1

1.1 O verview . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Dissertation Outline ........................... 7

2 ENERGY METHODS FOR SHAPE MODELING .............. 9

2.1 Symmetry-Seeking Deformable Models ................. 9
2.1.1 The deformable model ...................... 10
2.1.2 Extrinsic forces .......................... 13
2.2 Snake: The Active Contour Model .................... 14
2.2.1 Snake dynamics .......................... 15

3 FRONT PROPAGATION PROBLEM .................... 17

3.1 Lagrangian Formulation ......................... 17
3.2 Level Set Formulation .......................... 18

4 SHAPE RECOVERY WITH FRONT PROPAGATION ........... 22

4.1 Stopping Criterion: An Image-based Speed Term ........... 22
4.2 Extension of Image-based Speed Term ................. 24
4.3 Shape Recovery in 3D .......................... 27

5 NUMERICAL SOLUTION .......................... 29

5.1 Entropy-satisfying Upwind Numerical Scheme ............. 30
5.2 Narrow-band Updation Scheme ..................... 33
5.3 Numerical Issues in 3D .......................... 40

6 EXPERIMENTAL RESULTS ......................... 43

6.1 2D Shape Recovery Results ...... ....................... 44
6.2 Experimental Results in 3D ...... ....................... 50

V











7 CONCLUDING REMARKS .......................... 62

REFERENCES ................................... 64

BIOGRAPHICAL SKETCH ............................ 68











Abstract of Dissertation
Presented to the Graduate School of the University of Florida
in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy


A TOPOLOGY-INDEPENDENT SHAPE MODELING SCHEME By

Ravikanth Malladi

December 1993


Chairman: Dr. Baba C. Vemuri
Cochairman: Dr. James A. Sethian
Major Department: Computer and Information Sciences

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 dissertation 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-Jacobi" 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.

vii










The resulting equation of motion is solved by employing entropy-satisfying upwind finite difference schemes. We also introduce a new algorithm for rapid advancement of the front using what we call a narrow-band updation scheme. This leads to significant improvement in the time complexity of the shape recovery procedure in 2D. 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. We also demonstrate the recovery of 3D shapes.
















CHAPTER 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 model-based 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 be either rigid or nonrigid. In the former case, we have what are popularly known as correlation-based template matching techniques, while the later involves a dynamic model fitting process to the data. Shape models serve the purpose of an intermediate stage in object recognition tasks, since they provide a more stable and useful description than the original raw data. In this dissertation, we present a new dynamic approach to shape modeling which retains the most attractive features of existing methods, and overcomes their prominent limitations. For the rest of this dissertation 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. The shape models must aid in recovery of detailed structure from noisy data using only the weakest of the possible assumptions about the observed shape. Several variational shape reconstruction methods have been proposed and there is abundant literature on the same [5, 37, 42, 6, 50, 22]. Generalized spline models with continuity constraints are well suited for fulfilling the goals of shape recovery (see [6, 7, 41]). Generalized splines are the key ingredient of the dynamic shape modeling paradigm introduced by Terzopoulos et al., [45]. 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 [17, 46, 47, 49, 51, 10, 12]. 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 "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 dissertation, 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 [35].

The framework of energy minimization has also been used successfully in the past for extracting salient image contours - edges and lines. Kass et al. [171 used energyminimizing "snakes" that exhibit an affinity toward image features such as edges points and edge segments, while maintaining piecewise smoothness. The weights of the smoothness and image 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 surface models to

















PR
AB



(a) CT image (b) DSA image (c) Holed shapes


Figure 1.1. Test bed for our topology-independent shape modeling scheme.

image data by means of energy minimization [45]. The scheme seeks to design energy functionals whose local minima comprise the set of alternative solutions available to high-level processes. In the absence of a well-developed high-level 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 [9] defines an inflation force on the active contour. This 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.1(b). Moreover, despite a good











initialization, the active contour model, due to its arc-length 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 thin-plate-membrane-spline, in an adaptive environment wherein the material parameters controlling the relative strengths of elasticity and rigidity are adapted (see [36]). 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. Terzopoulos [45] has 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 "sprout" branches during the evolutionary process. Once the shape has been segmented from the image, its constituent part structure can be inferred using the algorithm described in [18].
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 [35]. Returning to











static domain, more recently molecular dynamics has been used to model surfaces of arbitrary topology [40]. Smoothness and continuity constraints are imposed by subjecting a 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 dissertation can be applied 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 shape of interest (see figure 1.1(c)), has the ability to split freely to represent each shape [25, 26]. 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 6.9). Finally, we apply our method to 3D images and recover the surface shapes of 3D objects (see figure (6.12)).

1.1 Ovelview

In this section we briefly outline our scheme to model complex shapes. This method is inspired by ideas first introduced in Osher and Sethian [30, 39] 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 [38] 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











[30, 39] achieve this by viewing the propagating surface as a specific level set of a higher-dimensional function. The equation of motion for this function is reminiscent of an initial valued Hamilton-Jacobi equation with a parabolic right-hand 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 [301, we adopt a global approach and view the (N - 1) dimensional moving surface as a level set of a timedependent function 0 of N space dimensions. The equations of motion written for this higher dimensional function are then amenable to stable entropy-satisfying numerical schemes designed to approximate hyperbolic conservation laws. Topological changes can be handled naturally in this approach, since a particular level set { = 0} of the function V) 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 an negative speed function from the image. Secondly, we have to construct an extension of this speed function to other level sets {1 = C} in order to give a consistent meaning to the image-based speed term at all points in the image (see figure 3.1). 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 [30, 38, 39], has been applied in the area of computer vision for shape characterization by Kimia et al. [18, 19], 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 boundary of 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 4 and recover complex shapes by propagating it along its gradient field. Shape characterization can be done once the object shape is extracted.

An important contribution of this dissertation over our previous work in [26] is the design and implementation of a faster numerical technique for solving the governing equation in 2D Our new algorithm exploits the fact that the front, which is a particular level set { = 0} of a higher dimensional function, can be advanced by updating the function 0 at a small set of points. This scheme is an alternative to updating the 4 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 { = 0}. Since the narrow-band 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 narrow-band techniques for interface propagation may be found in [1]. We also show that the surface shapes of objects can be recovered from 3D image data.

1.2 Dissertation Outline

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







8


for high-level processing. The remainder of this dissertation is organized as follows: chapter 2 introduces the curvature-dependent front propagation problem and establishes a link between Hamilton-Jacobi equations and a hyperbolic conservation law. In chapter 3 we explain our level set algorithm for shape recovery and in chapter 4 we outline the details of our narrow-band approach. Finally, in chapter 5 we present some experimental results of applying our method to some noisy and low contrast medical images in 2D as well as 3D. We close with a discussion of advantages of our approach and direction of future research in chapter 6.
















CHAPTER 2
ENERGY METHODS FOR SHAPE MODELING

In this chapter, we review the some of the literature on existing shape modeling methods in computer vision. Specifically, we provide the details of two "energybased" methods for shape modeling and reconstruction. The first is the symmetryseeking deformable model [451 and the second is a shape recovery scheme in 2D called active contour model or "snake" [17].

Existing shape modeling methods roughly fall under two categories: passive and active models. Models of shape such as generalized cylinders, introduced by Binford [4], and lumped-parameter family of shapes such as superquadric models [3, 32] are purely geometric, hence passive. Generalized cylinders are used to model elongated shapes with axial symmetry, while the superquadric shape models are well suited for object recognition tasks because one can express them compactly using a small set of parameters. On the contrary, active models are governed by the principles of Lagrangian dynamics and react to external forces in a way that is similar to real elastic objects. Generalized spline models with continuity constraints [7, 10, 12, 41, 44, 49, 51] are prime examples of the active shape modeling paradigm.

2.1 Symmetry-Seeking Deformable Models

The deformable model is a physically based modeling framework for shape and motion reconstruction of free-form flexible objects. In this framework the objects are modeled as elastically deformable bodies subject to continuum mechanical laws. Constraints are expressed as forces which deform the elastic models and propel them











through potentially complicated motions such that they satisfy the available constraints over time. There are two types of forces, intrinsic and extrinsic.

Intrinsic constraints reflect generically valid assumptions about natural objects while the extrinsic constraints reflect observations about the environment that can be extracted from sensory data. The model is embedded in a force field which encodes the available sensory information. The ambient force field will then mold the deformable model to make their 3D shape consistent with the observed data. A plethora of external force fields can be synthesized from the data using a wide variety of constraints. Contrary to the conventional models of 3D shape, the deformable models are active in that they react to extrinsic forces as one would expect real elastic objects to react to physical forces. This is because the deformable models are governed by the principles of elasticity theory as expressed through Lagrangian dynamics. Specifically, the continium mechanical equation a2v av S$(v) fv 21
P--+ --- = f (v) (2.1) �t I -at 6v

governs the nonrigid motion of the deformable model in response to a net extrinsic force f(v), where i is the mass density function of the deformable body and 7 is the viscosity function of the ambient medium.

2.1.1 The deformable model

We begin this subsection by providing an informal description of the deformable model. Consider a deformable sheet made of elastic material (a membrane-thin-plate hybrid). This sheet is rolled to form a tube. Next, a deformable spine made of similar material is passed down the length of the tube. At regularly spaced points along the spine, it is coupled to the tube with radially projecting forces so as to maintain the spine in approximate axial position within the tube. Additional forces are included that coerce the tube into a quasi-symmetric shape around the spine. Finally, extra











control is provided over the shape by introducing expansion/contraction forces radiating from the spine. The rigidities of the spine and the tube are independently controllable, and their natural rest metrics and curvatures can either be prescribed or modified dynamically. If the circumferential metric of the tube is set to zero, for instance, the tube will tend to contract around the spine, unless the other forces prevail. The model will shorten or lengthen as the longitudinal metrics of the tube and the spine are modified. In short, a variety of interesting behaviors can be obtained by adjusting the control variables designed into the model.

The deformable tube and the spine can be described by geometric mappings from material coordinate domains into Euclidean three-space W. The mapping is expressed by position vectors in space whose component functions are time varying. The spine is a deformable space curve defined by mapping a univariate material coordinate domain s E [0, 1] into R'

v : [0,1] x [0,x3] --+ R' such that v(s,t) = (X(s,t),Y(s,t),Z(s,t)) (2.2) The tube is constructed from a deformable sheet that is defined by mapping a bivariate material coordinate domain (u, v) C [0, 1]2 into R' V: [0,12 x [0,ooj - W3 such that v(u,v,t)= (X(u,v,t),Y(u,v,t),Z(u,v,t)) (2.3)


The tubular shape is created by specifying boundary conditions on two opposite edges of the sheet which effectively glue these edges together. The edges u = 0 and u = 1 are "seamed" together, letting v span the length of the tube. The required periodic boundary conditions are

v(,v,t)= v(1,v,t), (2.4) au 0v











The strain energy E' associated with the spine function v(s, t) is given by

I01 +WV 12 2V 12 ds, (2.6)
.6S(m = Wl - ._9S]_ 2 .$


where the vertical bars denote the Euclidean norms and the quantities Wl, w2 are both functions of (s,t). wi(s,t) determines the local tension along the spine and w2(s,t) determines its local rigidity. The weighting functions wI, w2 control the elastic properties of the spine. We can regulate the equilibrium shape of the spine by suitably defining the weighting functions.

The strain energy ET associated with the tube function v(u, v, t) is given by

ET(V 11W10I &V 12 +01 1 I2+2 Iav1
0 92 v a2
fol 1V 2 V 2 12V 1

2w O1 2W02 I 12 dudv. (2.7)

The weighting functions wl0 and w0l define the natural metric of the tube along each parameter curve. These functions locally control the tension of the deformable sheet (that constitutes the tube) along each material coordinate curve. Analogous expressions for W20(U, v), wn(u, v) and wo2(U, v) determining the natural curvature of the sheet can be written down. These functions locally control the rigidity of the sheet. The continuity of the deformable spine can be controlled via the functions w, and w2. Position discontinuity is introduced by setting wi(s,t) = w2(s,t) = 0 and tangent discontinuity by setting w2(s,t) = 0. Similarly, we can introduce discontinuities in position as well as tangents for the deformable tube.
The spine and the tube functions are coupled by first identifying that v - S, bringing into correspondence the spine coordinate with the coordinate along the length of the tube. The mapping function of the spine vS(v, t) is then distinguished from that of the tube vT(u, v, t) with superscripts S and T.











Summing all the forces into the equations of motion (2.1) associated with the spine and tube, the following dynamic system describing the motion of the deformable model can be obtained:
02vS Ovs es P (2.8) j-2 + t VS -VS
02vT VT bET 6pT
O-2 Olt + - v (2.9) where fS and fT are the coupling forces between the spine and tube that force the spine in a central position within the tube and coerce the tube to radial symmetry around the spine. PS(vs) and pT(vT) are the generalized potential functions associated with the spine and tube respectively. The extrinsic forces on the right hand side of equation are the variational derivatives of the generalized potential functions Ps and pT that are synthesized from the data constraints.

2.1.2 Extrinsic forces

In this subsection, we follow the argument in [49] and describe the extrinsic forces that are applied to the deformable model to obtain conformation to the range data. Basically, fitting a deformable model to range data is equivalent to applying an extrinsic force to the model until the model achieves maximal conformation with the data. The extrinsic forces, as mentioned earlier, can be defined as variational derivatives of potential functions. Techniques for generating suitable potential functions from monocular, binocular and dynamic image sequences are described in [45]. The extrinsic force on the tube due to the depth data is defined by the variational derivative of the potential function, namely

Pd(T lv - a ll (2.10)
2

where fQ specifies the data domain, vT denotes the vector function of the tube, and d denotes the data vector from the range finder. ad is the inverse of the variance aa











along the vector a defined from the data point to the point of influence on the model. The force associated with the above potential is given by f (VT) = d(VT - d) (2.11) Assuming no external forces on the spine i.e., the potential function P's = 0 for the spine. The dynamic equations that are solved for the surface fitting problem are then given by,
O2vT OVT 6gT (.2 a- '9t T - IT = ad(vT - d) + fT, (2.12) a2vS (vS w5�s
at2 +7-5- + vs2 fs (2.13) The differential equations (2.12) and (2.13) pose a nonlinear initial boundary value problem. A solution is obtained by discretizing in spatial and time domains. The discretization in spatial domain is carried out using standard finite central difference formulas on regular grid of nodes [20]. Since the extrinsic forces are obtained from a static field, the integration in time is performed until y, the viscous damping term, dissipates all the kinetic energy in the model thereby bringing it to a static equilibrium. A detailed discussion of various numerical procedures that can be employed for solving the above system of differential equation can be found in Terzopoulos et al. [431.

2.2 Snake: The Active Contour Model

Snakes are an example of a more general technique of matching a deformable model to an image by means of energy-minimization. The idea of using energyminimizing splines to extract image features of interest has been introduced by Kass et al. [17]. Snake is guided by external image forces that guide it toward features such as lines and edges.










The snake is a model of a deformable contour composed of abstract elastic materials. Two types of materials determine the snake configuration: string and rod. The former makes the snake resistant to stretching and the latter makes it resistant to bending. This deformable curve is then embedded in a potential functional synthesized from a variety of image constraints. Depending on the nature of the external potential field, the snake will extract different image feature. Unlike most other techniques for finding salient image features, the snake model is active. It is always minimizing its energy functional and therefore exhibits dynamic behavior.

2.2.1 Snake dynamics

The basic snake model is a controlled-continuity spline [41] under the influence of image forces and external constraint forces. The internal spline forces serve to impose a piecewise smoothness constraint and the external constraint forces are responsible for putting the snake near the desired local minimum.

Consider a deformable curve v(s, t) with parameters s and t (time), defined on given intervals 0 = [0..11, and T, respectively. Let this deformable curve be a function of two coordinate variables x and y with the same parameterization: v(s,t) = (x(s,t),y(s,t)) : s E Q,t E T. (2.14) The potential energy functional of the snake Eak,(V) is defined as

( j(Et + Ext)ds, (2.15) where Eit is the potential associated with external constraint forces. Ein denotes the internal potential functional of the snake and is given by:

Eint(V) Wl(S) I v. I' +w2(s) vs 12 ds. (2.16) The first-order term wi(s) I v, 12 is responsible for C' continuity in the snake and the second-order term w2(s)I v." 12 enforces Ca continuity. The weight w1(s) regulates











the tension of the snake, whereas w2(s) regulates its rigidity. Position and tangent discontinuities may be introduced along a snake by setting these weights to zero.

Shape recovery in 2D can be accomplished by defining a simple energy functional from a given image. If we set

Eext(v) = - I f VG,, * I(x,y) 12 (2.17) then the snake is attracted by the local minima of the potential, which corresponds to the local maxima of the image gradient. Here, G, * I(x, y) denotes the image function convolved with a Gaussian smoothing filter whose characteristic width is a. Potential functionals modeling other external constraints defined by the user may be added to the above expression.

Given the total potential energy of the snake (see equation (2.15)) for a specific initial position, the equilibrium configuration can be reached by minimizing it. Specifically, if the function v* represents a local minimum for Es,,ak, it satisfies the the following Euler-Lagrange equation:

*- a( v) (2.18)
fv O s ()V*) + j-j(W2(S)VS*3) = VEx4tV*). 2.8

In this formulation, each term appears as a force applied on the curve. Given the value of v(s, t = 0), the initial snake configuration, the above equation can be solved numerically by adopting direct or iterative finite difference schemes [17, 20].
















CHAPTER 3
FRONT PROPAGATION PROBLEM

In this chapter, we introduce the front propagation problem. As a starting point and motivation for the level set approach, we first discuss the interface motion in terms of a parameterization of the moving front. For details and an expository review, see Sethian [38, 39]. The numerical schemes based on discrete approximations to these equations are shown to have limitations. We then present the level set technique due to Osher and Sethian [30], which overcomes the limitations inherent in the former approach.

3.1 Lagrangian Formulation

Consider a closed curve moving in the plane, that is, let ^1(0) be a smooth, closed initial curve in Euclidean plane R2, and let -y(t) be the one-parameter 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 -y(t) by s, 0 < s < S. The curve is parameterized so that the interior is to the left of the direction of increasing s. With n and K(s, t) as the outward normal and curvature of the curve respectively, the equality n xt = F(K) holds. The equation of motion in terms of individual components of x (x, y) is


[x + y) 1 ( (xS )

yt = F Y.1X - XSSYJ]( . (3.1) 2~x + y.2)1 (X2 Y
(X3 Iy)/J (x �y')1/2J
This is called a "Lagrangian" representation because the physical coordinate system moves with the front.











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. However, there are several problems with this approach, as discussed in Sethian [381. 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 dimensions.

3.2 Level Set Formulation

As an alternative, the central idea in the level set approach of Osher and Sethian [301 is to represent the front -y(�) as the level set {1 = 0} of a function 4'. To motivate this approach, we consider the example of an expanding circle. Suppose the initial front / at t = 0 is a circle in the xy-plane (figure 3.1(a)). We imagine that the circle is the level set { = 0} of an initial surface z = 4'(x,y,t = 0) in R3 (see figure 3.1(b)). We can then match the one-parameter family of moving curves -Y(t) with a one-parameter family of moving surfaces in such a way that the level set {4= 0} always yields the moving front (figures 3.1(c) & 3.1(d)).

In the general case, let -y(O) be a closed, nonintersecting, (N - 1) dimensional hypersurface. Let 0 (x,t), x E RN, be the scalar function such that 4'(x,0)= +d(x), where d(x) is the signed distance from x to the hypersurface y(O). We use the plus










Y Z = V(XYt-O)



7(e) X * J'W W = Sd
Y x Y I~ 'W x,y,t)

(C) (d) -- - - - -- - - - -


x --- -- V2~e 4=0 D7to) w-:_-x

Figure 3.1. Level set formulation of equations of motion - (a) & (b) show the curve -Y and the surface t(x, y) at t = 0, and (c) & (d) show the curve y and the corresponding surface -0 (x, y) at time t. sign if x is outside -y(O) and minus sign if x is inside. Each level set of V) 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 V).
Consider the motion of some level set {f? = C}. Following the derivation in [30], let x(t) be the trajectory of a particle located on this level set, so V'(x(t),t) = C. (3.2) The particle speed Ox/at in the direction n normal to -y(t) is given by the speed function F. Thus,

a- - n -- F, (3.3)


where the normal vector n is given by n = VO/ I Vj 1. By the chain rule we get, ft
Ot + -V�k = 0 (3.4)










and substitution yields

O, + F I V�k 1= 0, (3.5)


with an initial condition O(x,0) = �d(x). We refer to equation (3.5) as the level set "Hamilton-Jacobi" formulation. Note that at any time, the moving front 1(t) is simply the level set { (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 {P = 0} need not be simply connected. The signed distance function 4(x, t) always remains a function, even if the level surface { = 0} corresponding to the front -y(t) changes topology, or forms sharp corners. The geometric and differential properties of -y(t) are captured in the function 0 and can be readily extracted. As an example, if x E R', the curvature is given by K = (Oy, - 20.zyik ~ (3.6) (IP2 + 2)3/23

This approach can also be easily extended to higher dimensions and appropriate expressions can be obtained for the mean curvature and the Gaussian curvature [30]. For example, the mean curvature at a point (x, y, z) on a moving surface can be computed using the following formula:

g + 7k2) + ?k(yy(0 + ?)2) + i).2(k'2 + 02)_ 20.yO.Oy_ 20.�O.O - 20y. Oye) (0 2 + V) + 0'2)3/2

By substituting F(K) = 1 - K as a typical speed function in equation (3.5), the equation of motion becomes


I,+ IVP =J K IV4 I (3.I-


(3.7)







21


Equation (3.7) resembles a Hamilton-Jacobi equation with viscosity, where "viscosity" refers to the second-order parabolic right-hand side. This equation can be solved using the stable, entropy-satisfying finite difference schemes, borrowed from the literature on hyperbolic conservation laws (see [301).
















CHAPTER 4
SHAPE RECOVERY WITH FRONT PROPAGATION

In this chapter, we describe how the level set formulation for the front propagation problem discussed in the previous chapter 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.

4.1 Stopping Criterion: An Image-based Speed Term

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 [9]. The second term Fa, 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 thin-plate-membrane splines [17 (see the figure (6.3)). We rewrite equation (3.7) by splitting the influence of F as ?kt + FA I VO I +Fc I VO J= 0. (4.1) 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 F1 to be

fikx, Y) = OW { VG, * I(x,y) 1.- A121 (4.2)


where M1 a-.d M11 are the maximum and minimmm values of the magnitude of image gradient I VG, * I(x, y) 1, (x, y) C Q. The expression G, * I denotes the image convolved with a Gaussian smoothing filter whose characteristic width is o,. Alternately, we could use a smoothed zero-crossing image to synthesize the negative speed function. The zero-crossing image is produced by detecting zero-crossings in the function V2G, * I, which is the original image convolved with a Laplacian-of-Gaussian filter whose characteristic width is a. The equation of motion with the addition of image-based speed is


+(t+(FA + F) V 1= 0. (4.3) FI is called an extension of F1 to points away from the boundary 'Y(t), i.e. at points

(xy) E (9 - -y(t)), and is equal to F on -(t). We shall return to the issue of











extension shortly. The value of F, 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 kz. The term k1, which is defined as
1
ki(x,y) = 1+ 1 VG, * I(x,y) 1' (4.4) 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 "steerable" filters [16]. The modified equation of motion is given by Ot - kI(FA + F) I V4 1= 0. (4.5)

4.2 Extension of Image-based Speed Term

We now come to an important juncture in our discussion. The image-based speed term, be it FI or kI, has meaning only on the boundary -y(t), i.e. on the level set {� = 0}. This follows from the fact that they were designed to force the propagating level set { = 0} to a complete stop in the neighborhood of an object boundary. However, the equation of motion (4.3) is written for the function 0, which is made up of infinitely many level curves. In other words, equations (4.3) & (4.5) 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











Y








C

SX

Figure 4.1. Huygen's principle construction sets, i.e. at every point (x, y) C fQ. Speed functions such as FG which are functions of geometric properties of the surface z = b(x, y), can be readily computed at any (x, y) E Ql. However, F1 is not such a function. It derives its meaning not from the geometry of b but from the configuration of the level set {V = 0 in the image plane. Thus, our goal is to construct an image-based speed function P, that is globally defined. We call it an extension of FI off the level set {0 = 0} because it extends the meaning of F1 to other level sets. Note that the level set {O/ = 0} lies in the image plane and therefore P, must equal F on {4' = 0}. The same argument applies to the coefficient kx.

If the level curves are moving with a constant speed, i.e. F- = 0, then at any time t, a typical level set {V = C}, C E R, is a distance C away from the level set {4 = 0} (see figure 4.1). 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


V












Q (x,y)








X

Figure 4.2. Extension of image-based speed terms to other level sets

function F = FA + FG is continuous (see [141). With the above remarks in mind we state the following:

Property 4.2.1 Any external (image-based) speed function that is used in the equation of motion written for the function V should not cause the level sets to collide and cross each other during the evolutionary process.

Recall that the function O(x, i) has been initialized to d(x), where d(x) is the signed distance from a point x E Q to the boundary -y(O). Since we cannot attribute any geometric meaning to the function F (ki) at points away from the level set {I = 01, we look for a meaning that is consistent with property (4.2.1). Therefore, the question to ask is: what is the value of Fj (or k1) at a point (x, y) lying on a level set {= C}? We answer this question in the following construction (see figure 4.2). Construction 4.2.1 The value of Fi ('ki) at a point P lying on a level set {= C} is exactly the value of F, (k1) at a point Q, such that point Q is a distance C away from P and lies on the level set {k = o}.

It is easy to see that 1i reduces to F1 on {' = 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 4.2.1, then it will also be consistent with the property 4.2.1. Thus, having ascribed the intent of pseudodifferential equations (4.3) & (4.5) in the context of shape modeling, we can use finite difference schemes to solve them numerically. Since 0 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.3 Shape Recovery in 3D

In this section, we show that the surface shapes of three-dimensional objects can be recovered using a scheme that is analogous to the one described above. We begin by writing the equation of motion for a moving surface; this represents front propagation in higher dimension. Let the surface be closed, nonintersecting, and the level set {k = 0} of a function 0 (x,t), x E R'. The function #(x,t = 0) = �d(x), where d(x) is the signed distance function from the point x to the initial surface; plus sign is used if x is outside the initial surface and minus sign if x is inside. Each level set of V moves along its gradient field with speed F(K). Following the derivation in chapter 2, we arrive at the equation of motion for the moving surface, namely �t + F(K) I V4 J= 0, (4.6) where K is either mean or gaussian curvature of the surface. Recall that the surface at any time t is simply the level set {V = 0}.

In this work, we only address the issue of recovering surface shapes from 3D images. The input data constitutes a set of MRI (Magnetic Resonance Imaging) or CT (Computed Tomography) images. A stack of such images represents an image function I(x, y, z). Our next step is to define an image-based speed term which will


V











cause the level set {?, = O} to stop near the 3D object boundary. To this end, we apply on the front the speed equal to the product of F(K) and a term ki. The term kI, which is defined as
1
ki(x,y,z) = (4.7) 1+ 1 VG, *I(x, y,z) I

has values closer to zero in the vicinity of object boundaries and values that axe closer to unity in regions of relatively constant image intensity. To attribute a consistent meaning to the image-based speed term at all the level surfaces, we construct an extension to k1 in a manner similar to the construction (4.2.1). Therefore, the modified equation of motion that we solve to extract complex surface shapes from 3D images is


Ot + kIF(K) I VV) J= (.


(4.8)
















CHAPTER 5
NUMERICAL SOLUTION

The equation (3.7) poses an initial valued problem. It is rewritten here as

ot + (,o2+ 2)1/2 = 6V' ]V (5.1) X Y (5.1)

with (x,y,t = 0) = distance from (x,y) to y(O). As shown in Sethian [381, for c > 0, the parabolic right-hand side diffuses sharp gradients and forces 7P to stay smooth at all values of t. For - = 0, the boundary moves with unit speed, and a corner must develop from smooth initial data. Once a corner develops, it is not dear how to propagate the front in the normal direction, since the derivative is not defined at the corner. A variety of "weak" 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 [38]: 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 [30, 39) are motivated by the fact that the entropy condition for propagating fronts is identical to the one for hyperbolic conservation laws, where stable, consistent, entropy-satisfying algorithms have a rich history.











5.1 Entropy-satisfying Upwind Numerical Scheme

In this section, almost without a change, we present the arguments discussed in Sethian [391. For complete details of the following scheme, we refer the reader to Osher and Sethian [30, 39].

First, consider the one-dimensional version of the level set equation, with constant normal velocity FA = 1. Then the equation (3.7) becomes a standard HamiltonJacobi equation

Ot + H(O.) = 0, (5.2) where H(�,) =-(, , and with a given initial value of q. Let u = 0.. Taking the derivative with respect to x, equation (5.2) becomes ut + [H(u)]X = 0, (5.3) where H(u) = -(u2)"/2. Equation (5.3) 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 (5.3) are integrated in an arbitrary interval [a, b] to produce

d 1bu(xt)dt = H[u(a,t)] - H[u(b,t)]. (5.4)


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, entropy-satisfying solution to equation (5.4)? The answer is found in the following definition:











Definition 5.1.1 Let u 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 71n+1 - u7 g(u,u - g(u - ,)
At Ax

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 entropy-satisfying weak solution, monotonicity is required, i.e., that u!'+ be an increasing function of the arguments u7_, u7, and un 1. The main fact is: A conservative, monotone scheme produces a solution that satisfies the entropy condition. Equation (5.5) is a scheme for the slope u, which must be converted into a scheme for 0 itself. First write equation (5.2) with a forward difference in time as,

+ = - AtH(u). (5.6) Since the numerical flux function g approximates H, the solution to equation (5.6) may be approximated by

0�n+1 = - Atg(ui_1/2, Ui+1/2) = '-Atg(D-Oi,D+�i), (5.7) where g is an appropriate numerical flux function and the standard definitions of the forward and the backward difference operators have been used, namely, D X -x -A x

D+ _ - (5.8) X - 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 Hamilton-Jacobi flux function given in [30]:

g(ui-ll2, Ui+1l2) =g( D; Oi,DS+ i)

f((max(D-, 0))2 + (min(D+� , 0))2). (5.9) 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 finite-sized 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) =-(u2)1/2, equation (5.7) reduces to + At{(max(D- ,0))2 + (min(D+�,0))2}/2. (5.10) This algorithm produces the correct entropy-satisfying weak solution to the propagating front problem.

The above discussion can be extended to problems in more than one space dimension (see Osher and Sethian [301). Recall that the original intent was to solve equations (4.3) and (4.5). In two dimensions, the scheme given in equation (5.10) is extended by differencing in each direction to produce the following numerical scheme for equation (4.3):
.+ - At[FA + (PI)iiI{(max(D;�,,,0))2 + (min(D+Vkjj, 0))2


+(max(D-j, 0))2 + (min(D+ ,jj, 0))211/2. (5.11) Similarly, the numerical scheme for equation (4.5) is,

= - AtFA(ki) ij{(max(D - ,j, 0))2 + (min(D+Oj, 0))2 (5.12)

+(max(D-O,j, 0))2 + (min(D+ ij,0))2}11/2- AtFGkc I\7" I











The second term FGkl I VV I is not approximated in the above equation; one may use a straightforward central difference approximation to this term.

5.2 Narrow-band Updation Scheme

In this section we first consider a straightforward approach to solve the equations (5.11) & (5.13) 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 { = 01.

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 {V = 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 image-based speed term on the zero set (equation (4.5)). An extension function is constructed to attribute a consistent meaning to the image-based speed term at points away from the zero set (see discussion in section 3). Therefore, a 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 image-based 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 { = 0}, and is closest to the point (iAx,jAy). The value of image-based speed term at the current

point is simply its value at the pcint q.











2. With the value of extended speed term (I),, and V, calculate 0'l using

the upwind, finite difference scheme given in equation (5.13).

3. Construct an approximation for the level set {0 = 0} from . This is

required to visualize the current position of the front in the image plane. A piecewise linear approximation for the front -y(t) is constructed as follows. Given

a cell C(i,j), if max(oij, Vbi+,,j, Oi,j+,, Oi+l,j+l) < 0 or min(Oij, Vbi+l,j, 0i,j+i,

4'O+i,j+l) > 0, then C(i,j) -y(t) and is ignored, else, the entry and exit points where 0b = 0 are found by linear interpolation. This provides two nodes on 'Y(t) and thus, one of the line segments which form the approximation to -(t). The collection of all such line segments constitutes the approximation to the level set {b = 0}, which is used for future evaluation of the image-based speed term

in the update equation (5.13).

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 image-based speed term. This is because at each grid point, we must search for the closest point lying on the level set { = 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 (5.13), 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 image-based speed (Il)i,,. Alternately, we could down-sample 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











Y














Figure 5.1. A narrow-band of width 6 around the level set {4= 0}.

embed our level set algorithm in a multiresolution framework. Instead, we improve the time complexity of our algorithm by updating the level set function , 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.1) the bold curve depicts the level set { = 0} and the shaded region around it is the narrow-band. The narrow-band is bounded on either side by two curves which are a distance 6 apart, i.e., the two curves are the level sets {O' = �6/2}. The value of 6 determines the number of grid points that fall within the narrow-band. Since, during a given time step the value of ?ij is not updated at points lying outside the narrow-band, the level sets {i >1 6/2 1} remain stationary. The zero set which lies inside moves until it collides with the boundary of the narrow-band. Which boundary the front collides with depends on whether it is propagating inward or outward; either which way, it cannot move past the narrow-band's boundary.











Therefore, as a consequence of our update strategy, the front can be moved through a maximum distance of b/2, either inward or outward, at which point we must rebuild an appropriate narrow band. We reinitialize the ' function by treating the current zero set configuration, i.e., {b = 0}, as the initial curve y(O) (see section 2). Note that the reinitialization procedure must account for the case when { = 0} changes topology. This procedure will restore the meaning of 0 function by correcting the inaccuracies introduced as a result of our update algorithm. Once a new 4 function is defined on the grid, we can create a new narrow-band 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 narrow-band, compute the extension

k1 of image-based speed term.

3. With the above value of extended speed term (kI )1,,, and , calculate 0!'t'

using the upwind, finite difference scheme given in equation (5.13).

4. Construct a polygonal approximation for the level set { 0 =0} from 4' +'. A

contour tracing procedure is used to obtain a polygonal approximation. Given a cell (ij) which contains -y(t), this procedure traces the contour by scanning the neighboring cells in order to find the next cell which contains -i(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 '-y(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 {4 = 0}. A polygonal approximation is required in step 2 for the evaluation of image-based speed term and more importantly, in step 6 for

reinitializing the V) function.

5. Increment m by one. If the value of rn equals 1, go to step 6, else, go to step 2.

6. Compute the value of signed distance function b by treating the polygonal
approximation of { = 0} as the initial contour -y(O). As mentioned earlier, a more general method of reinitialization is required when { 0 O} 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 -' was initialized and updated on a square grid. To update Oij at a boundary point (ij), it is necessary to specify the boundary conditions in order to numerically approximate the derivatives. We have chosen free-end 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 0 at points lying in the narrow-band, the issue of specifying boundary conditions for points lying on the edge of the band becomes pertinent. The focus of front propagation in the context of shape recovery is merely rapid advancement of the front. Therefore, at the expense of a little inaccuracy in the zero set configuration, we ignore the sensitive issue of specifying correct boundary conditions. In our implementation, derivatives are evaluated normally even if it involves a point that does not belong











in the narrow-band. On the other hand, in 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 (5.2), 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 (3.6). 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 0 function was recomputed once every (I =) 40 time steps. In figure 5.2(a), we show the initial curve along with the level sets {- = �i0.04}, where i C [1..5]. After 40 narrow-band updations (figure 5.2(b)), only the level sets { <1 0.08 1} move and the rest remain stationary. We note the inconsistency between the level sets lying on either side of the narrow-band, making the reinitialization step necessary in order to restore the meaning of the b function. Following the reinitialization step, another 40 update steps are applied (figure 5.2(c)), which "diffuses" 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 5.2(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.

We conclude this section by highlighting some advantages of the narrow-band updation scheme. The first is the obvious improvement in time complexity as a result of considering only a fraction of the total number of grid points for updation. The second is related to the issue of computing the value of image-based speed term at points away from the zero set, i.e., the extension. In our original method, the extension of image-based soeed term could be unreliable at points that are relatively far from the zero set. In addition, when the front is a circular one, it is not clear























(b) After 40 iterations


(c) After 80 iterations


(d) After 120 iterations


(e) After 160 iterations (f) After 200 iterations


Figure 5.2. Narrow-band updation algorithm applied to a star-shaped front propagating with speed F = -K. Calculations were done on a 64 x 64 grid with a time step At = 0.0003. b was recomputed after every 40 time iterations.


(a) Initialization











what value of image-based speed term to use at its center. Here, the fact that the center of a circle is equidistant from all the points lying on it causes an ambiguity in the choice of the shortest distance point on the circular front. Clearly, this effect is undesirable. The narrow-band updation algorithm does not require the image-based speed term's evaluation at far away points. Therefore, the extension is very accurate and reliable.

5.3 Numerical Issues in 3D

In this section, we provide some details of the numerical method based upon the Eulerian formulation [30] that is used to solve the equation (4.8). To begin, the discussion in the previous sections can be readily extended to 3D by differencing in each space direction to produce the following numerical scheme to approximate the surface motion:
I = 0'i, At(O)ij,k{(max(D-lk,J,k, 0))2 + (min(DOi,J,k, 0))2 (5.13)
i,j,k itj,k - tkiJkMa( 5.3

+(max(D-ij,, 0))2 + (min(D+O,j,k, 0))2

+(max(D - Oij,k, 0))2 + (min(D+ ki,j,k, 0))2}/2 - At(eK)(ki) I VO' I.

Note that the level surfaces are moving with speed F = k1(1 - -K), where K is the mean curvature.

In the context of of 3D shape recovery, we move the level surfaces by updating the V) function on a three dimensional grid according to the above update equation. The level surface {4' = 0}, which models the shape of interest, is forced to a stop in the vicinity of desired object boundaries by applying an image-based speed coefficient on it. In order to to compute the value of image-based speed coefficient at points away from the zero set, an extension is constructed. More specifically, the shape recovery procedure in 3D consists of the following steps:













Algorithm 3

1. At each grid point (%Ax, jAy, kAz), the extension of image-based speed coefficient k, is computed.

2. With the value of extended speed coefficient (k Ii,j,k, and lk~j,k, calculate,,+

using the upwind, finite difference scheme given in equation (5.14).

3. Construct an approximation for the level surface { 0= O} from 0'J'. This is

done in the following fashion. Across every two vertices of a given cell C(i, j, k), the value of b function is sampled. If the 0 function changes sign across a given edge, the exact point at which 0 = 0 is found by linear interpolation. Note that at each cell, the 4 function can change sign across a number of edges that lies between 3 and 6. The collection of all such points constitutes the approximation to the level surface {' = 0}, which is used for the evaluation of k, during the

next time step.

4. Replace n by n + 1 and return to step 1.

In the previous section, it has been established that the extension calculation is very expensive. This observation has motivated the narrow-band approach for front propagation. The situation in 3D is even worse due to an extra dimension. In addition, there is a large degree of repetition in the zero set due to the lack of an organized method in which to populate the zero set. By organized method we mean a 3D analog of a 2D contour tracing procedure. This limitation adds to the size of the zero set and subsequently to the time complexity of extension calculation. In 2D, we alleviated this situation by adopting a narrow-band update strategy, but implementing a similar scheme in 3D could prove overly cumbersome. Therefore. we exploit








42


the inherent data independence that is present in our grid calculation, and devise a parallel implementation to move the level surfaces. The parallel implementation was done on a KSR (Kendall Square Research) shared memory parallel machine. For a grid with N3 points, we have used N processors and realized a substantial improvement in the running time of our 3D shape recovery algorithm. Steps 1 & 2 of the above algorithm are done in parallel and step 3 is executed in serial mode.
















CHAPTER 6
EXPERIMENTAL RESULTS

In this chapter we present several shape recovery results that were obtained by applying the narrow-band level set algorithm to image data. Subsequently, we present some results in 3D. First, in 2D, given an image, our method requires the user to provide an initial contour -y(O). 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 'y(O), our algorithm does not require any further user interaction.

The initial value of the function �- i.e., O(x, 0) is computed from 'y(O). We first discretize the level set function V on the image plane and denote Oi,j as the value of 0 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 2k x 2k 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 '(0). The magnitude of tkj is set to this value. We use the plus sign if (i,j) is outside -1(0) and minus sign if (i,j) is inside. Once the value of Oij is computed at time t = 0 by following the above procedure, we use the update equations from the previous chapter to move the front.











6.1 2D Shape Recovery Results

We now present our shape recovery results in 2D. We first consider a 256 x 256 CT (computed tomography) image of an abdominal section shown in figure 6.1(a). Our intention is to recover the shape of the liver in this particular slice. The function 0 has been discretized on a 128 x 128 mesh, i.e., calculations are performed at every second pixel. In figure 6.2(a), we show the closed contour that the user places inside the desired shape at time t 0. The function 0 is then made to propagate in the normal direction with speed F k(-1.0 - 0.025K). We have employed the narrow-band updation algorithm to move the front. The time step size was set to At = 0.0005 and the 0 function was recomputed after every 50 time steps. Figure 6.1(b) shows the image-based speed term which is synthesized according to equation (4.4). This term when applied on the propagating front, acts as a stopping criterion, thereby causing the front to come to a rest in the neighborhood of a desired shape. Note that in figure 6.1(b), k1(x, y) values lying in the interval [0..1] have been mapped into the interval [0..255]. In figures 6.2(b) through 6.2(e) we depict the configuration of the level set {' = 0} at four intermediate time instants. The final result is achieved after 575 time iterations and is shown in figure 6.2(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 6.1(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 high 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 --K. The objective
























(a) Original image


Figure 6.1. Image-based speed term ki(x, y) = 1 with a = 3.25, synthesized from the CT image.

of our next experiment is to demonstrate smoothness control in the context of shape recovery. In figures 6.3(a) through 6.3(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 4(x, y) from these curves. The level sets are then made to propagate with speed F = ki(1.0-eK). First, as shown in figure 6.3(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 6.3(b) and 0.75 in figure 6.3(c). Clearly, with every increment in the value of -, the level set { = 0} attains a configuration that is relatively smoother. This is analogous to the controlled-continuity provided by thin-plate-membrane splines [41, 17].

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 branches and significant protrusions. In this experiment we compare the performance of our scheme with the active contour


(b) Speed term ki






















(b) After 100 iterations


(c) After 175 iterations


(e) After 450 iterations


(d) After 300 iterations


(f) After 575 iterations


Figure 6.2. Recovery of the liver shape from a CT image of an abdominal section. Narrow-band computation was done on a 128 x 128 grid - the front was made to propagate with speed F = ki(-1.o- 0.025K) and the time step At was set to
0.0005. k was recomputed once every 50 time steps.


(a) Initialization






















(a) c = 0.05 (b) E = 0.25 (c) - = 0.75


Figure 6.3. Smoothness control in shape recovery can be achieved by varying the curvature component in the speed F = ki(1.0 - eK). model to bring the limitations of the later into focus. We first attempt to reconstruct the complex arterial structure using a snake model with inflation forces [9]. In figures 6.4(a) through 6.4(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.4(c), 6.4(f), & 6.4(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.4(g)), the snake "snaps" back into a relatively "bumpless" configuration in figure 6.4(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 deformiation energies. Now,











we apply our level set algorithm to reconstruct the same shape. After initialization in figure 6.5(a), the front is made to propagate in the normal direction. We employ the narrow-band 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 "flows" into the branches and finally in 6.5(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 single instance of our shape model "sprouts" branches and recovers all the connected components of a given shape. In figures 6.6(a)-(i) we plot the other level sets to elucidate the process of extending the image-based speed function to points away from the zero set. 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 transformation to reconstruct the constituent shapes in an image. The image shown in figure 6.7(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 { = 0} first wraps itself tightly around the objects (see figures 6.7(c) & 6.7(d)) and subsequently splits into four separate closed curves (figure 6.7(e)). While the first three closed segments of {O = 0} recover the three distinct shapes, the one in the middle (see figure 6.7(e)), since it does not enclose any object, eventually disappears. Figure 6.7(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 6.8(a)-(i) we show the level sets { = +iC}, i E [-5.. + 5], with C = 0.05. The function 0 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 hand-written characters. The characters "A" and "B" in the figure (6.9) are examples of shapes with holes. Recovery of hand-written 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 [21, 8] or skeletons [29, 23]. 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 6.9(a), we show the initial contour which encloses all the characters. This contour is then made to propagate inward with a constant speed. Figures 6.9(b) shows an intermediate stage in the front's evolution and in 6.9(c), the front splits into four separate contours. The calculation comes to a halt when in figure 6.9(d), the level set { = 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 6.9(d) as an initial state, and propagate all the fronts inward by momentarily relaxing the image-based speed term. This causes the zero set to move into the character shapes as shown in figure 6.9(e). Finally in figure 6.9(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 128 x 128 grid and the time step At was set to 0.00025.











6.2 Experimental Results in 3D

In this section, we present some front propagation and shape recovery examples in 3D. The first two experiments demonstrate the basic surface motion scheme in 3D with a constant and curvature dependent speed. To visualize the level set { = 0}, we use a tetrahedronization algorithm. Given a function Oi,j,k on a grid, this algorithm first marks off all the points with negative V; values (here we assume that the points lying inside the surface have negative 0 values and the points lying outside have positive 7 values). A tetrahedronization procedure is then applied to these points. The triangular facets that lie on the boundary of the tetrahedron constitutes our approximation to the level surface { = 0}. Using these surface triangles the level set { -= 0}, which forms our 3D surface model, is rendered as a shaded shell.

In our first 3D experiment, we evolve a toroidal surface with constant speed F = 1. The computational box is a cubical one which extends from -0.75 to 0.75 in all three coordinate directions. Initial value of 4b is defined according to the following equation:


(xly, z) = 0.015 - {(x2 + y2)1/2 - 0.4}2 _ z2. (6.1)


The set of all points (x, y, z) at which the value of / = 0, lie on the surface of the torus (see figure 6.10(a)). We discretize the 0 on a 32 x 32 x 32 grid and follow a nested set of toroidal shapes by updating it on the grid according to the equation (5.14). Figures 6.10(b)-6.10(d) depict three intermediate stages wherein the inner radius gradually becomes smaller until in figure 6-10(e) the inner radius collapses to zero, normals collide, and topology changes. The surface at this stage looks like a dented sphere and will asymptotically attain a spherical configuration. In figure 6.10(f), the level set {4 = 0} intersects with the edge of the computational box causing the slicing of the toroidal surface.











In our next experiment, we collapse a "pinched" superquadric shape under its mean curvature. The initialization was done using the implicit function form of a superquadric, namely

i/(x,y,z) + +-(s)] -1.35, (6.2)

where C1, 62 > 0 are "squareness" parameters. By varying the values of c, and (2, several interesting shapes can be generated. Flat beveled shapes are produced when either c or E2 = 2 and pinched shapes are produced when either 61 or E2 > 2. Figure 6.11(a) shows the initial superquadric shape with sharp corners. This shape is made to move with speed F = -K, where K is the mean curvature. As shown in figures 6.11(b)- 6.11(d), the high curvature regions corresponding to the sharp corners on the surface smoothly diffuse until the superquadric collapses into an ellipsoid in figure 6.11(e). Note that an ellipsoid is a superquadric with c, = 2 = 1. Since the mean curvature on the surface of an ellipsoid is constant, the surface eventually collapses to a point without changing its shape (see the smaller ellipsoid in figure 6.11(f)). Calculations were done on a 32 x 32 x 32 grid with a time step At = 0.0001.
In our last experiment, we recover the shape of a flat superquadric using the level set front propagation scheme in 3D. Image data for this experiment consists of 32 images each with a particular cross section of the superquadric. The image-based speed term ki is computed from these images according to equation (4.7). A sphere, which is the level surface {V = 0} of a function ?P(x,y,z) = x2 + y2 + z2 -0.01, forms our initialization (see figure 6.12(a)). This initial surface is moved with speed F =ki by updating the value of 0b on a discrete grid according to the details given in algorithm 3. The initial surface expands smoothly in all directions until a portion of it collides with the superquadric boundary. At points with high gradient, the k, values are close to zero and cause the zero set to locally come to stop near the boundary








52


of the superquadric shape. This situation is depicted in figures 6.12(b)- 6.12(e), wherein the initial spherical shape transforms into a flat superquadric. Finally, in figure 6.12(f), all the points on our shape model are stopped, thereby recovering the entire shape of the flat superquadric. Calculations were done on a 32 x 32 x 32 with a time step At = 0.0025.























(a) Initialization 1 (b) 500 iterations


(d) Initialization 2 (e) 500 iterations


(g) Initialization 3 (h) 500 iterations


(f) 1000 iterations


(i) 1000 iterations


Figure 6.4. An tsuccessful attempt to reconstruct a complex shape with "significant" 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.


(c) 1000 iterations























(b) After 60 iterations


(c) After 123 iterations


(e) After 275 iterations


(d) After 200 iterations


(f) After 391 iterations


Figure 6.5. Reconstruction of a shape with "significant" protrusions: an arterial tree structure. Computation was done on a 64 x 64 grid with a time step At = 0.001. A narrow-band updation scheme was used with a band width of & = 0.075.


(a) Initialization
























(a) Initialization (b) 60 iterations


(c) 123 iterations


(d) 170 iterations


(e) 200 iterations


(f) 275 iterations


(g) 320 iterations (h) 350 iterations


(i) 390 iterations


Figure 6.6. Reconstruction of a shape with protrusions: figure shows the level sets
{ =+iC}, i E [-5.. + 5], and C = 0.05.
























(b) After 30 iterations


(c) After 60 iterations


(d) After 100 iterations


(e) After 125 iterations (f) After 140 iterations


Figure 6.7. 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.


(a) Initialization


























(a) Initialization (b) 60 iterations


(d) 100 iterations (e) 120 iterations


(f) 130 iterations


(g) 154 iterations


(h) 170 iterations


(i) 185 iterations


Figure 6.8. Topological split example: level sets shown are {V = �iC}, i [-5.. �5], with C = 0.05.


(c) 89 iterations








ANM
(a) Initialization


AN
(c) 760 iterations
5 B

AH
(e) Second stage


(b) 545 iterations
so
Anl
AN
(d) Outer boundaries
5SB

AH
(f) Inner boundaries


Figure 6.9. Shapes with holes: a two-stage scheme is used to arrive at a complete shape description of both simple shapes and shapes with holes. Computation was performed on 128 x 128 grid and the time step At was set to 0.00025.

























(b) After 10 iterations


(c) After 30 iterations


(e) After 90 iterations


(d) After 70 iterations


(f) After 130 iterations


Figure 6.10. Expanding torus (F = 1): change of topology is depicted in this example. Calculations were done on a 32 x 32 x 32 grid with a time step At = 0.0005.


(a) Initialization
























(b) After 10 iterations


(c) After 30 iterations


(d) After 50 iterations


(e) After 70 iterations (f) After 90 iterations



Figure 6.11. Motion of a "pinched" superquadric shape under its mean curvature. Calculations were done on a 32 x 32 x 32 grid with a time step At = 0.0001.


(a) Initialization

























(b) After 20 iterations


(c) After 40 iterations


(e) After 90 iterations


(d) After 70 iterations


(f) After 120 iterations


Figure 6.12. Shape recovery in 3D: a flat superquadric shape. Calculations were done on a 32 x 32 x 32 grid with a time step At = 0.0025.


(a) Initialization

















CHAPTER 7
CONCLUDING REMARKS

In this dissertation 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 and Sethian [30] to the problem of shape recovery. With this approach, complex shapes can be reconstructed. 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 Hamilton-Jacobi equation with a parabolic righthand side and are amenable to stable entropy-satisfying 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








63


complexity is realized by adopting the narrow-band update strategy. In the narrowband scheme, the front is advanced by updating the 0 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 focussed on extending this method to other application domains. We have implemented 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 chapter 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 narrow-band strategy has drastically reduced the computational cost in 2D, but implementing a similar scheme in 3D could prove cumbersome. One approach is to embed our level set algorithm in a multiresolution framework, where the loss of accuracy at rapidly-converging coarse grid computation is compensated by highly accurate slowly-converging fine grid calculation. Alternately, sophisticated parallel implementations can be devised. One such scheme is briefly discussed in chapter 4.















REFERENCES


[1] D. Adalsteinsson and J. A. Sethian, "A narrow-band approach to level set techniques for propagating interfaces," LBL Report, Oct. 1993.
[2] R. Bajcsy and S. Kova6i, "Multiresolution elastic matching," Computer Vision,
Graphics, and Image Processing, Vol. 46, pp. 1-21, 1989.
[3] R. Bajcsy and F. Solina, "Three dimensional object representation revisited," in
Proceedings of First International Conference on Computer Vision, pp. 231-240,
London, England, 1987.
[4] T. 0. Binford, "Visual Perception by computer," invited talk, IEEE Systems
and Control Conference, Miami, Florida.
[5] R. M. Bolle and B. C. Vemuri, "On three-dimensional surface reconstruction
methods," IEEE Trans. on Pattern Analysis and Machine Intelligence, Vol.
PAMI 13, No. 1, pp. 1-13, 1991.
[6] T.E. Boult and J.R. Kender, "Visual surface reconstruction using sparse depth
data," in Proc. IEEE Conf. on Computer Vision and Pattern Recognition, pp.
68-76, June 1986.
[7] A. Blake and A. Zisserman, Visual Reconstruction, MIT Press, Cambridge, MA.
[8] H. Blum, "A transformation for extracting new descriptors of shape," Models for
the Perception of Speech and Visual Form (W. Wathen-Dunn, ed.), Cambridge
MA: MIT Press, 1967.
[9] L. D. Cohen, "On Active Contour Models and Balloons," Computer Vision, Graphics, and Image Processing, Vol. 53, No. 2, pp. 211-218, March 1991.
[10] 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, Urbana, Illinois, June 1992.
[11] G. Dahlquist, and A. Bj6rck, Numerical Methods, Prentice-Hall, Englewood
Cliffs, NJ, 1974.
[12] 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 Recognition, pp. 467-472, Maui Hawaii, June
1991.
[13] M. P. Do Carmo, Differential Geometry of Curves and Surfaces, Prentice-Hall,
Englewood Cliffs, New Jersey, 1976.











[14] L. C. Evans and J. Spruck, "Motion of level sets by mean curvature. I," Journal
of Differential Geometry, Vol. 33, pp. 635-681, 1991.
[15] J. A. Fessler and A. Macovski, "Object-based 3D reconstruction of arterial trees
from magnetic resonance angiograms," IEEE Transactions on Medical Imaging,
Vol. 10, No. 1, pp. 25-39, March 1991.
[16] W. T. Freeman and E. H. Adelson, "Steerable filters for early vision, image
analysis, and wavelet decomposition," in Proceedings of ICCV, pp. 406-415,
Osaka, Japan, 1990.
[17] M. Kass, A. Witkin, and D. Terzopoulos, "Snakes: Active Contour Models,"
International Journal of Computer Vision, pp. 321-331, 1988.
[18] B. B. Kimia, A. Tenenbaum and S. W. Zucker, "Toward a computational theory
of shape," Lecture Notes in Computer Science, Vol. 427, pp. 402-407, 1990.
[19] B. B. Kimia, A. R. Tannenbaum, and S. W. Zucker, "Shapes, shocks, and deformations I: The components of shape and reaction-diffusion space," to appear
in International Journal of Computer Vision, 1992.
[20] L. Lapidus and G. Pinder, Numerical Solutions of Partial Differential Equations
in Science and Engineering, John Wiley & Sons, New York, 1988.
[21] D. T. Lee, "Medial Axis Transformation of a Planar Shape," IEEE Trans. Patt.
Anal. Machine Intell. PAMI-4, No. 4, pp. 363-369, 1982.
[22] D. Lee and T. Pavlidis, "One-dimensional regularization with discontinuities,"
IEEE Trans. on Pattern Analysis and Machine Intelligence, Vol. PAMI 10, pp.
822-829, 1986.
[23] F. Leymarie and M. D. Levine, "Simulating the grass-fire transform using an
active contour model," IEEE Trans. Patt. Anal. Machine Intell. PAMI-14 No.
1, pp. 56-75, 1992.
[24] R. Malladi, "Deformable models: Canonical parameters for surface representation and multiple view integration," M.S. Thesis, Department of computer and
information sciences, University of Florida, Gainesville, May 1991.
[25] R. Malladi, J. A. Sethian, and B. C. Vemuri, "A topology-independent shape
modeling scheme," in Proceedings of SPIE Conference on Geometric Methods in
Computer Vision II, San Diego, July 1993.
[26] 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
PAM-589, Univ. of California, Berkeley, August 1993.
[27] R. Malladi, J. A. Sethian, and B. C. Vemuri, "A fast level set based algorithm
for topology-independent shape modeling," UF-CIS Technical Report TR-93031, University of Florida, Gainesville, Sept. 1993.
[28] R. Malladi, J. A. Sethian, and B. C. Vemuri, "Evolutionary fronts for topologyindependent shape modeling and recovery," submitted to the Third European
Conference on Computer Vision, Stockholm, Sweden, Oct. 1993.











[29] N. Mayya and V. T. Rajan, "Voronoi Diagrams of Polygons - a robust and practical algorithm," Tech. Report UF-CIS TR-93-29, Univ. of Florida, Gainesville,
1993.
[30] S. Osher and J. A. Sethian, "Fronts propagating with curvature dependent speed:
Algorithms based on Hamilton-Jacobi formulation," Journal of Computational
Physics, Vol. 79, pp. 12-49, 1988.
[311 P. Parent and S.W. Zucker, "Trace inference, curvature consistency and curve
detection," in IEEE Trans. on Pattern Analysis and Machine Intelligence, Vol.
PAMI-11, pp. 823-839, 1989.
[32] A. P. Pentland, "Parts: Structured descriptions of shape," in Proceedings of
AAAI, Philadelphia, PA, pp. 695-701, 1986.
[33] P. Perona and J. Malik, "Scale-space and edge detection using anisotropic diffusion," IEEE Trans. on Pattern Analysis and Machine Intelligence, Vol. 12, No.
7, pp. 629-639, 1990.
[34] W. Press, B. Flannery, S. Teukolsky, and W. Vetterling, Numerical Recipes in
C, Cambridge University Press, Cambridge, 1988.
[35] R. Samadani, "Changes in connectivity in active contour models," Proceedings
of the Workshop on Visual Motion, pp. 337-343, Irvine, California, March 1989.
[36] R. Samadani, "Adaptive snakes: control of damping and material parameters,"
Proceedings of SPIE Conference on Geometric Methods in Computer Vision,
Vol. 1570, pp. 202-213, San Diego, California, July 1991.
[37] L.L. Schumaker, "Fitting Surfaces to Scattered data," in Approximation Theory
II, G.G. Lorentz, C.K. Chui, and L.L. Schumaker, (eds.), pp. 203-267, Academic
Press, New York, 1976.
[38] J. A. Sethian, "Curvature and the evolution of fronts," Commun. in Mathematical Physics, Vol. 101, pp. 487-499, 1985.
[39] J. A. Sethian, "Numerical algorithms for propagating interfaces: HamiltonJacobi equations and conservation laws," Journal of Differential Geometry, Vol.
31, pp. 131-161, 1990.
[40] R. Szeliski and D. Tonnesen, "Surface modeling with oriented particle systems,"
Computer Graphics SIGGRAPH, Vol. 26, No. 2, pp. 185-194, July 1992.
[41] D. Terzopoulos, "Regularization of inverse visual problems involving discontinuities," IEEE Trans. on Pattern Analysis and Machine Intelligence, Vol. PAMI
8, No. 2, pp. 413-424, 1986.
[42] D. Terzopoulos, "The computation of visible surface representations," IEEE
Trans. on Pattern Analysis and Machine Intelligence, Vol. PAMI 4, Vol. 10,
pp. 417-438, 1988.
[43] D. Terzopoulos, "Matching deformable models to images: Direct and iterative
solutions," Proc of Optical Society of America, Tropical Meeting on Machine Vision, Incline Village, NV, March 18-20, 1987.











[44] D. Terzopoulos, A. Witkin, and M. Kass, "Symmetry seeking models and 3D
object reconstruction," Int. Journal of Computer Vision, Vol. 1, No. 3, pp. 211221, 1987.

[45] D. Terzopoulos, A. Witkin, and M. Kass, "Constraints on deformable models:
Recovering 3D shape and nonrigid motion," Artificial Intelligence, 36, pp. 91123, 1988.
[46] 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. 724-725, Maui Hawaii, June
1991.
[47] B. C. Vemuri and R. Malladi, "Surface griding with intrinsic parameters," Pattern Recognition Letters, Vol. 13, No. 11, pp. 805-812, November 1992.
[48] B. C. Vemuri and R. Malladi, "Intrinsic parameters for surface representation
using deformable models," IEEE Trans. on Systems, Man, & Cybernetics, Vol.
23, No. 2, pp. 614-623, March/April 1993.
[49] B. C. Vemuri and R. Malladi, "Constructing intrinsic parameters with active
models for invariant surface reconstruction," IEEE Trans. on Pattern Analysis
and Machine Intelligence, Vol. 15, No. 7, pp. 668-681, July 1993.
[50] B. C. Vemuri, A. Mitiche, and J. K. Aggarwal, "Curvature-based representation
of objects from range data," Int. Journal of Image and Vision Computing, 4,
pp. 107-114, 1986.
[51] Y. F. Wang and J. F. Wang, "Surface reconstruction using deformable models
with interior and boundary constraints," in Proceedings of ICCV, pp. 300-303,
Osaka, Japan, 1990.
[52] D. M. Young, Iterative Solution of Large Linear Systems, Academic Press, New
York, 1971.














BIOGRAPHICAL SKETCH


Ravikanth Malladi was born in Vijayawada, India, in 1966. He received his B.Eng. degree with honors in electrical engineering and M.Sc. degree in physics from Birla Institute of Technology and Science, India, in 1988. In summer of 1985 he worked at the National Geophysical Research Institute, India. During spring of 1988 he held a research assistant position with the semiconductor physics group in Central Electrical and Electronics Research Institute, India, where he was involved in modeling the exposure and development of positive photoresist. He subsequently received the M.S. degree in 1991 and expects to receive the Ph.D. degree in computer and information sciences from the University of Florida, Gainesville, in 1993. Since the summer of 1990 he was a graduate research assistant at the University of Florida, and was affiliated to the Center for Computer Vision and Visualization. His research interests are focused on computational vision, shape modeling, surface reconstruction, computer graphics, image processing, and medical image interpretation.













I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy.


aba '-. Ve-muri, Chairman
Associate Professor of Computer and Information Sciences

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy.


James A. Sethian, Cochairman
Professor of Mathematics
University of California at Berkeley

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of hilosophy.



ofessor of Computer and
Information Sciences

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy.



Andrew F. Laine
Assistant Professor of Computer and Information Sciences













I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for t, degree of Doctor o Philosophy.


Padl A. Fishwick
Associate Professor of Computer and Information Sciences


I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy.


Murali Rao
Professor of Mathematics



This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degreor of Philosophy.

December 1993 ___ _ _ _,_ /,4 Winfred M. Phillips
Dean, College of Engineering



Karen A. Holbrook
Dean, Graduate School




Full Text

PAGE 1

A TOPOLOGY-INDEPENDENT SHAPE MODELING SCHEME By RAVIKANTH MALLADI A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1993

PAGE 2

Copyright 1993 by Ravikanth Malladi

PAGE 3

To my parents, Pardhasaradhi and Vasanthavalli Malladi, my sister Lakshmi, and my uncle, Chandrasekhar Nori, who have always encouraged my academic aspirations.

PAGE 4

ACKNOWLEDGEMENTS First, I would like to express my gratitude to my advisor, Dr. Baba C. Vemuri, for his continuous encouragement, for providing financial support, for numerous discussions on research topics, and guidance throughout the course of this work. I am grateful to Dr. James A. Sethian for consenting to serve on my supervisory committee, it was his involvement and insightful suggestions that made progress possible. I thank him for being generous with his time and providing me with more ideas than I could implement. I would also like to thank the other members of my committee, Drs. Gerhard X. Ritter, Andrew F. Laine, Paul A. Fishwick, and Murali Rao, for their help and constructive criticism. Special thanks are due to Randy Fischer for his help in acquainting me with the video equipment and the KSR machine. I owe a great deal to my friend Niranjan Mayya for his constant support and understanding throughout the past five years. An early acquaintance with Mark Burch has grown into a strong friendship which I shall always cherish. I thank my other “old” friends Amitabh Mungale, Pamela Rolfe, and Vidhya Gangar for making my life in Gainesville a very enjoyable one. More recently, in Roxana Dastur, Eman Anwar, Pankaj Agarwal, and Chitrank Sheshadri, I found excellent friends with whom I have enjoyed numerous lively conversations. Financial support for this work was provided in part by the Department of Radiology and the Department of Computer and Information Sciences at the University of Florida. IV

PAGE 5

TABLE OF CONTENTS ACKNOWLEDGEMENTS iv ABSTRACT vii CHAPTERS 1 INTRODUCTION 1 1.1 Overview 5 1.2 Dissertation Outline 7 2 ENERGY METHODS FOR SHAPE MODELING 9 2.1 Symmetry-Seeking Deformable Models 9 2.1.1 The deformable model 10 2.1.2 Extrinsic forces 13 2.2 Snake: The Active Contour Model 14 2.2.1 Snake dynamics 15 3 FRONT PROPAGATION PROBLEM 17 3.1 Lagrangian Formulation 17 3.2 Level Set Formulation 18 4 SHAPE RECOVERY WITH FRONT PROPAGATION 22 4.1 Stopping Criterion: An Image-based Speed Term 22 4.2 Extension of Image-based Speed Term 24 4.3 Shape Recovery in 3D 27 5 NUMERICAL SOLUTION 29 5.1 Entropy-satisfying Upwind Numerical Scheme 30 5.2 Narrow-band Updation Scheme 33 5.3 Numerical Issues in 3D 40 6 EXPERIMENTAL RESULTS 43 6.1 2D Shape Recovery Results 44 6.2 Experimental Results in 3D 50 v

PAGE 6

7 CONCLUDING REMARKS 62 REFERENCES 64 BIOGRAPHICAL SKETCH 68 VI

PAGE 7

Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy A TOPOLOGY-INDEPENDENT SHAPE MODELING SCHEME By Ravikanth Malladi December 1993 Chairman: Dr. Baba C. Vemuri Cochairman: Dr. James A. Sethian Major Department: Computer and Information Sciences 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 dissertation 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 interlace (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-Jacobi” 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. vii

PAGE 8

1 he resulting equation of motion is solved by employing entropy-satisfying upwind finite difference schemes. We also introduce a new algorithm for rapid advancement of the front using what we call a narrow-band updation scheme. This leads to significant improvement in the time complexity of the shape recovery procedure in 2D. 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. We also demonstrate the recovery of 3D shapes. viii

PAGE 9

CHAPTER 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 model-based 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 be either rigid or nonrigid. In the former case, we have what are popularly known as correlation-based template matching techniques, while the later involves a dynamic model fitting process to the data. Shape models serve the purpose of an intermediate stage in object recognition tasks, since they provide a more stable and useful description than the original raw data. In this dissertation, we present a new dynamic approach to shape modeling which retains the most attractive features of existing methods, and overcomes their prominent limitations. For the rest of this dissertation 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. The shape models must aid in recovery of detailed structure from noisy data using only the weakest of the possible assumptions about the observed shape. Several variational shape reconstruction methods have been proposed and there is abundant literature on the same [5, 37, 42, 6, 50, 22]. Generalized spline models with continuity constraints are well suited for fulfilling the goals of shape recovery (see [6, 7, 41]) Generalized splines are the key ingredient of the dynamic shape modeling paradigm introduced by Ierzopoulos et al ., [45]. Incorporating dynamics into shape modeling 1

PAGE 10

2 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 [17, 46, 47, 49, 51, 10, 12]. 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 “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 dissertation, 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 ot the shape initialization. Our method can also recover shapes whose topology changes over time, e.g., the cell boundary in cell division [35], The framework of energy minimization has also been used successfully in the past for extracting salient image contours edges and lines. Kass et al. [17] used energyminimizing snakes that exhibit an affinity toward image features such as edges points and edge segments, while maintaining piecewise smoothness. The weights of the smoothness and image 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 surface models to

PAGE 11

3 P A (a) CT image (b) DSA image (c) Holed shapes Figure 1.1. Test bed for our topology-independent shape modeling scheme. image data by means of energy minimization [45]. The scheme seeks to design energy functionals whose local minima comprise the set of alternative solutions available to high-level processes. In the absence of a well-developed high-level 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 [9] defines an inflation force on the active contour. This 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.1(b). Moreover, despite a good

PAGE 12

4 initialization, the active contour model, due to its arc-length 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 ID thin-plate-membrane-spline, in an adaptive environment wherein the material parameters controlling the relative strengths of elasticity and rigidity are adapted (see [36]). 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. Terzopoulos [45] has 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 “ sprout ” branches during the evolutionary process. Once the shape has been segmented from the image, its constituent part structure can be inferred using the algorithm described in [18]. 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 [35]. Returning to

PAGE 13

5 static domain, more recently molecular dynamics has been used to model surfaces of arbitrary topology [40]. Smoothness and continuity constraints are imposed by subjecting a 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 dissertation can be applied 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 shape of interest (see figure 1.1(c)), has the ability to split freely to represent each shape [25, 26]. 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 6.9). Finally, we apply our method to 3D images and recover the surface shapes of 3D objects (see figure (6.12)). 1.1 Overview In this section we briefly outline our scheme to model complex shapes. This method is inspired by ideas first introduced in Osher and Sethian [30, 39] 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 [38] 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

PAGE 14

6 [30, 39] achieve this by viewing the propagating surface as a specific level set of a higher-dimensional function. The equation of motion for this function is reminiscent of an initial valued HamiltonJacobi equation with a parabolic right-hand 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 I< is the curvature of the hypersurface. As in [30], we adopt a global approach and view the (TV — 1) dimensional moving surface as a level set of a timedependent function xp oi N space dimensions. The equations of motion written for this higher dimensional function are then amenable to stable entropy-satisfying numerical schemes designed to approximate hyperbolic conservation laws. Topological changes can be handled naturally in this approach, since a particular level set {xp = 0} of the function xp 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 lequired that we stop the hypersurface in the neighborhood of the desired shape. We do this by synthesizing an negative speed function from the image. Secondly, we have to construct an extension of this speed function to other level sets {xp = C} in order to give a consistent meaning to the image-based speed term at all points in the image (see figure 3.1). 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 [30, 38, 39], has been applied in the area of computer vision for shape characterization by Kimia et al. [18, 19], who unify many diverse aspects of shape by defining a continuum of shapes (reaction/diffusion space), which places shapes within

PAGE 15

7 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 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 if and recover complex shapes by propagating it along its gradient field. Shape characterization can be done once the object shape is extracted. An important contribution of this dissertation over our previous work in [26] is the design and implementation of a faster numerical technique for solving the governing equation in 2D Our new algorithm exploits the fact that the front, which is a particular level set {ip — 0} of a higher dimensional function, can be advanced by updating the function ip at a small set of points. This scheme is an alternative to updating the ip 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 {ip — 0}. Since the narrow-band 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 narrow-band techniques for interface propagation may be found in [1], We also show that the surface shapes of objects can be recovered from 3D image data. 1.2 Dissertation Outline 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

PAGE 16

8 for high-level processing. The remainder of this dissertation is organized as follows: chapter 2 introduces the curvature-dependent front propagation problem and establishes a link between HamiltonJacobi equations and a hyperbolic conservation law. In chapter 3 we explain our level set algorithm for shape recovery and in chapter 4 we outline the details of our narrow-band approach. Finally, in chapter 5 we present some experimental results of applying our method to some noisy and low contrast medical images in 2D as well as 3D. We close with a discussion of advantages of our approach and direction of future research in chapter 6.

PAGE 17

CHAPTER 2 ENERGY METHODS FOR SHAPE MODELING In this chapter, we review the some of the literature on existing shape modeling methods in computer vision. Specifically, we provide the details of two “energybased methods for shape modeling and reconstruction. The first is the symmetryseeking deformable model [45] and the second is a shape recovery scheme in 2D called active contour model or “snake” [17]. Existing shape modeling methods roughly fall under two categories: passive and active models. Models of shape such as generalized cylinders, introduced by Binford [4], an( ^ lumped-parameter family of shapes such as superquadric models [3, 32] are purely geometric, hence passive. Generalized cylinders are used to model elongated shapes with axial symmetry, while the superquadric shape models are well suited for object recognition tasks because one can express them compactly using a small set of parameters. On the contrary, active models are governed by the principles of Lagrangian dynamics and react to external forces in a way that is similar to real elastic objects. Generalized spline models with continuity constraints [7, 10, 12, 41, 44, 49, ol] are pume examples of the active shape modeling paradigm. 2J — Symmetry-Seeking Deformable MoHpN The deformable model is a physically based modeling framework for shape and motion reconstruction of free-form flexible objects. In this framework the objects are modeled as elastically deformable bodies subject to continuum mechanical laws. Constiaints are expressed as forces which deform the elastic models and propel them 9

PAGE 18

10 through potentially complicated motions such that they satisfy the available constraints over time. There are two types of forces, intrinsic and extrinsic. Intrinsic constraints reflect generically valid assumptions about natural objects while the extrinsic constraints reflect observations about the environment that can be extracted from sensory data. The model is embedded in a force field which encodes the available sensory information. The ambient force field will then mold the deformable model to make their 3D shape consistent with the observed data. A plethora of external force fields can be synthesized from the data using a wide variety of constraints. Contrary to the conventional models of 3D shape, the deformable models are active in that they react to extrinsic forces as one would expect real elastic objects to react to physical forces. This is because the deformable models are governed by the principles of elasticity theory as expressed through Lagrangian dynamics. Specifically, the continuum mechanical equation <9 2 v dt 2 4. ^(v) dt + 7~H7 = f( v ) <5v ( 2 . 1 ) governs the nonrigid motion of the deformable model in response to a net extrinsic force f(v), where n is the mass density function of the deformable body and 7 is the viscosity function of the ambient medium. 2.1.1 The deformable model We begin this subsection by providing an informal description of the deformable model. Consider a deformable sheet made of elastic material (a membrane-thin-plate hybrid). This sheet is rolled to form a tube. Next, a deformable spine made of similar material is passed down the length of the tube. At regularly spaced points along the spine, it is coupled to the tube with radially projecting forces so as to maintain the spine in approximate axial position within the tube. Additional forces are included that coerce the tube into a quasi-symmetric shape around the spine. Finally, extra

PAGE 19

11 control is provided over the shape by introducing expansion/contraction forces radiating from the spine. The rigidities of the spine and the tube are independently controllable, and their natural rest metrics and curvatures can either be prescribed or modified dynamically. If the circumferential metric of the tube is set to zero, for instance, the tube will tend to contract around the spine, unless the other forces prevail. The model will shorten or lengthen as the longitudinal metrics of the tube and the spine are modified. In short, a variety of interesting behaviors can be obtained by adjusting the control variables designed into the model. The deformable tube and the spine can be described by geometric mappings from material coordinate domains into Euclidean three-space !ft 3 . The mapping is expressed by position vectors in space whose component functions are time varying. The spine is a deformable space curve defined by mapping a univariate material coordinate domain s € [0, 1] into 3J 3 v : [0, 1] x [0, oo] 3f? 3 such that v(s, t) = (X(s, t), Y(s, t ), Z(s, t)) (2.2) The tube is constructed from a deformable sheet that is defined by mapping a bivariate material coordinate domain (u, v ) 6 [0, l] 2 into 3£ 3 v : [0, l] 2 x [0, oo] -> 5R 3 such that v(t x, u, t) = (X(u, v, t), Y(u , v, t), Z(u, v, t)) (2.3) The tubular shape is created by specifying boundary conditions on two opposite edges of the sheet which effectively glue these edges together. The edges u = 0 and u = 1 are “seamed” together, letting v span the length of the tube. The required periodic boundary conditions are v (0, v,t) = v(l, v, t), &v . dv . du ih, ' (2.4) (2.5)

PAGE 20

12 The strain energy £ s associated with the spine function v(s,t) is given by where the vertical bars denote the Euclidean norms and the quantities w 1 ,w 2 are both functions of (s,t). uq(s,i) determines the local tension along the spine and w 2 (s,t) determines its local rigidity. The weighting functions uq, w 2 control the elastic properties of the spine. We can regulate the equilibrium shape of the spine by suitably defining the weighting functions. The strain energy £ associated with the tube function v(u,u,t) is given by The weighting functions uqo and u>oi define the natural metric of the tube along each parameter curve. These functions locally control the tension of the deformable sheet (that constitutes the tube) along each material coordinate curve. Analogous expressions for w 20 (u, v), w n (u, v) and w 02 (u,v) determining the natural curvature of the sheet can be written down. These functions locally control the rigidity of the sheet. The continuity of the deformable spine can be controlled via the functions uq and w 2 . Position discontinuity is introduced by setting «q( s ,£) = w 2 (s,t) — 0 and tangent discontinuity by setting w 2 (s,t) = 0. Similarly, we can introduce discontinuities in position as well as tangents for the deformable tube. The spine and the tube functions are coupled by first identifying that v = s, bringing into correspondence the spine coordinate with the coordinate along the length of the tube. The mapping function of the spine v s (v,t) is then distinguished from that of the tube v T (u,v,t) with superscripts S and T. (2.6) 2wn I I 2 +^02 (2.7)

PAGE 21

13 Summing all the forces into the equations of motion (2.1) associated with the spine and tube, the following dynamic system describing the motion of the deformable model can be obtained: d 2 v s dv s 8S s 8V S ' dt 2 + t 3( 8w s 8\ s d 2 v T dv T 8f T 8V t dt 2 + 7 m + 8v T ~ ~ 8v T where f s and f T are the coupling forces between the + f S , (2.8) + f T , (2.9) spine and tube that force the spine in a central position within the tube and coerce the tube to radial symmetry around the spine. V s (v s ) and V T {y T ) are the generalized potential functions associated with the spine and tube respectively. The extrinsic forces on the right hand side of equation are the variational derivatives of the generalized potential functions V and V T that are synthesized from the data constraints. 2.1.2 Extrinsic forces In this subsection, we follow the argument in [49] and describe the extrinsic forces that are applied to the deformable model to obtain conformation to the range data. Basically, fitting a deformable model to range data is equivalent to applying an extrinsic force to the model until the model achieves maximal conformation with the data. The extrinsic forces, as mentioned earlier, can be defined as variational derivatives of potential functions. Techniques for generating suitable potential functions from monocular, binocular and dynamic image sequences are described in [45], The extrinsic force on the tube due to the depth data is defined by the variational derivative of the potential function, namely ^J(v T ) = ^5I a «*ll vT d l| 2 » (2.10) z n where Q specifies the data domain, v T denotes the vector function of the tube, and d denotes the data vector from the range finder. a d is the inverse of the variance <7 a

PAGE 22

14 along the vector a defined from the data point to the point of influence on the model. The force associated with the above potential is given by fJ(v T ) = a d (v r d) ( 2 . 11 ) Assuming no external forces on the spine i.e., the potential function V s = 0 for the spine. The dynamic equations that are solved for the surface fitting problem are then given by, d 2 v T dv T 88 T = ad(v T d) + f T , d 2 v s dv s 88 s '‘‘si 5 ' + 1 ~dt + W = f s . ( 2 . 12 ) (2.13) The differential equations (2.12) and (2.13) pose a nonlinear initial boundary value problem. A solution is obtained by discretizing in spatial and time domains. The discretization in spatial domain is carried out using standard finite central difference formulas on regular grid of nodes [20], Since the extrinsic forces are obtained from a static field, the integration in time is performed until qq the viscous damping term, dissipates all the kinetic energy in the model thereby bringing it to a static equilibrium. A detailed discussion of various numerical procedures that can be employed for solving the above system of differential equation can be found in Terzopoulos et al. [43], 2.2 Snake: The Active Contour Model Snakes are an example of a more general technique of matching a deformable model to an image by means of energy-minimization. The idea of using energyminimizing splines to extract image features of interest has been introduced by Kass et al. [17]. Snake is guided by external image forces that guide it toward features such as lines and edges.

PAGE 23

15 The snake is a model of a deformable contour composed of abstract elastic materials. Two types of materials determine the snake configuration: string and rod. The former makes the snake resistant to stretching and the latter makes it resistant to bending. This deformable curve is then embedded in a potential functional synthesized from a variety of image constraints. Depending on the nature of the external potential field, the snake will extract different image feature. Unlike most other techniques for finding salient image features, the snake model is active. It is always minimizing its energy functional and therefore exhibits dynamic behavior. 2.2.1 Snake dynamics The basic snake model is a controlledcontinuity spline [41] under the influence of image forces and external constraint forces. The internal spline forces serve to impose a piecewise smoothness constraint and the external constraint forces are responsible for putting the snake near the desired local minimum. Consider a deformable curve v(s,f) with parameters s and t (time), defined on given intervals 0 = [0..1], and T, respectively. Let this deformable curve be a function of two coordinate variables x and y with the same parameterization: v (a, i) = {x(s, t), y(s , t)) : s € ft, < eT. (2.14) The potential energy functional of the snake E S nake{ v) is defined as E snake ( v) = f ( E{ n t -)E ex t)ds , (2.15) where E ext is the potential associated with external constraint forces. Ei n % denotes the internal potential functional of the snake and is given by: Ei nt(v) = f Wx{s) I v, I 2 +u; 2 (s) I v ss I 2 ds. (2.16) J The first-order term w^s) | v s | 2 is responsible for C° continuity in the snake and the second-order term w 2 (s) | v ss | 2 enforces C l continuity. The weight uq(.s) regulates

PAGE 24

16 the tension of the snake, whereas w 2 {s) regulates its rigidity. Position and tangent discontinuities may be introduced along a snake by setting these weights to zero. Shape recovery in 2D can be accomplished by defining a simple energy functional from a given image. If we set E ext (v) = [ | VG> * I(x,y) | 2 , (2.17) J 17 then the snake is attracted by the local minima of the potential, which corresponds to the local maxima of the image gradient. Here, G a * I(x,y) denotes the image function convolved with a Gaussian smoothing filter whose characteristic width is a. Potential functionals modeling other external constraints defined by the user may be added to the above expression. Given the total potential energy of the snake (see equation (2.15)) for a specific initial position, the equilibrium configuration can be reached by minimizing it. Specifically, if the function v* represents a local minimum for E sna k e , it satisfies the the following EulerLagrange equation: d d 2 7V * ~ = -VLexi(v’). (2.18) In this formulation, each term appears as a force applied on the curve. Given the value of v(s, t — 0), the initial snake configuration, the above equation can be solved numerically by adopting direct or iterative finite difference schemes [17, 20].

PAGE 25

CHAPTER 3 FRONT PROPAGATION PROBLEM In this chapter, we introduce the front propagation problem. As a starting point and motivation for the level set approach, we first discuss the interface motion in terms of a parameterization of the moving front. For details and an expository review, see Sethian [38, 39]. The numerical schemes based on discrete approximations to these equations are shown to have limitations. We then present the level set technique due to Osher and Sethian [30], which overcomes the limitations inherent in the former Consider a closed curve moving in the plane, that is, let 7 (0) be a smooth, closed initial curve in Euclidean plane R 2 , and let 7 (<) be the one-parameter 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,/) be the position vector which parameterizes 7 (t) by s, 0 < s < S. The curve is parameterized so that the interior is to the left of the direction of increasing s. With n and K{s, t ) as the outward normal and curvature of the curve respectively, the equality n • x t = F(K) holds, d he equation of motion in terms of individual components of x = (x,y) is This is called a “Lagrangian” representation because the physical coordinate system moves with the front. approach. 3T Lagrangian Formulation (3.1) 17

PAGE 26

18 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,*), 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. However, there are several problems with this approach, as discussed in Sethian [38]. 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 dimensions. 3.2 Level Set Formulation As an alternative, the central idea in the level set approach of Osher and Sethian [30] is to represent the front 7 (Z) as the level set {ij, = 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 zy-plane (figure 3.1(a)). We imagine that the circle is the level set {0 = 0 } of an initial surface z = = 0 ) in R 3 (see figure 3.1(b)). We can then match the one-parameter family of moving curves 7 (i) with a one-parameter family of moving surfaces in such a way that the level set {*/> = 0 } always yields the moving front (figures 3 . 1 (c) & 3 . 1 (d)). In the general case, let 7 (0) be a closed, nonintersecting, (N — 1 ) dimensional hypersurface. Let x £ R N , be the scalar function such that 0(x,O) = ±d(x), where d(x) is the signed distance from x to the hypersurface 7 (0). We use the plus

PAGE 27

19 Figure 3.1. Level set formulation of equations of motion (a) & (b) show the curve 7 and the surface ip{x, y) at t — 0, and (c) & (d) show the curve 7 and the corresponding surface ip(x,y) at time t. sign if x is outside 7(0) and minus sign if x is inside. Each level set of ip flows along its gradient field with speed F(K). The gradient S/ip{x,t) is normal to the (N — 1) dimensional level set passing through x. Now, we derive the equation of motion for function ip. Consider the motion of some level set {ip = C}. Following the derivation in [30], let x(t) be the trajectory of a particle located on this level set, so V’(x(t),f) = C. (3.2) The particle speed dx/dt in the direction n normal to 7 (t) is given by the speed function F. Thus, dx ~dt' n = ( 3 3 ) where the normal vector n is given by n = Vipf | \7ip |. By the chain rule we get, . dx + -57 • Vt/^ = 0 ot (3.4)

PAGE 28

20 and substitution yields i/>t + F | V0 |= 0, ( 3 . 5 ) with an initial condition 0(x,O) = ±d(x). We refer to equation (3.5) as the level set “HamiltonJacobi” formulation. Note that at any time, the moving front 7 (f) is simply the level set {0(x,f) = 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 m the front can be handled naturally by exploiting the property that the level surface {if, = 0} need not be simply connected. The signed distance function *(*,*) alwa * s remains a function, even if the level surface {0 = 0 } corresponding to the front 7 (f) changes topology, or forms sharp corners. The geometric and differential properties of 7 (i) are captured in the function $ and can be readily extracted. As an example, if x 6 i? 2 , the curvature is given by UT _ (M l ~ tyxMxy + i>xAl) (V> 2 + V^y)3/ 2 ' ( 3 ‘ 6 ) This approach can also be easily extended to higher dimensions and appropriate expressions can be obtained for the mean curvature and the Gaussian curvature [30]. For example, the mean curvature at a point (x,y,z) on a moving surface can be computed using the following formula: K = + ^vvi ^l + $ 1 ) + + V> 2 ) 2i\) xy ^ x ij, y 2 ip yz ip y ip z ( 0 2 + ^ 2 + 0 2 ) 3/2 " By substituting F(K) = l-eK as a typical speed function in equation (3.5), the equation of motion becomes A+ I V0 |= el< | | (3.7)

PAGE 29

21 Equation (3.7) resembles a Hamilton-Jacobi equation with viscosity, where “viscosity” refers to the second-order parabolic right-hand side. This equation can be solved using the stable, entropy-satisfying finite difference schemes, borrowed from the literature on hyperbolic conservation laws (see [30]).

PAGE 30

CHAPTER 4 SHAPE RECOVERY WITH FRONT PROPAGATION In this chapter, we describe how the level set formulation for the front propagation problem discussed in the previous chapter 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. — Stopping Criteri on: An Image-based Speed Term 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 + Fq. 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 F a depending on its sign and is analogous to the inflation 22

PAGE 31

23 force defined in [9]. The second term Fa, 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 thin-plate-membrane splines [17] (see the figure (6.3)). We rewrite equation (3.7) by splitting the influence of F as + Fa | V0 | +Fq | Vip |= 0. (4.1) First consider the case when the front moves with a constant speed, i.e. F = F A . 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 Fj to be F/(x,y) = Fa (Mi M 2 ) VGV * I(x,y) j -M 2 } , (4.2) where M x a id M 2 are the maximum and minimum values of the magnitude of image gradient | VG<, * I(x,y) |, (x,y) £ ft. The expression G a * I denotes the image convolved with a Gaussian smoothing filter whose characteristic width is a. Alternately, we could use a smoothed zero-crossing image to synthesize the negative speed function. The zero-crossing image is produced by detecting zerocrossings in the function V 2 G a * /, which is the original image convolved with a Laplacian-of-Gaussian filter whose characteristic width is a. The equation of motion with the addition of image-based speed is Vt + {Fa + Fj) j V?/’ |= 0. (4.3) I * s ca ^ e d an extension of b j to points away from the boundary 7(f), i.e. at points [ x iij) £ (0 — 7 (0), an d is equal to Fj on 7 (f). We shall return to the issue of

PAGE 32

24 extension shortly. The value of F 7 lies in the range [-i^,0] as the value of image giadient varies between Mi and M 2 . 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. F g # 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 = F A + F G with a quantity kj. The term kj , which is defined as 1+ | VG a *I(x,y) [’ ( 4 * 4 ) 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 “steerable” filters [16]. The modified equation of motion is given by 0 « + kj(F A + F g ) I V0 ]= 0. (4_5) M2 — Extension of I mage-based Spe^d Term We now come to an important juncture in our discussion. The image-based speed term, be it F r or k h has meaning only on the boundary 7 (f), i.e. on the level set {0 = 0}. This follows from the fact that they were designed to force the propagating level set {0 = 0} to a complete stop in the neighborhood of an object boundary. However, the equation of motion (4.3) is written for the function 0, which is made up of infinitely many level curves. In other words, equations (4.3) k (4.5) control the evolution of a family of level sets. Therefore, it is imperative that the net speed used m the evolution equation has a consistent physical meaning for all the level

PAGE 33

25 y Figure 4.1. Huygen’s principle construction sets, i.e. at every point (x,y) € H. Speed functions such as F G which are functions of geometric properties of the surface z = iP(x,y), can be readily computed at any (^by) £ H. However, Fj is not such a function. It derives its meaning not from the geometry of ip but from the configuration of the level set {ip = 0} in the image plane. Thus, our goal is to construct an image-based speed function Fj that is globally defined. We call it an extension of Fj off the level set {ip = 0} because it extends the meaning of Fj to other level sets. Note that the level set {0 = 0} lies in the image plane and therefore Fj must equal Fj on {0 = 0}. The same argument applies to the coefficient kj. If the level curves are moving with a constant speed, i.e. F a = 0, then at any time f, a typical level set {ip = C}, C G R, is a distance C away from the level set { = 0} (see figure 4.1). 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 F g £ 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

PAGE 34

26 Figure 4.2. Extension of image-based speed terms to other level sets function F = F A + F G is continuous (see [14]). With the above remarks in mind we state the following: Pr operty 4.2.1 Any external (image-based) speed function that is used in the equation of motion written for the function ii> should not cause the level sets to collide and cross each other during the evolutionary process. Recall that the function 0(x,<) has been initialized to d(x), where d(x) is the signed distance from a point x e fUo the boundary 7 (0). Since we cannot attribute any geometric meaning to the function Fj (kf) at points away from the level set {^ = 0}, we look for a meaning that is consistent with property (4.2.1). Therefore, the question to ask is: what is the value of F, (or k f ) at a point (x,y) lying on a level set {if = Cyi We answer this question in the following construction (see figure 4.2). Construct io n 4.2.1 The value of F, (kj) at a point P lying on a level set {if} = C } zs exactly the value of F I (kj) at a point. Q, such that point Q is a distance C away from P and lies on the level set {if = 0}. It is easy to see that Fj reduces to F, on {if = 0}. We use the same construction to determine the value ol ki at a point (x,y) lying on some level set {if = C }.

PAGE 35

27 Note that if the definition of a speed function adheres to construction 4.2.1, then it will also be consistent with the property 4.2.1. Thus, having ascribed the intent of pseudodifferential equations (4.3) & (4.5) in the context of shape modeling, we can use finite difference schemes to solve them numerically. Since ip 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.3 Shape Recovery in 3D In this section, we show that the surface shapes of three-dimensional objects can be recovered using a scheme that is analogous to the one described above. We begin by writing the equation of motion for a moving surface; this represents front propagation in higher dimension. Let the surface be closed, nonintersecting, and the level set {ip = 0} of a function 0(x, t), x £ R 3 . The function ip(x,t = 0) = ±d(x), where d(x) is the signed distance function from the point x to the initial surface; plus sign is used if x is outside the initial surface and minus sign if x is inside. Each level set of ip moves along its gradient field with speed Following the derivation in chapter 2, we arrive at the equation of motion for the moving surface, namely ipt + F(K) | S7ip |= 0, (4.6) where I< is either mean or gaussian curvature of the surface. Recall that the surface at any time t is simply the level set {ip — 0}. In this work, we only address the issue of recovering surface shapes from 3D images. The input data constitutes a set of MRI (Magnetic Resonance Imaging) or CT (Computed Tomography) images. A stack of such images represents an image function 7 (x,?/, 2;). Our next step is to define an image-based speed term which will

PAGE 36

28 cause the level set {0 = 0) to stop near the 3D object boundary. To this end, we apply on the front the speed equal to the product of F(K) and a term £/. The term ki, which is defined as to unity in regions of relatively constant image intensity. To attribute a consistent meaning to the image-based speed term at all the level surfaces, we construct an extension to kj in a manner similar to the construction (4.2.1). Therefore, the modified equation of motion that we solve to extract complex surface shapes from 3D images is (4.7) has values closer to zero in the vicinity of object boundaries and values that are closer i’t + kjF(K ) | V?/> |= 0. (4.8)

PAGE 37

CHAPTER 5 NUMERICAL SOLUTION The equation (3.7) poses an initial valued problem. It is rewritten here as with xf>(x,y,t = 0) = ± distance from (x,y) to 7(0). As shown in Sethian [38], for t > 0, the parabolic right-hand side diffuses sharp gradients and forces to stay smooth at all values of t. For £ = 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 “weak” 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 [38]: 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 [30, 39] are motivated by the fact that the entropy condition for propagating fronts is identical to the one for hyperbolic conservation laws, where stable, consistent, entropy-satisfying algorithms have a rich history. (5.1) 29

PAGE 38

30 Sil — Entropy-satisf ying Upwind Numerical Scheme In this section, almost without a change, we present the arguments discussed in Sethian [39]. For complete details of the following scheme, we refer the reader to Osher and Sethian [30, 39 ], First, consider the one-dimensional version of the level set equation, with constant normal velocity F A = 1 . Then the equation (3.7) becomes a standard HamiltonJacobi equation t + H(4: r) = 0, (5.2) where H( x ) = — (<^|) 1/2 , and with a given initial value of . Let u = 4> x . Taking the derivative with respect to x, equation ( 5 . 2 ) becomes u t + [H(u)] x = 0 , ( 5 . 3 ) where H(u ) = — (u 2 ) 1/2 . Equation (5.3) 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 (5.3) are integrated in an arbitrary interval [a, b] to produce d r b u{x , t)dt = H[u{a, <)] H[u(b, t)]. ( 5 . 4 ) 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, entropy-satisfying solution to equation (5.4)? The answer is found in the following definition:

PAGE 39

31 Definition 5.1.1 Let u" 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 itself. First write equation (5.2) with a forward difference in time as, A tH(u). ( 5 . 6 ) Since the numerical flux function g approximates H, the solution to equation (5.6) may be approximated by = ? &tg(ui_ 1/2 ,u i+1/2 ) = W ~ ^tg(D-i,D+i), ( 5 . 7 ) where g is an appropriate numerical flux function and the standard definitions of the forward and the backward difference operators have been used, namely, « DUi = (5.8)

PAGE 40

32 Finally, an appropriate numerical flux function g is required. In the special case where H(u) may be written as a function of u\ i.e., H{u) = f{u 2 ) for some function /, one can use the Hamilton-Jacobi flux function given in [30]: 9{ u i-i/2, Ui+ 1 / 2 ) = g{D~ ,-, D+ i) = /((max(/};^,0)) 2 + (min(£»+^,0)) 2 ). (5.9) 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 finite-sized 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 F A = I, so that f(u ) = — (u 2 ) 1 / 2 , equation (5.7) reduces to $* +1 = ? ~ At{(max(Z)J
PAGE 41

33 The second term F G k I | Vi/> | is not approximated in the above equation; one may use a straightforward central difference approximation to this term. A2 — Narrow-band Updation Scheme In this section we first consider a straightforward approach to solve the equations (5.11) & (5.13) 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 {ip = 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 {V> = 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 image-based speed term on the zero set (equation (4.5)). An extension function is constructed to attribute a consistent meaning to the image-based speed term at points away from the zero set (see discussion in section 3). Therefore, a straightforward algorithm consists of advancing from one time step to the next as follows: Algorithm 1 1. At each grid point (iAxJAy), where Ax and Ay are step sizes in either coordinate directions, the extension of image-based speed term is computed. This is done m accordance with the construction described in previous section; i.e., by searching for a point q which lies on the level set {ip = 0}, and is closest to the point (iAxJAy). The value of image-based speed term at the current point is simply its value at the pcint q.

PAGE 42

34 2. With the value of extended speed term and 0£, calculate using the upwind, finite difference scheme given in equation (5.13). 3. Construct an approximation for the level set {0 = 0} from 0’ 1 + 1 . This is required to visualize the current position of the front in the image plane. A piecewise linear approximation for the front 7 (f) is constructed as follows. Given a cell C(i,j ), if max(0 ! j,0, +l j) 0,j + i ) 0 i+lj . +1 ) < 0 or min(0 t>J , 0 t+1 j, 0,j +1 , V’.+ij+i) > 0, then C(i,j) (£ 7 (f) and is ignored, else, the entry and exit points where 0 = 0 are found by linear interpolation. This provides two nodes on 7 (f) 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 = 0}, which is used for future evaluation of the image-based speed term in the update equation (5.13). 4. Replace n by n -f1 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 image-based speed term. This is because at each grid point, we must search for the closest point lying on the level set {0 = 0}. Moreover, if F g = 0, then the stability requirement for the explicit method for solving our level set equation is At = O(Az). For the full equation (5.13), the stability requirement is At = 0( Ax 2 ). 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 image-based speed (k^. Alternately, we could down-sample 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

PAGE 43

X Figure 5.1. A narrow-band of width 8 around the level set {ip = 0}. embed our level set algorithm in a multiresolution framework. Instead, we improve the time complexity of our algoritnm by updating the level set function ip 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.1) the bold curve depicts the level set {if, = 0} and the shaded region around it is the narrow-band. The narrow-band is bounded on either side by two curves which are a distance 8 apart, i.e., the two curves are the level sets {ip = ±8/2}. The value of 8 determines the number of grid points that fall within the narrow-band. Since, during a given time step the value of ipij is not updated at points lying outside the narrow-band, the level sets {V> >| 8/2 |} remain stationary. The zero set which lies inside moves until it collides with the boundary of the narrow-band. 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.

PAGE 44

36 Therefore, as a consequence of our update strategy, the front can be moved through a maximum distance of <5/2, either inward or outward, at which point we must rebuild an appropriate narrow band. We reinitialize the xp function by treating the cunent zero set configuration, i.e., {xp = 0 }, as the initial curve 7 ( 0 ) (see section 2). Note that the reinitialization procedure must account for the case when {xp = 0} changes topology. This procedure will restore the meaning of xp function by correcting the inaccuracies introduced as a result of our update algorithm. Once a new xp function is defined on the grid, we can create a new narrow-band around the zero set, and go through another set of, say /, iterations in time to move the front ahead by a distance equal to 8/2. The value of l is set to the number of time steps required to move the front by a distance roughly equal to 8/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 narrow-band, compute the extension ki of image-based speed term. 3. With the above value of extended speed term and xp^j, calculate V >™ +1 using the upwind, finite difference scheme given in equation (5.13). 4. Construct a polygonal approximation for the level set {xp = 0 } from xp™* 1 . A contour tracing procedure is used to obtain a polygonal approximation. Given a cell (*, 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 (f). Once such a cell is found, the process is repeated until the contour closes on itself.

PAGE 45

37 The set of nodes visited during this tracing process constitutes the polygonal approximation to 7 (<). 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 {0 = 0 }. A polygonal approximation is required in step 2 for the evaluation of image-based speed term and more importantly, in step 6 for reinitializing the 0 function. 5 . Increment m by one. If the value of m equals /, go to step 6, else, go to step 2. 6. Compute the value of signed distance function 0 by treating the polygonal pproximation of {0 — 0 } as the initial contour 7(0). As mentioned earlier, a more general method of reinitialization is required when {0 = 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 0 was initialized and updated on a square grid. To update 0 tJ at a boundary point it is necessary to specify the boundary conditions in order to numerically approximate the derivatives. We have chosen free-end 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 0 at points lying m the narrow-band, the issue of specifying boundary conditions for points lying on the edge of the band becomes pertinent. The focus of front propagation in the context of shape recovery is merely rapid advancement of the front. Therefore, at the expense of a little inaccuracy in the zero set configuration, we ignore the sensitive issue of specifying correct boundary conditions. In our implementation, derivatives are evaluated normally even if it involves a point that does not belong

PAGE 46

38 in the narrow-band. On the other hand, in 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 (5.2), we show the result of applying narrowband updation algorithm to a star shaped front propagating with speed F = —A”, where A is the curvature as in equation (3.6). 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 ip function was iccomputed once every (/ =) 40 time steps. In figure 5.2(a), we show the initial curve along with the level sets {ip = ±f0.04), where i e [1..5]. After 40 narrow-band updations (figure 5.2(b)), only the level sets {ip <| 0.08 |} move and the rest remain stationary. We note the inconsistency between the level sets lying on either side of the narrow-band, making the reinitialization step necessary in order to restore the meaning of the ip function. Following the reinitialization step, another 40 update steps are applied (figure 5.2(c)), which “diffuses” 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 5.2(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. We conclude this section by highlighting some advantages of the narrow-band updation scheme. The first is the obvious improvement in time complexity as a result of considering only a fraction of the total number of grid points for updation. The second is related to the issue of computing the value of image-based speed term at points away from the zero set, i.e., the extension. In our original method, the extension of image-based speed term could be unreliable at points that are relatively fai fiom the zero set. In addition, when the front is a circular one, it is not clear

PAGE 47

39 (a) Initialization (c) After 80 iterations (e) After 160 iterations (b) After 40 iterations (d) After 120 iterations (f) After 200 iterations Figure 5.2. Narrow-band updation algorithm applied to a star-shaped front gating with speed F = -K. Calculations were done on a 64 x 64 grid with step At = 0.0003. 0 was recomputed after every 40 time iterations. propaa time

PAGE 48

40 what value of image-based speed term to use at its center. Here, the fact that the center of a circle is equidistant from all the points lying on it causes an ambiguity in the choice of the shortest distance point on the circular front. Clearly, this effect is undesirable. The narrow-band updation algorithm does not require the image-based speed term s evaluation at far away points. Therefore, the extension is very accurate and reliable. 5.3 Numerical Issues in 3D In this section, we provide some details of the numerical method based upon the Eulerian formulation [30] that is used to solve the equation (4.8). To begin, the discussion in the previous sections can be readily extended to 3D by differencing in each space direction to produce the following numerical scheme to approximate the surface motion: = ~ A< (^)., J ,fc{(max(D;V’. J ,^0)) 2 + (min(D+^ tJl fc, 0)) 2 (5.13) +(max(D“ V’tj.fc, 0)) 2 + (mm(D+^0)) 2 +(max(Dj tpij'k, 0)) 2 + (min(D+^\j,fc,0)) 2 } 1/2 A t{eK)(kj) | V0 | . Note that the level surfaces are moving with speed F = jfc 7 ( 1 e K), where I< is the mean curvature. In the context of of 3D shape recovery, we move the level surfaces by updating the function on a three dimensional grid according to the above update equation. The level surface {tp = 0}, which models the shape of interest, is forced to a stop in the vicinity of desired object boundaries by applying an image-based speed coefficient on it. In order to to compute the value of image-based speed coefficient at points away from the zero set, an extension is constructed. More specifically, the shape recovery procedure in 3D consists of the following steps:

PAGE 49

41 Algorithm 3 1. At each grid point. (iAx,j Ay, kAz), the extension of image-based speed coefficient ki is computed. 2. With the value of extended speed coefficient (k/)ij,k, and tp^ jk , calculate ipffl using the upwind, finite difference scheme given in equation (5.14). 3. Construct an approximation for the level surface {ip = 0} from ipff l k . This is done in the following fashion. Across every two vertices of a given cell C{i,j, k), the value of ip function is sampled. If the ip function changes sign across a given edge, the exact point at which ip — 0 is found by linear interpolation. Note that at each cell, the ip function can change sign across a number of edges that lies between 3 and 6. The collection of all such points constitutes the approximation to the level surface {0 = 0}, which is used for the evaluation of k r during the next time step. 4. Replace n by n + 1 and return to step 1. In the previous section, it has been established that the extension calculation is very expensive. This observation has motivated the narrow-band approach for front propagation. The situation in 3D is even worse due to an extra dimension. In addition, there is a large degree of repetition in the zero set due to the lack of an organized method in which to populate the zero set. By organized method we mean a 3D analog of a 2D contour tracing procedure. This limitation adds to the size of the zero set and subsequently to the time complexity of extension calculation. In 2D, we alleviated this situation by adopting a narrow-band update strategy, but implementing a similar scheme in 3D could prove overly cumbersome. Therefore, we exploit

PAGE 50

42 the inherent data independence that is present in our grid calculation, and devise a parallel implementation to move the level surfaces. The parallel implementation was done on a I\SR (Kendall Square Research) shared memory parallel machine. For a grid with N 3 points, we have used N processors and realized a substantial improvement in the running time of our 3D shape recovery algorithm. Steps 1 & 2 of the above algorithm are done in parallel and step 3 is executed in serial mode.

PAGE 51

CHAPTER 6 EXPERIMENTAL RESULTS In this chapter we present several shape recovery results that were obtained by applying the narrowband level set algorithm to image data. Subsequently, we present some results in 3 D. Eirst, in 2 D, 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 plans 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 ip i.e., ^(x, 0 ) is computed from 7(0). We first discretize the level set function ip on the image plane and denote ip itj as the value of ip 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 2 k x 2 h 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 ipij 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 ipij is computed at time t = 0 by following the above procedure, we use the update equations from the previous chapter to move the front. 43

PAGE 52

44 6.1 2D Shape Recovery Results We now present our shape recovery results in 2D. We first consider a 256 x 256 CT (computed tomography) image of an abdominal section shown in figure 6.1(a). Our intention is to recover the shape of the liver in this particular slice. The function ip has been discretized on a 128 x 128 mesh, i.e., calculations are performed at every second pixel. In figure 6.2(a), we show the closed contour that the user places inside the desired shape at time t = 0. The function ip is then made to propagate in the normal direction with speed F — A:/ ( — 1 .0 — 0.025 K). We have employed the narrow-band updation algorithm to move the front. The time step size was set to At = 0.0005 and the ip function was recomputed after every 50 time steps. Figure 6.1(b) shows the image-based speed term which is synthesized according to equation (4.4). This term when applied on the propagating front, acts as a stopping criterion, thereby causing the front to come to a rest in the neighborhood of a desired shape. Note that in figure 6.1(b), ki{x,y) values lying in the interval [0..1] have been mapped into the interval [0..255]. In figures 6.2(b) through 6.2(e) we depict the configuration of the level set {ip = 0} at four intermediate time instants. The final result is achieved after 575 time iterations and is shown in figure 6.2(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 6.1(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 high 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

PAGE 53

45 (b) Speed term kj (a) Original image Figure 6.1. Image-based speed term */(x,y) = with a = 3.25, synthesized from the CT image. of our next experiment is to demonstrate smoothness control in the context of shape recovery. In figures 6.3(a) through 6.3(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 fj>(x, y ) from these curves. The level sets are then made to propagate with speed F = Ar/ ( 1 .0 — eK ). First, as shown in figure 6.3(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 6.3(b) and 0.75 in figure 6.3(c). Clearly, with every increment in the value of e, the level set {ij> = 0} attains a configuration that is relatively smoother. This is analogous to the controlledcontinuity provided by thin-plate-membrane splines [41, 17]. 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 branches and significant protrusions. In this experiment we compare the performance of our scheme with the active contour

PAGE 54

(a) Initialization (b) After 100 iterations (c) After 175 iterations (d) After 300 iterations (e) After 450 iterations (f) After 575 iterations Figure 6.2. Recovery of the liver shape from a CT image of an abdominal section. Narrow-band computation was done on a 128 x 128 grid the front was made to propagate with speed F = *,(-1.0 0.025/t) and the time step At was set to 0.0005. iJj was recomputed once every 50 time steps.

PAGE 55

47 (a) e 0.05 (b) e = 0.25 (c) e = 0.75 1' igure 6.3. Smoothness control in shape recovery can be achieved by varying the curvature component in the speed F = &/(1.0 el<). model to bring the limitations of the later into focus. We first attempt to reconstruct the complex arterial structure using a snake model with inflation forces [9]. In figures 6.4(a) through 6.4(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.4(c), 6.4(f), & 6.4(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.4(g)), the snake “snaps” back into a relatively bumpless configuration in figure 6.4(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,

PAGE 56

A 48 we apply our level set algorithm to reconstruct the same shape. After initialization in figure 6.5(a), the front is made to propagate in the normal direction. We employ the narrow-band updation scheme with a band width of 8 = 0.075 to move the front. It can be seen that in subsequent frames the front literally “flows” into the branches and finally in 6.5(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 single instance of our shape model “sprouts” branches and recovers all the connected components of a given shape. In figures 6.6(a)-(i) we plot the other level sets to elucidate the process of extending the image-based speed function to points away from the zero set. 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 transformation to reconstruct the constituent shapes in an image. The image shown in figure 6.7(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 {if = 0} first wraps itself tightly around the objects (see figures 6.7(c) & 6.7(d)) and subsequently splits into four separate closed curves (figure 6.7(e)). While the first three closed segments of {if = 0} recover the three distinct shapes, the one in the middle (see figure 6.7(e)), since it does not enclose any object, eventually disappears. Figure 6.7(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

PAGE 57

49 figures 6.8(a)-(i) we show the level sets {ip = ±iC}, i € [—5.. + 5], with C = 0.05. The function ip 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 hand-written characters. The characters “A” and “B” in the figure (6.9) are examples of shapes with holes. Recovery of hand-written 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 [21, 8] or skeletons [29, 23]. 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 6.9(a), we show the initial contour which encloses all the characters. This contour is then made to propagate inward with a constant speed. Figures 6.9(b) shows an intermediate stage in the front’s evolution and in 6.9(c), the front splits into four separate contours. The calculation comes to a halt when in figure 6.9(d), the level set {ip = 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 6.9(d) as an initial state, and propagate all the fronts inward by momentarily relaxing the image-based speed term. This causes the zero set to move into the character shapes as shown in figure 6.9(e). Finally in figure 6.9(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 128 x 128 grid and the time step At was set to 0.00025.

PAGE 58

50 6.2 Experimental Results in 3D In this section, we present some front propagation and shape recovery examples in 3D. The first two experiments demonstrate the basic surface motion scheme in 3D with a constant and curvature dependent speed. To visualize the level set {ip = 0}, we use a tetrahedronization algorithm. Given a function ipij^ on a grid, this algorithm first marks off all the points with negative ?/; values (here we assume that the points lying inside the surface have negative ip values and the points lying outside have positive ip values). A tetrahedronization procedure is then applied to these points. The triangular facets that lie on the boundary of the tetrahedron constitutes our approximation to the level surface {ip = 0}. Using these surface triangles the level set {ip — 0}, which forms our 3D surface model, is rendered as a shaded shell. In our first 3D experiment, we evolve a toroidal surface with constant speed F = 1. The computational box is a cubical one which extends from -0.75 to 0.75 in all three coordinate directions. Initial value of ip is defined according to the following equation: y , z) = 0-015 {(a: 2 + y 2 ) 1/2 0.4} 2 z 2 . (6.1) The set of all points (x,y,z) at which the value of ip = 0, lie on the surface of the torus (see figure 6.10(a)). We discretize the ip on a 32 x 32 x 32 grid and follow a nested set of toroidal shapes by updating it on the grid according to the equation (5.14). Figures 6.10(b)-6. 10(d) depict three intermediate stages wherein the inner radius gradually becomes smaller until in figure 6 10(e) the inner radius collapses to zero, normals collide, and topology changes. The surface at this stage looks like a dented sphere and will asymptotically attain a spherical configuration. In figure 6.10(f), the level set {ip = 0} intersects with the edge of the computational box causing the slicing of the toroidal surface.

PAGE 59

51 In our next experiment, we collapse a “pinched” superquadric shape under its mean curvature. The initialization was done using the implicit function form of a superquadric, namely where e x ,e 2 > 0 are “squareness” parameters. By varying the values of e x and e 2 , several interesting shapes can be generated. Flat beveled shapes are produced when either e x or e 2 = 2 and pinched shapes are produced when either e x or e 2 > 2. Figure 6.11(a) shows the initial superquadric shape with sharp corners. This shape is made to move with speed F = —A , where K is the mean curvature. As shown in figures 6.11(b)6.11(d), the high curvature regions corresponding to the sharp corners on the suiface smoothly diffuse until the superquadric collapses into an ellipsoid in figure 6.11(e). Note that an ellipsoid is a superquadric with e x = e 2 = 1. Since the mean curvature on the surface of an ellipsoid is constant, the surface eventually collapses to a point without changing its shape (see the smaller ellipsoid in figure 6.11(f)). Calculations were done on a 32 x 32 x 32 grid with a time step At = 0.0001. In our last experiment, we recover the shape of a flat superquadric using the level set front propagation scheme in 3D. Image data for this experiment consists of 32 images each with a particular cross section of the superquadric. The image-based speed term ki is computed from these images according to equation (4.7). A sphere, which is the level surface — 0} of a function ip(x,y,z) = x 2 + y 2 -f z 2 — 0.01, forms our initialization (see figure 6.12(a)). This initial surface is moved with speed F = ki by updating the value of if) on a discrete grid according to the details given in algorithm 3. The initial surface expands smoothly in all directions until a portion of it collides with the superquadric boundary. At points with high gradient, the k t values are close to zero and cause the zero set to locally come to stop near the boundary

PAGE 60

52 of the superquadric shape. This situation is depicted in figures 6.12(b)6.12(e), wherein the initial spherical shape transforms into a flat superquadric. Finally, in figure 6.12(f), all the points on our shape model are stopped, thereby recovering the entire shape of the flat superquadric. Calculations were done on a 32 x 32 x 32 with a time step At = 0.0025.

PAGE 61

53 (a) Initialization 1 (d) Initialization 2 (g) Initialization 3 (b) 500 iterations (h) 500 iterations (c) 1000 iterations (f) 1000 iterations (i) 1000 iterations Figure 6.4. An unsuccessful attempt to reconstruct a complex shape with “significant 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.

PAGE 62

(a) Initialization (b) After 60 iterations (c) After 123 iterations (d) After 200 iterations (e) After 275 iterations (f) After 391 iterations igure 6.5. Reconstruction of a shape with “significant” protrusions: an arterial tree stiucture. Computation was done on a 64 x 64 grid with a time step At = 0.001. A narrow-band updation scheme was used with a band width of 8 = 0.075.

PAGE 63

55 (a) Initialization (b) 60 iterations (c) 123 iterations (d) 170 iterations (e) 200 iterations (g) 320 iterations (h) 350 iterations (f) 275 iterations (i) 390 iterations Figure 6.6. Reconstruction of a shape with protrusions: figure shows the level sets W = ±iC}, t e [-5.. + 5], and C = 0.05.

PAGE 64

56 (a) Initialization (c) After 60 iterations (e) After 125 iterations (b) After 30 iterations (d) After 100 iterations (f) After 140 iterations Figure 6.7. Topological split: a single instance of tl i model splits into three was done on a 64 x 64

PAGE 65

57 (a) Initialization (d) 100 iterations (g) 154 iterations (b) 60 iterations (e) 120 iterations (h) 170 iterations (c) 89 iterations (f) 130 iterations (i) 185 iterations Figure 6.8. Topological split example: level sets shown are {xb = ±iC\ i G 1-5.. + 51 with C = 0.05.

PAGE 66

58 (b) 545 iterations (a) Initialization (c) 760 iterations (d) Outer boundaries Figuie 6.9. Shapes with holes: a two-stage scheme is used to arrive at a complete shape description of both simple shapes and shapes with holes. Computation was performed on 128 x 128 grid and the time step At was set to 0.00025.

PAGE 67

(a) Initialization (c) After 30 iterations (e) After 90 iterations 59 (b) After 10 iterations (d) After 70 iterations (f) After 130 iterations Figure 6.10. Expanding torus ( F — 1): change of topology is depicted in this example. Calculations were done on a 32 x 32 x 32 grid with a time step At = 0.0005.

PAGE 68

(a) Initialization (b) After 10 iterations (c) After 30 iterations (d) After 50 iterations (e) After 70 iterations (f) After 90 iterations Figure 6.11. Motion of a “pinched” superquadric shape under its mean curvatun Calculat.ons were done on a 32 x 32 x 32 grid with a time step At = 0.000L

PAGE 69

61 (a) Initialization (b) After 20 iterations (c) After 40 iterations (d) After 70 iterations (e) After 90 iterations (f) After 120 iterations Fl gure 6. 12. Shape recovery in 3D: a flat superquadric shape. Calculations were done on a 32 x 62 x 32 grid with a time step A/ = 0.0025.

PAGE 70

CHAPTER 7 CONCLUDING REMARKS In this dissertation 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 and Sethian [30] to the problem of shape recovery. With this approach, complex shapes can be reconstructed. 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 extendibihty 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 entropy-satisfying 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

PAGE 71

63 complexity is realized by adopting the narrow-band, update strategy. In the narrowband scheme, the front is advanced by updating the i\) 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 focussed on extending this method to other application domains. We have implemented 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 chapter 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 narrow-band strategy has drastically reduced the computational cost in 2D, but implementing a similar scheme in 3D could prove cumbersome. One approach is to embed our level set algorithm m a multiresolution framework, where the loss of accuracy at rapidly-converging coarse grid computation is compensated by highly accurate slowly-converging fine grid calculation. Alternately, sophisticated parallel implementations can be devised. One such scheme is briefly discussed in chapter 4.

PAGE 72

REFERENCES [1] D. Adalsteinsson and J. A. Sethian, “A narrow-band approach to level set techniques for propagating interfaces,” LBL Report, Oct. 1993. [2] R. Bajcsy and S Kovacic, “Multiresolution elastic matching,” Computer Vision Graphics, and Image Processing , Vol. 46, pp. 1-21, 1989. [3] R. Bajcsy and F. Solina, “Three dimensional object representation revisited,” in l roceedings of First International Conference on Computer Vision , pp 231-240 London, England, 1987. ’ [4] T. 0 Binford, “Visual Perception by computer,” invited talk, IEEE Systems and Control Conference, Miami, Florida. ® R ’ ^ G ‘ Vemuri > “ 0n three-dimensional surface reconstruction J EEE Trans on Pattern Analysis and Machine Intelligence, Vol. LAM1 13, No. 1, pp. 1-13, 1991. [6] T.E. Boult and TR. Render, “Visual surface reconstruction using sparse depth 68-76 LnKM98™ E ^ C ° mputer Vision and Pattern Recognition, pp. [7] A. Blake and A. Zisserman, Visual Reconstruction , MIT Press, Cambridge, MA. [8] H. Blum, “A transformation for extracting new descriptors of shape,” Models for M A ^MIT^Press°^l 967 CC/l ^ W ' Wathen ~ Duni b ed -), Cambridge [9] L. D Cohen “On Active Contour Models and Balloons,” Computer Vision Graphics, and Image Processing , Vol. 53, No. 2, pp. 211-218, March 1991. [iO] L I). Cohen and I Cohen, ^“Deformable models for 3D medical images using hnite elements and balloons, in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition , Urbana, Illinois, June 1992. ^ ci i fR a NJ U 1974 and A ' Bj6rck ’ Nurnerical Methods, Prentice-Hall, Englewood [12] 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 Recognition, pp. 467-472, Maui Hawaii, June [13] M. P. Do Carmo, Differential Geometry of Curves and Surfaces, Prentice-Hall Englewood Cliffs, New Jersey, 1976. 64

PAGE 73

65 [14] L. C. Evans and J. Spruck, “Motion of level sets by mean curvature. I,” Journal of Differential Geometry , Vol. 33, pp. 635-681, 1991. [15] J. A. Fessler and A. Macovski, “Object-based 3D reconstruction of arterial trees Irom magnetic resonance angiograms,” IEEE Transactions on Medical Imaging Vol. 10, No. 1, pp. 25-39, March 1991. [16] W. T. Freeman and E. IF Adelson, “Steerable filters for early vision, image analysis, and wavelet decomposition,” in Proceedings of ICCV, pp 406-415 Osaka, Japan, 1990. [17] M. Kass, A. Witkin, and D. Terzopoulos, “Snakes: Active Contour Models,” International Journal of Computer Vision , pp. 321-331, 1988. [Id] B. B. Kimia, A. Tenenbaum and S. W. Zucker, “Toward a computational theory of shape,” Lecture Notes in Computer Science , Vol. 427, pp. 402-407, 1990. [19] B. B. Kimia, A. R. Tannenbaum, and S. W. Zucker, “Shapes, shocks, and deformations I: The components of shape and reaction-diffusion space,” to appear in International Journal of Computer Vision , 1992. [20] L. Lapidus and G. Pinder, Numerical Solutions of Partial Differential Equations in Science and Engineering , John Wiley k Sons, New York, 1988. [21] D. T. Lee, “Medial Axis Transformation of a Planar Shape,” IEEE Trans. Patt Anal. Machine Intell. PAMI-4, No. 4, pp. 363-369, 1982. [22] D. Lee and T. Pavlidis, “One-dimensional regularization with discontinuities,” IEEE Trans, on Pattern Analysis and Machine Intelligence , Vol. PAMI 10 dd S90_QOO 1 nee J ’ ’ rF[23j F. Leymaiie and M. D. Levine, “Simulating the grass-fire transform using an active contour model,” IEEE Trans. Patt. Anal. Machine Intell. PAMI14 No 1, pp. 56-75, 1992. [24] R. Malladi, “Deformable models: Canonical parameters for surface representation and multiple view integration,” M.S. Thesis, Department of computer and information sciences, University of Florida, Gainesville, May 1991. [25] R. Malladi, J. A. Sethian, and B. C. Vemuri, “A topology-independent shape modeling scheme, in Proceedings of SPIE Conference on Geometric Methods in Computer Vision If San Diego, July 1993. [26] R. Malladi, J. A. Sethian, and B. C. Vemuri, “Shape modeling with front gation: A level set approach,” Center for Pure and Applied Mathematics, PAM-589, Umv. of California, Berkeley, August 1993. propaReport [27] R. Malladi, J. A. Sethian, and B. C. Vemuri, “A fast level set based algorithm for topology-independent shape modeling,” UF-CIS Technical Report TR-93031, University of Florida, Gainesville, Sept. 1993. [2b] R. Malladi, J. A. Sethian, and B. C. Vemuri, “Evolutionary fronts for topologyindependent shape modeling and recovery,” submitted to the Third European Conference on Computer Vision , Stockholm, Sweden, Oct. 1993.

PAGE 74

66 [29] N. Mayya and V. T. Rajan, “Voronoi Diagrams of Polygons a robust and prac1993 a g ° nthm ’ Tech ' Report UF CIS TR-93-29, Univ. of Florida, Gainesville, [oO] S. Osher and J. A Sethian, “Fronts propagating with curvature dependent speed: f0rmUlati ° n ’” Journal °f Computational ^ R , Pa . rent „ and S -W7 ;; icker > “Trace inference, curvature consistency and curve PAMI-ll’ pp n m E m a im° n Pattern Analysis and Machine Intelligence , Vol. ^ ? AT P pn ; Par ' S . : Structured descriptions of shape,” in Proceedings of AAAI, Philadelphia, PA, pp. 695-701, 1986. J [33] P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffu7° n ’ aoQ^i ra ^on Pattern Analysis and Machine Intelligence , Vol. 12, No. *? pp* 'boy, iyyu. [34] W. Press, B. Flannery, S. Teukolsky, and W. Vetterling, Numerical Recipes in G, Cambridge University Press, Cambridge, 1988. [35] R Samadani “Changes in connectivity in active contour models,” Proceedings of the Workshop on Visual Motion , pp. 337-343, Irvine, California, March 1989. [36] R. Samadani, “Adaptive snakes: control of damping and material parameters,” Proceedings of SPIE Conference on Geometric Methods in Computer Vision, Vol. 1570, pp. 202-213, San Diego, California, July 1991. [37] L.L. Schumaker, Fitting Surfaces to Scattered data,” in Approximation Theory P^sfN^ r %rk,^976 ChU1 ’ ^ L ' L ‘ SdlUmaker ’ ( eds ')> pp ' 203-267, Academic ^ !i^Tr tUre the evolution of fronts,” Commun. in Mathematical Physics, Vol. 101, pp. 487-499, 1985. ^ LV 6thi f n ’ “ Nu , mericaI algorithms for propagating interfaces: Hamilton31 pp m lei" U)90 COnSerVatlOn aWS ’” J ° Urnal °f Di ff erentia l Geometry, Vol. [40] R. Szeliski and D. Tonnesen, “Surface modeling with oriented particle systems ” Computer Graphics SIGGRAPH, Vol. 26, No. 2, pp. 185-194, July 1992. [41] ?• aRegula n iZati ° n i in , verse visual problems involving discontinu8 No 2 pp 413 n 424° n i986^ ern Ana ySlS and Machine Intelligence , Vol. PAMI [42] D. Terzopoulos, “The computation of visible surface representations,” IEEE “nvoo Analysis and Machine Intelligence , Vol. PAMI 4, Vol. 10, pp. 41 (— 4oo, 1988. [43] D Terzopoulos, “Matching deformable models to images: Direct and iterative olutions, Proc of Optical Society of America, Tropical Meeting on Machine Vision, Incline Village, NV, March 18-20, 1987.

PAGE 75

67 [44] D. Terzopoulos, A. Witkin, and M. Kass, “Symmetry seeking models and 3D object reconstruction,” Int. Journal of Computer Vision, Vol. 1, No. 3, pp. 211221, 1987. 1 [45] D. Terzopoulos, A. Witkin, and M. Kass, “Constraints on deformable models: 3D Shape and nonri S id motion,” Artificial Intelligence , 36, pp. 91lzo, 1988. [46] B. C. Vemuri and R. Malladi, Deformable models: Canonical parameters for surface representation and multiple view integration,” Proc. IEEE Conference 1991 mpUtCr Vlsi ° n and Pattern Recognition , pp. 724-725, Maui Hawaii, June [47] B. C. Vemuri and R. Malladi, “Surface griding with intrinsic parameters,” Pattern Recognition Letters , Vol. 13, No. 11, pp. 805-812, November 1992. [48] B. C. Vemuri and R. Malladi, Intrinsic parameters for surface representation using deformable models,” IEEE Trans, on Systems, Man, & Cybernetics , Vol. 23, No. 2, pp. 614-623, March/April 1993. [49] B. C. Vemuri and R. Malladi, “Constructing intrinsic parameters with active models for invariant surface reconstruction,” IEEE Trans, on Pattern Analysis and Machine Intelligence , Vol. 15, No. 7, pp. 668-681, July 1993. [oO] B C. Vemuri, A. Mitiche, and J. K. Aggarwal, ” Curvature-based representation pp°107-n4 1986 anSe data ’” InL Journal of Image and Vision Computing , 4, [51] Y. F. Wang and J. F. Wang, “Surface reconstruction using deformable models with interior and boundary constraints,” in Proceedings of ICCV pp 300-303 Osaka, Japan, 1990. ’ [52] D. M. Young, Iterative Solution of Large Linear Systems , Academic Press, New Ynrk 1 U71

PAGE 76

BIOGRAPHICAL SKETCH Ravikanth Malladi was born in Vijayawada, India, in 1966. He received his B.Eng. degree with honors in electrical engineering and M.Sc. degree in physics from Birla Institute of Technology and Science, India, in 1988. In summer of 1985 he worked at the National Geophysical Research Institute, India. During spring of 1988 he held a research assistant position with the semiconductor physics group in Central Electrical and Electronics Research Institute, India, where he was involved in modeling the exposure and development of positive photoresist. He subsequently received the M.S. degree in 1991 and expects to receive the Ph.D. degree in computer and information sciences from the University of Florida, Gainesville, in 1993. Since the summer of 1990 he was a graduate research assistant at the University of Florida, and was affiliated to the Center for Computer Vision and Visualization. His research interests are focused on computational vision, shape modeling, surface reconstruction, computer graphics, image processing, and medical image interpretation. 68

PAGE 77

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. VP Q\x\ i c v i BabaTCTVemuri, Chairman Associate Professor of Computer and Information Sciences I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. James A. Sethian, Cochairman Professor of Mathematics University of California at Berkeley I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. ^rrofes <2L ard X. Ritter rofessor of Computer and Information Sciences I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Andrew F. Laine Assistant Professor of Computer and Information Sciences

PAGE 78

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for tjie degree of Doctor of Philosophy. A ^4 Pail 1 A. Fishwick Associate Professor of Computer and Information Sciences I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Kwmi ft 1 Jl Murali Rao Professor of Mathematics This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree^ofTJo^tor of Philosophy. December 1993 Z 7 Winfred M. Phillips Dean, College of Engineering Karen A. Holbrook Dean, Graduate School