
Citation 
 Permanent Link:
 http://ufdc.ufl.edu/AA00034778/00001
Material Information
 Title:
 Snake pedals : active geometric models for shape modeling and recovery
 Creator:
 Guo, Yanlin, 1971
 Publication Date:
 1998
 Language:
 English
 Physical Description:
 xii, 144 leaves : ill. ; 29 cm.
Subjects
 Subjects / Keywords:
 Coordinate systems ( jstor )
Curvature ( jstor ) Ellipses ( jstor ) Geometric shapes ( jstor ) Parametric models ( jstor ) Pedal points ( jstor ) Snakes ( jstor ) Stiffness matrix ( jstor ) Three dimensional modeling ( jstor ) Topology ( jstor ) Computer vision  Mathematical models ( fast ) Dissertations, Academic  Electrical and Computer Engineering  UF ( lcsh ) Electrical and Computer Engineering thesis, Ph. D ( lcsh ) Image processing  Digital techniques ( fast )
 Genre:
 bibliography ( marcgt )
nonfiction ( marcgt )
Notes
 Thesis:
 Thesis (Ph. D.)University of Florida, 1998.
 Bibliography:
 Includes bibliographical references (leaves 137143).
 Additional Physical Form:
 Also available online.
 General Note:
 Typescript.
 General Note:
 Vita.
 Statement of Responsibility:
 by Yanlin Guo.
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 nonprofit 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:
 029542136 ( ALEPH )
41284335 ( OCLC )

Downloads 
This item has the following downloads:

Full Text 
SNAKE PEDALS: ACTIVE GEOMETRIC MODELS
FOR SHAPE MODELING AND RECOVERY
By
YANLIN GUO
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
1998
To My Motherland
I wish you wealth, strength, and happiness.
ACKNOWLEDGMENTS
I would like to express my sincerest gratitude to my advisor and committee chairman, Dr. Baba C. Vemuri, for his invaluable guidance, constant support and encouragement, in every step of my education and research at the University of Florida. I have greatly benefited by his stimulating approach to research, complete dedication to science, and relentless pursuit of perfection.
My special thanks go to Dr. Christina M. Leonard for generously agreeing to serve as an external member of my dissertation committee, and for sharing her expertise in visual analysis of the human brain, and for providing many valuable comments on the applications of my work to automatic analysis of brain images as well as providing me with large amount of data. Many thanks to my dissertation committee members, Dr. Thomas E. Bullock, Dr. John G. Harris, and Dr. Stanley Su, for their willingness to serve on my committee and their interest and important comments and suggestions on various aspects of this dissertation.
I would like to thank Dr. Hsu, of the Department of Aerospace Engineering, Mechanics and Engineering Science, for his insightful discussion on finite difference methods.
Many thanks to former student ShangHong Lai (now in Siemens Corp. Research) for his stimulating discussions and ideas. I also would like to thank fellow students from the Computer Vision, Graphics and Medical Imaging lab Chhandomay
iii
Mandal, Yi Cao, Fengting Chen, Li Chen, Shuangying Huang, Li Wang, Jundong Liu, and Arun Srinivasan for the many discussions, ideas, advice, and friendship.
Finally, I would like to thank my parents and sisters for their endless love and support, and my husband, for sharing the happiness and struggle I experienced, without which I could never have succeeded.
The research was supported in part by the NIH under the grant number R10lLM05944 and the Whitaker Foundation.
iv
TABLE OF CONTENTS
ACKNOWLEDGMENTS..................................
LIST OF FIGURES......................................
viii
ABSTRACT......................................
CHAPTERS
1 INTRODUCTION.............
1.1 Problem Statement........
1.2 Literature Review.........
1.2.1 Geometric & Deformable
1
Models with Fixed Topology . ...
1.2.2 Topologically Adaptable Models .. .. .. . 1.3 Contributions.... .. .. .. .. .. .. .. .. .. .
1.3.1 Compact Representation of the Model ...
1.3.2
1.3.3
1.4 Thesis
Automatic Topological Change Handling Scheme. Fast Numerical Methods.................
Organization.........................
2 MATHEMATICAL FOUNDATIONS OF ACTIVE MODELS.........
2.1 Snakes: Energy Minimizing Active Contours................
2.2 Active Surface Models....... .. .. .. .. .. .. .. .. .. .. .. .
3 SNAKE PEDALS: GEOMETRIC MODELS WITH PHYSICSBASED CONTR O L . . . . . . . . . . . . . . . . . . . .
3.1 Definition of Regular Pedal Curves......................
3.2 Snake Pedals: The Proposed Model......................
3.3 Modified Snake Pedal Model...........................
3.4 Handling Changes in Topology.........................
3.5 Pedal Surfaces....................................
3.5.1 Uniform Tessellation of the Model..................
3.6 Summary: Salient Features of the Model.... .. .. .. .. .. .. ..
2
3
4
6
6
8
9
9
10
12
14 18
20
20 23 28 30 31
34 38
v
. . xi
4 PROBABILISTIC FRAMEWORK AND MODEL LEARNING........
4.1 4.2 4.3 4.4
Prior Model in the Bayesian Framework... .. .. .. .. .. .. ..
Sensor Model in the Bayesian Framework.. .. .. .. .. .. ....
Maximum a Posterior........ .. .. .. .. .. .. .. .. .. .. ..
Supervised Learning Scheme... .. .. .. .. .. .. .. .. .. .. .
5 NUMERICAL SOLUTIONS FOR MODEL FITTING WITH SNAKE PEDAL MODELS............................
5.1 Related Work....................................
5.2 Our Fast Numerical Algorithms........................
5.2.1 Preconditioned Nonlinear Conjugate Gradient Algorithm ..
5.2.2
5.2.3
5.2.4
Computation of the Approximate Hessian . ... LM Method for Global Parameter Estimation ADI Method for Local Parameter Estimation
6 HYBRID GEOMETRIC ACTIVE MODELS.............
6.1 Overview of Hybrid Geometric Active Models.........
6.2 Curve Evolution Theory......................
6.3 Traditional LevelSet Approach... .. .. .. .. .. ....
6.3.1 Embedding Mechanism... .. .. .. .. .. .. ..
6.3.2 Equation of Motion.... .. .. .. .. .. .. .. ..
6.3.3 Advantages of the Levelset Approach. .. .. .. .
6.4 Our Model: Hybrid Geometric Active Contours. .. .. .
6.4.1 Big Picture.... .. .. .. .. .. .. .. .. .. ..
6.4.2
6.4.3
6.4.4
6.4.5
Relation Between the Snake and Snake Pedal Curve Evolution PDE for the Snake and Snake Pedal Evolution. .. .. .. .. Relation Between Vol andV02... .. .. .. .. .. .. .. ..
PDE for the Higher Dimensional Function 01.. .. .. .. ..
7 NUMERICAL SOLUTIONS FOR SHAPE RECOVERY WITH HYBRID GEOMETRIC ACTIVE CONTOURS.. .. .. .. .. .. ..
7.1 Numerical Algorithm for Standard PDEbased Curve Evolution ...
7.2 Our Numerical Approach........ .. .. .. .. .. .. .. .. .. .. .
8 SHAPE RECOVERY EXPERIMENTS WITH SNAKE PEDALS .. .. .
8.1 8.2 8.3
Shape Recovery with Snake Pedals in 2D... .. .. .. .. .. .. ...
Shape Recovery with Snake Pedals in 3D... .. .. .. .. .. .. .. .
Fitting Results with the Hybrid Geometric Active Model .. .. .. .
40
42 45 46 46
50
51 52 53 56 58 62
70
72 73 77 77 79 81 82 82 86 88 91 95
97
98
102
105
105 108 115
vi
9 CONCLUSIONS AND FUTURE DIRECTIONS . ..... 120)
9.1 Salient Features of the Model .. .. .. ... ... ... ... ... ..120
9.1.1 Geometric Model with Physicsbased Control .. .. .. .....121
9.1.2 Compact Representation .. .. .. .. ... ... ... ......121
9.1.3 Automatic Handling of Topological Changes. .. .. .. ....121
9.2 Future Directions. .. .. .. ... ... ... ... ... ... ... ..122
9.2.1 Invariant Representation .. .. .. .. ... ... ... ......122
9.2.2 Introduce Dynamics in the Representation. .. .. .. ... ..122
9.2.3 Extend the hybrid geometric active model to 3D .. .. .. ....122
APPENDICES
A DERIVATION OF STIFFNESS MATRIX K .. .. .. .. ... ... .....124
A.1 Bilinear Quadrilateral Elements .. .. .. .. ... ... .... .....124
A.2 Derivation of the Local Stiffness Matrix .. .. .. .. ... ... .....126
A.3 Kronecker Tensor Product Representation of the Assembled Stiffness
Matrix. .. .. .. .. ... .... ... ... ... ... ... ... ..127
B RELATIONSHIP OF ke..0. .... .... .... .... .... ... ..135
REFERENCES. .. .. .. .. ... .... ... ... ... ... ... ... ..137
BIOGRAPHICAL SKETCH. .. .. .. .. ... ... ... ... .... .....144
vii
LIST OF FIGURES
2.1 Active surface being pulled by a spring force showing the effect of
various w, and W2 settings: (a) W1 =0.9, W2 =0.1, (b)wl W 0.5,
and (C) W1 =0.01,W2 =0.99 .. .. .. .. ... ... ... ... .....19
3.1 f is on the pedal curve of a with respect to the pedal point p.. .. ...20
3.2 Examples of pedal curves of an ellipse for different pedal point positions. Pedal points are shown by a dot in each case .. .. .. .. .. ..22
3.3 (a) The process of generating a snake pedal with an ellipse as a generator. (b) "snake pedal" controlled by the snake using an ellipse generator. 24
3.4 Radialbased correspondence association scheme: each point ae(tj) on
the generator is associated with a point pi on the snake. .. .. .. ...25 3.5 Examples of "snake pedals" using an ellipse generator and a snake. .28
3.6 Examples of the modified "snake pedals" using an ellipse generator
and a snake .. .. .. .. .. .. .... ... ... ... ... ... ....29
3.7 Comparisons of original (red) and modified (blue) "snake pedals"
using the same generators and same snakes in both case. .. .. .. ...30
3.8 Examples of topology change. (a) depicts a single snake pedal curve
(solid) obtained from an ellipse generator (dashed) and a single pedal point. (b) same as (a) except two distinct pedal points located outside the generator are used. (c) same as (a) except several pedal points on
a snake are used. (d) same as (c) except a discontinuous snake is used. 31
3.9 Examples of pedal surfaces with fixed pedal points located inside ellipsoidal generator surfaces. .. .. .. .. ... ... ... .... ....33
3.10 Examples of shaping snake pedal surfaces with snakes. Each figure
was generated by applying interactive mousebased forces to the snake. 34 3.11 Uniform tessellation result: (a) f (u) '~ u relation, (b) f (v) ', v relation. 38 4.1 Supervised learning to incorporate information into the prior model . 48
viii
5.1 Big picture of numerical schemes for model fitting.... .. .. .. . .. .
5.2 Model tessellation in the model coordinates .. .. .. .. .. ... ...64
6.1 Level set formulation of equations of motion  (a) and (b) show the
curve P and the surface O(x, y) at t =0, (c) and (d) show the curve P and the corresponding surface O(x, y) at time t. (e) and (f) show the curve P and surface O(x, y) at time t + A~t. The curve P has changed
topology in going from (c) to (e) .. .. .. .. .. .... ... ... ..78
6.2 Relation between P(t), Pe(t), 01, and 042: P (t) and Pe(t) are levelsets
of the higher dimensional surface 01 and 02 respectively. They are
related by the pedaling operation .. .. .. .. .. ... ... ... ....84
6.3 The hybrid geometric active model can capture the "global orientation" of clustered objects .. .. .. .. ... ... ... ... ... ....85
6.4 Three representative pointpairs on the snake and snake pedal .. .. ...92
7.1 Curve Evolution: The curve moves with unit speed in (a), and moves
with curvature dependent speed in (b). The solid lines are the initial
curves and the dashed lines are the curves after evolution. .. .. ....99
8.1 Fitting result (c) of a snake pedal to a caudate nucleus with initialization shown in (a) and intermediate stage shown in (b). .. .. .. ....106
8.2 Fitting result (c) of a snake pedal to a cerebellum with initialization
shown in (a) and intermediate stage shown in (b). .. .. .. .. .....107
8.3 Fitting result (c) of a snake pedal (in green) to a hippocampus with
initialization shown in (a) and intermediate stage shown in (b). The data points (in red) are placed by the expert along the hippocampus
boundary. .. .. .. .. .. .... ... ... ... ... ... ... ..107
8.4 Fitting results of snake pedal surfaces to a hippocampus (d), a gyrus
(h), and a cerebellum (1) with initializations shown in (b), (f), and (j), the intermediate stages are shown in (c), (g), and (k), respectively.
(a), (e) and (i) are slices of MRI scans for the hippocampus, gyrus,
and cerebellum, respectively. .. .. .. ... ... ... ... ... ...10
8.5 Fitting results of snake pedal surfaces to a bend (c) and twist (f) shape.
Initializations are shown in (a) and (d) followed by the intermediate
stages (b) and (e) respectively. .. .. .. .. .. ... ... ... ....113
8.6 Fitting results of snake pedal model to spock(c) and an anvil (f). Initializations are shown in (a) and (d) followed by the intermediate stages
of fitting in (b) and (e) respectively. .. .. .. .. .. ... ... ....114
ix
59
8.7 Hybrid geometric active contour model fit examples: (a) is model initializations (snake pedal in green and generator in red), (b) and (c)
are intermediate stages of evolution and final fits are shown in (d). . . 115
8.8 Hybrid geometric active contour model fit examples: (a) is model initializations (snake pedal in green and generator in red), (b) and (c)
are intermediate stages of evolution and final fits are shown in (d). . . 116
8.9 Hybrid geometric active contour model fit examples: Left to right,
model initialization (snake pedal in green and generator in red), intermediate stages of evolution and final fit .. .. .. .. ... ... ....117
8.10 Topological change examples with hybrid geometric active contour
models: Left to right, model initialization (snake pedal in yellow or green and generator in red), intermediate stages of evolution and final
fit .. .. .. .. .. ... ... ... ... ... ... .... ... ....119
A. 1 Bilinear quadrilateral element, where (u, v) are material coordinates
and (w, V)) are the reference coordinates .. .. .. .. .. .... .....125
A.2 A small 3 x 4 mesh in (u, v) space .. .. .. .. ... ... ... ....128
A.3 Membrane molecules .. .. .. .. ... ... .... ... ... ....129
A.4 Assembled membrane molecular representation. .. .. .. .. ... ..130
x
Abstract of Dissertation
Presented to the Graduate School of the University of Florida
iu Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy
SNAKE PEDALS: ACTIVE GEOMETRIC MODELS FOR SHAPE MODELING AND RECOVERY
By
Yanlin Guo
August 1998
Chairman: Dr. Baba C. Vemuri Major Department: Electrical and Computer Engineering
Modeling shapes, be it in 2D or 3D), is a fundamental constitutent of computer vision and computer graphics. Shape modeling techniques have been widely used in shape synthesis, shape reconstruction, recognition, motion tracking, and so on. In this dissertation, we introduce a compact and versatile geometric shape modeling scheme that can model a large class of shapes and is amenable to stable and efficient numerical implementations. Geometric models are traditionally well suited for representing global shapes but not the local detail. However, in this thesis we propose a powerful geometric shape modeling scheme which allows for the representation of global shapes
xi
with local detail and permits model shaping as well as topological changes via physicsbased control. The proposed geometric models are a blend of geometric and physicsbased models and are in spirit "similar" to the now popular deformable superquadric models but differ from them considerably in the expressiveness and numerical stability leading to greater applicability.
The proposed modeling scheme consists of representing shapes by pedal curves and surfacespedal curves/surfaces are the loci of the foot of perpendiculars to the tangents of a fixed curve/surface from a fixed point called the pedal point. By varying the location of the pedal point, one can synthesize a large class of shapes which exhibit both local and global deformations. We introduce physicsbased control for shaping these geometric models by letting the control point vary and use a dynamic spline, that is, a snake, to represent the position of this varying pedal point. The model dubbed as a "snake pedal" allows for interactive manipulation via forces applied to the snake. We extend the model to automatically cope with topological changes by introducing the concept of hybrid geometric active contour/surface models wherein the traditional snake in the snake pedal is replaced with a geometric active contour which is realized via a levelset implementation.
Efficient numerical algorithms for shape recovery from image data using the proposed snake pedal model are presented. These algorithms involve novel mathematical derivations that lead to efficient numerical solutions of the model fitting problem. We demonstrate the applicability of this modeling scheme via shape synthesis examples and shape estimation examples from image data.
xii
CHAPTER 1
INTRODUCTION
Shape modeling is a fundamental constituent of computer vision and computer graphics. In computer graphics, it can be used to synthesize shapes, object motion and object interaction for the purpose of computer aided design and computer animation. In computer vision, it is very crucial in 2D/3D shape representation, shape reconstruction, quantitative model extraction from image data for analysis, visualization, recognition and motion tracking. Recently, shape modeling techniques have been widely used in medical image analysis applications. Twodimensional and threedimensional shape models have been used in Xray, computer tomography
(CT), magnetic resonance (MR), and ultrasound images to segment, visualize, track, and quantify a variety of anatomic structures including brain, heart, cerebral, lungs, liver and skull.
Two primary tasks of importance in computer vision are shape recognition and shape recovery from image data. The shape recovery task imposes the requirement that the shape model be capable of representing fine detail and hence requires that the model have a large number of degrees of freedom. Shape recognition, on the other hand, requires that the model representation be compact in terms of the number of parameters so as to achieve easy storage of a large database of models and easy matching for retrieval. These two requirements are potentially conflicting. Therefore,
1
2
designing shape models that can possibly "satisfy" these requirements has been the focus of many researchers in the recent years in the computer vision community [3, 11, 14, 54, 78].
1.1 Problem Statement
There are numerous shape modeling techniques in the computer vision and computer graphics literature. Most modeling schemes fall into two major categories, namely, active models and passive models. Active models or deformable, physicsbased models are distributed parameter models. Typical examples are generalized spline models such as snakes [31]. Active models have broad geometric coverage and hence are good for shape recovery and reconstruction. On the other hand, they are not well suited for recognition since they require a large number of parameters for their representation. Passive models or geometric models are lumped parameter models. Typical examples are spheres, cylinders and prisms [4, 67]. Passive models can represent shapes with very few parameters and they are good for recognition but not for representing detailed shapes. In order to fill in the gap between the requirements of reconstruction and recognition, a class of hybrid models have been proposed. For example, Terzopotulos and Metaxas [78] proposed a hybrid model which combines the descriptive power of geometric and physicsbased models via superposition of a spline on a conventional superellipsoid. Hybrid models are a comnp~romise between the two extremes and occupy both ends of the spectrum in shape description. One end of this spectrum is occupied by lumped parameter models and the other by distributed parameter models. However, most hybrid models introduce
3
a lot of nonlinearities into the model representation to describe global deformations such as bending, tapering, and twisting. These nonlinearities increase the complexity in the model representation and affect the numerical stability of the algorithms used in fitting the model to image data. Moreover, the hybrid models can not handle shapes of unknown topologies without introducing additional nonlinearities into the modeling via blending schemes [14] or other special handling schemes.
In this thesis, we develop a compact and versatile shape modeling scheme with physicsbased control for shape design and shape recovery from image data. This modeling technique allows for the representation of global shapes with local detail and allows model shaping via physicsbased control. The local/global descriptive power of the model satisfies the conflicting requirements of both shape reconstruction and shape recognition. In addition, topological changes can be handled naturally and automatically via a levelset implementation in the proposed model.
In the remainder of this chapter, we will first review related shape modeling schemes in literature, with an emphasis on the advantages and disadvantages of each scheme, and then present the overview of our model.
1.2 Literature Review
There are numerous shape modeling schemes in computer vision and computer graphics literature.
4
1.2.1 Geometric & Deformable Models with Fixed Topology
Motived by the generalized cylinder idea, Kass et al. [31] proposed a deformable cylinder model constructed from generalized splines. They developed force field techniques for fitting their model to monocular, binocular, and dynamic image data. The distributive nature of this deformable model enhances its descriptive power and allows the representation of natural objects with asymmetries and fine detail. However, the generalized spline representation of the model does not explicitly provide a description of object shapes in terms of a few parameters.
Terzopoulos and Metaxas [781 developed a hybrid dynamic model a deformable superquadric model which combines the descriptive power of geometric and physicsbased models via superposition of a membrane spline on a core superquadric. Vemuri and Radisavljecvic [81] developed a multiresolution, hybrid modeling scheme which represents the displacement function of a deformable superquadric model in an orthonormal wavelet basis. This multiresolution basis provides the model with the ability to continuously transform from local to global shape deformations thereby allowing a continuum of shape models to be created and to be presented by relatively fewer parameters.
Pentland and Williams [541 introduced a novel shape modeling technique based on the modal representation, which can be used to represent global shapes with local detail. The number of parameters required for this representation depends on the level of detail needed.
5
Cohen and Cohen [11] and Han et al. [26] proposed a hybrid hyperquadric model that adds an exponential term to a hyperquadric primitive. The addition of the exponential term to the hyperquadric equation allows for local control of the shape generated. An arbitrary number of concavities can be generated with this modeling scheme. This modeling scheme differs from most of the schemes discussed above in that this scheme is basically a geometric model but it allows for global as well as local deformations. The model representation is quite compact in terms of the number of parameters needed to describe the shape. However, the computation of local concavities is nontrivial and human interaction is needed in the model fitting process.
Another modeling scheme described by several authors in [3, 27, 30, 59] involves superimposing a global volumetric deformationthe free form deformation (FFD)on the core superquadric. Both global and local deformations can be represented by this modeling scheme, however, it is difficult to describe the shape of an object with complex deformations and fine detail with this modeling scheme.
Recently, O'Donnell et al. [50] developed a novel solid shape modeling scheme. This model includes builtin offsets from a core base model. Local deformation over the scaledoffset model is allowed via displacement vectors attached to the default model. Model fitting to data is achieved in a dynamic framework as in Terzopoulos and Metaxas [78]. This modeling scheme is well sited for tasks such as cardiac motion recovery fromi tagged MIR images.
6
1.2.2 Topologically Adaptable Models
Surfaces of arbitrary topology were modeled using particle systems to model surfaces in Szeliski and Tonnesen [72]. Particles can be added and deleted dlynamnically to enlarge and trimi the surface.
An interesting and topologically adaptable deformable model was developed by DeCarlo and Metaxas [14] via the use of parameterized blending functions. Their scheme can determine the number and locations of the holes in an object with unknown topology and accordingly change the topology by use of appropriate error of fit criteria.
Numerous techniques based on the concept of levelsets have been reported in the literature [9, 10, 46, 58, 82]. These techniques allow a deformable contour not only to represent long tubelike shapes or shapes with bifurcations, but also to dynamically sense and change its topology. Another modeling scheme called "topologically adaptable snakes" was introduced by Mclnerney and Terzopoulos [48]. This novel version of snakes allows for automatic change in topology during shape recovery. A common problem associated with these models is that they do not possess a mechanism to characterize the global/core shape of an object.
1.3 Contributions
In this thesis, we propose a powerful geometric shape modeling scheme which permits the representation of global shapes with local detail using only a "~semiglobal" characterization and relatively simple and efficient numerical methods. The
7
proposed modeling scheme consists of representing shapes by pedal curves and surfacespedal curves/surfaces are the loci of the foot of perpendiculars to the tangents of a fixed curve/surface from a fixed point called the pedal point. By varying the location of the pedal point, one can synthesize a large class of shapes which exhibit both local and global deformations. We introduce physicsbased control for shaping these geometric models by letting the pedal point vary and using a dynamic spline, namely, the snake, to represent the position of this varying pedal point. This process of generating models is referred to as the "pedaling operation." The model, dubbed as a "snake pedal," allows for interactive manipulation via forces applied to the snake. In this model, automatic changes in topology can be realized by embedding the snake in a levelset based evolution.
The key contributions of this thesis are summarized below in the form of salient features of the developed modeling scheme, the snake pedal, and the associated numerical algorithms:
1. The snake pedal model is inherently geometric but can exhibit both local and
global deformations as would a physicsbased model.
2. The model can exhibit large global deformations such as bending, twisting, and
tapering without explicitly incorporating parameters for the same.
3. The modeling scheme allows (for the first time) the concept of a global shape
to be incorporated into the curve/surface evolution framework.
4. The model permits a levelset implementation and thereby seamlessly allows
for topological changes in shape during evolution.
5. Fast numerical algorithms for estimating the model fits to 2D/3D data sets.
6. A hierarchical neural network based learning scheme for incorporating prior
knowledge about shapes of interest from the application domain. We elaborate these contributions in the following sections.
1.3.1 Compact Representation of the Model
The snake pedal model is inherently geometric but can exhibit both local and global deformations as would a physicsbased model. Our modeling scheme is somewhat similar to the hybrid modeling schemes mentioned in Section 1.2.1 in that both techniques combine the descriptive power of physicsbased models and geometric models. The difference lies in that the hybrid models discussed earlier have to explicitly include the global parameters to incorporate global deformations such as bending, tapering, and twisting to achieve the required deformations, which causes additional nonlinearities in their representation, while our model does not. This compact representation of global deformations by our model can avoid the introduction of additional nonlinearities into the model. Therefore, numerical efficiency and stability can be largely improved in the model fitting to the image data.
This compact representation comes from the mechanism of creating shapes in the snake pedal model. The pedaling operation is a nonlinear operation. When applying the pedaling operation to geometric primitives, one can produce global
9
deformations of the geometric primitives by performing the nonlinear transformation. This nonlinear operation is only applied in the model building phase, and it does not affect the model fitting phase, as we will discuss in Chapter 5. Hence no additional nonlinearities are added to the model fitting process.
1.3.2 Automatic Topological Change Handling Scheme
The proposed snake pedal model allows for the estimation of shapes of arbitrary topology without resorting to the blending or other special operations as required in Szeliski and Tonnesen [72]. The topology changes are achieved naturally by embedding the modeling scheme into a levelset framework [46]. We introduce the novel concept of a global/core model into the now popular partial differential equation (PDE) based curve/surface evolution framework. This global/core model allows us to express any shape as a combination of a global shape such as ellipsoid and sphere, and a variable offset with respect to this global shape and an evolving curve/surface. Because of the presence of the global shape, the snake pedal model can be regarded as a unique hybrid geometric active model.
1.3.3 Fast Numerical Methods
Novel and fast numerical methods are proposed to achieve stable local optimal solutions in the model fitting process. These numerical methods are not a regurgitation of techniques that are described in a standard numerical analysis text book, or for that matter, the book on numerical recipes in C [84]. In all the methods that are discussed here, novel derivations are presented which lead to development of software modules that can either be easily embedded into or used with standard
10
numerical software such as the conjugate gradient (CG) code and the alternating direction implicit (ADI) code.
1.4 Thesis Organization
The remainder of the thesis is organized as follows:
Chapter 2 contains a review of the basic mathematical foundations of classical active models, especially the snake formulations in order to illustrate the physicsbased control that is also pertinent to the snake pedal model.
In Chapter 3, we present the definitions of pedal curves and surfaces and show examples of the same in 2D and 3D), followed by the development of the proposed model, namely, the "snake pedal" with a demonstration of the expressive power of the modeling scheme via illustrative examples of shape synthesis. In addition, we also discuss the ways of automatically handling change of topology.
In Chapter 4, we cast the model fitting process in a probabilistic framework, and present the advantages of the probabilistic framework over the deterministic approach. We will devise a novel supervising learning scheme to incorporate the prior knowledge of object shapes into the shape estimation from image data.
In Chapter 5, we present efficient and stable numerical algorithms for fitting the snake pedals to real image data. The model fitting process can be posed as a nonlinear optimization problem and can be dissected into global and local shape parameter estimation. We present efficient numerical solutions for both tasks, and then combine the global and local procedures appropriately to solve the overall optimization problem.
11
In Chapter 6, we extend the snake pedal model to cope with topological changes by introducing the concept of hybrid geometric active models. We first review the traditional geometric active curve evolution theory and its levelset implementation. Following this, we present our hybrid geometric active model evolution and its levelset implementation. The hybrid geometric active contour evolution involves the evolution of both snake and snake pedal curves /surfaces. However, by employing the insightful relationship between tie two curve/surface evolutions, we solve these evolutions simultaneously by only solving for one of them and thereby dramatically reduce the computational and storage requirements when applied to shape recovery tasks.
In Chapter 7, wve present the numerical techniques used for computing a stable solution to the hybrid geometric active model evolution equation. The equation of motion for the model can be posed as a special type of HamiltonJacobi partial differential equation (PDE). Based on the traditional levelset approach, the motion equation of the our model can be efficiently solved by employing entropysatisfying upwind finite difference plus minmod finite difference schemes.
In Chapter 8, we present the experimental results and the applications of the model in shape modeling and recovery and conclude in Chapter 9.
Appendix A gives the detail derivation of the stiffness matrix, which corresponds to the internal smoothness energy imposed on the snake in the snake pedal model. Appendix B contains the derivation of the relationship between the curvature of the snake pedal and that of the snake at a specific point.
CHAPTER 2
MATHEMATICAL FOUNDATIONS OF ACTIVE MODELS
In the proposed snake pedal model, we use an active curve/surface to control the model shapes and deformiations, this active curve/surface is an integratedl andl nonsejparable part of the nmodeling scheme. In tis chapter. wve wvill review the mnathemnatical foundations of active modlels in order to obtain a better understandingof' the behavior of the controlling curves/surfaces in the snake pedal model.
Active or deformable curve and surface models gained popularity after they were proposed for use in computer vision [79] and computer graphics [77] in the mid 1980s. Iii computer vision, deformable curve and surface models were proposed as solutions to illposed inverse problems such as edge detection and surface reconstruction [55, 75]. The theory of continuous deformable models was introduced by Terzopoulos in a Lagranigian dynamic setting [76], based on deformation energies in the form of (cont rolled continuity) generalized splines [75]. The controlledcontinuity spline is a generalization of Tikhonov and Arsenin stabilizer [80], and can formally be regarded as regularizing the illposed problems [56], by recasting them as wellposed functional minimization problems.
Active or deformable model can represent broad classes of shapes by employing geometric representations that involve many degrees of freedom, such as splines. The degrees of freedom are generally not permitted to evolve independently, but are
12
13
governed by physical principles that bestow intuitively meaningful behavior. The use of physical principles, especially elasticity theory, generally within a Lagrangian setting, give rise to the name "deformable models" or "active models." The physical interpretation views deformable models as elastic bodies which respond naturally to applied forces and constraints. Typically, deformation energy functions are defined in terms of the geometric degrees of freedom for the deformable models as for physically elastic objects. The energy of the model grows monotonically as it deforms away from a specified natural or "rest" shape and often includes terms that constrain the smoothness or symmetry of the model. In the Lagrangian setting, the elastic forces internal to the model contribute to the deformation energy, while the external potential energy functions are defined in terms of the data of interest to which the model is to be fitted. These potential energies give rise to external forces which deform the model such that it maximally fits the data.
From a theoretical point of view, the concept of active or deformable models reflects the integrated fluency of geometry, physics, and approximation theory. Geometry serves to represent object shape, physics imposes constraints on how the shape may vary over space and time, and the theory of optimal approximation provides the formal support to the mechanisms for fitting the models to measured data.
The active/deformable model that has attracted the most attention in the computer vision and computer graphics community is popularly known as "snakes" [5, 31].
14
Snakes or "active contour models" represent a special case of the general multidimensional deformable model theory. We will review their formulation in this chapter in order to illustrate the basic mathematical foundation of the controlling curves/surfaces in our proposed shape modeling scheme, namely, the snake pedal model.
2.1 Snakes: Energy Minimizing Active Contours
Snakes are energy minimizing splines guided by internal and external forces. Internal forces act as smnoothiness constraints, while external forces pull the snake toward the desired features such as edges and image boundaries. By initializing the snake around the target object, it can find the desired object contours by minimizing a userdefined energy function. Because the way the contours slither while minimizing their energy, we call them "snakes." The ability of the snakes to converge to the outline of the object despite the intensity gradient inside the object makes them very useful in many visual tasks, such as image segmentation, boundary detection, motion tracking, and stereo matching.
Geometrically, a snake consists of a set of discrete points, effectively connected by straight lines. Each point has a position, given by the coordinates in the plane (x, y), and a snake is completely specified by the number and coordinates of these points, namely, p(s) =(x(s), y(s)) with the domain parameter s E [0, 1]. When the snake is a closed curve, we apply the periodic boundary condition, p(0) =p(l), on the snake.
A snake is controlled by minimizing a functional which converts highlevel contour information, such as curvature and discontinuity, and lowlevelArsenin image
15
information, such as edge gradients into an energy represented by
E (s) = E (p(s) ) ds f {Ej i t(p (s) ) + Ee.. t(p (s)) ds (2.1)
where, Ei, and Eext are the internal and external energies associated with the snake. The shape and position of the snake corresponds to the minimum of this functional. The first termn of the above functional
E,, (p(.v)) I (s) Op 2 + W2 (S) & 2(2.2) O s Os 2
is the internal deformation energy of the snake. It characterizes the deformation f'an elastic contour, w1l(s) and w2(s) are nonnegative parameters controlling the "tension" and "rigidity" of this elastic contour.
The second termn of (2.1) consists of the external potential that guides the snake toward desired features. One( forin of this potential is a poutenltial function (lefimne oil the inlage plane. whose local inimlia coincide with intensity extreima. edlges. and( other image features of interest. For example, if we want to attract thme contour toward the edges of a gray level image I(x, y), the potential takes the form
Eert(p(s)) =jf Po(p(s)) ds, (2.3)
with
Po(p(s)) =Po(x, y) =c JVG, * I(x, y)I ,(24
(2.4)
16
where c controls the magnitude of the potential, V is the gradient operator, and G,* I denotes the image convolved with a Gaussian smoothing filter, and a is the width of the filter which controls the spatial extent of the local minima of Ee,,t.
Another form of the external energy potential is the user defined external constraint potential such as interactive springs, anchored springs, and volcanos [31]. For example, the snake can be pulled in the direction of the mouse cursor location (XMI Yin) by using a spring potential Eext =c[ (x  xm,)' + (y  y,,,)'], where c controls the stiffness of the spring. Note that the points on the snake that are affected by the spring can be restricted to a small section of the snake closest to the mouse.
The combination of external image potentials and external constraint potentials allows snakes to extract and represent a broad spectrum of shapes. Furthermore, external constraint potentials are an effective, flexible means from which highlevel control mechanisms can guide shape recovery, forming a sound basis for fully automatic image analysis.
Using the calculus of variations, we obtain the contours p(s) that minimize the energy E. They must satisfy the EulerLagrange equation
19S s1) T(2 ) + VPO(p(s)) =0. (2.5)
This vectorvalued partial differential equation expresses the balance of internal and external forces when the contour rests at equilibrium. The first two terms represent the internal stretching and bending forces respectively, while the third term represents the external forces that guide the snake to the image data.
17
In order to compute numerically a minimum energy solution, it is necessary to discretize the energy (2.1). By using the finite element method [87] or the finite difference method [70], the discrete form of the energy E(p(s)) can be written as:
E (p) = 1 TKp + VPo(p), (2.6)
where K is the stiffness matrix, and E(p) is the discrete version of the external potential. Differentiating (2.6) with respect to p results in the minimum energy solution, which is equivalent to solving the set of nonlinear algebraic equations:
Kp + f (:, y) =0, (2.7)
where f (x, y) =(OMo/Dx, OPo/Dy) is the generalized external force vector.
It is possible to introduce dynamics into the above static formulation of the snake. This may be done by introducing a dissipation energy functional which dissipates the kinetic energy during the snakes motion. The motion equation poses an initial boundary value problem which can be solved using efficient numerical techniques leading to a real time response from the snake, to applied forces. The initial value problem can be solved using explicit Euler integration. Let Y denote the Euler (time) step size, then the equations for time integration may be written as
P'=(K + Y) I(YPt  f Wt, Yt)) .
(2.8)
18
For compactness in representation, we can replace the snake with a Bsnake which is a snake constructed from Bsplines that have builtin smoothness. In this case, we simply replace p(s) by E'% CiBi(s), where C. are the control points of the chosen B3sp~line basis Bi and m are the number of control points. For more details on Bsnakes, we refer the reader to Marc et al. [47].
  Active Surface Models
A act ive/delormiable surface IS rel')r''eete(l b a vecto rvaluied paramietric represent atlonl p( a c) = (x(ii c).ja tl ) '(u7 c))T where (u, v) are domain parameters, and the parametric domain is the unit square [0, 1].
The deformation energy associated with the deformable surface, which simulates the thinplatemembrane material surface, is given by
E Op =o''~l+(P)2 2 2 o2 2
E21~t(p) W ]wi) + v) + 'W 2(522 + 2( ) 2 ) )du dv. (2.9)
Where Ei, canl be regarded as a controlledcontinuity spline defined by Terzopoulos [75]. As described before, and wi and W2 control the tension and rigidity of the spline respectively. Increasing wl tends to decrease the surface area of the material, while increasing W2 tends to make it smoother but less flexible. Fig. 2.1 demonstrates the effect of the value changes Of wi/W2 to the material.
Using variational principles, we can obtain the EulerLagrange equation for (2.9) and its corresponding stiffness matrix K. For the details of (2.9), we refer the reader to the Appendix A and Terzopoulos and Metaxas [78].
19
(a)(b
(C)
Figure 2.1. Active surface being pulled by a spring force showing the effect of various aw1 and W2 settings: (a) w, =O09,W2 =0.1, (b)'w, W 0.5, and (c) wi
001, W2 =0.99.
('14APTE R 3
SNAKE PEDALS: GEOMETRIC MODELS WITH PHYSICSBASED CONTROL
In this chapter we will first present definitions for the regular pedal curves/surfaces and illustrate the ideas via synthesized examples. \Ve will then present our proposed modeling scheme snakee pedal" model. Snake pedals are inherently geometric modeIs, but they can express local retaill, as well as cope with topology changes of shapes automatically. We will illustrate the power of this modeling scheme via some synthesized examples in this chapter as well as some model fitting examples in Chapter
3.1 Definition of Regular Pedal Curves
As shown in Fig. 3.1, let a be a p~lanlar curve, the pedal curve of a is defined as the locus of Points Onl the loot f of the p~erp~endicular from a fixed point p called the control point to a variable tangent of a. More specifically, let the planar curve
panar curve 7 pedal curve
pedal point p
Figure 3. 1. f is on the pedal curve of ae with respect to the pedal point p.
20
21
a be a lparameterized curve with domain parameter t, let ~3 be the pedal curve of a with respect to the fixed pedal point p, and let a(t) =g, 13 (t) =f at a particular value of parameter t, the projection of a (t)  p in the direction Ja' (t) must be '3
(t)  p, where a' (t) is the tangent of plane curve a (t), j : R2 ___ R2 is a linear map given by
J(X, Y) (X, Y) (3.1)
Note that applying the operation J on a vector (X, y) in R?2 Euclidean space results in a new vector (y,x.), and J1 can be geometrically interpreted as a rotation by7/ in a counterclockwise direction. \ecan thus formally define a regular pedal curve as follows [23]:
Definition :The pedal curve of a regular curve a :(c, d) __4 RJ2 with respect to a fixed (control) point p C R?2 is given by
pedal (t) = + (a JM ') J a(t Ja '(t). (3.2)
Where (c, d) is the parameter domain for t. The definitions of regular pedal curves and surfaces presented here may be found in any book on elementary differential geometry text [23, 38, 42]. Here we adopt the definition by Gray [23].
Some examples depicting the pedal curves of an ellipse for different positions of the pedal point are shown in Fig. 3.2. Usually the global or core shape depicts an overall size and orientation of an object, and local deformations depict fine detail in the local positions of an object. It is clear from Fig. 3.2 that the pedal curve is capable of exhibiting both global and local deformations. For example, the global
22
(a) (b) (C) (d)
Figure 3.2. Examples of pedal curves of an ellipse for different pedal point positions. Pedal points are shown by a (lot in each case.
shape of' a pedal curve of an ellipse is approximated by a new ellipse, whose long axis is almost the same as that of the original ellipse, while the short axis is much longer than that of the original ellipse. The new ellipse is pinched in the direction of the short axis to produce the pedal curve. The location and extent of this pinching operation are dlet erminedl bv the position of the pedal point. As shown in Fig. 3.2, if' the pedal point is located in the center of the original ellipse (Fig. 3.2 (a)), the p~iniching occurs around the (enter part of the pedal curve. As the pedal point moves to anl upper position (Fig. 3.2 (b)), the pinching position in the pedal curve moves up accordingly. Similarly, as the pedal point moves to the right of the center of the ellipse (Fig. 3.2 (c)), more pinching occurs on the right side in comparison to the left of the origin. This pinching operation constitutes the local deformation of the pedal curve, but the type of the local deformation that can be exhibited by a pedal curve is by no means limited to the pinchingtype or cusptype deformations. More complex local deformations can be realized in the pedal curve as shown later subsequently.
23
In summary, by moving the position of the pedal point, it is possible to synthesize a variety of global and local deformations in regular pedal curves/surfaces. The original curve ci(t) will henceforth be referred to as the generator for the pedal curve 03 (t) and the process of generating a pedal curve will be referred to as the pedaling operation, as depicted in Eq. (3.2). Note that the word "generator" here represents a curve in 2D or surface in 3D), which differs from its ordinary meaning in electricity industry.
.3.2 Snake Pedals: The Proposed Model
From the discussion in Section 3.1, we observed that when the pedal curves are generated from ellipses with respect to single fixed pedal points, the shapes and deformations exhibited by regular pedal curves are quite limited. To overcome this limitation, we can let the pedal points be different for each point of the generator when applying the pedaling operation to the generator, and thereby synthesize more general shapes of pedal curves/surfaces. We can let these pedal points be specified by another curve p(t), which can be represented by a active or deformable model (described in Chapter 2), that is, a standard snake or a Bsnake [31, 47] and then apply the pedaling operation to each point on the generator a(ti) with respect to the corresponding pedal point p(ti). Here we use p(t) to represent a parametric curve, p(ti) represents a particular point on the curve p(t) when t = t. Note that in this generalization of regular pedal curve definition, for different points on the generator, we apply the pedaling operation with respect to different control points, unlike in the regular pedal curves where only a single fixed pedal point is used.
24
Y a(
// snake edal  s nake pedal
snake
(t. enerator snake
W~)
x
Figure 3.3. (a) The process of'generating a snake pedal with anl ellipse as a generator.
(b) "snake pedal" controlled by the snake using an ellipse generator.
In this generalization process, we can employ a proper association mechanism to establish a correspondence between the points on the generator a(t) and the points on the curve p(t). This association scheine should keep) the continuitY property of the curve p(t) and be easY to impuIlemnut . In accordance with these two reqluiren wents, we developed a simplle correspondlence association mnechiauisrr as shown in Fig. 3.4. In this mechanism, for each point a(t2) on the generator , we draw a straight line from the center of the generator 0 to the point a(ti), and denote by p(ti) the intersection of the line 0Q(ti) and the curve p(t). When applying the pedaling operation to point a(ti), we carry out this operation with respect to point p(ti). The pedaling operation is performed for each of the points a (t3), J= 1, 2, . .., Al, on the generator resulting in a curve that we call the "snake pedal." Since 0ae(tj) and 0p(ti) are on the same radial line of the generator, we refer to this association scheme as the
25
P L(t4)
/P(ti) P 43) 42)
cL~j~fP~t~)P(t2)
snake
generator
Figure 3.4. Radialbased correspondence association scheme: each point ae(ti) on the generator is associated with a point pi on the snake.
radialbased association technique. It is obvious in this scheme that we use the same parameterization for the generator a(t) and curve p(t) with the same domain parameter t. It can be proved that snake pedal has the same continuity property as that of the curve p(t), given that the generator is differentiable, although this may not always be the case.
The generator can be either a parameterized or an implicit function representing a curve. In this chapter, we will focus on parameterized curves and in most applications, we use an ellipse as the generator in 2D. If the generator is an ellipse, we can
26
represent it in a parametric form by
a~) Cos 0 Sini ][ a, cos t ] T ml
sz 0cos 0 a2Sin t M
where a,, a2 are aspect ratio parameters and we use g, (a,, a2)T to denote these shape parameters, 0 is the rotation angle between the world coordinates and the object centered coordinates which is denoted by go = (0)', and g, (I, M2)T is the centroid of the generator in the world coordinates. We collect all the generator parameters into a vector g =(g', g', g')T (a,, a2, 0, I1, m2)T, and the vector g is referred to as the global shape parameter vector.
For the ellipse, the term Jci'(t) in the Eq. (3.2) of the pedal curve can be expressed as
Ja (t) II.a ~ i 0a O O (3.4)
172] [a, sint Cos 0+a2 COS tSinO 0 Note that 1jee'11 1 2 a~ 2S Z*12t+a~ 2cos2t, which is independent of the rotation angle 0. We will see that J&Z 2 is a very important quantity in the model fitting presented in later Chapters.
The curve p(t) on which the pedal points are located can be of any kind, which is not limited to be to a parameterized function. By using the correspondence association scheme described earlier, we can sample the curve at the same rate as that
27
of the generator so that every point on the generator a (t) can find a unique corresponding point on the snake p(t). The pedaling operation can then be applied to every point on the generator with respect to its corresponding point on the snake. We use a snake/Bsnake as this controlling curve and get a sequence of position vectors, pi1 (Xi, YJ), by sampling the snake curve. We collect all the snake position vectors in another vector, p =(p', p', ..., p' )", where M is the number of sampling points. Vector p is referred to as the local shape parameter vector. Vector p is the discretized version of continuous curve p(t).
In our implementation, since the pedal points are located onl a snake or Bsnake and the generalized pedal curve is created by applying the pedaling operation to the generator with respect to the points on the snake, we dub this generalized pedal curve as a snake pedal. With the radialbased correspondence associating scheme, a snake pedal is completely determined by the global shape parameter, namely, the generator parameter g, and the local shape parameter, namely, the snake position vector p. Fig. 3.3 demonstrates the process of creating the snake pedal curve.
With the correspondence association scheme, the generalized pedaling operation can create much broader classes of shapes with more deformations. We show some examples of snake pedal curves generated using snakes and an ellipse as the generator in Fig. 3.5. Note the variety of local deformations that can be generated using this modeling technique. We remind the reader that the snake pedal itself is a geometric model and that it is not directly responsive to the application of external forces unlike the standard snake models described in Chapter 2. Also, it is easy to synthesize
(a)
(d)
Figure 3.5. Examples
(b)
/
(c)
(e) Mf
of "snake pedals" using an ellipse generator and a snake.
snake pedals whose topology is different from that of the generator, as will be seen ill examples depicted subsequently. The uniqueness of our modeling scheme lies in the ability to generate snake like local deformations in a nonreactive geometric shape, namely, the pedal curve/surface.
3.3 Modified Snake Pedal Model
The pedal curve definition can be modified slightly by subtracting the second term from the first term in equation 3.2,
pedal(t) = Jp (at) 2 Ja ,(t). (3.5)
28
29
(b
(d)(e(f
Figure 3.6. Examples of the modified "snake pedals" using an ellipse generator and a snake.
This modification allows us to synthesize a larger class of shapes which are quite different from the shapes produced using Eq. (3.2). The key feature of using this modification is that the snake pedal curve allows for more local deformations including shrinkage and expansion of the snake pedal. The shrinkage and expansion behavior is controlled by the location of the snake. When the snake is located inside the generator, the snake pedal exhibits shrinkage behavior, while expansion behavior appears when the snake is outside the generator.
Fig. 3.6 depicts examples of shapes generated using this modification. Note the amount of local deformations allowed in the snake pedal curves using this modification. In Fig. 3.7, we compare the original and "modified" snake pedals using the same generators, same snakes and same correspondence associating schemes. It is
30
(a) (b)
Figure 3.7. Comparisons of original (red) and modified (blue) "snake pedals" using the same generators and same snakes in both case.
obvious that the "modified" snake pedals exhibit larger local deformations in both Fig.3.7 (a) and (b).
3.4 Handling Changes in Topology
Topology changes may be achieved by applying the pedaling operation to a single generator and several pedal points (or snakes). Fig. 3.8 depicts examples of the topology change realized in 2D. When the snake is located inside the generator, a single closed snake pedal curve is created (Fig. 3.8 (a)), when the snake is split into two parts with both parts located outside the generator, two separated closed snake pedal curves are obtained (Fig.3.8 (b)). In the second case, only part of the generator is needed to produce the two snake pedal curves. Further study shows that any number of closed snake pedals can be generated from one generator by using this procedure, if the position of snakes and parts of the generator are properly chosen. When using the correspondence associating scheme discussed previously, the number of closed snake pedal contours are determined by the number of closed snakes. More generally, the topological property of the snake pedal is determined by the topology of the controlling snake. Therefore, we can accomplish the automatic topological change
(b)
31
0 sk
0"
(d)
Figure 3.8. Examples of topology change. (a) depicts a single snake pedal curve (solid) obtained from an ellipse generator (dashed) and a single pedal point. (b) same as (a) except two distinct pedal points located outside the generator are used.
(c) same as (a) except several pedal points on a snake are used. (d) same as (c) except a discontinuous snake is used.
for the snake pedal model if the topology of the snake can be updated automatically in the data fitting applications. In Chapter 6, we will introduce a novel levelset based approach to realize the automatic topological change for the snake pedal model.
3.5 Pedal Surfaces
A regular pedal surface is the surface analog of the pedal curve. It is the locus of the points on the foot of the perpendicular from a fixed pedal point to a variable tangent plane of the surface. More precisely, it is defined as follows [23].
Definition :The pedal surface of a parameterized surface or a patch a U  R: with respect to a point p c R' i's defined by
peda~u, ) = , +(*u, V)  p) . au, X a,,
ped lau v) pvi 2 ' )aU x a,. (3.6)
32
Where U is an open subset of R2, ce(u, v) is a parameterized generator surface with parameters u and v and ct and ce, are tangent vectors in the u, v directions. For an ellipsoid generator, we have
a, cos ucos v MIn
a~~)R a2 cos usin v + Mn2 (3.7)
(:3 s'01 V T113
Where (I. I u( and a3 are aspect ratio parameters which can be collectedl inl gs
(o 1. ()T. \Ne IlSe gc = (Ml, Mn2, Mn3)T denote the centroid of the generator in the w( )rl(1 coordinate,. R. is thle rot ation inat rix between the world coordinates iio the object centered coordlinates. \\e( use (juat erions [6:3] to represent thle rotat ion inat rix 1 ecaiise the rotat ion mnatrix is easier to update t hroiigh the (jlaternionup 1Jdate anld remains orthogonal in the mean time. A quaternion [s, v] defines a rotation of the model from its reference position through an angle 0 =2 cos s around an axis aligned with vector v (v ',~)T .The rotation matrix corresponds to VS. v] is
1 2 + V2 ) 2(vIV2  SV3) 2(vl'c3 + S02) R 2(VIV2 + sv3) 1  2 (V2 + V2) 2(V2v3  SVI) (3.8)
2(VlV3  sv2 2(v2V3 + SV1) 1  2(v 1 2~ We collect the quaternion parameters into a new vector go (s, q T) =(S,vl, v2,V v)T. The ellipsoid generator is then completely represented by g =(gT' gT' gT)T. As in 2D case, g is referred to as the global shape parameter vector.
33
(a) (b)
(d)(e(f
Figure 3.9. Examples of pedal surfaces with fixed pedal points located inside ellipsoidal generator surfaces.
Fig. 3.9 depicts examples of pedal surfaces synthesized using an ellipsoidal generator surface and a fixed pedal point located inside the generator surface in each example. Note that change in topology is achieved by simply choosing an appropriate location of the pedal point as shown in Fig. 3.9. However, this type of topology change is different from what wve described in Fig. 3.8.
As in the 2D case, we can let the p~edal point vary for each point, on the generator surface. Thus we have the snake pedal surfaces in 3D whose shape can be controlled by snakes which are either curves or surfaces in 3D. Fig. 3.10 depicts some examples of snake pedal surfaces in 3D. They are synthesized by using an ellipsoid as a generator and a snake curve in 3D which lies inside the generator. The curve is shown outside of the generator adjacent to the snake pedal model for the clarity of visualization. Note that the snake pedal shape can be manipulated by applying external forces to the snake as illustrated in the Fig. 3.10 (a)(d).
34
(a)
(c)
(b)
(di)
Figure 3. 10. Examples of' shaping snake pedal surfaces with snakes. Each figure was generated by applying interactive mousebased forces to the snake.
For the 3D snake pedal examples obtained by ellipsoids and snake surfaces, we refer the reader to the next chapters.
3.5.1 Uniform Tessellation of the Model
From the definition (3.6) of the pedal surface, the points on the pedal surface must be located in the tangent of the generator for every point on the generator specified by (u, v). For an ellipsoid generator, the evenly spaced domain parameters u and 'v are often referred to as the "latitude" and "longitude." It should be noted that the snake pedal model exhibits nonuniformly tessellated surface in spite of the domain parameters u and v being uniformly tessellated /sampled. This is obvious from the nonlinear function defining the pedal surface in Eq. (3.6). For example, the tangent of the generator undergoes a rapid change (in orientation) as u  7/2, and
K
35
a slow change (in orientation) as u + 0. Therefore, for an evenly sampled generator and snake, when we use the association correspondence mechanism to produce a snake pedal surface, the snake pedal model will not be uniformly tessellated.
This uneven tessellation phenomenon is not noticeable when the aspect ratio parameters a,, a2 and a3 of the ellipsoid generator do not differ much from each other. But when the difference of the aspect ratio parameters increases, this nonuniformity becomes obvious, and leads to the distortion of the deformation properties of the snake pedal surfaces. In addition, numerical instabilities may be caused in the computation of derivatives for nonuniformly tessellated snake pedal. Hence a uniform tessellation scheme for pedal surfaces must be employed to overcome this problem. It is a natural way to first construct uniformly spaced tangent directions in the parameter direction u, then apply the same technique to achieve uniform spacing along the parameter direction v. We now present our scheme for the construction of uniform tessellation of snake pedal surfaces.
When constructing the uniform tessellation of an ellipsoid generator in the parameter direction u, we actually produce a uniformity in the tangent direction for an ellipse along a fixed longitude (v is fixed in this case). Since we do not need to consider the rotation and translation of the ellipse when constructing the uniform tessellation, for an arbitrary point (xO, YO) on an ellipse, the standard parametric representation is:
Io =0 a, cos u
(3.9)
yo a2 Sziu.
36
The straight line equation for the tangent of this ellipse on the point (xO, YO) is
X X O YY
2 2  (3.10)
a1 a2
Substituting Eq. (3.9) into (3.10), we obtain
a2COS U a2
Y _ x + .(3.11)
a, siriu sinu
The slope of this tangent line is: a2 Cos U which corresponds to an angle of 0 a, szn U'
tan'( cot ua). The angle 9 is not uniformly spaced when the variable u is evenly spaced, which results in a nonuniform tessellation of the model. The objective of the uniform tessellation is to modify the standard ellipse parametric formula so as to construct a uniformly tessellated 0 when u is evenly spaced. We change the standard ellipse equation (3.9) to: 1o = ,Cos f(u)
(3.12)
where f (a) is the function to lbe solved. The angle of the tangent directionl of' thle ellipse accordingly becomes 0 =tan'1 ( cot f (u)) To make this angle uniformly tessellated, we must set
d f 1a' a2Cof U C,(.3
(111 (11
37
where C is a constant. The boundary condition is:
{U= fu) 2' (3.14)
Basic algebraic manipulation of Eq. (3.13) results in the following ordinary differential equation (ODE) for f (u)
fP(u) C a ( in 2f(u) +a iC0.2 f (U)) (3.15)
a 1 a1I
where f'(u) is the derivative of f (u) with respect to u. This ordinary differential equation can be solved numerically. We use a third order RungeKutta method. The result is shown in Fig. 3.11 (a). It can be seen that when u is small (u  0), f (u) changes much faster than u, which tends to make up the fact that the tangent direction of the ellipse changes slowly in this region. While when u * ir/2. f (1) changes much slower than u, which is in an attempt to slow down the rapid change in tangent direction in this area. After changing u to f (u) in the ellipse expression, the change in tangent directions will become even. Carrying out the same procedure on the parameter v, we obtain the f (v) " v relationship as shown in Fig. 3.11 (b). Note that the retessellation scheme needs to be done whenever the parameters of the snake pedal are subject to large changes.
38
f(u) f(v)
o 0
U V
(a) (b)
Figure 3. 11. Uniform tessellation result: (a) f (u) " u relation, (b) f (v) "'v relation.
3.6 Summary: Salient Features of the Model
Compared with the hybrid models mention in Section 1.2.1, the snake pedal models are bestowed with several unique features which make them very powerful and attractive for use in shape synthesis and recovery applications. A list of these features are:
1. 'The snake pedlal model is inherently geoinet vic l)llt can exhibit bothI local and
global (leforinatiolls as would a physicsbased ilodlel.
2. The model can exhibit large global deformations such as bending, twisting, and
tapering without explicitly incorporating parameters for the same.
3. The modeling scheme allows (for the first time) the concept of a global shape
to be incorporated into the curve/surface evolution framework.
39
4. The model permits a levelset implementation and thereby seamlessly allows
for topological changes in shape during evolution.
CHAPTER 4
PROBABILISTIC FRAMEWORK AND MODEL LEARNING
In Chapter 3, we described our snake pedal model, which is fully characterized by the shape parameter vector q =(gT, pT)T with g, p being the global and local shape parameter vectors respectively. In many applications, the general shape, location and orientation of objects is often known a priori and this knowledge may be incorporated in the form of initial conditions, data constraints, model shape parameter constraints, or into the model fitting procedures. For example, in most medical images, the boundaries of the anatomical structures of interest are seldom well delineated from the background, either because the image quality is poor or the shape themselves gradually blend into the surroundings. In many medical imaging applications, it is necessary to use prior knowledge about the shape of interest by learning the shape from a human expert or from a shape atlas. Proper learning schemes must be developed to incorporate the human expert knowledge or the knowledge from an atlas. In addition, for automatic image interpretation, it is necessary to have a model that not only describes the size, shape, location and orientation of the target shape but that also permits expected variations in these characteristics.
To meet above requirements, a number of researchers have cast the model fitting process in a probabilistic framework to include prior knowledge of object shape by incorporating prior probability distributions on the shape variables to be estimated
40
41
[29, 69, 81, 86]. In Staib and Duncan [69], a deformable contour model was used on 2D echocardiograms and MR images to extract the left ventricle (LV) of the heart and the corpus callosum of the brain respectively. This closed contour model is parameterized using an elliptic Fourier decomposition and a priori shape information is included as a spatial probability expressed through the likelihood of each model parameter. The model parameter probability distributions are derived from a set of example object boundaries and serve to bias the contour model towards expected or more likely shapes.
Szekely et al. [71] have also developed Fourier parameterized models. In addition, they have added elasticity to their models to create "Fourier snakes" inl 2D and elastically deformable Fourier surface models in 3D. By using the Fourier paraineterization followed by a statistical analysis of a training set, they define mean organ models and their eigendeformations. An elastic fit of the mean model in the subspace of eigenmodes restricts possible deformations and finds an optimal match between the model surface and boundary candidates.
Cootes et al. [13] and Hill et al. [28] present a statistically based technique for building deformable shape templates and use these models to segment various anatomical structures from 2D and 3D) medical images. The statistical parameterization provides global shape constraints and allows the model to deform only in ways implied by1 the training set. These shape models represent anatomical structures bY sets of landmark points which are p~lacedl in the same way~ on an object boundary in each input image. For example, to extract the LV from echocardiograms, thley
42
choose points around the ventricle boundary, the nearby edge of the right ventricle, and the top of the left atrium. These points can be connected to form a deformable contour. By examining the statistics of training sets of handlabeled medical images, and using principal component analysis, a shape model is derived that described the average positions and the major modes of variation of the points. New shapes are generated using the mean shape and a weighted sum of the major modes of variation. Shape boundaries are then segmented using this "point distribution model" by examining a region around each model point to calculate the displacement required to mnove it towards the botindarv. these displacements are tlhen tused to update the shape p~aramfeter weights.
In this thesis. we cast the model fitting process in a Bayesian framework and incorporate prior distributions on the vector of shape parameters that are being estimated. By using a probabilistic approach rather than a deterministic approach, we can develop probabilistic models of data acquisition which closely match the characteristics of real image sensors, and more importantly, we can determine the uncertainty in our estimate, which can be used (if necessary) when integrating multiple sources of data/ information [73]. We now present the details of the Bayesian probabilistic framnework.
4.1 Prior Model in the Bayesian Framework
As discussed earlier, in many applications such as medical image analysis, the boundaries of anatomical structures are usually "fuzzy" because of the partial volume averaging effects and other effects. We must incorporate the prior information about
43
the shape of interest into the shape representation. Therefore, in the prior model, in addition to generic smoothness assumptions, we should also incorporate into the prior model the information about the possible shapes that an observed object can assume. We express the probability of a given configuration q through a Gibbs distribution [20], which assigns a low probability to configurations with large energy and vise versa. The configuration corresponding to , = 0 represents the rest state of the model. In the Gibbs distribution, the internal deformation energy, which is similar to the internal deformation energy in generalized splines described in Chapter 2, is transfered into a probabilistic distribution by:
1
where Z,, is a niormnalizing factor and the energy F,, (.) governs the dlefornmat ion of' the model away froiti its niat ural shape. There are two parts contained in this dleformation energy, one part corresponds to the local shape deformation, and the other part corresponds to the global deformation. For the local shape deformation, if we denote the rest shape of the local representation (snake positions) by p, which can be obtained by the scheme discussed in Section 4.4, the internal energy in a continuous form is given by
E(JIO 10  O~(((P  p))2 + (O(P  iP))2 ) +(42
(92 (P  1P) 2 (P P)2 a2+ P)
k~ 0,,2 ( 0o? ( 1 ) du dv,
44
with the discrete version
E(Iocal) !(p  p)K(p  iP), (4.3)
P 2
where K is the stiffness matrix for the membrane plus thin plate deformation energy imposed on the snake. Note that Eq. (4.2) and (2.9) have the similar form, if we set 15 0, Eq. (4.2) goes to (2.9). Eq. (2.9) is the special case of (4.2) where the local rest shape is zero.
For the global shape representation, if no additional information is made about the shape being modeled, the underlying global shape of the model, the generator, can take the mean value of the shapes that has been learned. By extending the above internal (leformat iOn eneI(rgy, we can include the global shape parameters of the model in our prior statistics in the similar way. The extendled energy cani be expressed as:
E E(11l) + j'(fihbol)
P  pp
 (p) K (pp
+ (Jg5 g )T C1 (g.s_ g) (4.4)
(gog)T C1' (go _ o)
1 ~(gcgc)T C1 (gc _)
Where g, go, and g, are the rest state of the aspect ratio vector, rotation vector, and translation vector of the generator discussed in Chapter 3. Cj'1, CO', and C,' represent the covariance matrices which are accumulated in the iterative training phase described in Section 4.4.
45
4.2 Sensor Model in the Bayesian Framework
The data acquisition (sensor) model in the Bayesian framework is based on linear measurements with Gaussian noise and given by:
where ZD is a constant, D is the data, and ED(') can be either the edgebased potential energy synthesized from image data, or the potential energy in the springs attached from 3D data points to the surface of the snake pedal, constraining it to conform to the observed data. These two types of sensor model energies are expressed as:
F>2 c(Di  xj)2 spring energy,
ED(D!II) (4.6)
[C Ic (G, * VI(x))112 image potential. In the first form of energy ED() c represents t he uncert ainty of thle sensor measureinents and( xi is the closest model p~oinit from the given data Di. In the second formn of' energy ED('), VI(x) is the gradient of the image intensity I(x), G,. is the Gaussian filter mask with width o,, and * represents the convolution operation.
Taking the partial derivative of ED with respect to q yields the external force:
OED (7
Oq
46
The external force pulls the model toward the image boundary or range data, thus molding the model into the desired shapes. This process is known as the model fitting process.
4.3 Maximum a Posterior
In the probabilistic framework, combining the prior model (4.1) arid sensor model (4.5) by' using Bayes rule, wve obtain the posterior clistribiit 10o1
p~qD) (Dlq)p(q) I
p1q(D) = ' p1(E(q)), (4.8)
where Z is a constant, andl the total energy is:
E (q) = E,(q) + ED(D . q). (4.9)
In the Bayesin framework, the objective is to compute the Maximum a Posterior, (MAP) estimate, that is, the value of q that maximizes the conditional probability p(qJD). This estimation provides the same result as finding the minimum energy configuration of the physically based model. However, compared with the deterministic model, the probabilistic framework has several advantages as described earlier.
4.4 Suiervised Learning Scheme
Once we have established the Bayesian statistical modeling scheme, we propose a supervised learning scheme to incorporate the information about the shape of iterest to the model. In this learning scheme, we divide thme shape recovery process ito a
47
training phase and an operational phase. In the training phase, prior knowledge about the shapes to be recovered from the image data is collected via a supervised learning technique. A training sequence of several iterations provides samples of correctly reconstructed shapes for specific classes of shapes, for example, the caudate nucleus and the hippocampus from the human brain. We use Pg denote the semiglobal part our model, which may correspond to the coarse level wavelet coefficients of the local shape parameter vector p. The parameter sets q(i) corresponding to these shapes represent the training sequence for the global part g and semiglobal part Pg of our prior model as shown in Fig. (4.1). With each iteration i, we improve estimates of the mean vector El and covariance matrix C through a simple equation for accumulating measurements as indicated in Fig. (4.1). Since we update only the statistics of the global parameters of the model, the local deformation remains constrained only by smoothness assumptions captured in the stiffness matrix K. Therefore, most of the flexibility of the locally deformable surface remains unchanged.
In the operational phase, the mean and covariance of the model state qj are used in initializing the model for fitting the model to a data set not included in the training phase. The main difference between training and operational phase is that, in the former, the model initialization for the fitting is far from the final desired fit because no prior information is available. Most existing schemes in literature that use a priori knowledge about anatomical shapes of interest manually construct their prior models [69, 73] which can be a very tedious process. The model fitting scheme that was described in this thesis in conjunction with the supervised learning mentioned
48
El(,)TCI(Surface fitting by maximizing p(q) (q  (i))Cli(q  Elp(qlD(i')) ,where current prior model p(q) is used q(i*)  final result
Figure 4.1 Supervised learning to incorporate information into the prior model
Iterative improvment of the mean and covariance of the global model Ej(i) = q(i*) + 1)
M2(i) +q jq'j) + 1) C(i) =M2(i) _(~Ti
49
earlier and discussed in detail by Vemuri and Radisavljevic [81] can be used for fast automatic construction of prior models. We will devote the next chapter to present the fast numerical algorithms for automatic model construction via fitting.
CHAPTER 3
NUMERICAL SOLUTIONS FOR MODEL FITTING WITH SNAKE PEDAL MODELS
In this chapter we will develop numerical techniques for fitting our snake pedal model to a sparse set of data points in 2D/3D as well as to image data for shape recovery. The objective of model fitting is to mold or sculpt the snake pedal model so that the model passes through the set of data points or conforms to the desired image boundaries as closely as possible, while maintaining a certain degree of smoothness. This fitting process is equivalent to finding the minimum of the total energy (4.9) of the model described in Chapter 4, which consists of two terms, namely, the internal strain energy Ep(.) corresponding to the snake deformation energy and the external energy ED(.). The energy minimization is posed as a nonlinear optimization problem and can be solved numerically. There are two types of parameters which need to be determined in this optimization process. One is the global shape parameters, namely, generator parameters g =(a, b, 0, inl, M2)T in 2D or g =(a,, a2, a3, V1, V21 V2) 8, I, Mn2, M3 )T in 3D); Another is the local shape parameters, that is, snake position p ={fpT} 1 with M being the number of sampling points on the snake. We collect the global and local shape parameters in a vector q =(gT, pT)T. In the model fitting, we need to find the optimum configuration of q so that the total energy in Eq. (4.9) is minimized. The estimation of both global
50
51
and local parameters in the model fitting makes the minimization problem highly nonlinear.
5.1 Related Work
A primary approach to solve energy minimization problems is the calculus of variations. In this approach, finite difference or finite element method is used to solve the EulerLagrange equations derived from the energy functional. This method can be regarded as gradient descent along the potential energy surface of the model. Vemuri et al. [35] used a gradient descent (GD) method that iteratively updates the shape parameter q,
q k+1 ~qk A aE (q) (5.1)
aq
where q k+1 and q k are the status of the shape parameter vector q at iteration k and k + 1 respectively. A is a diagonal matrix containing step sizes and consequently effecting the speed of convergence. Partial derivatives of the total energy ( q can aq
be represented by the sum of the partial derivatives of the prior model and data energies defined in Chapter 4. While the iterative equation (5.1) may eventually converge to the correct estimate, in practice it may be unacceptably slow. In addition, gradient descent does not guarantee finding the global minimum if the energy surface is not convex and the initial contour model is far from the target.
Amini et al. [2] use dynamic programming (DP) to carry out a more extensive search for a global minima. Although DP provides numerical stability, the stabilitN is achieved at a high cost in computational complexity. Williams and Shah [85]
52
proposed an alternative greedy algorithm to DP which drastically cuts numerical costs, however, it does not guarantee numerical stability.
Blake and Zisserman [6] propose an graduated nonconvexity algorithm, which is based on iterative approximation of the energy functional by a convex function. It can bridge low ridges in the potential field and eliminate termination in local valleys.
Poon et al. [57] and Grzeszczuk and Levin [25] minimize the energy of active models using simulated annealing (SA) method which is known to give global solutions and allows time incorporation of non differentiable constraints, however, deterininist ic algorithmis usuallY preferred due to their faster convergence.
5.2 Our Fast Numerical Algorithms
In the following sections, we present two algorithms for solving the energy minlimizmng problem efficiently. The first one involves a Preconditioned Nonlinear Conjugate Gradient (PNCG) method for estimating the global and local parameters of the model. The preconditioner we use for this case is the diagonal Hessian for the global as well as local parameters. In the second numerical scheme, we first solve for the global parameters using the LevenbergMarquardt (LM) method in the outer loop and employ time alternating directions implicit (ADI) scheme in the inner loop to solve for the local parameters. The first scheme is less sophisticated and easier to program, the second scheme is more sophisticated and relatively difficult to program. The next sections will be devoted to the presentation of the details of each of these schemes.
53
5.2.1 Preconditioned Nonlinear Conjugate Gradient Algorithm
In general, the Conjugate Gradient (CG) algorithm is well suited for solving quadratic optimization problems whereas, the nonlinear version can be used for finding the locally optimal solution of a nonlinear nonconvex function [661. In this work, we develop a Preconditioned Nonlinear Conjugate Gradient (PNCG) algorithm and apply it to the energy minimization discussed in the previous section.
Given any nvariate continuous function f (x) =f (X1, X2, ... , X,) with gradient f'(X) = [ f(X) ___f(X) ... Opff], we can use Nonlinear Conjugate Gradient (NCG) algorithm to minimize it. The residual ri f'(xi) indicates how far we are from the correct value that minimizes f (x), where Z denotes the number of iterations. We use do, d, ._..... dn1 to denote a set of orthogonal search directions, ai, and 01i the step sizes, and we choose a preconditioner M that approximates f"(x) and has the property that M1r is easy to compute for any vector r [36]. Let imaz be the maximum number or iterations and e be the residue threshold for the function f (x), the outline of the Preconditioned Nonlinear Conjugate Gradient method is as follows:
1. Let i =, ro  f'(xo).
2. Calculate a preconditioner M0 ;zz f"(xo).
3. do Mo 1 ro.
4. Find a, that minimize f (xi + aidi).
5. xi+1 ==xi aii
6. Calculate the preconditioner Mi+1 ; "x+)
54
7. rj~ H1i+)
8. or1r~ ~1rf
r T 1''ri
~ =max ( ij1 i+1rj+' ?'r1)
9. dj~j rjï¿½1 + f3j 1dj.
10. While (i' < ~mTor IfI < e) I= i+ 1, go to step 4.
In step 8, using the first expression for fij~ leads to the Fle tcher Reeves method, and using the second formula leads to the PolakRibi~re method. In step 4, the liiction f (x + (Id) is a)proximlately nminimnized by setting n~  11 Whee dTfIfd he,
f"(x) is the Hess iain atrix of the continuous function f (x).
24 a2a2f
H .. ~xd (5.2)
In our case, it cam be written as
02 E 02 E 02 E
H  dq2ql 0q2q2 Oq2Oqn (5.3)
____ 02 E 02 E
T&q1,q 6Jq1 d ..
55
Since it is expensive to compute the Hessian at each iteration, to perform an exact line search without computing f"(x), the Secant method is used leading to the following approximation
[f' (x)]Td(54
[f'(x + od)]Td  [ff(x)]Td (4
Typically, we will choose an arbitrary small nonzero number for 0' on the first Secant method iteration; on subsequent iterations we will choose x + od to be the value of x from the previous Secant method iteration. In other words, if we let ai denote the value of G~ calculated during Secant iteration i, then ci~l al.
In NCG, r(i) =f'(xi) and GramSchmidt conjugation of the residuals is used to compute the search directions. As for computing Q3, two choices are given in the above algorithm. The FletcherReeves technique converges for a good initial guess close to the desired optimum. PolakRibi~re scheme often converges much faster however, can cycle infinitely without converging in rare cases. Convergence in this case is guaranteed by choosing (0 = max (/3R, 0). In our implementation, we use the PolakRibere scheme.
When it's too expensive to compute the full Hessian matrix, it is reasonable to use the diagonal of the Hessian as a preconditioner. However, if the initial guess is too far from a local optimum, the diagonal elements may not all be positive. Since a preconditioner should be positive definite, nonpositive diagonal entries will violate this property and hence can not be allowed. A conservative approach in this case is to use the identity matrix as a preconditioner (M =I), that is, not to precondition.
56
In the following section, we describe a seminumerical approximation for computing the Hessian of the energy function.
5.2.2 Computation of the Approximate Hessian
In this section, we will introduce a seminumerical method for approximating the Hessian of the energy function by dissecting the structure of the Hessian matrix. In the Hessian matrix formula, it can be seen that the elements in one row of Hessian matrix, say row i can be obtained by taking the partial derivative with respect to gi of the vector (,9E)T, where gi is the ith entry of the vector g. Therefore, the ith row of the Hessian has the form:
a [ OEOE aEl a [ OF 1(5)
09gi O91 O992 (g, [ 'g I' (55
where n is total number of global parameters representing the snake pedal model. In 3D), the subblock corresponding to the snake pedal global aspect ratio parameters
g, [a, a2 a3 ]T has the following form:
O (Ox 0 O9x 0O
M o9 M a(a) Ox a) 0 Ox
57
19X 09x ax a9x
2 [aEaa E ,(.6)
Xax ax axla
where XI, X2, X3 are the world coordinates in which the model is described.
Similarly, the part corresponding to go [vi v2 A3 81T is,
a aE Qa QEaE ]
DVI [7V IaD2 T73J
a (ax __9 ax a (ax a a9x'
V [a'Val aV ax aV ax aV ax a a
QEQ Ea ax a ax a ax a ax
0(x 9 ax a ax a ax
ax ax ax ax a9x
arVaa~~ a [a 1a2 aV a l
a9 a9X _9_ a ax
only~ ~ ~ ~ ~ ~ ~ ~~~~F th folwn lmns edt ecmutd i(yVa2 )aVV3 PS~ an ou ~ bre ht the smlrpoeu.Sicmp eutat onlpey ilgl entrmesl dteesaed
if only those diagonal elements are required to be computed. From Eq. (4.6) in Section 4, we can obtain the analytical form for computing [aE aF QE and ~IDX aX2 DX3
58
[02E 02E 2El for and w
_T _j 3 owever, there is no analytical form fra and g we 9
use the forward difference to compute these partial derivatives numerically. Our experiments show that this seminumerical method is very efficient and convenient for the computation of the diagonal of Hessian matrix.
5.2.3 LM Method for Global Parameter Estimation
Although the Preconditioned Nonlinear Conjugate Gradient (PNCG) method performs quite well in the global parameter minimization, however, we found that the Levenberg Marquardt (LM) method has larger convergence range and converges faster than the pure PNCG algorithm. Therefore, in our second model fitting scheme, we utilize the LM algorithm to estimate the global shape parameters. In the local shape parameter estimation, by representing the stiffness matrix in a Schur complement formula [53], we use the alternating direction implicit method in the third model fitting scheme to achieve a 0(N) convergence rate in the local parameter estimation, where N is the number of nodes in the finite element discretization of the snake positions. Fig. (5.1) shows the framework of the second and third model fitting schemes [37, 85]. We now first introduce the Levenberg Marquardt method in the global shape parameter minimization.
The LM method essentially combines the Newton method and the gradient descent (GD) method in solving the nonlinear function minimization problem. When the current solution is far from the true solution, the LM method adopts the gradient descent~ (GD) algorithm. On the other hand, when the current solution is close to
59
While (Energy > E)
One iteration of
solving global parameters
Solve local parameters
completely , A
(a)
Figure 5.1. Big picture of numerical schemes for model fitting.
the true solution, the nonlinear function can be approximated by a quadratic form and the LM method switches to the Newton method.
NVe now adopt the LMI method to our model fitting problem. Assume that wve only needle to fit the snake pedal model to a set of 3D) data points, then the external energy used here becomes the springbased deformation energy ED (D, g, p) introduced in Section 4. More specifically, to fit the snake pedal to a given data set of 7n points {Dj}j'j in 3D, we minimize the sum of the squared distances from each data point Di to the corresponding closest point xi on the pedal curve/surface. Let h be the function that returns the closest point on the snake pedal from a given data point Di, that is, xi hi (g, pi, ti). Let {pi}I', be a set of M discrete points along the snake p(t), {x}My be the corresponding points on the snake pedal curve. Let D ={fDT}1T1 and p {pT}M1, the quantity to be minimized for the model fitting
60
can be written as
2i=1
where xi hi (g, pi, ti),I and the magnitude of parameter c is related to the uncertainty of data measurements.
In LM, the gradient of ED(D3, g, p) with respect to generator parameters g
(9192,) =(aj,a2,a3,V,v2,v3,S,m,m2,m3 )T in 3D), and g =(g12) (a,a2,0, M1,M2 )T in 2D) has components:
OED n Zc (Di  hig i i)ah g i i (5.8)
Ogk 10.9k
If we assume that there is a onetoone correspondence between the points on the generator and the pedal snake then the partial derivative OEacan be written as:
OED_ c (Dk  hk(g, Pk7tk)) Ak 9k, (5.9)
09k 19
Taking an additional partial derivative gives ak( which is needed in the LM iteration. For the derivative of the energy required in the LM algorithm, we use the chain rule to get L a= ( Fg E O 9h)* More specifically, if we denote each component hi of hi as hi (hil, hi2 )T, then in 2D)
OEc [Dii hi i i r
(5.10)
61
If we define:
Ok ED 02 ED
19gk ' Ok 091e
(5.11)
using the Newton method to minimize the energy function gives
(5.12)
On the other hand, using the gradient descent method gives
corist. x 31 = il,
Ace,,
(5.13)
where we choose the constant to be 1~ in Eq. (5. 13). This choice of step size in Eq. (5.13) leads to a preconditioned gradient descent method with a diagonal Hessian as the preconditioner. Furthermore, we define
a ajj (I + A),
(5.14)
aik =aJ/ j k),
(5.15)
and combining Eq. (5.12) and (5.13) yields
Z a'16g, k
(5.16)
When A is very large, Eq. (5.16) becomes identical to (5.13); and as A approaches to zero, Eq. (5.16) becomes the same as Eq. (5.12).
E ak169k = Ok 
62
From above discussion, given an initial guess for the set of parameters g, the LM method for the global parameter estimation is described as follows:
1. Compute ED (D, g, p),
2. Pick a suitable value for A, for example, A =0.001,
3. (t) Solve thle linear Eq. (5.16) for 6g and evaluate ED (D, g + 6g, p).
4. If ED(D, g + 6g, p) > ED(D, g, p), increase A by a factor of 10 (or any other
substantial factor) and go back to (t),
5. If ED(D, g +Sg, p) < ED(D, p, q), decrease A by a factor of 10, update the
trial solution g * g + 6g, and go back to (M)I until the energy ED (D, g, p) is
less than a predefined threshold.
After the estimation of the global parameters, we solve the local shape parameters in the inner 1001) of the paranieter estimation using the Preconditioned Conjugate Gradient method or the alternating direc;tionlal implicit method. .5.2.4 ADI Mlethod for Local Parameter Estimation
When solving for the local parameters in the inner loop for fixed global parameters, the quantity to be minimized has a quadratic form and its minimization leads to solving a linear system, namely, Kx =b, where K is the stiffness matrix, x is the local parameter vector, and b = is the external force. Note that we dp
uise x to represent the local shape paramieter vector and~ b to represent the external force, (they were rep~resentedl by p and f in previous chapters.) which is the standard
63
method to represent unknowns and righthand side in linear algebra. Because of the special structure of the matrix K in this linear system, we can use the alternating direction implicit (ADI) method in the local shape parameter estimation. We can prove that the ADI takes only 0(N) time to converge with soft data constraints where N is the number of nodes in the rectangular finite element discretization of the snake in our model. Note that in soft data constraints, the model is not required to interpolate the data but is expected to approximate it.
We now present the structure of the stiffness matrix K. In the snake pedal model, the controlling snake p is a controlled continuity spline under the influence of image forces or other external constraint forces. In our model, the spline deformation energy imposed on the snake is a membrane plus thin plate energy given in continuous form by
a (() 2 () 2 W2 P 2 2 )2+(0a2P)2)
E, (p) f1 + +p W2 2) +2( 2 ) du dv, (5.17)
Oi OI v ODy2O
where w, and W2 are the tension and rigidity factors respectively.
The external energy used can be the same as the energy ED (D, q, p) discussed in Section 5.2.3. In this case, the external force impose on the snake is given by
b OpD, E Oh p c [Dii  hil Di2  hi2] . h< (5.18)
64
South Pole
Total number of nodes:
nV
Denote : N = n Xn
NotU Pl
(a)
Figure 5.2. Model tessellation in the model coordinates.
w h e r e h hA in 2n  n
OpUA [ Lzf1f2 ?1f2 (5.19)
with nj, ni2 defined as before.
Notice that 9ED1,9hj is the external force in a standard snake, but in our case, the external force is 9ED/aPi, which is related to the force in a standard snake via a Jacobian matrix Ohi/O9p, in Eq. (5.19).
We then discretize the model in model coordinates using finite elements. Geometrically, the snake pedal surfaces are closed surfaces in space whose intrinsic coordinates u =:(u, v) are defined on a domain Q. The positions of points on the model in an inertial frame of reference are given by a vectorvalued function pedal(u) =(XI(U),X2(U),X3(U)). The tessellation of u =(u,v) in intrinsic coordinates is shown in Fig. 5.2. We use a combination of quadrilateral and triangular
65
elements in the discretization. In the inner quadrilateral mesh, we have (n x ni) =N nodes and 2 additional nodes are used in representing the north and south pole leading to a total of (N + 2) nodes. Due to the presence of the north and south poles, after the discretization of the model in intrinsic coordinates u using finite elements, the stiffness matrix K C Rt(N 2)(N 2) in our problem can be written as a (2 x 2) block matrix:
K = 21 (5.20)
where K11 C fNxN contains coefficients which corresponds to the inner quadrilateral mesh and K22 E R 212 is a diagonal block matrix containing the coefficients corresponding to the north and south poles, and K12 =K'1 contains the coefficients of interaction between the poles and other nodes.
To solve the linear system Kx =b for the local parameter estimation, we employ the following Schur complement formula [53]
K'b Ki ilK2_KlKl1lllSlb, (5.21)
where
S =K22  K21K711K12. 5.2
(5.22)
66
From Eq. (5.21), it is evident that to obtain the solution of K'b, we need the solution of several Klbs. Further analysis (see Appendix A) shows that the component K1, of the stiffness matrix K can be expressed as the sum of two Kronecker products:
KI1 = wj(AoB BoA)+W2A0DA
=AoCiCoA, (5.23)
where C =wiB + (W2/2)A with A and B being (ni x ni) tridiagonal and symmetric positive definite (SPD) matrices of the form
1i2 1
1 2 1 1 4 1
A E RX7L; B *. * ..
1 2 1 1 4 1
1 1 2
(5.24)
The Kronecker tensor product or Kronecker product of C and A, denoted by C 0 A, is the matrix of order N x N (where N =n x ni) and given by
67
C21A C22A ... C2nA
COA= (5.25)
When computing with Kronecker products, it is computationally convenient for us to represent vectors using matrices. When solving (C 09 A)x, we represent the vector x of order n' by the mnatrix X ={XIk I of order n x n defined by
1ki,= Xk ,() (5.26)
It can be shown that this representation can lead to efficient procedures for computing (C 0 A)x and solving (C 0 A)x =b, respectively.
LEMMA 1
Let C { Ckl}1. A ={fUkl I and1( X ={ik! I be matrices of order n x 7?. then the 'n x n. matrix representation of, 42 x 1 vector (C 0 A)x is given by (C 0 A)x 4(C(AX)T)T.
Using this representation, the linear system Kjjx =b can be written as:
CXA T + AXCT =D, (.7
(5.27)
where D, X are the n x n representation of the n' vector b and x. Since A and C are also SPD matrices, Eq. (5.27) can simply be written as
CXA +AXC =D. (5.28)
After some simple linear algebra operations, it follows that above equation is equivalent to
A'X +XA' = D', (5.29)
with A' :A'C and D' =A1DA'.
Eq. (5.29) is the standard Lyapunov matrix equation [43]. The ADI method can be applied to solve this matrix equation and has the following steps in each iteration
(A'+ pjI) X, D'  Xj,1(A' pI) (5.30)
Xj (A' + pj1) D'  (A'  pjI) Xj .1 (5.31)
For each iteration in ADI, we need to solve two equations, namely, Eq. (5.30) and (5.31), and therefore Xjiis introduced as an intermediate variable. Since A' is a SPD matrix, we can use Cholesky decomposition for the tridiagonal SPD matrix (A' pjl) to compute the update solution. This only takes 0(N) time per iteration. The parameters p2 are called ADI parameters and J is the number of iterations needed. They are chosen according to the method described by Lai and Vemutri [36].
69
In order to fit the snake pedal curve to salient features in an image, we use the following equation as the energy to be minimized:
E = H1(G, * VI)112 11(G, * VI(x, y)11'. (5.32)
Where (x, y) is on the snake pedal curve, G,. is a Gaussian with standard derivation aT, andl * dlenotes the convolution operator. More sophisticated boundary finding energies can be synthesized and we refer the reader to Kass ct U1. [311. By miniinuzing the above energy, the nodel is inadle to fit to points ofigh gradlient, in the image. The procedure to minimize above energy is similar to that of fitting of snake pedal (llrves/slirfaces to dlata points.
CHAPTER 6
HYBRID GEOMETRIC ACTIVE MODELS
Since the inception of snakes in the vision/graphics community by Kass et al. [31] and Terzopoulos et al. [76], these elastically deformable contours /surfaces have been widely used for a variety of applications, such as edge detection, segmentation, motion tracking, and shape modeling. The classical approach to object shape recovery using the snakes is based on deforming an initial configuration of the snake P0 towards the boundary of the object to be detected. As described in Chapter 2, snakes are energy minimizing splines, their deformation under constraints is achieved iby minimizing a functional, that can be regarded as the bending energy stored in a thin flexible beam/rod or a stretchable string subject to loading.
There are several problems associated with this approach, such as initializat ions, convergence to the local minimum, and the specification of the physical or elasticity parameters. Moreover, this energy model requires that the topology of the shape to be estimated be known a prior. However, in many applications, such as motion tracking, it is not always possible to specify the topology of an object prior to its recovery. For example, during the cell mitosis, the close contours of cell boundary change connectivity or split, therefore undergoing a topological transformation. Hence it is necessary to develop a model which has the ability to split and merge freely during the
70
7 1
shape recovery without making the a prior assumption about the object's topology in these applications.
Recently, Caselles et al. [91 and Malladi et al. [45] proposed a novel geometric
model of active contours dubbed geodesic/geometric active contours or geodesic/geometric snakes. These models are based on the the theory of curve evolution and geometric flows. Automatic changes in topology can be handled in a natural way in this modeling technique when we implement the curve evolution using the levelset embedding schemes. In essence, this model is and implicit surfaceevolution model where the surface propagates with a velocity which is directly related to the curvature of the evolving curve which is embedded in the surface. Unlike the classical energy models, these models are not the result of minimizing an energy functional. In this framework, boundary detection can be considered equivalent, to finding a curve of minimal weighted length. More recently, several researchers independently demonstrated the unity of the curve evolution approaches and the classical energy minimization methods [10, 32, 64].
The challenge in the geometric active contour evolution is to devise numerical schemes for the equations of the evolving curve, which will accurately approximate solution that can be very unstable due to formation of shocks or discontinuities. Osher and Sethian [52] developed special numerical techniques based on upwind finite difference schemes leading to a very stable solution. They introduced levelset techniques to seamlessly handle any topological changes of the curve/surface that may
721
occur during evolution. The equation of motion for the embedded curve is reminiscent of an initial value "Hamilton Jacobi" equation with a parabolic righthand side and is closely related to a viscous hyperbolic conservation law. For details on the theory of curve/surface evolution and its levelset implementation, we refer the reader to the relative papers [9, 10, 33, 46, 51, 45, 65, 74].
6.1 Overview of Hybrid Geometric Active Models
In many computer vision applications, characterizing the global shape of an object is very useful, for example, in object recognition. Traditional geometric active contour/surface models do not possess the capability to characterize the global shape of an object. They are however useful in characterizing the local shape details in an object. In this section, we introduce a novel concept of a global/core model into the PDEbased curve evolution framework by embedding the snake pedal model into a levelset framework. Instead of characterizing an object boundary by the position of every point on the boundary, our proposed model, referred to as the hybrid geometric active model, now describes an object as a combination of a global/core shape, for example, circle, superellipse, and so on. and a variable offset defined with respect to the global shape. The variable offsets are controlled by this global shape and an evolving curvethe controlling snake in the snake pedal model. Augmentation of the global shape/core in the model representation is very useful in object recognition and image indexing applications. In these applications, the global description, such as center, size, and orientation, of a group of objects or structures needs to be retrieved initially which is then followed by a retrieval of the local description of each
73
object (if necessary). For the model to recover the object boundary, we introduce a robust and efficient numerical method which consist of a global plus local shape estimation technique in the model fitting. We use the Levenberg Marquardt (LM) described in Chapter 5 for estimating the global shape and a combination of upwind and minmod finite difference schemes in a levelset framework for estimating the local shape. The LM method has a broad convergence range, converges to the local optimum very quickly, and therefore is well suited for the global shape recovery. The hybrid geometric active contour/surface model retains all the advantages of traditional geometric active models (for example. topology change, ability to model complex geometries and~ amenability to stable numerical implement at ion) and has the added ability/ advant age of being able to compactly represent global shape of all object.
\Ve first present the theoretical foundation of the traditional curve evolution in Section 6.2, followed by the levelset implementation framework of the PDEbased curve evolution in Section 6.3. We then introduce our hybrid geometric active models in Section 6.4. The numerical issues of the model fitting process will be discussed in Chapter 7, followed by the implementation results in Chapter 8.
6.2 Curve Evolution Theory
The mathematical foundation of the geodesic or geometric active contour model is based on the Euclidean curve shortening theory, which guides the curve evolution in Euclidean space [15]. Let p =(p.,py)T denote an arbitrary point on the geometric active contour, the family of planar curves flow according to the geometric heat
14
equation or curve shortening equation
t kN, (6.1)
where k is the curvature and N is the inward unit normal. When the curve evolves according to (6.1), the Euclidean perimeter shrinks as quickly as possible. We demonstrate this concept in the following.
Let P(O) be a smooth, closed initial curve in Euclidean plane R'2, let P(t) be the oneparameter family of curves generated by moving P(O). Let p(s', t) be the position vector which parameterizes P(t) by s', where 0 < s' < 1 and t is the time factor. For closed curves, we assume that p(O, t) =p(l, t) for any time t and a similar condition is posted on the first derivatives. We define the length functional:
C(t) =JII P ds', (6.2)
where .denotes the standard 2normi. Taking the first variation [19] andl using integration by parts in Eq. (6.2), we obtain ap 02p
~as'H
(6.3)
Op1
a I jSj> ds',
as'as
75
where <, > is the inner product of two vectors. Note that the arc length ds can be expressed by
ds flp' ds', (6.4)
using the definition of curvature k and normal N [7, 8, 49], the integral in (6.3) becomes
<2(L 1Lt < kN > ds. (6.5)
Setting equation (6.5) to zero, we obtain the direction in which 12(t) decreases most rapidly as
= kN.
at
It has been demonstrated [16, 17, 18, 24, 601 that when simple closed curves evolve according to (6.1) without developing singularities, they will converge to "round" points. High curvature regions will be smoothed out during the evolution. We can easily understand this phenomenon by further examining Eq. (6.1). In (6.1), the planar curve evolves in the direction of its inward normal, and the rate of evolution is proportional to its local curvature. In locations where the curve has a larger curvature, the rate of change is large, thus the curve evolves faster in the high curvature points. For the same reason, the curve evolves relatively slower at the low curvature points. Therefore, the portions of the curve in high curvature regions undergo more changes toward the rest state, while the portions in the low curvature regions remain smooth until all portions of curve reach the rest statea circle.
76
Since the curvature k determines the rate of change in (6.1), it is regarded as the speed function of the curve in the evolution. More general speed functions can be allowed in the curve evolution. Let F(k) denote the speed function, the curve then evolves according to
 =F (k) N. (6.6)
In general, F(k) can be split into two terms:
F(k) F4j + F;, (6.7)
where the first terni F4 is referred to as the advection term and is independent of the geometry of' the curve. The curve uniformly expands or contracts with speed FA depending on its sign, the behavior of the curve is analogous to that of balloon model [12] with an inflation/deflation force. The second term FG is referred to as the diffusion term and is dependent on the geometry of the curve, such as the local curvature. The diffusion term smoothes out the high curvature regions of the curve and has the same regularization effect on the curve as the internal deformation energy term in the generalized splines described in Chapter 2.
When implementing the curve evolution equations such as (6.6), there are a number of numerical considerations. One numerical approach is first taking the above Lagrangian description of the problem, producing equations of motion for the position vector p(s', t), then discretizing the parameterization with a set of discrete markers lying on the moving curve. These discrete markers are updated in time
77
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 [61]. First, small errors in the computed marker positions are tremendously amplified by the curvature term, and calculations are unstable unless an extremely small step is employed. Secondly, singularities may develop (when FA :A 0) in the propagating curve even with a smooth initial curve. A typical example of this phenomenon is when the smoothing curvature (viscous) term is absent, namely, F(k) =FA in (6.6), and an entropy condition must be imposed to extract the correct weak solution. Thirdly, it is difficult to handle the topological changes if the evolving curve breaks and merges. Finally, significant bookkeeping problems occur in the extension of this technique to three dimensions.
6.3 Traditional LevelSet Amproach
Oslier and Sethian [52, 61] proposed an algorithm based on the theory of viscosity solutions to PDEs . The algorithm provides a stable numerical solution for curve/surface evolution. In their approach, the curve is first embedded in a two dimensional surface, the equations of motion are solved using finite difference techniques and hyperbolic conservation laws.
6.3.1 Embedding Mechanism
The embedding step can be carried out by representing the curve P(t) by the zero level set of a smooth and Lipschitz continuous function 0 : 2X [0, T) > RJ. As an illustration of this approach, we consider the example of an expanding circle. Let the initial curve P at t =0 be a circle in the xyplane (Fig. 6.1 (a)). We imagine
Iy
(a)
x
P(0
y
(C)
x
P; (t)
Z= (P (x,y,t=O)
(b)    (P C
y
P(O) Level Set (~
x
iZ=( (x,y,t)
(d) . P=C
   y
 IP(t) Level Set 4)=0
x
y
(e)
P (t+A t)

x
Z = 0 (X,y, t+ At)
P=C
iP(tit) = Level Set 4)=0
x
Figure 6.1. Level set formulation of equations of motion  (a) and (b) show the curve P and the surface O(x, y) at t = 0, (c) and (d) show the curve P and the corresponding surface O(x, y) at time t. (e) and (f) show the curve P and surface O(x, y) at time t + At. The curve P has changed topology in going from (c) to (e).
78
79
that the circle is the level set {0 =O0 of an initial surface z = (x, y, t =0) in R' (Fig. 6.1(b)). We can then match a oneparameter family of moving curves 2(t) with a oneparameter family of moving surfaces in such a way that the level set to = 0} always yields the moving curve (Fig. 6.1 (c) and (d)).
6.3.2 Equation of Motion
In the general case, let P(0) be a closed, nonintersecting, N  1 dimensional hypersurface. Let 0(p, t), p E RN1, be a scalar function such that 0(p, 0) = d(p), where d(p) is the signed distance from p to the hypersuface P(0). Assume that 0 is negative in the interior and positive in the exterior of the zero level set.
Let~ P(0) lbe a planar curve, represented as the zero level set of 6(p, t), that is.
,p(t) C SJR2 : 0 (p, t) =0. (6.8)
By differentiating (6.8) with respect to t we obtain:
V (')(p, t)  O + Ot (p, 0t) 0. (6.9)
Note that the gradient V70(p, t) is normal to the N  1 dimensional level set passing through p in general. When N =3, for the zero level set, the following relation holds:
O  N, (6.10)
ITO5
80
where N is the curve normal as defined before. In this equation, the left hand side is in terms of the surface 4), while the right hand side is in terms of a property of the curve P. Combining Eq. (6.1), (6.10), and (6.9) yields:
Ot=F(k)11V4)fl. (6.11)
The curve P, evolving by Eq. (6.6), is obtained by the taking the zero level set of the function 4), which evolves by (6.11). This scheme is referred to as Eulerian formulation for curve evolution, because it is written in terms of a fixed coordinate system. Also, for an implicitly defined curve, the curvature can be computed as M~organ [49]
k = V VO(6.12) or explicitly as
k X4)  20).0yo4) _ Y (6.13)
(o2~ 43/2
Eq. (6. 11), together with the initial condition, 4)(p, 0) d(p), is the governing equation for the curve evolution in a levelset framework. Since this equation resembles the HamiltonJacobi equation, it is referred to as the "HamiltonJacobi" formulation [52].
To attract the active contour to the desired image features, the speed function in (6. 11) is multiplied by a quantity:
g~x~y) 1 + IIVG, * 11(6.14
81
where I is the greyscale image and G, is the Gaussian smoothing filter with width o,and n =1 or 2 depends on the application [44, 83]. The function g(x, y) depends on the given image and is used as a "stopping term." For example, when the curve is close to an edge, g(x, y) becomes very small and retards the curve evolution. A modification of (6.11) based on the use of Riemanian metric was reported in Caselles et al. and Kichenassamy et al. [10, 32], it involves the minimization of a weighted Euclidean length and takes the form:
Ot=g(x, y) F(k) IVOI5 + VO  Vg, (6.15)
Note that for g as in (6.14), Vg behaves like a "doublet" near an edge. The effect of Vg is to attract the evolving contour as it approaches an edge, and to push the contour back if it passes the edge. The inclusion of Vq5. Vg term allows the use of larger inflation forces, resulting in features being extracted in relatively fewer time steps.
6.3.3 Advantages of the Levelset Approach
By using the level set approach discussed above, topological changes can be handled naturally. This advantage lies in the fact that in the curve evolution, the higher dimensional function 0(p, t) always remains a function, while its level surface {O = 0}, which corresponds to the evolving curve p(t), needs not to be simply connected. The moving curve can change its topology or form sharp corners easily. As shown in Fig. 6.1 (f), 0(p, t) remains a simple continuous saddlelike function. while
82
its level set {0 0}, the moving curve, splits into two separated closed contours, thus undergoes a topological change from its initial configuration shown in Fig. 6.1 (a). Another advantage of the level set approach lies in the fact that the geometric and differential properties of p(t) can be computed from 0. For example, the curvature of p(t) can be computed in terms of the derivatives of 0 expressed in Eq. (6.13). In addition, the level set approach can be easily extended to higher dimensions and appropriate expressions can be obtained for the mean curvature and the Gaussian curvature. Moreover, the governing equation (6.11) is reminiscent of an initial value HamiltonJacobi equation with viscosity, and can be solved using the stable, entropysatisfying finite difference schemes in accordance with the hyperbolic conservation laws [40, 52, 62].
6.4 Our Model: Hybrid Geometric Active Contours
As described in the former sections, the level set approach provides an elegant way to handle topological changes with active models, however, traditional PDEbased approaches do not possess a mechanism to characterize the global shape of an object. In this section, we describe how the level set formulation of the curve evolution can be applied to our snake pedal model to realize topological changes (when necessary) as well as capture the global shape representation of an object.
6.4.1 Big Picture
In our approach, to build some amount of regularization on the snake pedal, we impose the regularization constraint via Euclidean arclength minimization on the snake pedal. This would lead to the standard geodesic active contour. Let us
S:3
consider a snake pedal denoted by Pe(t), with its position denoted by p, =eI, Pe2}, then the standard curve evolution formula for the snake pedal can be written as
aPe
at=F(k,)Ne, (6.16)
where ke and Ne are the curvature and normal of the snake pedal respectively. We would like to know how the "snake" would evolve if the snake pedal were evolving as a function of its local curvature. Indeed, we will not evolve the snake pedal curves directly, instead, we will first solve the "snake" position under the constraints that the arclength of the snake pedal is minimized, and the pedal curves can then be determined by the pedaling operation defined in Chapter 3, given the position of the generator. The problem therefore can be solved by the following procedure:
1. Derive the curve evolution equation for the "snake" P(t) by minimizing the
arclength of the snake pedal Pe(t);
2. Embed the evolving curve P(t) in a higher dimensional surface 01 and formulate
the equation of motion for 01
3. Solve the equation of motion for 01 using proper numnerical techniques;
4. Determine the snake pedal curve from the evolving snake via thle pedaling
operation, given the generator.
We should remind the reader that in thle above procedure, the snake pedal curve also evolves, but it evolves as a standard geometric active contour embedded
z
x
2 (x,y,t) = Level Set
0=0
t) =Level Set 2~ = 0
(a)
Figure 6.2. Relation between 2(t), 2e(t), Pi, and 02: 2P(t) and 2P,(t) are levelsets of the higher dimensional surface 01 and 02 respectively. They are related by the pedaling operation.
in another higher dimensional surface 02. The equation of motion for q02 is the same as the standard one as shown in (6.11), but we do not need to solve this equation directly for function 02 0k2 can be regarded as being implicit in the procedure, we do not obtain the levelset Of 02, namely, the snake pedal curve Pe(t), from evaluating 0f2 and determining the points 02 =0, instead, we first evaluate 01 and determine its zero set, then we apply the pedaling operation to obtain the snake pedal Pe(t). In addition, even though we never solve 02, 02 is very useful in the derivation of the governing PDE for the function 01. The development of the relationship between the gradient of 01 and 02 is a crucial step, we will derive this relationship later in Section 6.4.3. Note that the "snake" model discussed here is not a classical energyminimizing model any, more, since its deformation is not obtained by minimizing some
84
Figure 6.3. The hybrid geometric active model can capture the "global orientation"~ of clustered objects.
internal deformation energy, but we will still uise the name "snake" to represent the controlling active contours in the snake pedal model. Fig. 6.2 depicts this important relationship between P (t), Pe(t), 01, and 02 *
Note that an important feature of this modeling technique is the incorporation of a global parameterized shape, namely, the generator, into the curve evolution framework. This global parameterized shape can be very useful in applications involving recognition of clustered objects. For example, in Fig. (6.3), there are three clusters of objects in 2D. Each cluster has its center, size and orientation, which can be represented by an ellipse. The snake pedal model can capture this "global orientation" of the cluster via the global part of the model, namely, the generator, while recovering every single object in the cluster. Since a geometric active contour is used to represent the controlling "snake" in the model, and also the model can capture the
86
"global orientation" of a group of objects, this snake pedal model is referred to as a hybrid geometric active contour model. Detailed discussion on the related concepts will be given in the following sections, starting with the derivation of the relationship between the snake and the snake pedal evolution.
6.4.2 Relation Between the Snake and Snake Pedal Curve Evolution
Consider a snake p =(PI, P2)T and a generator a =(a,, a2 )T with normal
a (nl,n12)T, where p and a are related by the radiusbased correspondence associating scheme described in Chapter 3, we obtain the corresponding snake pedal Pe =(Pei, Pe2 )T via the following pedaling operation:
ap*
(6.17)
here we use the modified pedaling operation given in Chapter 3 for the purpose of explanation. Rewriting the Eq. (6.17) in component form, we obtain
[Pei
P e2
I
I
+ 1 + 2 fl2122
iiLA!2 Fe Pt 1i a n+a2 n2 nI (6.18)
2 Jn + n2 2 11+ 2 n
or simpllly,
P, JA p  b
(6.19)
where
I+4 n 2 JA1~22
n n
1+ 2 , n,+ n2
(6.20)
871
and
n2 + n2 [2
Taking the inverse Of JA, denoted by JB, we have 2+2n2 n
JB = A  12 72n2n 2ll+n2 (6.22)
We should remind the reader that by using the association scheme described in Chapter 3, the forms Of JB and b do not change over time, but their contents do change. But, for the following reason, we can still treat .JB as a constant matrix and b as a constant vector with respect to time when evolving p(t) or pe(t). First, the generator does not evolve over time, we only estimate the generator parameters iteratively; Second, the changes Of JB and b are very small compared with that of p or p,. Therefore, Eq. (6.19) can be explicitly written as a function of time:
Pe (t) =JA P (t)  b. (6.23)
Taking partial derivative on both sides of Eq. (6.23) with respect to time t yields
&pe (t) p MA 6.4
at JAat (.4
or
a9p(t) &pe (t)
at = Bat (6.25)
00
Eq. (6.24) and (6.25) reveal a very crucial relationship between the snake and snake pedal evolution: The evolution of two curves are linearly related by a Jacobian matrix. We remind the reader that the relationship between the snake and snake pedal is not linear, but their evolutions over time are related by a linear transformation! This linear transformation between the evolution of the snake and the evolution of the snake pedal largely simplifies our further derivation and computation, as seen in the following sections.
6.4.3 PDE for the Snake and Snake Pedal Evolution
With the relation between the snake and the snake pedal curve evolution given by (6.25), we now derive the equation of motion for the higher dimensional surface 01, in which the snake is embedded. As discussed in Section 6.3.2, the zero level set of the higher dimensional surface 01 is represented by Eq. (6.8),
Differentiating Eq. (6.8) with respect to t yields
0 +Vol1 O0 (6.27)
at at
Similarly, for the higher dimensional function q52, we have the following formula
_0_+V0 __ OP=0 (6.28)

Full Text 
PAGE 1
SNAKE PEDALS: ACTIVE GEOMETRIC MODELS FOR SHAPE MODELING AND RECOVERY .By YANLIN GUO 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 1998
PAGE 2
To My Motherland I wish you wealth, strength, and happiness
PAGE 3
ACKNOWLEDGMENTS I would like to express my sincerest gratitude to my advisor and committee chairman, Dr. Baba C. Vemuri, for his invaluable guidance, constant support and encouragement, in every step of my education and research at the University of Florida. I have greatly benefited by his stimulating approach to research, complete dedication to science, and relentless pursuit of perfection. My special thanks go to Dr. Christina M. Leonard for generously agreeing to serve as an external member of my dissertation committee, and for sharing her expertise in visual analysis of the human brain, and for providing many valuable comments on the applications of my work to automatic analysis of brain images as well as providing me with large amount of data. Many thanks to my dissertation committee members, Dr. Thomas E. Bullock, Dr. John G. Harris, and Dr. Stanley Su, for their willingness to serve on my committee and their interest and important comments and suggestions on various aspects of this dissertation. I would like to thank Dr. Hsu, of the Department of Aerospace Engineering, Mechanics and Engineering Science, for his insightful discussion on finite difference methods. Many thanks to former student ShangHong Lai (now in Siemens Corp. Research) for his stimulating discussions and ideas. I also would like to thank fellow students from the Computer Vision, Graphics and Medical Imaging lab Chhandomay iii
PAGE 4
Mandal, Yi Cao, Fengting Chen, Li Chen, Shuangying Huang, Li Wang, Jundong Liu, and Arun Srinivasan for the many discussions, ideas, advice, and friendship. Finally, I would like to thank my parents and sisters for their endless love and support, and my husband, for sharing the happiness and struggle I experienced, without which I could never have succeeded. The research was supported in part by the NIH under the grant number ROlLM05944 and the Whitaker Foundation. iv
PAGE 5
TABLE OF CONTENTS ACKNOWLEDGMENTS iii LIST OF FIGURES viii ABSTRACT xi CHAPTERS 1 INTRODUCTION 1 1.1 Problem Statement 2 1.2 Literature Review 3 1.2.1 Geometric & Deformable Models with Fixed Topology .... 4 1.2.2 Topologically Adaptable Models 6 1.3 Contributions 6 1.3.1 Compact Representation of the Model 8 1.3.2 Automatic Topological Change Handling Scheme 9 1.3.3 Fast Numerical Methods 9 1.4 Thesis Organization 10 2 MATHEMATICAL FOUNDATIONS OF ACTIVE MODELS 12 2.1 Snakes: Energy Minimizing Active Contours 14 2.2 Active Surface Models 18 3 SNAKE PEDALS: GEOMETRIC MODELS WITH PHYSICSBASED CONTROL 20 3.1 Definition of Regular Pedal Curves 20 3.2 Snake Pedals: The Proposed Model 23 3.3 Modified Snake Pedal Model 28 3.4 Handling Changes in Topology 30 3.5 Pedal Surfaces 31 3.5.1 Uniform Tessellation of the Model 34 3.6 Summary: Salient Features of the Model 38 V
PAGE 6
4 PROBABILISTIC FRAMEWORK AND MODEL LEARNING 40 4.1 Prior Model in the Bayesian Framework 42 4.2 Sensor Model in the Bayesian Framework 45 4.3 Maximum a Posterior 46 4.4 Supervised Learning Scheme 46 5 NUMERICAL SOLUTIONS FOR MODEL FITTING WITH SNAKE PEDAL MODELS 50 5.1 Related Work 51 5.2 Our Fast Numerical Algorithms 52 5.2.1 Preconditioned Nonlinear Conjugate Gradient Algorithm ... 53 5.2.2 Computation of the Approximate Hessian 56 5.2.3 LM Method for Global Parameter Estimation 58 5.2.4 ADI Method for Local Parameter Estimation 62 6 HYBRID GEOMETRIC ACTIVE MODELS 70 6.1 Overview of Hybrid Geometric Active Models 72 6.2 Curve Evolution Theory 73 6.3 Traditional LevelSet Approach 77 6.3.1 Embedding Mechanism 77 6.3.2 Equation of Motion 79 6.3.3 Advantages of the Levelset Approach 81 6.4 Our Model: Hybrid Geometric Active Contours 82 6.4.1 Big Picture 82 6.4.2 Relation Between the Snake and Snake Pedal Curve Evolution 86 6.4.3 PDF for the Snake and Snake Pedal Evolution 88 6.4.4 Relation Between and V02 91 6.4.5 PDE for the Higher Dimensional Function 95 7 NUMERICAL SOLUTIONS FOR SHAPE RECOVERY WITH HYBRID GEOMETRIC ACTIVE CONTOURS 97 7.1 Numerical Algorithm for Standard PDFbased Curve Evolution ... 98 7.2 Our Numerical Approach 102 8 SHAPE RECOVERY EXPERIMENTS WITH SNAKE PEDALS 105 8.1 Shape Recovery with Snake Pedals in 2D 105 8.2 Shape Recovery with Snake Pedals in 3D 108 8.3 Fitting Results with the Hybrid Geometric Active Model 115 vi
PAGE 7
9 CONCLUSIONS AND FUTURE DIRECTIONS 120 9.1 Salient Features of the Model 120 9.1.1 Geometric Model with Physicsbased Control 121 9.1.2 Compact Representation 121 9.1.3 Automatic Handling of Topological Changes 121 9.2 Future Directions 122 9.2.1 Invariant Representation 122 9.2.2 Introduce Dynamics in the Representation 122 9.2.3 Extend the hybrid geometric active model to 3D ....... . 122 APPENDICES A DERIVATION OF STIFFNESS MATRIX K 124 A.l Bilinear Quadrilateral Elements 124 A. 2 Derivation of the Local Stiffness Matrix 126 A. 3 Kronecker Tensor Product Representation of the Assembled Stiffness Matrix 127 B RELATIONSHIP OF ~ 135 REFERENCES 137 BIOGRAPHICAL SKETCH 144 vii
PAGE 8
Â» ;I", ? . LIST OF FIGURES 2.1 Active surface being pulled by a spring force showing the effect of various wi and W2 settings: (a) Wi = 0.9, = 0.1, (b)u;i = W2 = 0.5, and (c) wi = 0.01, u;2 = 0.99 19 3.1 f is on the pedal curve of a with respect to the pedal point p 20 3.2 Examples of pedal curves of an ellipse for different pedal point positions. Pedal points are shown by a dot in each case 22 3.3 (a) The process of generating a snake pedal with an ellipse as a generator, (b) "snake pedal" controlled by the snake using an ellipse generator. 24 3.4 Radialbased correspondence association scheme: each point cx{ti) on the generator is associated with a point p, on the snake 25 3.5 Examples of "snake pedals" using an ellipse generator and a snake. . 28 3.6 Examples of the modified "snake pedals" using an ellipse generator and a snake 29 3.7 Comparisons of original (red) and modified (blue) "snake pedals" using the same generators and same snakes in both case 30 3.8 Examples of topology change, (a) depicts a single snake pedal curve (solid) obtained from an ellipse generator (dashed) and a single pedal point, (b) same as (a) except two distinct pedal points located outside the generator are used, (c) same as (a) except several pedal points on a snake are used, (d) same as (c) except a discontinuous snake is used. 31 3.9 Examples of pedal surfaces with fixed pedal points located inside ellipsoidal generator surfaces 33 3.10 Examples of shaping snake pedal surfaces with snakes. Each figure was generated by applying interactive mousebased forces to the snake. 34 3.11 Uniform tessellation result: (a) f{u) ~ u relation, (b) f{v) ~ v relation. 38 4.1 Supervised learning to incorporate information into the prior model . 48 viii
PAGE 9
5.1 Big picture of numerical schemes for model fitting 59 5.2 Model tessellation in the model coordinates 64 6.1 Level set formulation of equations of motion (a) and (b) show the curve V and the surface ^{x, y) att = 0, (c) and (d) show the curve V and the corresponding surface (l){x, y) at time t. (e) and (f) show the curve V and surface y) at time t + Ai. The curve V has changed topology in going from (c) to (e) 78 6.2 Relation between 'P{t),Ve{t), 0i, and (/)2: V{t) and Ve{t) are levelsets of the higher dimensional surface (pi and 02 respectively. They are related by the pedaling operation 84 6.3 The hybrid geometric active model can capture the "global orientation" of clustered objects 85 6.4 Three representative pointpairs on the snake and snake pedal 92 7.1 Curve Evolution: The curve moves with unit speed in (a), and moves with curvature dependent speed in (b). The solid lines are the initial curves and the dashed lines are the curves after evolution 99 8.1 Fitting result (c) of a snake pedal to a caudate nucleus with initialization shown in (a) and intermediate stage shown in (b) 106 8.2 Fitting result (c) of a snake pedal to a cerebellum with initialization shown in (a) and intermediate stage shown in (b) 107 8.3 Fitting result (c) of a snake pedal (in green) to a hippocampus with initialization shown in (a) and intermediate stage shown in (b). The data points (in red) are placed by the expert along the hippocampus boundary 107 8.4 Fitting results of snake pedal surfaces to a hippocampus (d), a gyrus (h), and a cerebellum (1) with initializations shown in (b), (f), and (j), the intermediate stages are shown in (c), (g), and (k), respectively, (a), (e) and (i) are slices of MRI scans for the hippocampus, gyrus, and cerebellum, respectively 109 8.5 Fitting results of snake pedal surfaces to a bend (c) and twist (f) shape. Initializations are shown in (a) and (d) followed by the intermediate stages (b) and (e) respectively 113 8.6 Fitting results of snake pedal model to spock(c) and an anvil (f). Initializations are shown in (a) and (d) followed by the intermediate stages of fitting in (b) and (e) respectively 114 ix
PAGE 10
8.7 Hybrid geometric active contour model fit examples: (a) is model initializations (snake pedal in green and generator in red), (b) and (c) are intermediate stages of evolution and final fits are shown in (d). . . 115 8.8 Hybrid geometric active contour model fit examples: (a) is model initializations (snake pedal in green and generator in red), (b) and (c) are intermediate stages of evolution and final fits are shown in (d). . . 116 8.9 Hybrid geometric active contour model fit examples: Left to right, model initialization (snake pedal in green and generator in red), intermediate stages of evolution and final fit 117 8.10 Topological change examples with hybrid geometric active contour models: Left to right, model initialization (snake pedal in yellow or green and generator in red), intermediate stages of evolution and final fit 119 A.l Bilinear quadrilateral element, where {u,v) are material coordinates and (a;, il>) are the reference coordinates 125 A. 2 A small 3x4 mesh in {u,v) space 128 A. 3 Membrane molecules 129 A.4 Assembled membrane molecular representation 130
PAGE 11
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 SNAKE PEDALS: ACTIVE GEOMETRIC MODELS FOR SHAPE MODELING AND RECOVERY By Yanlin Guo August 1998 Chairman: Dr. Baba C. Vemuri Major Department: Electrical and Computer Engineering Modeling shapes, be it in 2D or 3D, is a fundamental constitutent of computer vision and computer graphics. Shape modeling techniques have been widely used in shape synthesis, shape reconstruction, recognition, motion tracking, and so on. In this dissertation, we introduce a compact and versatile geometric shape modeling scheme that can model a large class of shapes and is amenable to stable and efficient numerical implementations. Geometric models are traditionally well suited for representing global shapes but not the local detail. However, in this thesis we propose a powerful geometric shape modeling scheme which allows for the representation of global shapes xi
PAGE 12
with local detail and permits model shaping as well as topological changes via physicsbased control. The proposed geometric models are a blend of geometric and physicsbased models and are in spirit "similar" to the now popular deformable superquadric models but differ from them considerably in the expressiveness and numerical stability leading to greater applicability. . . , The proposed modeling scheme consists of representing shapes by pedal curves and surfaces Â— pedal curves/surfaces are the loci of the foot of perpendiculars to the tangents of a fixed curve/surface from a fixed point called the pedal point. By varying the location of the pedal point, one can synthesize a large class of shapes which exhibit both local and global deformations. We introduce physicsbased control for shaping these geometric models by letting the control point vary and use a dynamic spline, that is, a snake, to represent the position of this varying pedal point. The model dubbed as a "snake pedal" allows for interactive manipulation via forces applied to the snake. We extend the model to automatically cope with topological changes by introducing the concept of hybrid geometric active contour /surface models wherein the traditional snake in the snake pedal is replaced with a geometric active contour which is realized via a levelset implementation. Efficient numerical algorithms for shape recovery from image data using the proposed snake pedal model are presented. These algorithms involve novel mathematical derivations that lead to efficient numerical solutions of the model fitting problem. We demonstrate the applicability of this modeling scheme via shape synthesis examples and shape estimation examples from image data. xii
PAGE 13
CHAPTER 1 INTRODUCTION Shape modeling is a fundamental constituent of computer vision and computer graphics. In computer graphics, it can be used to synthesize shapes, object motion and object interaction for the purpose of computer aided design and computer animation. In computer vision, it is very crucial in 2D/3D shape representation, shape reconstruction, quantitative model extraction from image data for analysis, visualization, recognition and motion tracking. Recently, shape modeling techniques have been widely used in medical image analysis applications. Twodimensional and threedimensional shape models have been used in Xray, computer tomography (CT), magnetic resonance (MR), and ultrasound images to segment, visualize, track, and quantify a variety of anatomic structures including brain, heart, cerebral, lungs, liver and skull. Two primary tasks of importance in computer vision are shape recognition and shape recovery from image data. The shape recovery task imposes the requirement that the shape model be capable of representing fine detail and hence requires that the model have a large number of degrees of freedom. Shape recognition, on the other hand, requires that the model representation be compact in terms of the number of parameters so as to achieve easy storage of a large database of models and easy matching for retrieval. These two requirements are potentially conflicting. Therefore, 1
PAGE 14
2 designing shape models that can possibly "satisfy" these requirements has been the focus of many researchers in the recent years in the computer vision community [3, 11, 14, 54, 78]. 1.1 Problem Statement There are numerous shape modeUng techniques in the computer vision and computer graphics literature. Most modeling schemes fall into two major categories, namely, active models and passive models. Active models or deformable, physicsbased models are distributed parameter models. Typical examples are generalized spline models such as snakes [31]. Active models have broad geometric coverage and hence are good for shape recovery and reconstruction. On the other hand, they are not well suited for recognition since they require a large number of parameters for their representation. Passive models or geometric models are lumped parameter models. Typical examples are spheres, cylinders and prisms [4, 67]. Passive models can represent shapes with very few parameters and they are good for recognition but not for representing detailed shapes. In order to fill in the gap between the requirements of reconstruction and recognition, a class of hybrid models have been proposed. For example, Terzopoulos and Metaxas [78] proposed a hybrid model which combines the descriptive power of geometric and physicsbased models via superposition of a spline on a conventional superellipsoid. Hybrid models are a compromise between the two extremes and occupy both ends of the spectrum in shape description. One end of this spectrum is occupied by lumped parameter models and the other by distributed parameter models. However, most hybrid models introduce
PAGE 15
3 a lot of nonlinearities into the model representation to describe global deformations such as bending, tapering, and twisting. These nonlinearities increase the complexity in the model representation and affect the numerical stability of the algorithms used in fitting the model to image data. Moreover, the hybrid models can not handle shapes of unknown topologies without introducing additional nonlinearities into the modeling via blending schemes [14] or other special handling schemes. In this thesis, we develop a compact and versatile shape modeling scheme with physicsbased control for shape design and shape recovery from image data. This modeling technique allows for the representation of global shapes with local detail and allows model shaping via physicsbased control. The local/global descriptive power of the model satisfies the conflicting requirements of both shape reconstruction and shape recognition. In addition, topological changes can be handled naturally and automatically via a levelset implementation in the proposed model. In the remainder of this chapter, we will first review related shape modeling schemes in literature, with an emphasis on the advantages and disadvantages of each scheme, and then present the overview of our model. 1.2 Literature Review There are numerous shape modeling schemes in computer vision and computer graphics literature.
PAGE 16
4 1.2.1 Geometric fc Deformable Models with Fixed Topology Motived by the generalized cylinder idea, Kass et al. [31] proposed a deformable cylinder model constructed from generalized splines. They developed force field techniques for fitting their model to monocular, binocular, and dynamic image data. The distributive nature of this deformable model enhances its descriptive power and allows the representation of natural objects with asymmetries and fine detail. However, the generalized spline representation of the model does not explicitly provide a description of object shapes in terms of a few parameters. Terzopoulos and Metaxas [78] developed a hybrid dynamic model Â— a deformable superquadric model which combines the descriptive power of geometric and physicsbased models via superposition of a membrane spline on a core superquadric. Vemuri and Radisavljecvic [81] developed a multiresolution, hybrid modeling scheme which represents the displacement function of a deformable superquadric model in an orthonormal wavelet basis. This multiresolution basis provides the model with the ability to continuously transform from local to global shape deformations thereby allowing a continuum of shape models to be created and to be presented by relatively fewer parameters. Pentland and Williams [54] introduced a novel shape modeling technique based on the modal representation, which can be used to represent global shapes with local detail. The number of parameters required for this representation depends on the level of detail needed.
PAGE 17
5 Cohen and Cohen [11] and Han et al. [26] proposed a hybrid hyperquadric model that adds an exponential term to a hyperquadric primitive. The addition of the exponential term to the hyperquadric equation allows for local control of the shape generated. An arbitrary number of concavities can be generated with this modeling scheme. This modeling scheme differs from most of the schemes discussed above in that this scheme is basically a geometric model but it allows for global as well as local deformations. The model representation is quite compact in terms of the number of parameters needed to describe the shape. However, the computation of local concavities is nontrivial and human interaction is needed in the model fitting process. Another modeling scheme described by several authors in [3, 27, 30, 59] involves superimposing a global volumetric deformation Â— the free form deformation (FFD) Â— on the core superquadric. Both global and local deformations can be represented by this modeling scheme, however, it is difficult to describe the shape of an object with complex deformations and fine detail with this modeling scheme. Recently, O'Donnell et al. [50] developed a novel solid shape modeling scheme. This model includes builtin offsets from a core base model. Local deformation over the scaledoffset model is allowed via displacement vectors attached to the default model. Model fitting to data is achieved in a dynamic framework as in Terzopoulos and Metaxas [78]. This modeling scheme is well suited for tasks such as cardiac motion recovery from tagged MR images.
PAGE 18
6 1.2.2 Topologically Adaptable Models Surfaces of arbitrary topology were modeled using particle systems to model surfaces in Szeliski and Tonnesen [72]. Particles can be added and deleted dynamically to enlarge and trim the surface. An interesting and topologically adaptable deformable model was developed by DeCarlo and Metaxas [14] via the use of parameterized blending functions. Their scheme can determine the number and locations of the holes in an object with unknown topology and accordingly change the topology by use of appropriate error of fit criteria. Numerous techniques based on the concept of levelsets have been reported in the literature [9, 10, 46, 58, 82]. These techniques allow a deformable contour not only to represent long tubelike shapes or shapes with bifurcations, but also to dynamically sense and change its topology. Another modeling scheme called "topologically adaptable snakes" was introduced by Mclnerney and Terzopoulos [48]. This novel version of snakes allows for automatic change in topology during shape recovery. A common problem associated with these models is that they do not possess a mechanism to characterize the global/core shape of an object. 1.3 Contributions In this thesis, we propose a powerful geometric shape modeling scheme which permits the representation of global shapes with local detail using only a "semiglobal" characterization and relatively simple and efficient numerical methods. The
PAGE 19
proposed modeling scheme consists of representing shapes by pedal curves and surfacespedal curves/surfaces are the loci of the foot of perpendiculars to the tangents of a fixed curve/surface from a fixed point called the pedal point. By varying the location of the pedal point, one can synthesize a large class of shapes which exhibit both local and global deformations. We introduce physicsbased control for shaping these geometric models by letting the pedal point vary and using a dynamic spline, namely, the snake, to represent the position of this varying pedal point. This process of generating models is referred to as the "pedaling operation." The model, dubbed as a "snake pedal," allows for interactive manipulation via forces applied to the snake. In this model, automatic changes in topology can be realized by embedding the snake in a levelset based evolution. Â• Â• \ The key contributions of this thesis are summarized below in the form of salient features of the developed modeling scheme, the snake pedal, and the associated numerical algorithms: . ; ' 1. The snake pedal model is inherently geometric but can exhibit both local and global deformations as would a physicsbased model. 2. The model can exhibit large global deformations such as bending, twisting, and tapering without explicitly incorporating parameters for the same. 3. The modeling scheme allows (for the first time) the concept of a global shape to be incorporated into the curve/surface evolution framework.
PAGE 20
8 4. The model permits a levelset implementation and thereby seamlessly allows for topological changes in shape during evolution. 5. Fast numerical algorithms for estimating the model fits to 2D/3D data sets. 6. A hierarchical neural network based learning scheme for incorporating prior knowledge about shapes of interest from the application domain. We elaborate these contributions in the following sections. 1.3.1 Compact Representation of the Model The snake pedal model is inherently geometric but can exhibit both local and global deformations as would a physicsbased model. Our modeling scheme is somewhat similar to the hybrid modeling schemes mentioned in Section 1.2.1 in that both techniques combine the descriptive power of physicsbased models and geometric models. The difference lies in that the hybrid models discussed earlier have to explicitly include the global parameters to incorporate global deformations such as bending, tapering, and twisting to achieve the required deformations, which causes additional nonlinearities in their representation, while our model does not. This compact representation of global deformations by our model can avoid the introduction of additional nonlinearities into the model. Therefore, numerical efficiency and stability can be largely improved in the model fitting to the image data. This compact representation comes from the mechanism of creating shapes in the snake pedal model. The pedaling operation is a nonlinear operation. When applying the pedaling operation to geometric primitives, one can produce global
PAGE 21
9 deformations of the geometric primitives by performing the nonlinear transformation. This nonlinear operation is only applied in the model building phase, and it does not affect the model fitting phase, as we will discuss in Chapter 5. Hence no additional nonlinearities are added to the model fitting process. 1.3.2 Automatic Topological Change Handling Scheme The proposed snake pedal model allows for the estimation of shapes of arbitrary topology without resorting to the blending or other special operations as required in Szeliski and Tonnesen [72]. The topology changes are achieved naturally by embedding the modeling scheme into a levelset framework [46]. We introduce the novel concept of a global/core model into the now popular partial differential equation (PDE) based curve/surface evolution framework. This global/core model allows us to express any shape as a combination of a global shape such as ellipsoid and sphere, and a variable offset with respect to this global shape and an evolving curve/surface. Because of the presence of the global shape, the snake pedal model can be regarded as a unique hybrid geometric active model. 1.3.3 Fast Numerical Methods Novel and fast numerical methods are proposed to achieve stable local optimal solutions in the model fitting process. These numerical methods are not a regurgitation of techniques that are described in a standard numerical analysis text book, or for that matter, the book on numerical recipes in C [84]. In all the methods that are discussed here, novel derivations are presented which lead to development of software modules that can either be easily embedded into or used with standard
PAGE 22
10 numerical software such as the conjugate gradient (CG) code and the alternating direction implicit (ADI) code. 1.4 Thesis Organization The remainder of the thesis is organized as follows: Chapter 2 contains a review of the basic mathematical foundations of classical active models, especially the snake formulations in order to illustrate the physicsbased control that is also pertinent to the snake pedal model. In Chapter 3, we present the definitions of pedal curves and surfaces and show examples of the same in 2D and 3D, followed by the development of the proposed model, namely, the "snake pedal" with a demonstration of the expressive power of the modeling scheme via illustrative examples of shape synthesis. In addition, we also discuss the ways of automatically handling change of topology. In Chapter 4, we cast the model fitting process in a probabilistic framework, and present the advantages of the probabilistic framework over the deterministic approach. We will devise a novel supervising learning scheme to incorporate the prior knowledge of object shapes into the shape estimation from image data. In Chapter 5, we present efficient and stable numerical algorithms for fitting the snake pedals to real image data. The model fitting process can be posed as a nonlinear optimization problem and can be dissected into global and local shape parameter estimation. We present efficient numerical solutions for both tasks, and then combine the global and local procedures appropriately to solve the overall optimization problem.
PAGE 23
11 In Chapter 6, we extend the snake pedal model to cope with topological changes by introducing the concept of hybrid geometric active models. We first review the traditional geometric active curve evolution theory and its levelset implementation. Following this, we present our hybrid geometric active model evolution and its levelset implementation. The hybrid geometric active contour evolution involves the evolution of both snake and snake pedal curves/surfaces. However, by employing the insightful relationship between the two curve/surface evolutions, we solve these evolutions simultaneously by only solving for one of them and thereby dramatically reduce the computational and storage requirements when applied to shape recovery tasks. In Chapter 7, we present the numerical techniques used for computing a stable solution to the hybrid geometric active model evolution equation. The equation of motion for the model can be posed as a special type of HamiltonJacobi partial differential equation (PDE). Based on the traditional levelset approach, the motion equation of the our model can be efficiently solved by employing entropysatisfying upwind finite difference plus minmod finite difference schemes. In Chapter 8, we present the experimental results and the applications of the model in shape modeling and recovery and conclude in Chapter 9. Appendix A gives the detail derivation of the stiffness matrix, which corresponds to the internal smoothness energy imposed on the snake in the snake pedal model. Appendix B contains the derivation of the relationship between the curvature of the snake pedal and that of the snake at a specific point.
PAGE 24
CHAPTER 2 MATHEMATICAL FOUNDATIONS OF ACTIVE MODELS In the proposed snake pedal model, we use an active curve/surface to control the model shapes and deformations, this active curve/surface is an integrated and nonseparable part of the modeling scheme. In this chapter, we will review the mathematical foundations of active models in order to obtain a better understanding of the behavior of the controlling curves/surfaces in the snake pedal model. Active or deformable curve and surface models gained popularity after they were proposed for use in computer vision [79] and computer graphics [77] in the mid 1980s. In computer vision, deformable curve and surface models were proposed as solutions to illposed inverse problems such as edge detection and surface reconstruction [55, 75]. The theory of continuous deformable models was introduced by Terzopoulos in a Lagrangian dynamic setting [76], based on deformation energies in the form of (controlledcontinuity) generalized splines [75]. The controlledcontinuity spline is a generalization of Tikhonov and Arsenin stabilizer [80], and can formally be regarded as regularizing the illposed problems [56], by recasting them as wellposed functional minimization problems. Active or deformable model can represent broad classes of shapes by employing geometric representations that involve many degrees of freedom, such as splines. The degrees of freedom are generally not permitted to evolve independently, but are 12
PAGE 25
13 governed by physical principles that bestow intuitively meaningful behavior. The use of physical principles, especially elasticity theory, generally within a Lagrangian setting, give rise to the name "deformable models" or "active models." The physical interpretation views deformable models as elastic bodies which respond naturally to applied forces and constraints. Typically, deformation energy functions are defined in terms of the geometric degrees of freedom for the deformable models as for physically elastic objects. The energy of the model grows monotonically as it deforms away from a specified natural or "rest" shape and often includes terms that constrain the smoothness or symmetry of the model. In the Lagrangian setting, the elastic forces internal to the model contribute to the deformation energy, while the external potential energy functions are defined in terms of the data of interest to which the model is to be fitted. These potential energies give rise to external forces which deform the model such that it maximally fits the data. From a theoretical point of view, the concept of active or deformable models reflects the integrated fluency of geometry, physics, and approximation theory. Geometry serves to represent object shape, physics imposes constraints on how the shape may vary over space and time, and the theory of optimal approximation provides the formal support to the mechanisms for fitting the models to measured data. The active/deformable model that has attracted the most attention in the computer vision and computer graphics community is popularly known as "snakes" [5, 31].
PAGE 26
14 Snakes or "active contour models" represent a special case of the general multidimensional deformable model theory. We will review their formulation in this chapter in order to illustrate the basic mathematical foundation of the controlling curves/surfaces in our proposed shape modeling scheme, namely, the snake pedal model. * 2.1 Snakes: Energy Minimizing Active Contours Snakes are energy minimizing splines guided by internal and external forces. Internal forces act as smoothness constraints, while external forces pull the snake toward the desired features such as edges and image boundaries. By initializing the snake around the target object, it can find the desired object contours by minimizing a userdefined energy function. Because the way the contours slither while minimizing their energy, we call them "snakes." The ability of the snakes to converge to the outline of the object despite the intensity gradient inside the object makes them very useful in many visual tasks, such as image segmentation, boundary detection, motion tracking, and stereo matching. Geometrically, a snake consists of a set of discrete points, effectively connected by straight lines. Each point has a position, given by the coordinates in the plane {x,y), and a snake is completely specified by the number and coordinates of these points, namely, p(s) = {x{s),y{s)) with the domain parameter s e [0, 1]. When the snake is a closed curve, we apply the periodic boundary condition, p(0) = p(l), on the snake. A snake is controlled by minimizing a functional which converts highlevel contour information, such as curvature and discontinuity, and lowlevel Arsenin image
PAGE 27
15 information, such as edge gradients into an energy represented by E{s) = r E{p{s)) ds = [\eUp{s)) + EeMs))} ds (2.1) Jo Jo where, Eint and Ee^t are the internal and external energies associated with the snake. The shape and position of the snake corresponds to the minimum of this functional. The first term of the above functional EinMs)) = Wi{s) dp ds d'^p ds' (2.2) is the internal deformation energy of the snake. It characterizes the deformation of an elastic contour, Wi{s) and W2{s) are nonnegative parameters controlling the "tension" and "rigidity" of this elastic contour. The second term of (2.1) consists of the external potential that guides the snake toward desired features. One form of this potential is a potential function defined on the image plane, whose local minima coincide with intensity extrema, edges, and other image features of interest. For example, if we want to attract the contour toward the edges of a gray level image I{x, y), the potential takes the form EexMs)) = Po{p{s)) ds, Jo (2.3) with Po{p{s)) = Po{x, y) = c I VG, * I{x, y)\ , (2.4)
PAGE 28
16 where c controls the magnitude of the potential, V is the gradient operator, and Ga * I denotes the image convolved with a Gaussian smoothing filter, and a is the width of the filter which controls the spatial extent of the local minima of EextAnother form of the external energy potential is the user defined external constraint potential such as interactive springs, anchored springs, and volcanos [31]. For example, the snake can be pulled in the direction of the mouse cursor location i^m, Urn) by usiug a spring potential Ee^t = c[{x Â— + {y Â— J/m)^], where c controls the stiffness of the spring. Note that the points on the snake that are affected by the spring can be restricted to a small section of the snake closest to the mouse. The combination of external image potentials and external constraint potentials allows snakes to extract and represent a broad spectrum of shapes. Furthermore, external constraint potentials are an effective, flexible means from which highlevel control mechanisms can guide shape recovery, forming a sound basis for fully automatic image analysis. Using the calculus of variations, we obtain the contours p(s) that minimize the energy E. They must satisfy the EulerLagrange equation This vectorvalued partial differential equation expresses the balance of internal and external forces when the contour rests at equilibrium. The first two terms represent the internal stretching and bending forces respectively, while the third term represents the external forces that guide the snake to the image data.
PAGE 29
17 In order to compute numerically a minimum energy solution, it is necessary to discretize the energy (2.1). By using the finite element method [87] or the finite difference method [70], the discrete form of the energy E{p{s)) can be written as: where K is the stiffness matrix, and E{p) is the discrete version of the external potential. Differentiating (2.6) with respect to p results in the minimum energy solution, which is equivalent to solving the set of nonlinear algebraic equations: where f{x,y) = {dPo/dx,dPo/dy) is the generalized external force vector. It is possible to introduce dynamics into the above static formulation of the snake. This may be done by introducing a dissipation energy functional which dissipates the kinetic energy during the snakes motion. The motion equation poses an initial boundary value problem which can be solved using efficient numerical techniques leading to a real time response from the snake, to applied forces. The initial value problem can be solved using explicit Euler integration. Let 7 denote the Euler (time) step size, then the equations for time integration may be written as Â£;(p) = ^p^Kp + VPo(p) (2.6) Kp + f (x, y) = 0, (2.7) P (K + 7l)i(7p*f(x*,y')). (2.8)
PAGE 30
18 For compactness in representation, we can replace the snake with a Bsnake which is a snake constructed from Bsplines that have builtin smoothness. In this case, we simply replace p(s) by XI^o C'iBi(s), where are the control points of the chosen Bspline basis and m are the number of control points. For more details on Bsnakes, we refer the reader to Marc et al. [47]. 2.2 Active Surface Models A active/deformable surface is represented by a vectorvalued parametric representation p{u,v) = {x{u,v),y{u,v), z{u,v))^, where {u,v) are domain parameters, and the parametric domain is the unit square [0, 1]. The deformation energy associated with the deformable surface, which simulates the thinplatemembrane material surface, is given by eup) = I^Mfj ^ (^?) + M(^r + 2(^/ + i'^n * d.. (2.9) Where Eint can be regarded as a controUedcontinuity spline defined by Terzopoulos [75]. As described before, and wi and W2 control the tension and rigidity of the spline respectively. Increasing wi tends to decrease the surface area of the material, while increasing W2 tends to make it smoother but less flexible. Fig. 2.1 demonstrates the effect of the value changes of Wi/w2 to the material. Using variational principles, we can obtain the EulerLagrange equation for (2.9) and its corresponding stiflfness matrix K. For the details of (2.9), we refer the reader to the Appendix A and Terzopoulos and Metaxas [78].
PAGE 31
19 (a) (b) (c) Figure 2.1. Active surface being pulled by a spring force showing the effect of various wi and W2 settings: (a) Wi = 0.9,1^2 = 0.1, {h)wi = W2 = 0.5, and (c) Wi = 0.01, = 0.99.
PAGE 32
CHAPTER 3 SNAKE PEDALS: GEOMETRIC MODELS WITH PHYSICSBASED CONTROL In this chapter we will first present definitions for the regular pedal curves/surfaces and illustrate the ideas via synthesized examples. We will then present our proposed modeling scheme Â— "snake pedal" model. Snake pedals are inherently geometric models, but they can express local detail, as well as cope with topology changes of shapes automatically. We will illustrate the power of this modeling scheme via some synthesized examples in this chapter as well as some model fitting examples in Chapter 8. 3.1 Definition of Regular Pedal Curves As shown in Fig. 3.1, let a be a planar curve, the pedal curve of a is defined as the locus of points on the foot f of the perpendicular from a fixed point p called the control point to a variable tangent of a. More specifically, let the planar curve pedal point p Figure 3.1. f is on the pedal curve of a with respect to the pedal point p. 20
PAGE 33
21 a be a parameterized curve with domain parameter t, let /3 be the pedal curve of at with respect to the fixed pedal point p, and let a{t) = g, /3 (t) = f at a particular value of parameter t, the projection of a (t) Â— p in the direction Ja (t) must be P (t) Â— p, where a' {t) is the tangent of plane curve a (i), J : Â— y is a linear map given by J{x,y) = {x,y) (3.1) Note that applying the operation J on a vector (x, y) in 5?^ Euclidean space results in a new vector (Â— y, x), and J can be geometrically interpreted as a rotation by n/2 in a counterclockwise direction. We can thus formally define a regular pedal curve as follows [23]: ; ; .' . Definition : The pedal curve of a regular curve a : (c, d) Â— > with respect to a fixed ( control) point p G 9?^ is given by , . (a(t) p) J a ' (t) ^ ,,. , ^ pedal(i) = p + II Ja '(t) ^^'^^ Where (c, d) is the parameter domain for t. The definitions of regular pedal curves and surfaces presented here may be found in any book on elementary differential geometry text [23, 38, 42]. Here we adopt the definition by Gray [23]. Some examples depicting the pedal curves of an ellipse for different positions of the pedal point are shown in Fig. 3.2. Usually the global or core shape depicts an overall size and orientation of an object, and local deformations depict fine detail in the local positions of an object. It is clear from Fig. 3.2 that the pedal curve is capable of exhibiting both global and local deformations. For example, the global
PAGE 34
22 (a) (b) (c) (d) Figure 3.2. Examples of pedal curves of an ellipse for different pedal point positions. Pedal points are shown by a dot in each case. shape of a pedal curve of an ellipse is approximated by a new ellipse, whose long axis is almost the same as that of the original ellipse, while the short axis is much longer than that of the original ellipse. The new ellipse is pinched in the direction of the short axis to produce the pedal curve. The location and extent of this pinching operation are determined by the position of the pedal point. As shown in Fig. 3.2, if the pedal point is located in the center of the original ellipse (Fig. 3.2 (a)), the pinching occurs around the center part of the pedal curve. As the pedal point moves to an upper position (Fig. 3.2 (b)), the pinching position in the pedal curve moves up accordingly. Similarly, as the pedal point moves to the right of the center of the ellipse (Fig. 3.2 (c)), more pinching occurs on the right side in comparison to the left of the origin. This pinching operation constitutes the local deformation of the pedal curve, but the type of the local deformation that can be exhibited by a pedal curve is by no means limited to the pinchingtype or cusptype deformations. More complex local deformations can be realized in the pedal curve as shown later subsequently.
PAGE 35
23 . \ t '.> *;' , ' In summary, by moving the position of the pedal point, it is possible to synthesize a variety of global and local deformations in regular pedal curves/surfaces. The original curve a(f) will henceforth be referred to as the generator for the pedal curve j3 (t) and the process of generating a pedal curve will be referred to as the pedaling operation, as depicted in Eq. (3.2). Note that the word "generator" here represents a curve in 2D or surface in 3D, which differs from its ordinary meaning in electricity industry. 3.2 Snake Pedals: The Proposed Model From the discussion in Section 3.1, we observed that when the pedal curves are generated from ellipses with respect to single fixed pedal points, the shapes and deformations exhibited by regular pedal curves are quite limited. To overcome this limitation, we can let the pedal points be different for each point of the generator when applying the pedaling operation to the generator, and thereby synthesize more general shapes of pedal curves/surfaces. We can let these pedal points be specified by another curve p{t), which can be represented by a active or deformable model (described in Chapter 2), that is, a standard snake or a Bsnake [31, 47] and then apply the pedaling operation to each point on the generator a{ti) with respect to the corresponding pedal point p{ti). Here we use p{t) to represent a parametric curve, p{ti) represents a particular point on the curve p(^) when t = U. Note that in this generalization of regular pedal curve definition, for different points on the generator, we apply the pedaling operation with respect to different control points, unlike in the regular pedal curves where only a single fixed pedal point is used.
PAGE 36
24 snake pedal snake (a) (b) Figure 3.3. (a) The process of generating a snake pedal with an ellipse as a generator, (b) "snake pedal" controlled by the snake using an ellipse generator. In this generalization process, we can employ a proper association mechanism to establish a correspondence between the points on the generator a(^) and the points on the curve p(^). This association scheme should keep the continuity property of the curve p(t) and be easy to implement. In accordance with these two requirements, we developed a simple correspondence association mechanism as shown in Fig. 3.4. In this mechanism, for each point (x{ti) on the generator , we draw a straight line from the center of the generator O to the point a{ti), and denote by p{ti) the intersection of the line Oa{ti) and the curve p{t). When applying the pedaling operation to point a{ti), we carry out this operation with respect to point p{ti). The pedaling operation is performed for each of the points a{tj),j = 1, 2, M, on the generator resulting in a curve that we call the "snake pedal." Since Oa{ti) and OpiU) are on the same radial line of the generator, we refer to this association scheme as the
PAGE 37
25 0C(t4) a(ti)/ a(tj) "\ a(t2) a(ti) snake generator (a) Figure 3.4. Radialbased correspondence association scheme: each point oi{ti) on the generator is associated with a point on the snake. radialbased association technique. It is obvious in this scheme that we use the same parameterization for the generator a(t) and curve p(t) with the same domain parameter t. It can be proved that snake pedal has the same continuity property as that of the curve p(t), given that the generator is differentiable, although this may not always be the case. The generator can be either a parameterized or an implicit function representing a curve. In this chapter, we will focus on parameterized curves and in most applications, we use an ellipse as the generator in 2D. If the generator is an ellipse, we can
PAGE 38
26 represent it in a parametric form by cos 6 sin 0 Â—sin 9 cos 9 ai cos t a2 sin t mi m2 (3.3) where ai,a2 are aspect ratio parameters and we use = (01,02)^ to denote these shape parameters, 9 is the rotation angle between the world coordinates and the object centered coordinates which is denoted by ge Â— (9)'^, and gc = (mi, 7712)^ is the centroid of the generator in the world coordinates. We collect all the generator parameters into a vector g = (gf , g^", gj)"^ = (ci, a2,9, mi, 7712)^, and the vector g is referred to as the global shape parameter vector. For the ellipse, the term Ja'{t) in the Eq. (3.2) of the pedal curve can be expressed as Ja'{t) = ni n2 ai sin t sin 9 Â— a2 cos t cos 9 ai sin t cos 9 + a2 cos t sin 9 (3.4) Note that  Ja'lp = n^+nl = sin'^t+al cos'^t, which is independent of the rotation angle 9. We will see that Ja:'p is a very important quantity in the model fitting presented in later Chapters. The curve p{t) on which the pedal points are located can be of any kind, which is not limited to be to a parameterized function. By using the correspondence association scheme described earlier, we can sample the curve at the same rate as that
PAGE 39
. 27 of the generator so that every point on the generator a [t) can find a unique corresponding point on the snake p{t). The pedaUng operation can then be applied to every point on the generator with respect to its corresponding point on the snake. We use a snake/Bsnake as this controlling curve and get a sequence of position vectors, Pi = {xi,yiY , by sampUng the snake curve. We collect all the snake position vectors in another vector, p = (pj", pj, Pm)*^; where M is the number of sampling points. Vector p is referred to as the local shape parameter vector. Vector p is the discretized version of continuous curve p(t). In our implementation, since the pedal points are located on a snake or Bsnake and the generalized pedal curve is created by applying the pedaling operation to the generator with respect to the points on the snake, we dub this generalized pedal curve as a snake pedal. With the radialbased correspondence associating scheme, a snake pedal is completely determined by the global shape parameter, namely, the generator parameter g, and the local shape parameter, namely, the snake position vector p. Fig. 3.3 demonstrates the process of creating the snake pedal curve. With the correspondence association scheme, the generalized pedaling operation can create much broader classes of shapes with more deformations. We show some examples of snake pedal curves generated using snakes and an ellipse as the generator in Fig. 3.5. Note the variety of local deformations that can be generated using this modeling technique. We remind the reader that the snake pedal itself is a geometric model and that it is not directly responsive to the application of external forces unlike the standard snake models described in Chapter 2. Also, it is easy to synthesize
PAGE 40
28 (a) (b) (c) . , (d) (e) (f) Figure 3.5. Examples of "snake pedals" using an ellipse generator and a snake. snake pedals whose topology is different from that of the generator, as will be seen in examples depicted subsequently. The uniqueness of our modeling scheme lies in the ability to generate snake like local deformations in a nonreactive geometric shape, namely, the pedal curve/surface. 3.3 Modified Snake Pedal Model The pedal curve definition can be modified slightly by subtracting the second term from the first term in equation 3.2, pedal(i) = p ^AJ^>_Â±lj^ (3.5)
PAGE 41
29 (d) (e) (f) Figure 3.6. Examples of the modified "snake pedals" using an ellipse generator and a snake. This modification allows us to synthesize a larger class of shapes which are quite different from the shapes produced using Eq. (3.2). The key feature of using this modification is that the snake pedal curve allows for more local deformations including shrinkage and expansion of the snake pedal. The shrinkage and expansion behavior is controlled by the location of the snake. When the snake is located inside the generator, the snake pedal exhibits shrinkage behavior, while expansion behavior appears when the snake is outside the generator. Fig. 3.6 depicts examples of shapes generated using this modification. Note the amount of local deformations allowed in the snake pedal curves using this modification. In Fig. 3.7, we compare the original and "modified" snake pedals using the same generators, same snakes and same correspondence associating schemes. It is
PAGE 42
30 Figure 3.7. Comparisons of original (red) and modified (blue) "snake pedals" using the same generators and same snakes in both case. obvious that the "modified" snake pedals exhibit larger local deformations in both Fig.3.7 (a) and (b). 3.4 Handling Changes in Topology , Topology changes may be achieved by applying the pedaling operation to a single generator and several pedal points (or snakes). Fig. 3.8 depicts examples of the topology change realized in 2D. When the snake is located inside the generator, a single closed snake pedal curve is created (Fig. 3.8 (a)), when the snake is split into two parts with both parts located outside the generator, two separated closed snake pedal curves are obtained (Fig.3.8 (b)). In the second case, only part of the generator is needed to produce the two snake pedal curves. Further study shows that any number of closed snake pedals can be generated from one generator by using this procedure, if the position of snakes and parts of the generator are properly chosen. When using the correspondence associating scheme discussed previously, the number of closed snake pedal contours are determined by the number of closed snakes. More generally, the topological property of the snake pedal is determined by the topology of the controlling snake. Therefore, we can accomplish the automatic topological change
PAGE 43
31 (a) (b) (c) (d) Figure 3.8. Examples of topology change, (a) depicts a single snake pedal curve (solid) obtained from an ellipse generator (dashed) and a single pedal point, (b) same as (a) except two distinct pedal points located outside the generator are used, (c) same as (a) except several pedal points on a snake are used, (d) same as (c) except a discontinuous snake is used. for the snake pedal model if the topology of the snake can be updated automatically in the data fitting applications. In Chapter 6, we will introduce a novel levelset based approach to realize the automatic topological change for the snake pedal model. 3.5 Pedal Surfaces A regular pedal surface is the surface analog of the pedal curve. It is the locus of the points on the foot of the perpendicular from a fixed pedal point to a variable tangent plane of the surface. More precisely, it is defined as follows [23]. Definition : The pedal surface of a parameterized surface or a patch ct : U Â— > 3?^ with respect to a point p e is defined by
PAGE 44
32 Where U is an open subset of 3?^, cx{u,v) is a parameterized generator surface with parameters u and v and ttu and aÂ„ are tangent vectors in the u, v directions. For an ellipsoid generator, we have Qi cos U COS V mi a{u, v) = H 02 COS u sin V + m2 aa sin v (3.7) Where ai , 02 and 03 are aspect ratio parameters which can be collected in gs = (oi, 02, as)^. We use gc = (mi, 7712, ma) ^ denote the centroid of the generator in the world coordinates. R is the rotation matrix between the world coordinates and the object centered coordinates. We use quaternions [63] to represent the rotation matrix because the rotation matrix is easier to update through the quaternion update and remains orthogonal in the mean time. A quaternion [s, v] defines a rotation of the model from its reference position through an angle 9 = 2 cos~^s around an axis aligned with vector v = {vi,V2,V3)'^ . The rotation matrix corresponds to [s,v] is 1 2{vl + vj) 2{viV2 SV3) 2{ViV3 + SV2) = 2{viV2 + SV3) 1 2{vj + vj) 2{v2V3 svi) (38) 2{ViV3 SV2 2{v2Vs + SVi) 1 2{vj + vl) We collect the quaternion parameters into a new vector g^ = (s, q^) = (s, Vi,V2,V3) The ellipsoid generator is then completely represented by g = (gj', gj, gj)^. As in 2D case, g is referred to as the global shape parameter vector.
PAGE 45
33 (d) (e) (f) 't ' N 1 Figure 3.9. Examples of pedal surfaces with fixed pedal points located inside ellipsoidal generator surfaces. Fig. 3.9 depicts examples of pedal surfaces synthesized using an ellipsoidal generator surface and a fixed pedal point located inside the generator surface in each example. Note that change in topology is achieved by simply choosing an appropriate location of the pedal point as shown in Fig. 3.9. However, this type of topology change is diflferent from what we described in Fig. 3.8. As in the 2D case, we can let the pedal point vary for each point on the generator surface. Thus we have the snake pedal surfaces in 3D whose shape can be controlled by snakes which are either curves or surfaces in 3D. Fig. 3.10 depicts some examples of snake pedal surfaces in 3D. They are synthesized by using an ellipsoid as a generator and a snake curve in 3D which lies inside the generator. The curve is shown outside of the generator adjacent to the snake pedal model for the clarity of visualization. Note that the snake pedal shape can be manipulated by applying external forces to the snake as illustrated in the Fig. 3.10 (a)(d).
PAGE 46
34 (c) (d) Figure 3.10. Examples of shaping snake pedal surfaces with snakes. Each figure was generated by applying interactive mousebased forces to the snake. For the 3D snake pedal examples obtained by ellipsoids and snake surfaces, we refer the reader to the next chapters. 3.5.1 Uniform Tessellation of the Model From the definition (3.6) of the pedal surface, the points on the pedal surface must be located in the tangent of the generator for every point on the generator specified by {u,v). For an ellipsoid generator, the evenly spaced domain parameters u and V are often referred to as the "latitude" and "longitude." It should be noted that the snake pedal model exhibits nonuniformly tessellated surface in spite of the domain parameters u and v being uniformly tessellated/sampled. This is obvious from the nonlinear function defining the pedal surface in Eq. (3.6). For example, the tangent of the generator undergoes a rapid change (in orientation) as u n/2, and
PAGE 47
35 a slow change (in orientation) as u 0. Therefore, for an evenly sampled generator and snake, when we use the association correspondence mechanism to produce a snake pedal surface, the snake pedal model will not be uniformly tessellated. This uneven tessellation phenomenon is not noticeable when the aspect ratio parameters Oi , 02 and 03 of the ellipsoid generator do not differ much from each other. But when the difference of the aspect ratio parameters increases, this nonuniformity becomes obvious, and leads to the distortion of the deformation properties of the snake pedal surfaces. In addition, numerical instabilities may be caused in the computation of derivatives for nonuniformly tessellated snake pedal. Hence a uniform tessellation scheme for pedal surfaces must be employed to overcome this problem. It is a natural way to first construct uniformly spaced tangent directions in the parameter direction u, then apply the same technique to achieve uniform spacing along the parameter direction v. We now present our scheme for the construction of uniform tessellation of snake pedal surfaces. When constructing the uniform tessellation of an ellipsoid generator in the parameter direction u, we actually produce a uniformity in the tangent direction for an ellipse along a fixed longitude {v is fixed in this case). Since we do not need to consider the rotation and translation of the ellipse when constructing the uniform tessellation, for an arbitrary point (xo,yo) on an ellipse, the standard parametric representation is: V. (3.9)
PAGE 48
36 The straight line equation for the tangent of this elUpse on the point (xq, yo) is ^ + ^ = 1, (3.10) Substituting Eq. (3.9) into (3.10), we obtain a2 COSU a2 in ^^\ y = : Â— x + Â— . (3.11) ai sinu smu The slope of this tangent line is: Â— ^ , which corresponds to an angle of 0 = tan~^{ ^ cotu). The angle 9 is not uniformly spaced when the variable u is evenly spaced, which results in a nonuniform tessellation of the model. The objective of the uniform tessellation is to modify the standard ellipse parametric formula so as to construct a uniformly tessellated 9 when u is evenly spaced. We change the standard ellipse equation (3.9) to : xq = ai cos f{u) (3.12) yo = a2sinf{u), where f{u) is the function to be solved. The angle of the tangent direction of the ellipse accordingly becomes 9 = tan~^{ ^ cot f{u) ). To make this angle uniformly tessellated, we must set . ^ ^{tan\"^cotf{u))} = C, (3.13)
PAGE 49
37 where C is a constant. The boundary condition is fin) 2 ' (3.14) < u = + 2 /(Â«) Basic algebraic manipulation of Eq. (3.13) results in the following ordinary differential equation (ODE) for f{u) where f'{u) is the derivative of f{u) with respect to u. This ordinary differential equation can be solved numerically. We use a third order RungeKutta method. The result is shown in Fig. 3.11 (a). It can be seen that when u is small {u Â— > 0), f{u) changes much faster than u, which tends to make up the fact that the tangent direction of the ellipse changes slowly in this region. While when u Â— 7r/2, f{u) changes much slower than u, which is in an attempt to slow down the rapid change in tangent direction in this area. After changing u to f{u) in the ellipse expression, the change in tangent directions will become even. Carrying out the same procedure on the parameter v, we obtain the f{v) ~ v relationship as shown in Fig. 3.11 (b). Note that the retessellation scheme needs to be done whenever the parameters of the snake pedal are subject to large changes. f'{u) = C {'^szn'f{u) + ^cos'f{u)), (3.15)
PAGE 50
38 (a) (b) ' " r Figure 3.11. Uniform tessellation result: (a) f{u) ~ u relation, (b) f{v) ~ v relation. 3.6 Summary: Salient Features of the Model Compared with the hybrid models mention in Section 1.2.1, the snake pedal models are bestowed with several unique features which make them very powerful and attractive for use in shape synthesis and recovery applications. A list of these features are: 1. The snake pedal model is inherently geometric but can exhibit both local and global deformations as would a physicsbased model. 2. The model can exhibit large global deformations such as bending, twisting, and tapering without explicitly incorporating parameters for the same. 3. The modeling scheme allows (for the first time) the concept of a global shape to be incorporated into the curve/surface evolution framework.
PAGE 51
4. The model permits a levelset implementation and thereby seamlessly allows for topological changes in shape during evolution.
PAGE 52
CHAPTER 4 PROBABILISTIC FRAMEWORK AND MODEL LEARNING In Chapter 3, we described our snake pedal model, which is fully characterized by the shape parameter vector q = (g^,p^)^ with g,p being the global and local shape parameter vectors respectively. In many applications, the general shape, location and orientation of objects is often known a priori and this knowledge may be incorporated in the form of initial conditions, data constraints, model shape parameter constraints, or into the model fitting procedures. For example, in most medical images, the boundaries of the anatomical structures of interest are seldom well delineated from the background, either because the image quality is poor or the shape themselves gradually blend into the surroundings. In many medical imaging applications, it is necessary to use prior knowledge about the shape of interest by learning the shape from a human expert or from a shape atlas. Proper learning schemes must be developed to incorporate the human expert knowledge or the knowledge from an atlas. In addition, for automatic image interpretation, it is necessary to have a model that not only describes the size, shape, location and orientation of the target shape but that also permits expected variations in these characteristics. To meet above requirements, a number of researchers have cast the model fitting process in a probabilistic framework to include prior knowledge of object shape by incorporating prior probability distributions on the shape variables to be estimated 40
PAGE 53
' ' 41 [29, 69, 81, 86]. In Staib and Duncan [69], a deformable contour model was used on 2D echocardiograms and MR images to extract the left ventricle (LV) of the heart and the corpus callosum of the brain respectively. This closed contour model is parameterized using an elliptic Fourier decomposition and a priori shape information is included as a spatial probability expressed through the likelihood of each model parameter. The model parameter probability distributions are derived from a set of example object boundaries and serve to bias the contour model towards expected or more likely shapes. Szekely et al. [71] have also developed Fourier parameterized models. In addition, they have added elasticity to their models to create "Fourier snakes" in 2D and elastically deformable Fourier surface models in 3D. By using the Fourier parameterization followed by a statistical analysis of a training set, they define mean organ models and their eigendeformations. An elastic fit of the mean model in the subspace of eigenmodes restricts possible deformations and finds an optimal match between the model surface and boundary candidates. Cootes et al. [13] and Hill et al. [28] present a statistically based technique for building deformable shape templates and use these models to segment various anatomical structures from 2D and 3D medical images. The statistical parameterization provides global shape constraints and allows the model to deform only in ways implied by the training set. These shape models represent anatomical structures by sets of landmark points which are placed in the same way on an object boundary in each input image. For example, to extract the LV from echocardiograms, they
PAGE 54
42 choose points around the ventricle boundary, the nearby edge of the right ventricle, and the top of the left atrium. These points can be connected to form a deformable contour. By examining the statistics of training sets of handlabeled medical images, and using principal component analysis, a shape model is derived that described the average positions and the major modes of variation of the points. New shapes are generated using the mean shape and a weighted sum of the major modes of variation. Shape boundaries are then segmented using this "point distribution model" by examining a region around each model point to calculate the displacement required to move it towards the boundary, these displacements are then used to update the shape parameter weights. In this thesis, we cast the model fitting process in a Bayesian framework and incorporate prior distributions on the vector of shape parameters that are being estimated. By using a probabilistic approach rather than a deterministic approach, we can develop probabilistic models of data acquisition which closely match the characteristics of real image sensors, and more importantly, we can determine the uncertainty in our estimate, which can be used (if necessary) when integrating multiple sources of data/information [73]. We now present the details of the Bayesian probabilistic framework. f 4.1 Prior Model in the Bavesian Framework As discussed earlier, in many applications such as medical image analysis, the boundaries of anatomical structures are usually "fuzzy" because of the partial volume averaging effects and other effects. We must incorporate the prior information about
PAGE 55
43 the shape of interest into the shape representation. Therefore, in the prior model, in addition to generic smoothness assumptions, we should also incorporate into the prior model the information about the possible shapes that an observed object can assume. We express the probability of a given configuration q through a Gibbs distribution [20], which assigns a low probability to configurations with large energy and vise versa. The configuration corresponding to Ep = 0 represents the rest state of the model. In the Gibbs distribution, the internal deformation energy, which is similar to the internal deformation energy in generalized splines described in Chapter 2, is transfered into a probabilistic distribution by: p(q) = ^exp{Ep{q)), (4.1) where Zp is a normalizing factor and the energy Ep{) governs the deformation of the model away from its natural shape. There are two parts contained in this deformation energy, one part corresponds to the local shape deformation, and the other part corresponds to the global deformation. For the local shape deformation, if we denote the rest shape of the local representation (snake positions) by p, which can be obtained by the scheme discussed in Section 4.4, the internal energy in a continuous form is given by
PAGE 56
44 with the discrete version l(p_p)K(p_p), (4.3) where K is the stiffness matrix for the membrane plus thin plate deformation energy imposed on the snake. Note that Eq. (4.2) and (2.9) have the similar form, if we set p = 0, Eq. (4.2) goes to (2.9). Eq. (2.9) is the special case of (4.2) where the local rest shape is zero. For the global shape representation, if no additional information is made about the shape being modeled, the underlying global shape of the model, the generator, can take the mean value of the shapes that has been learned. By extending the above internal deformation energy, we can include the global shape parameters of the model in our prior statistics in the similar way. The extended energy can be expressed as: ^ _ ^{local) __ ^(global) = lippr K (pp) + l(gsgs)^ C;' (gsgs) (44) + l(gcgc)^ (gcgc). Where gs, ge, and gc are the rest state of the aspect ratio vector, rotation vector, and translation vector of the generator discussed in Chapter 3. Cs''\C0~^, and Cc"^ represent the covariance matrices which are accumulated in the iterative training phase described in Section 4.4.
PAGE 57
45 4.2 Sensor Model in the Bayesian Framework The data acquisition (sensor) model in the Bayesian framework is based on linear measurements with Gaussian noise and given by: p{B/q) = ^exp{ED{D,q)), (4.5) where Zd is a constant, D is the data, and jE'd() can be either the edgebased potential energy synthesized from image data, or the potential energy in the springs attached from 3D data points to the surface of the snake pedal, constraining it to conform to the observed data. These two types of sensor model energies are expressed as: {h Et c{Di Xi)2 spring energy, (4.6) Â— c(Gct * V/(x))p image potential. In the first form of energy Ed{), c represents the uncertainty of the sensor measurements and Xj is the closest model point from the given data Di. In the second form of energy E'd()) ^^(x) is the gradient of the image intensity /(x), is the Gaussian filter mask with width a, and * represents the convolution operation. Taking the partial derivative of with respect to q yields the external force:
PAGE 58
46 The external force pulls the model toward the image boundary or range data, thus molding the model into the desired shapes. This process is known as the model fitting process. 4.3 Maximum a Posterior In the probabilistic framework, combining the prior model (4.1) and sensor model (4.5) by using Bayes rule, we obtain the posterior distribution ^(^ID) . Pj^m^ ^ i.xp(Â£(q)), (4.8) where Z is a constant, and the total energy is: E{ci) = Ep{q) + ED{B,q). (4.9) In the Bayesian framework, the objective is to compute the Maximum a Posterior (MAP) estimate, that is, the value of q that maximizes the conditional probability p(qD). This estimation provides the same result as finding the minimum energy configuration of the physically based model. However, compared with the deterministic model, the probabilistic framework has several advantages as described earlier. 4.4 Supervised Learning Scheme Once we have established the Bayesian statistical modeling scheme, we propose a supervised learning scheme to incorporate the information about the shape of interest to the model. In this learning scheme, we divide the shape recovery process into a
PAGE 59
47 training phase and an operational phase. In the training phase, prior knowledge about the shapes to be recovered from the image data is collected via a supervised learning technique. A training sequence of several iterations provides samples of correctly reconstructed shapes for specific classes of shapes, for example, the caudate nucleus and the hippocampus from the human brain. We use Pg denote the semiglobal part our model, which may correspond to the coarse level wavelet coefficients of the local shape parameter vector p. The parameter sets q(z) corresponding to these shapes represent the training sequence for the global part g and semiglobal part Pg of our prior model as shown in Fig. (4.1). With each iteration i, we improve estimates of the mean vector q and covariance matrix C through a simple equation for accumulating measurements as indicated in Fig. (4.1). Since we update only the statistics of the global parameters of the model, the local deformation remains constrained only by smoothness assumptions captured in the stiffness matrix K. Therefore, most of the flexibility of the locally deformable surface remains unchanged. In the operational phase, the mean and covariance of the model state q are used in initializing the model for fitting the model to a data set not included in the training phase. The main difference between training and operational phase is that, in the former, the model initialization for the fitting is far from the final desired fit because no prior information is available. Most existing schemes in literature that use a priori knowledge about anatomical shapes of interest manually construct their prior models [69, 73] which can be a very tedious process. The model fitting scheme that was described in this thesis in conjunction with the supervised learning mentioned
PAGE 60
Set of training data samples D(z) obtained from sensor Surface fitting by maximizing p(qD(i)) , where current prior model Pi^) is used q(0 final result Iterative improvment of the mean and covariance of the global model q{i) = iq(i) + i^q(z 1) = fqWq^W + ^m2(7 1) Cii) = m,{z) q(z)q^(i) Figure 4.1 Supervised learning to incorporate information into the prior model
PAGE 61
49 earlier and discussed in detail by Vemuri and Radisavljevic [81] can be used for fast automatic construction of prior models. We will devote the next chapter to present the fast numerical algorithms for automatic model construction via fitting.
PAGE 62
CHAPTER 5 NUMERICAL SOLUTIONS FOR MODEL FITTING WITH SNAKE PEDAL MODELS In this chapter we will develop numerical techniques for fitting our snake pedal model to a sparse set of data points in 2D/3D as well as to image data for shape recovery. The objective of model fitting is to mold or sculpt the snake pedal model so that the model passes through the set of data points or conforms to the desired image boundaries as closely as possible, while maintaining a certain degree of smoothness. This fitting process is equivalent to finding the minimum of the total energy (4.9) of the model described in Chapter 4, which consists of two terms, namely, the internal strain energy Ep{) corresponding to the snake deformation energy and the external energy Ed{). The energy minimization is posed as a nonlinear optimization problem and can be solved numerically. There are two types of parameters which need to be determined in this optimization process. One is the global shape parameters, namely, generator parameters g = (a, 6, 6*, mi, 7712)^ in 2D or g = (ai, a2, 03, Vi, V2, V2, s, mi, m2, in 3D; Another is the local shape parameters, that is, snake position p = {pf }^i with M being the number of sampling points on the snake. We collect the global and local shape parameters in a vector q = (g^,P^)^In the model fitting, we need to find the optimum configuration of q so that the total energy in Eq. (4.9) is minimized. The estimation of both global 50
PAGE 63
51 and local parameters in the model fitting makes the minimization problem highly nonlinear. 5.1 Related Work A primary approach to solve energy minimization problems is the calculus of variations. In this approach, finite diff'erence or finite element method is used to solve the EulerLagrange equations derived from the energy functional. This method can be regarded as gradient descent along the potential energy surface of the model. Vemuri et al. [35] used a gradient descent (GD) method that iteratively updates the shape parameter q, ^''^ = (5.1) where q''^^ and q*^ are the status of the shape parameter vector q at iteration k and k + 1 respectively. A is a diagonal matrix containing step sizes and consequently effecting the speed of convergence. Partial derivatives of the total energy ^^^^ can be represented by the sum of the partial derivatives of the prior model and data energies defined in Chapter 4. While the iterative equation (5.1) may eventually converge to the correct estimate, in practice it may be unacceptably slow. In addition, gradient descent does not guarantee finding the global minimum if the energy surface is not convex and the initial contour model is far from the target. Amini et al. [2] use dynamic programming (DP) to carry out a more extensive search for a global minima. Although DP provides numerical stability, the stability is achieved at a high cost in computational complexity. Williams and Shah [85]
PAGE 64
52 proposed an alternative greedy algorithm to DP which drastically cuts numerical costs, however, it does not guarantee numerical stability. Blake and Zisserman [6] propose an graduated nonconvexity algorithm, which is based on iterative approximation of the energy functional by a convex function. It can bridge low ridges in the potential field and eliminate termination in local valleys. Poon et al. [57] and Grzeszczuk and Levin [25] minimize the energy of active models using simulated annealing (SA) method which is known to give global solutions and allows the incorporation of nondifferentiable constraints, however, deterministic algorithms usually preferred due to their faster convergence. 5.2 Our Fast Numerical Algorithms In the following sections, we present two algorithms for solving the energy minimizing problem efficiently. The first one involves a Preconditioned Nonlinear Conjugate Gradient (PNCG) method for estimating the global and local parameters of the model. The preconditioner we use for this case is the diagonal Hessian for the global as well as local parameters. In the second numerical scheme, we first solve for the global parameters using the LevenbergMarquardt (LM) method in the outer loop and employ the alternating direction implicit (ADI) scheme in the inner loop to solve for the local parameters. The first scheme is less sophisticated and easier to program, the second scheme is more sophisticated and relatively difficult to program. The next sections will be devoted to the presentation of the details of each of these schemes.
PAGE 65
53 5.2.1 Preconditioned Nonlinear Conjugate Gradient Algorithm In general, the Conjugate Gradient (CG) algorithm is well suited for solving quadratic optimization problems whereas, the nonlinear version can be used for finding the locally optimal solution of a nonlinear nonconvex function [66]. In this work, we develop a Preconditioned Nonlinear Conjugate Gradient (PNCG) algorithm and apply it to the energy minimization discussed in the previous section. Given any nvariate continuous function /(x) = f{xi,X2,...,Xn) with gradient /'(x) = [af^/(x) 5f;/(x) ... af;;/(x)]^, we can use Nonlinear Conjugate Gradient (NCG) algorithm to minimize it. The residual r, = Â— /'(x,) indicates how far we are from the correct value that minimizes /(x), where i denotes the number of iterations. We use do,di, ,dn_i to denote a set of orthogonal search directions, ai and ^i the step sizes, and we choose a preconditioner M that approximates /"(x) and has the property that M~^r is easy to compute for any vector r [36]. Let imax be the maximum number or iterations and e be the residue threshold for the function /(x), the outline of the Preconditioned Nonlinear Conjugate Gradient method is as follows: 1. Let2 = 0,ro = /'(xo). 2. Calculate a preconditioner Mq ~ /"(xq). 3. do = Mo'ro. 4. Find ai that minimize /(xj Ittidj). 5. Xi+i = Xi + Oidi. 6. Calculate the preconditioner Mi+i ^ /"(xi+i).
PAGE 66
54 7. ri+i = /'(xi+i). = max ( ^ rTM.^ri ' 9. di+i = Ti+i + A+ldi. 10. While [i < imax or \f \ < e) i = i+1, go to step 4. :] In step 8, using the first expression for leads to the FletcherReeves method, and using the second formula leads to the PolakRibiere method. In step 4, the function /(x + ad) is approximately minimized by setting a = Â—^rjw^Where, /"(x) is the Hessian matrix oi the continuous function /(x). H OXiOXi Â£7X2 ax 1 dxjidxi dxnO: In our case, it can be written as OX1OX2 OX2OX2 'X2 dxi d. ^X X2OXJ1 (5.2) H d'^E d'^E dqidqi dqidq2 dqidqn d'E d^E d^E dq2dqi dq2dq2 dq2dqn d'E d'^E d^E (5.3)
PAGE 67
55 Since it is expensive to compute the Hessian at each iteration, to perform an exact line search without computing /"(x), the Secant method is used leading to the following approximation [f (x)]^d ^[/'(x + ad)Fd[/'(x)rd ^^^^ Typically, we will choose an arbitrary small nonzero number for a on the first Secant method iteration; on subsequent iterations we will choose x + ad to be the value of x from the previous Secant method iteration. In other words, if we let a; denote the value of a calculated during Secant iteration i, then (Ti+i = Â— aj. In NCG, r(z) = /'(x,) and GramSchmidt conjugation of the residuals is used to compute the search directions. As for computing /?, two choices are given in the above algorithm. The FletcherReeves technique converges for a good initial guess close to the desired optimum. PolakRibiere scheme often converges much faster however, can cycle infinitely without converging in rare cases. Convergence in this case is guaranteed by choosing P = max{j3^^, 0). In our implementation, we use the PolakRibere scheme. When it's too expensive to compute the full Hessian matrix, it is reasonable to use the diagonal of the Hessian as a preconditioner. However, if the initial guess is too far from a local optimum, the diagonal elements may not all be positive. Since a preconditioner should be positive definite, nonpositive diagonal entries will violate this property and hence can not be allowed. A conservative approach in this case is to use the identity matrix as a preconditioner (M = I), that is, not to precondition.
PAGE 68
56 In the following section, we describe a seminumerical approximation for computing the Hessian of the energy function. 5.2.2 Computation of the Approximate Hessian In this section, we will introduce a seminumerical method for approximating the Hessian of the energy function by dissecting the structure of the Hessian matrix. In the Hessian matrix formula, it can be seen that the elements in one row of Hessian matrix, say row i can be obtained by taking the partial derivative with respect to gi of the vector ()^, where gi is the ith entry of the vector g. Therefore, the ith row of the Hessian has the form: where n is total number of global parameters representing the snake pedal model. In 3D, the subblock corresponding to the snake pedal global aspect ratio parameters gs = [ai 02 aa]^ has the following form: d dE dE dE d_ r % [ 5g (5.5) % L dg2 "' dgn d dE dE dE da\ I da\ da2 da^ dxi dx2 dxs dE dE dE
PAGE 69
57 d'^E d'^E d^E dx\ dxl dx] dai dx2 dal dx^ da: ai J dx[ dxi ~5a2 8x2 3x2 dx] 0X2 8x2, ^3:3 da[ da2 dai , (5.6) where xi,X2,X3 are the world coordinates in which the model is described. Similarly, the part corresponding to ge = [vl v2 v3 s]^ is, d dvi dE dE dE dE dvi Bv^ 5u3 ~Us dE dE 8E dxi dx2 dxs /(^) /(l^) /{^) dvi ^ dvi ' dvi ^ dv2 ' dvi ^ dvz ' dvi ^ ds ' d I dx2 \ dvi ^ dvi ' d (dx2_\ d / dx2 \ dvi ^ dv2 ' dvi ^ dvs ' ov\ ^ ds ' dvi ^ dv\ ' dvi ^ dv2 ' dvi ^ dv^ ' dv\ ^ ds ' + d'^E d'^E d'^E dxf dx'2 dxi dx[ dvi dx2 dvi 3x2, dvi dx\ dx\ dx[ dx[ dvi dv2 dvz 3s 3X2 3x2 3x2 3x2 3vi 3v2 3vl 3s dxrj 93:3 dxs dx^ dvi 3v2 3v3 ds We can get dE dE dE 3E 3vi 3v2 3v3 ~5s 3 3E 3E 3E 3E 3vi 3v2 3v3 ~ds^ , by carrying out the similar procedure. Since we require only the diagonal entries of the Hessian, only the following elements need to be computed, ^(^), ^(Â§)> ^(f^) ^(^). Observe that the computational complexity will be tremendously decreased if only those diagonal elements are required to be computed. From Eq. (4.6) in Section 4, we can obtain the analytical form for computing 3E 3E 3E 3xi 3x2 3x3 and
PAGE 70
58 However, there is no analytical form for ^ and Qg.Qg. : use the forward difference to compute these partial derivatives numerically. Our experiments show that this seminumerical method is very efficient and convenient for the computation of the diagonal of Hessian matrix. 5.2.3 LM Method for Global Parameter Estimation Although the Preconditioned Nonlinear Conjugate Gradient (PNCG) method performs quite well in the global parameter minimization, however, we found that the LevenbergMarquardt (LM) method has larger convergence range and converges faster than the pure PNCG algorithm. Therefore, in our second model fitting scheme, we utilize the LM algorithm to estimate the global shape parameters. In the local shape parameter estimation, by representing the stiffness matrix in a Schur complement formula [53], we use the alternating direction implicit method in the third model fitting scheme to achieve a 0{N) convergence rate in the local parameter estimation, where N is the number of nodes in the finite element discretization of the snake positions. Fig. (5.1) shows the framework of the second and third model fitting schemes [37, 85]. We now first introduce the LevenbergMarquardt method in the global shape parameter minimization. The LM method essentially combines the Newton method and the gradient descent (GD) method in solving the nonlinear function minimization problem. When the current solution is far from the true solution, the LM method adopts the gradient descent (GD) algorithm. On the other hand, when the current solution is close to d'^E d^E d^E 'dxf 'd^ 'dx^
PAGE 71
59 While ( Energy > e ) One iteration of Levenberg Mai'quardt Mcthcxl (Ref : Numerical Recipes) solving global parameters Solve local parameters Aiternating Directional Implicit Method completely (Ref : Lai & Vemuri " 97 LAA) (a) Figure 5.1. Big picture of numerical schemes for model fitting. the true solution, the nonlinear function can be approximated by a quadratic form and the LM method switches to the Newton method. We now adopt the LM method to our model fitting problem. Assume that we only need to fit the snake pedal model to a set of 3D data points, then the external energy used here becomes the springbased deformation energy Ed(D, g, p) introduced in Section 4. More specifically, to fit the snake pedal to a given data set of m points {D;}^^ in 3D, we minimize the sum of the squared distances from each data point Dj to the corresponding closest point Xi on the pedal curve/surface. Let h be the function that returns the closest point on the snake pedal from a given data point D,, that is, x^ = hi(g,pi,ij). Let {pi}^^ be a set of M discrete points along the snake p{t), {xi}fii be the corresponding points on the snake pedal curve. Let D = {Df and p = {pj}fii, the quantity to be minimized for the model fitting
PAGE 72
60 can be written as m ED(D,g,p) = ;^ c(D,x,)^ (5.7) where Xj = hi(g, p,, tj), and the magnitude of parameter c is related to the uncertainty of data measurements. In LM, the gradient of ED(D,g,p) with respect to generator parameters g = {91,92, V = {aua2,a3,vi,V2,V3,s, mi, 7712,1713)'^ in 3D, and g = {91,92, V = (ai, 02, 9, nil, 1712)^ in 2D has components: dE D dgk c(Dihi(g,pi,i,)) 1=1 dhi{g,pi,ti) dgk (5.8) If we assume that there is a onetoone correspondence between the points on the generator and the pedal snake then the partial derivative can be written as: c (Dk hk(g,Pk,4))dgk d9k (5.9) Taking an additional partial derivative gives Qg^^ which is needed in the LM iteration. For the derivative of the energy required in the LM algorithm, we use the chain rule to get = (%f )(Â§)More specifically, if we denote each component hj of h as hi = {hii,hi2Y , then in 2D as Dii hii Di2 h i2 dhj] dhu dhji dh^ dhu da db dU drrli 81712 dhn dhg dhg dhi2 dh,,2 da db 36 dvrii dml . (5.10)
PAGE 73
61 If we define Pk = ^ > (^kl= ^ , 5.11 using the Newton method to minimize the energy function gives J2^ki5gk = Pk(5.12) On the other hand, using the gradient descent method gives Sgi = const, x /3i = (5.13) where we choose the constant to be ^ in Eq. (5.13). This choice of step size in Eq. (5.13) leads to a preconditioned gradient descent method with a diagonal Hessian as the preconditioner. Furthermore, we define 4 = Â«J^(1 + ^)' (514) "jfc = "jfc (Jt^^), (5.15) and combining Eq. (5.12) and (5.13) yields ^oi'f.i6gi = pk, (5.16) When A is very large, Eq. (5.16) becomes identical to (5.13); and as A approaches to zero, Eq. (5.16) becomes the same as Eq. (5.12).
PAGE 74
62 From above discussion, given an initial guess for the set of parameters g, the LM method for the global parameter estimation is described as follows: 1. Compute Ed(D, g, p), 2. Pick a suitable value for A, for example, A = 0.001, 3. (t) Solve the linear Eq. (5.16) for Sg and evaluate ED(D,g + 5g,p), 4. If Ed(D, g + 6g, p) > Ed(D, g, p), increase A by a factor of 10 (or any other substantial factor) and go back to (f), 5. If ED(D,g + (^g,p) < ED(D,p,q), decrease A by a factor of 10, update the trial solution g < Â— g + Sg, and go back to (f), until the energy Ed(D, g, p) is less than a predefined threshold . After the estimation of the global parameters, we solve the local shape parameters in the inner loop of the parameter estimation using the Preconditioned Conjugate Gradient method or the alternating directional implicit method. 5.2.4 API Method for Local Parameter Estimation When solving for the local parameters in the inner loop for fixed global parameters, the quantity to be minimized has a quadratic form and its minimization leads to solving a linear system, namely, Kx = b, where K is the stiffness matrix, X is the local parameter vector, and b = is the external force. Note that we use X to represent the local shape parameter vector and b to represent the external force, (they were represented by p and f in previous chapters.) which is the standard
PAGE 75
63 method to represent unknowns and righthand side in linear algebra. Because of the special structure of the matrix K in this linear system, we can use the alternating direction implicit (ADI) method in the local shape parameter estimation. We can prove that the ADI takes only 0{N) time to converge with soft data constraints where N is the number of nodes in the rectangular finite element discretization of the snake in our model. Note that in soft data constraints, the model is not required to interpolate the data but is expected to approximate it. We now present the structure of the stiffness matrix K. In the snake pedal model, the controlling snake p is a controlled continuity spline under the influence of image forces or other external constraint forces. In our model, the spline deformation energy imposed on the snake is a membrane plus thin plate energy given in continuous form by where Wi and W2 are the tension and rigidity factors respectively. The external energy used can be the same as the energy Ed(D, q, p) discussed in Section 5.2.3. In this case, the external force impose on the snake is given by ^,(P) = /Â».(()^ + ()^) + Mi^? ^ 2i^/ + (0)'^)
PAGE 76
64 South Pole Total number of nodes n X n + 2 Denote : N = n X n North Pole (a) Figure 5.2. Model tessellation in the model coordinates. where dhi2 2 I 2 "1 i Â— r> 17 n( + n2 (5.19) with ni,n2 defined as before. Notice that dEo/dhi is the external force in a standard snake, but in our case, the external force is dEo/dpi, which is related to the force in a standard snake via a Jacobian matrix dhi/dpi in Eq. (5.19). We then discretize the model in model coordinates using finite elements. Geometrically, the snake pedal surfaces are closed surfaces in space whose intrinsic coordinates u = {u, v) are defined on a domain 1]. The positions of points on the model in an inertial frame of reference are given by a vectorvalued function pedal(u) = (a;i(u),Z2(u),X3(u)). The tessellation of u = {u,v) in intrinsic coordinates is shown in Fig. 5.2. We use a combination of quadrilateral and triangular
PAGE 77
65 elements in the discretization. In the inner quadrilateral mesh, we have {nx n) = N nodes and 2 additional nodes are used in representing the north and south pole leading to a total of {N + 2) nodes. Due to the presence of the north and south poles, after the discretization of the model in intrinsic coordinates u using finite elements, the stiffness matrix K E i?(^+2)(Ar+2) j^ ^^j. problem can be written as a (2 x 2) block matrix: K = Kii Ki2 K21 K22 (5.20) where Kn G fi^^^ contains coefficients which corresponds to the inner quadrilateral mesh and K22 Â£ R^^^ is a diagonal block matrix containing the coefficients corresponding to the north and south poles, and K12 = Kji contains the coefficients of interaction between the poles and other nodes. To solve the linear system Kx = b for the local parameter estimation, we employ the following Schur complement formula [53] K^b Kr/ + Kr/Ki2siK2iKr/ Kr/Ki2si (5.21) where S Â— K22 Â— K21K1/K12. (5.22)
PAGE 78
66 From Eq. (5.21), it is evident that to obtain the solution of K~^b, we need the solution of several K^i bs. Further analysis (see Appendix A) shows that the component Kii of the stiffness matrix K can be expressed as the sum of two Kronecker products: Kii = u;i(A(8)B + B(8) A) + i(;2A(8) A = A(8)C + CÂ®A, (5.23) where C = + {w2/2)A with A and B being (n x n) tridiagonal and symmetric positive definite (SPD) matrices of the form 1 1 2 1 1 2 1 1 4 1 A = 1 2 1 1 4 1 1 1 2 (5.24) The Kronecker tensor product or Kronecker product of C and A, denoted by C (8) A, is the matrix of order N x N (where N = n x n) and given by
PAGE 79
67 ciiA C12A ClnA C21A C22A C2nA C0A = (5.25) CjiiA. CJ12A. ... CfijjA When computing with Kronecker products, it is computationally convenient for us to represent vectors using matrices. When solving (C (E> A)x, we represent the vector X of order Â•n? by the matrix X = {xki} of order n x n defined by It can be shown that this representation can lead to efficient procedures for computing (C (g) A)x and solving (C (g) A)x = b, respectively. LEMMA 1 Let C = {cki}, A = {cfc/} and X = {xki} be matrices of order n x n, then the n X n matrix representation of x 1 vector (C (g) A)x is given by Xkl Â— Xk+n{ll)(5.26) (CO A)x (C(AX)T)T. Using this representation, the linear system Knx = b can be written as: CXA^ + AXC^ = D, (5.27)
PAGE 80
. 68 where D,X are the n x n representation of the n'^ vector b and x. Since A and C are also SPD matrices, Eq. (5.27) can simply be written as CXA + AXC = D. (5.28) After some simple linear algebra operations, it follows that above equation is equivalent to A'X + XA' = D', (5.29) with A' = A^C and D' = A'^DA^ Eq. (5.29) is the standard Lyapunov matrix equation [43]. The ADI method can be applied to solve this matrix equation and has the following steps in each iteration J = 1,2,..., J, (A' + p,I)X^._i = D'X,_i(A'p,I) (5.30) X,(A'+p,I) = D'(A'p,I)X^_,. (5.31) For each iteration in ADI, we need to solve two equations, namely, Eq. (5.30) and (5.31), and therefore X:! is introduced as an intermediate variable. Since A' is a SPD matrix, we can use Cholesky decomposition for the tridiagonal SPD matrix (A'+Pjl) to compute the update solution. This only takes 0{N) time per iteration. The parameters pj are called ADI parameters and J is the number of iterations needed. They are chosen according to the method described by Lai and Vemuri [36].
PAGE 81
69 In order to fit the snake pedal curve to salient features in an image, we use the following equation as the energy to be minimized: E = (G,*V/)r = (G,*V/(x,2/)p. (5.32) , Where (x, y) is on the snake pedal curve, is a Gaussian with standard derivation (7, and * denotes the convolution operator. More sophisticated boundary finding energies can be synthesized and we refer the reader to Kass et al. [31]. By minimizing the above energy, the model is made to fit to points of high gradient in the image. The procedure to minimize above energy is similar to that of fitting of snake pedal curves/surfaces to data points.
PAGE 82
CHAPTER 6 HYBRID GEOMETRIC ACTIVE MODELS Since the inception of snakes in the vision/graphics community by Kass et al. [31] and Terzopoulos et al. [76], these elastically deformable contours/surfaces have been widely used for a variety of applications, such as edge detection, segmentation, motion tracking, and shape modeling. The classical approach to object shape recovery using the snakes is based on deforming an initial configuration of the snake Vo towards the boundary of the object to be detected. As described in Chapter 2, snakes are energyminimizing splines, their deformation under constraints is achieved by minimizing a functional, that can be regarded as the bending energy stored in a thin flexible beam/rod or a stretchable string subject to loading. There are several problems associated with this approach, such as initializations, convergence to the local minimum, and the specification of the physical or elasticity parameters. Moreover, this energy model requires that the topology of the shape to be estimated be known a prior. However, in many applications, such as motion tracking, it is not always possible to specify the topology of an object prior to its recovery. For example, during the cell mitosis, the close contours of cell boundary change connectivity or split, therefore undergoing a topological transformation. Hence it is necessary to develop a model which has the ability to split and merge freely during the 70
PAGE 83
71 shape recovery without making the a prior assumption about the object's topology in these applications. Recently, Caselles et al. [9] and Malladi et al. [45] proposed a novel geometric model of active contours dubbed geodesic/ geometric active contours or geodesic/ geometric snakes. These models are based on the the theory of curve evolution and geometric flows. Automatic changes in topology can be handled in a natural way in this modeling technique when we implement the curve evolution using the levelset embedding schemes. In essence, this model is and implicit surfaceevolution model where the surface propagates with a velocity which is directly related to the curvature of the evolving curve which is embedded in the surface. Unlike the classical energy models, these models are not the result of minimizing an energy functional. In this framework, boundary detection can be considered equivalent to finding a curve of minimal weighted length. More recently, several researchers independently demonstrated the unity of the curve evolution approaches and the classical energy minimization methods [10, 32, 64]. The challenge in the geometric active contour evolution is to devise numerical schemes for the equations of the evolving curve, which will accurately approximate solution that can be very unstable due to formation of shocks or discontinuities. Osher and Sethian [52] developed special numerical techniques based on upwind finite difference schemes leading to a very stable solution. They introduced levelset techniques to seamlessly handle any topological changes of the curve/surface that may
PAGE 84
72 occur during evolution. The equation of motion for the embedded curve is reminiscent of an initial value "HamiltonJacobi" equation with a parabolic righthand side and is closely related to a viscous hyperbolic conservation law. For details on the theory of curve/surface evolution and its levelset implementation, we refer the reader to the relative papers [9, 10, 33, 46, 51, 45, 65, 74]. 6.1 Overview of Hybrid Geometric Active Models In many computer vision applications, characterizing the global shape of an object is very useful, for example, in object recognition. Traditional geometric active contour/surface models do not possess the capability to characterize the global shape of an object. They are however useful in characterizing the local shape details in an object. In this section, we introduce a novel concept of a global/core model into the PDEbased curve evolution framework by embedding the snake pedal model into a levelset framework. Instead of characterizing an object boundary by the position of every point on the boundary, our proposed model, referred to as the hybrid geometric active model, now describes an object as a combination of a global/core shape, for example, circle, superellipse, and so on. and a variable offset defined with respect to the global shape. The variable offsets are controlled by this global shape and an evolving curve Â— the controlling snake in the snake pedal model. Augmentation of the global shape/core in the model representation is very useful in object recognition and image indexing applications. In these applications, the global description, such as center, size, and orientation, of a group of objects or structures needs to be retrieved initially which is then followed by a retrieval of the local description of each
PAGE 85
73 object (if necessary). For the model to recover the object boundary, we introduce a robust and efficient numerical method which consist of a global plus local shape estimation technique in the model fitting. We use the LevenbergMarquardt (LM) described in Chapter 5 for estimating the global shape and a combination of upwind and minmod finite difference schemes in a levelset framework for estimating the local shape. The LM method has a broad convergence range, converges to the local optimum very quickly, and therefore is well suited for the global shape recovery. The hybrid geometric active contour/surface model retains all the advantages of traditional geometric active models (for example, topology change, ability to model complex geometries and amenability to stable numerical implementation) and has the added ability/ advantage of being able to compactly represent global shape of an object. We first present the theoretical foundation of the traditional curve evolution in Section 6.2, followed by the levelset implementation framework of the PDEbased curve evolution in Section 6.3. We then introduce our hybrid geometric active models in Section 6.4. The numerical issues of the model fitting process will be discussed in Chapter 7, followed by the implementation results in Chapter 8. 6.2 Curve Evolution Theory The mathematical foundation of the geodesic or geometric active contour model is based on the Euclidean curve shortening theory, which guides the curve evolution in Euclidean space [15]. Let p = {px.PyY denote an arbitrary point on the geometric active contour, the family of planar curves flow according to the geometric heat
PAGE 86
74 equation or curve shortening equation where k is the curvature and N is the inward unit normal. When the curve evolves according to (6.1), the Euclidean perimeter shrinks as quickly as possible. We demonstrate this concept in the following. Let 7^(0) be a smooth, closed initial curve in Euclidean plane 3fJ^, let V{t) be the oneparameter family of curves generated by moving 'P(O). Let p{s\t) be the position vector which parameterizes V{t) by s', where 0 < s' < 1 and t is the time factor. For closed curves, we assume that p(0, t) = p{l,t) for any time t and a similar condition is posed on the first derivatives. We define the length functional: (6.2) where {. denotes the standard 2norm. Taking the first variation [19] and using integration by parts in Eq. (6.2), we obtain dp d'^p (6.3) as'" 1 d dp ds' 11011 > ds',
PAGE 87
75 where <, > is the inner product of two vectors. Note that the arc length ds can be expressed by using the definition of curvature k and normal N [7, 8, 49], the integral in (6.3) becomes Setting equation (6.5) to zero, we obtain the direction in which C{t) decreases most rapidly as It has been demonstrated [16, 17, 18, 24, 60] that when simple closed curves evolve according to (6.1) without developing singularities, they will converge to "round" points. High curvature regions will be smoothed out during the evolution. We can easily understand this phenomenon by further examining Eq. (6.1). In (6.1), the planar curve evolves in the direction of its inward normal, and the rate of evolution is proportional to its local curvature. In locations where the curve has a larger curvature, the rate of change is large, thus the curve evolves faster in the high curvature points. For the same reason, the curve evolves relatively slower at the low curvature points. Therefore, the portions of the curve in high curvature regions undergo more changes toward the rest state, while the portions in the low curvature regions remain smooth until all portions of curve reach the rest state Â— a circle. ds = Â— 7II ds', (6.4) (6.5)
PAGE 88
76 Since the curvature k determines the rate of change in (6.1), it is regarded as the speed function of the curve in the evolution. More general speed functions can be allowed in the curve evolution. Let F{k) denote the speed function, the curve then evolves according to ^ = nm(6.6) In general, F{k) can be split into two terms: F{k) = FA + FG, (6.7) where the first term Fa is referred to as the advection term and is independent of the geometry of the curve. The curve uniformly expands or contracts with speed Fa depending on its sign, the behavior of the curve is analogous to that of balloon model [12] with an inflation/deflation force. The second term Fq is referred to as the diffusion term and is dependent on the geometry of the curve, such as the local curvature. The diffusion term smoothes out the high curvature regions of the curve and has the same regularization effect on the curve as the internal deformation energy term in the generalized splines described in Chapter 2. When implementing the curve evolution equations such as (6.6), there are a number of numerical considerations. One numerical approach is first taking the above Lagrangian description of the problem, producing equations of motion for the position vector p{s',t), then discretizing the parameterization with a set of discrete markers lying on the moving curve. These discrete markers are updated in time
PAGE 89
77 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 [61]. First, small errors in the computed marker positions are tremendously amplified by the curvature term, and calculations are unstable unless an extremely small step is employed. Secondly, singularities may develop (when Fa 7^ 0) in the propagating curve even with a smooth initial curve. A typical example of this phenomenon is when the smoothing curvature (viscous) term is absent, namely, F{k) = F4 in (6.6), and an entropy condition must be imposed to extract the correct weak solution. Thirdly, it is difficult to handle the topological changes if the evolving curve breaks and merges. Finally, significant bookkeeping problems occur in the extension of this technique to three dimensions. 6.3 Traditional LevelSet Approach Osher and Sethian [52, 61] proposed an algorithm based on the theory of viscosity solutions to PDEs . The algorithm provides a stable numerical solution for curve/surface evolution. In their approach, the curve is first embedded in a two dimensional surface, the equations of motion are solved using finite diflference techniques and hyperbolic conservation laws. 6.3.1 Embedding Mechanism The embedding step can be carried out by representing the curve V{t) by the zero level set of a smooth and Lipschitz continuous function ^ : x [0, r) > As an illustration of this approach, we consider the example of an expanding circle. Let the initial curve P at ii = 0 be a circle in the xyplane (Fig. 6.1 (a)). We imagine
PAGE 90
78 X Figure 6.1. Level set formulation of equations of motion (a) and (b) show the curve V and the surface (i){x,y) at i = 0, (c) and (d) show the curve V and the corresponding surface 0(a;, y) at time t. (e) and (f ) show the curve V and surface
PAGE 91
79 that the circle is the level set {0 = 0} of an initial surface z = (f){x,y,t = 0) in dt^ (Fig. 6.1(b)). We can then match a oneparameter family of moving curves Vit) with a oneparameter family of moving surfaces in such a way that the level set = 0} always yields the moving curve (Fig. 6.1 (c) and (d)). 6.3.2 Equation of Motion In the general case, let 7^(0) be a closed, nonintersecting, N Â— 1 dimensional hypersurface. Let 0(p,i), p 6 be a scalar function such that 0(p,O) = Â±d{p), where d{p) is the signed distance from p to the hypersuface 'P(O). Assume that 0 is negative in the interior and positive in the exterior of the zero level set. Let 7^(0) be a planar curve, represented as the zero level set of (j){p,t), that is, V{t) e : 0(p,t) = 0, (6.8) By differentiating (6.8) with respect to t we obtain: V(j>{p,t)^ + MP,t)=0. (6.9) Note that the gradient V(j){p, t) is normal to the N 1 dimensional level set passing through p in general. When A^^ = 3, for the zero level set, the following relation holds:
PAGE 92
80 where N is the curve normal as defined before. In this equation, the left hand side is in terms of the surface (f), while the right hand side is in terms of a property of the curve V. Combining Eq. (6.1), (6.10), and (6.9) yields: (j>t = F{k)\\V^\\. (6.11) The curve V, evolving by Eq. (6.6), is obtained by the taking the zero level set of the function 0, which evolves by (6.11). This scheme is referred to as Eulerian formulation for curve evolution, because it is written in terms of a fixed coordinate system. Also, for an implicitly defined curve, the curvature can be computed as Morgan [49] or explicitly as Eq. (6.11), together with the initial condition, (f){p,0) = Â±d{p), is the governing equation for the curve evolution in a levelset framework. Since this equation resembles the HamiltonJacobi equation, it is referred to as the "HamiltonJacobi" formulation [52]. To attract the active contour to the desired image features, the speed function in (6.11) is multiplied by a quantity:
PAGE 93
81 where / is the greyscale image and Ga is the Gaussian smoothing filter with width (7, and n = 1 or 2 depends on the application [44, 83]. The function g{x,y) depends on the given image and is used as a "stopping term." For example, when the curve is close to an edge, g{x,y) becomes very small and retards the curve evolution. A modification of (6.11) based on the use of Riemanian metric was reported in Caselles et al. and Kichenassamy et al. [10, 32], it involves the minimization of a weighted Euclidean length and takes the form: (l>t = g{x, y) F{k) II + V0 Â• V^, (6.15) Note that for g as in (6.14), Vg behaves like a "doublet" near an edge. The effect of Vg is to attract the evolving contour as it approaches an edge, and to push the contour back if it passes the edge. The inclusion of Vg term allows the use of larger inflation forces, resulting in features being extracted in relatively fewer time steps. 6.3.3 Advantages of the Levelset Approach By using the level set approach discussed above, topological changes can be handled naturally. This advantage lies in the fact that in the curve evolution, the higher dimensional function 0(p, t) always remains a function, while its level surface {0 = 0}, which corresponds to the evolving curve p{t), needs not to be simply connected. The moving curve can change its topology or form sharp corners easily. As shown in Fig. 6.1 (f), (f>{p,t) remains a simple continuous saddlelike function, while
PAGE 94
82 its level set {0 = 0}, the moving curve, splits into two separated closed contours, thus undergoes a topological change from its initial configuration shown in Fig. 6.1 (a). Another advantage of the level set approach lies in the fact that the geometric and differential properties of p{t) can be computed from 0. For example, the curvature of p{t) can be computed in terms of the derivatives of 0 expressed in Eq. (6.13). In addition, the level set approach can be easily extended to higher dimensions and appropriate expressions can be obtained for the mean curvature and the Gaussian curvature. Moreover, the governing equation (6.11) is reminiscent of an initial value HamiltonJacobi equation with viscosity, and can be solved using the stable, entropysatisfying finite difference schemes in accordance with the hyperbolic conservation laws [40, 52, 62]. , 6.4 Our Model: Hybrid Geometric Active Contours As described in the former sections, the level set approach provides an elegant way to handle topological changes with active models, however, traditional PDFbased approaches do not possess a mechanism to characterize the global shape of an object. In this section, we describe how the level set formulation of the curve evolution can be applied to our snake pedal model to realize topological changes (when necessary) as well as capture the global shape representation of an object. 6.4.1 Big Picture In our approach, to build some amount of regularization on the snake pedal, we impose the regularization constraint via Euclidean arclength minimization on the snake pedal. This would lead to the standard geodesic active contour. Let us
PAGE 95
83 consider a snake pedal denoted by Ve{t), with its position denoted by Pe = {pei,Pe2}) then the standard curve evolution formula for the snake pedal can be written as ^ = F{k,)N,, (6.16) where ke and Ng are the curvature and normal of the snake pedal respectively. We would like to know how the "snake" would evolve if the snake pedal were evolving as a function of its local curvature. Indeed, we will not evolve the snake pedal curves directly, instead, we will first solve the "snake" position under the constraints that the arclength of the snake pedal is minimized, and the pedal curves can then be determined by the pedaling operation defined in Chapter 3, given the position of the generator. The problem therefore can be solved by the following procedure: 1. Derive the curve evolution equation for the "snake" 'P{t) by minimizing the arclength of the snake pedal 'Pe(t); 2. Embed the evolving curve Vit) in a higher dimensional surface 0i and formulate the equation of motion for
PAGE 96
84 z P e (t) = Level Set ^^^=0 ;?(t) = Level Set
PAGE 97
85 of clustered obj internal deformation energy, but we will still use the name "snake" to represent the controlling active contours in the snake pedal model. Fig. 6.2 depicts this important relationship between 'P{t),Ve{t),(j)i, and (j)2. Note that an important feature of this modeling technique is the incorporation of a global parameterized shape, namely, the generator, into the curve evolution framework. This global parameterized shape can be very useful in applications involving recognition of clustered objects. For example, in Fig. (6.3), there are three clusters of objects in 2D. Each cluster has its center, size and orientation, which can be represented by an ellipse. The snake pedal model can capture this "global orientation" of the cluster via the global part of the model, namely, the generator, while recovering every single object in the cluster. Since a geometric active contour is used to represent the controlling "snake" in the model, and also the model can capture the
PAGE 98
86 "global orientation" of a group of objects, this snake pedal model is referred to as a hybrid geometric active contour model. Detailed discussion on the related concepts will be given in the following sections, starting with the derivation of the relationship between the snake and the snake pedal evolution. 6.4.2 Relation Between the Snake and Snake Pedal Curve Evolution Consider a snake p = (pi,P2)^ and a generator a Â— {ai,a2)^ with normal = (rii,n2)^, where p and a are related by the radiusbased correspondence associating scheme described in Chapter 3, we obtain the corresponding snake pedal Pe = {Pei,Pe2)^ via the following pedaling operation: here we use the modified pedaling operation given in Chapter 3 for the purpose of explanation. Rewriting the Eq. (6.17) in component form, we obtain Pe (6.17) (6.18) or simply. Pe = J4 P b, (6.19) where 3a 1 + (6.20) 7l{+ 7I2
PAGE 99
87 and b = a 1 Â• ni + Qo Â• 712 (6.21) Taking the inverse of J^, denoted by J^, we have Â— 2 (6.22) We should remind the reader that by using the association scheme described in Chapter 3, the forms of and b do not change over time, but their contents do change. But, for the following reason, we can still treat as a constant matrix and b as a constant vector with respect to time when evolving p(t) or Pe(t). First, the generator does not evolve over time, we only estimate the generator parameters iteratively; Second, the changes of and b are very small compared with that of p or Pe. Therefore, Eq. (6.19) can be explicitly written as a function of time: Taking partial derivative on both sides of Eq. (6.23) with respect to time t yields Pe(t) = P(^) b. (6.23) dt dp{t) dt (6.24) or dpit) dPejt) dt (6.25) dt
PAGE 100
88 Eq. (6.24) and (6.25) reveal a very crucial relationship between the snake and snake pedal evolution: The evolution of two curves are linearly related by a Jacobian matrix. We remind the reader that the relationship between the snake and snake pedal is not linear, but their evolutions over time are related by a linear transformation! This linear transformation between the evolution of the snake and the evolution of the snake pedal largely simplifies our further derivation and computation, as seen in the following sections. 6.4.3 PDE for the Snake and Snake Pedal Evolution With the relation between the snake and the snake pedal curve evolution given by (6.25), we now derive the equation of motion for the higher dimensional surface 01, in which the snake is embedded. As discussed in Section 6.3.2, the zero level set of the higher dimensional surface 0i is represented by Eq. (6.8), {V{t) e 3f?' : {p,t) = 0} (6.26) Differentiating Eq. (6.8) with respect to t yields + V(/)i ~ = 0. (6.27) at Similarly, for the higher dimensional function (j)2, we have the following formula
PAGE 101
89 From the discussion in the previous section, we have d or (6.29) As described earher, we want to minimize the arclength of the snake pedal, which amounts to imposing a smoothness constraint on the snake pedal curve, hence we require that dt (6.30) where F{ke) is the speed function, which depends on the curvature of the snake pedal kg. Ne is the unit normal of the snake pedal. Similar to the relation between the normal of the snake and its higher dimensional function (/), : N = m^^, for the iiv0iir zero level set of 02, the following relation holds: IIV02 (6.31) Combining Eq. (6.28), (6.30) and (6.31), we obtain d(t>2 dt = F{k,)\\Vci>2\ (6.32)
PAGE 102
90 Eq. (6.32) is the standard levelset evolution. For the governing equation of the snake curve, we substitute Eq. (6.29), (6.31), and (6.30) into Eq. (6.27) to obtain ^ + V01 Â• {JBF{ke)Ne) = 0 (6.33) Substituting Ne in Eq. (6.33) with Eq. (6.31) yields the following expression For simplicity, Eq. (6.34) can be rewritten as Eq. (6.32) and (6.35) are the equations of motion for 0i and (^2, respectively. We rewrite them together as follows: (6.36) Eq. (6.36) represents a coupled equation for the evolution of i and (f)2. One approach to solve the PDE system (6.36) is to use a combination of central differences and the finite upwind diflFerence method [52] to solve the first equation in every iteration, and then solve the second equation using a similar method subsequently. We will discuss the upwind finite difference method in Chapter 7. This approach is straightforward,
PAGE 103
91 but needs the memory storage for both 0i and 02 at every grid point. Note that in Eq. (6.36), only the gradients of 01 and 02 are involved on the right hand sides of both equations, we propose a more elegant approach to solve the PDE system of (6.36) simultaneously with much less storage requirement, by employing the intrinsic relation between V0i and V02 as discussed in the next section. 6.4.4 Relation Between V0i and Vc!>9 To discover the relation between V0i and V02, at first, we arbitrarily choose representative points on the snake. Let p = (x, yY denote the coordinates of point a. The other two representative points h and c are chosen by shifting point a in x direction and in y direction, such that p = {x+/S.x, yY and p = (x, y+/S.yY are their coordinates, respectively. Let a', 6', and c' be the corresponding points on the snake pedal that are obtained by applying the pedaling operation with respect to point a,6, and C, with coordinates p^ = {Pex^PeyV , Pe = {Pex,PeyY and p^ = {Pex,PeyV^ respectively. More specifically, (6.37) Pex Â— Pex{x + Ax,y) < (6.38) Pey Â— Pey{x + Ax,y), and (6.39)
PAGE 104
92 (x,y+Ay) (x,y) (x+Ax,y) a, b, c : points on the snake a', b', c' : points on the snake pedal (a) Figure 6.4. Three representative pointpairs on the snake and snake pedal. The higher dimensional function values at points a, b, and c are given by 0i, 4>i and 4>i respectively. We construct the higher dimensional surfaces (/>i and ^2 in such a way that the zero level sets of (pi and 02 represent the snake and the corresponding snake pedal curve, respectively. The same relation holds for the other level sets of (^1 and (j)2. Therefore, the higher dimensional function values of (f)2 at position a', b' and c' are 0i, ^1, and 4>i, respectively. Since (^, ^)'^ is the gradient of
PAGE 105
93 (6.40) and from a' to c' the second relation holds: df2 dx (PexPex) + ^{PeyPey) =
PAGE 106
94 Since Ax = 0 from a' to d , PexPex = Pey Pey Ay. (6.45) Substituting (6.41), (6.43), and (6.45) in to (6.40) yields: ' ^ . dp^Ax + ^ Â• ^Ax = I ox ox oy ox OX dd dy (6.46) Eliminating Ax and Ay on both sides of Eq. (6.46) results in dy dx dx ex , 902 . dPey ^ 90^ dy^ ~By 'By ox ox d^.dp, ox oy (6.47) Writing in the matrix form, we obtain 901 dpex dx OX d(l>i dpex L dy J L dy dPey dx 902 dx dpey dy J 902 L dy J (6.48) or simply, V01 = 3a V02, (6.49) or V02 = 3 b y(f>u (6.50)
PAGE 107
95 where the Jacobian between the evolution of the snake and snake pedal curve and Jfi = J^^ are defined in Eq. (6.20) and Eq. (6.22), respectively. After a simple yet tedious derivation, we obtain a surprisingly simple relation between the V0i and V(?!>2. At first glance, this relation is seemingly contradictory to between the gradients of 0i and (j)2, the relation should be V02 = J>iV^i, rather than * t * The development of the relationship between V(/Â»i and V02 is a crucial step in the derivation of the evolution of (pi , which will be given in the next section. 6.4.5 PDE for the Higher Dimensional Function Since our objective is to obtain the equation of motion for the higher dimensional surface (pi, we substitute (f)2 with Jb 0i in Eq. (6.35) to get Note that we use the symmetric property of Jb in the above derivation from the second step to the third step. The representation of kg in terms of (pi is quite complicated, we refer the reader to Appendix B for these details. Eq. (6.51) is the equation of motion for (pi and constitutes the primary equation in our application. When using the snake pedal model for recovery of the shape of our imagination. One would think that since 5^ = J^^l^, if there is any Jacobian (6.51)
PAGE 108
96 interest from image data, we need to solve the more general equation of motion for the higher dimensional function (^i. From Eq. (6.15), by using a similar approach, we obtain the more general equation of motion: ^ = y) F{K) II Jb V<^i II + Vc^i Â• Vp, (6.52) where g{x, y) is an image feature based function and is used to stop the curve evolution when the contour is close to the desired edges. V
PAGE 109
CHAPTER 7 NUMERICAL SOLUTIONS FOR SHAPE RECOVERY WITH HYBRID GEOMETRIC ACTIVE CONTOURS In Chapter 6, we discussed that how the snake pedal may be used to recover the object boundaries from an image using a novel global plus local shape estimation procedure. For the global shape estimation, we employ the LevenbergMarquardt (LM) method; while for the local shape estimation, we adopt a modified levelset method. We have discussed the LM method in Chapter 5, we will now discuss the modified levelset framework for local shape recovery. In Chapter 6, we illustrated that the solution for the snake pedal evolution can be achieved by first solving the snake evolution, which is embedded in a higher dimensional surface (^i, then applying the pedaling operation on the snake. Therefore, developing reliable and efficient numerical algorithms for solving the governing equation of that is, Eq. (6.51), becomes one of the main tasks in applying the snake pedal model for extracting shapes of interest from image data. The governing equation of motion for 0i, in the simplest form is given in Eq. (6.51), we rewrite this in the following for the purpose of discussion, ^ = F(A;e) JbV
PAGE 110
98 where Js is the Jacobian between the snake and snake pedal evolution. In our implementation, we assume that is a constant matrix with respect to time based on the reasons stated in Section 6.4.2. kg is the curvature of the snake pedal, and F{.) is the speed function. To discuss the numerical behavior of (7.1), we will first discuss the solution for Eq. (7.2) is equivalent to (7.1) if Jb is an identity matrix and kg is equal to k, the curvature of the snake. Eq. (7.2) is the standard levelset equation and its solution has been widely discussed in literature [52, 62]. We will show that the solution technique for (7.2) with minor modifications can be apphed to solve (7.1). We will then present some shape recovery results by applying the snake pedal model to image data in Chapter 8. 7.1 Numerical Algorithm for Standard PDEbased Curve Evolution In Chapter 6, we discussed the eff'ect of the advection term Fa and the diffusion term Fg in the speed function F{k). Fig. 7.1 depicts the curve evolution with two different speed functions. In Fig. 7.1 (a), the curve moves with unit speed, while in Fig. 7.1 (b), the speed depends on the curvature of the curve. To illustrate the numerical behavior of (7.2), we substitute F{k) = 1 eA; as a typical speed function and the equation of motion becomes F{k) \\Vcf> (7.2) dt {Fa + Fg) \\\/cf>,\\ = {lek) dt (7.3)
PAGE 111
99 Constant speed Curvature dependence F(k)=l F(k)=lÂ£k (a) (b) Figure 7.1. Curve Evolution: The curve moves with unit speed in (a), and moves with curvature dependent speed in (b). The soUd Unes are the initial curves and the dashed lines are the curves after evolution. where we use Fa = 1 and Fq = ek. Eq.(7.3) poses an initial value problem: 0u + (<^L + 0U^/^ = eV(^), (7.4) with the initial condition: (l)i{^,y,i^O)^Â±d{x,y), (7.5) where d{x, y) is the distance from {x, y) to the initial snake curve 'P(O). As discussed in Sethian [62], with the presence of the diffusion term Fq (curvature dependent term), for e > 0, the parabolic righthand side diffuses sharp gradients and force to stay smooth for all values of t. When there is no diffusion term F^, e = 0, the
PAGE 112
100 curve moves with unit speed, and a corner may develop from smooth initial data. Once a corner develops, it is not clear how to propagate the curve in the normal direction, since the derivative is not defined there. A variety of "weak" solutions have been proposed to evolve the curve beyond the formation of singularities. Of all such weak solutions, one is interested in the one that is the limit of smooth solutions as e > 0. Another way to characterize this weak solution is via the following "entropy condition" [52]: if the curve is viewed as the boundary of a propagating flame, then once a particle is burnt, it stays burnt. Hence, approximations to the spatial derivatives are sought that do not artificially smooth sharp corners but pick out the correct entropy solution beyond the occurrence of singularities. Careful adherence to this stipulation produces the Huygens principle construction. In fact, the schemes given by Osher and Sethian [52, 62] are motivated by the fact that the entropy condition for evolving curves is identical to the one for hyperbolic conservation laws, where stable, consistent, entropysatisfying algorithms have a rich history [22, 33, 40, 41]. Let's consider the curve evolution in one dimensional space with constant normal velocity F = = 1, e = 0. Eq. (7.1) becomes a standard HamiltonJacobi equation 0u + i^(0ix) = 0, (7.6) with H{(j)ix) = {(klxY'"^ and a given initial value of (f)^. Let u
PAGE 113
101 with H{u) Â— {u^yi'^. Eq. (7.7) is a scalar hyperbolic conservation law in one space variable. Numerous research results reported in literature [52, 61] have demonstrated that if we use the standard forward difference, central difference, or backward difference mechanisms to discretize this equation, the solutions can develop discontinuous jumps, known as shocks, even with smooth initial data. As in Osher and Sethian [52], we use a conservative, monotone scheme called upvirind finite diff"erence method to discretize Eq.(7.7). The salient point to be noted is that this conservative scheme produces a solution that satisfies the entropy condition. An upwind scheme automatically differences in the outwardflowing direction during the curve evolution, thus naturally accommodates for the evolution behavior. We now introduce the upwind finite difference scheme to solve Eq. (7.6). Writing Eq. (7.6) using a forward difference in time yields: lt' = K^tH{ul (7.8) where (l)^^ denotes the value of 0i at grid location i Ax and time n At, with Ax and At being space and time intervals in the finite difference. In the upwind finite difference, we use an appropriate numerical flux function G to approximate H rit' = AiG(u,_i/2,n,+i/2) (7.9) = 4>l^tG{Dcl>,,,Dtci>,,),
PAGE 114
102 where and are standard forward and backward difference operators: C;^, = Â«^t2 (7.10) Dth. = (7.11) We can use the HamiltonJacobi fiux function given in Osher and Sethian [52] as the flux function G. When H{u) can be written as a function of u^, that is, H{u) Â— f{u^) for some function /, the flux function is: G(Mi_i/2,u,+i/2) = G{D(j)u,D+(f)u) (7.12) = f{{max{Dcl>u, 0)f + {min{D+u, O))^). This conservative monotone scheme is an upwind method in that it differences in the direction of propagating characteristics. In the case when Fa = 1, /(u^) = (w^)^/^, Eq. (7.9) reduces to = 1 At{ {max{D^u, 0))' + {min{Dthi, 0))' Y'\ (7.13) This algorithm produces the correct entropysatisfying weak solution to the curve evolution problem. ' 7.2 Our Numerical Approach The above discussion can be extended to problems in higher dimensions [45, 52]. Recall that the original intent was to solve (7.3). In two dimensions, the scheme given
PAGE 115
103 in Eq. (7.13) is extended by differencing in each direction to produce the following scheme for Eq. (7.3) x{(maa;(D0i,(ij),O))2 + (mm(D+u = \J a{x, y)(f>l^ + b{x, y)(j)'iy + c{x, y)(f)ix(t>iy, (7.15) where a{x,y),b{x,y) and c{x,y) are determined by the entries of 3b and do not change over time. and (f^jy can still be approximated by the upwind finite difference scheme discussed earlier. But for the iy term, we use the minmod finite difference approximation discussed by Kimmel et al. [34]. The minmod
PAGE 116
104 finite derivative is defined as: sign{a)min{\a\, \b\) if a6 > 0 minmod{a,b} = I (716) 0 otherwise. We use this definition to approximate (j)ix4>iy by \x(t>\y\x=iAx,y=i^y = minmod{D^ (pi^^ij), D' (j)^^^,^^) minmod{D^(j)^^(i^j),Dy(j)x,{i,j)), (7.17) where D^,D~,D^, and D~ are as defined before. Combining these finite differences yields a first order numerical scheme for solving (7.15), the advection term in Eq. (7.1): [aij((maa;(LÂ»(^i,(,j),0))2 + (mzn(L)+0i,(,j), 0))^) +6ij((maa;(D;0i,(i,,), 0))^ + (mm(L>+
PAGE 117
CHAPTER 8 SHAPE RECOVERY EXPERIMENTS WITH SNAKE PEDALS This chapter presents 2D and 3D shape recovery examples using the snake pedal model developed in this thesis. Examples include shape recovery from a variety of 2D/3D data sets, including range data, MRI and CT image data. The examples shown depict a wide range of shapes with complex geometries and topologies. These experiments serve to demonstrate the power and versatility of the snake pedal model for shape recovery and shape synthesis. 8.1 Shape Recovery with Snake Pedals in 2D In this section, we present a set of three shape recovery experiments in 2D. All the experiments are performed on MRI brain scans. In all the experiments, the spatial resolution of the images is 256 x 256 and a region of interest is identified interactively. In the first experiment, we show shape estimation of a caudate nucleus, a cortical structure in the human brain, from a slice of a brain MRI. The snake pedal initialization is shown in Fig. 8.1 (a), while (b) and (c) show the intermediate stage and final model fit respectively. This fitting is achieved at interactive rates on an SGI02. Fig. 8.2 depicts a similar experiment carried out on a cerebellum structure in the human brain, from a slice of brain MRI. Fig. 8.2 (a) (c) show the snake i 105
PAGE 118
106 (a) (b) (c) Figure 8.1. Fitting result (c) of a snake pedal to a caudate nucleus with initialization shown in (a) and intermediate stage shown in (b). pedal initialization, intermediate stage, and final fits respectively. Fig. 8.3 (a) (c) depict a 2D fitting example for a hippocampus in the human brain. Fig. 8.3 is the initialization, (b) and (c) are the intermediate stage and the final fit respectively. In the caudate nucleus and the cerebellum examples, we use both image gradient potential as well as springbased potential from the data points as the external force to guide the model to conform to the desired image features. In the hippocampus example, we only use the springbased data point force, since the image boundary is not well defined. The points (shown in red) are placed by an expert neuroscientist. Note that the model initialization is quite far from the final true shape in each example. Also, the model fits are quite accurate qualitatively (for visual inspection). The accurate extraction of the structures in the human brain plays an important role for the computer aided medical diagnosis. For example, when a person is afflicted with epilepsy, the volume of the hippocampus in the person's brain changes significantly from the volume of the hippocampus in the normal brain [21, 39].
PAGE 119
107 (a) (b) (c) Figure 8.3. Fitting result (c) of a snake pedal (in green) to a hippocampus with initialization shown in (a) and intermediate stage shown in (b). The data points (in red) are placed by the expert along the hippocampus boundary.
PAGE 120
108 8.2 Shape Recovery with Snake Pedals in 3D In this section, we present a set of seven experiments, three of these experiments demonstrate the shape reconstruction with the snake pedal model for medical image analysis applications. The third and fourth examples are designed to illustrate the global deformations of the model without explicitly including the bending, tapering, and twisting parameters in the model representation. In these two examples, the model is fitted to 3D synthetic data from a bent shape (hair pin) and a twisted shape, respectively. The sixth and seventh examples demonstrate model fitting to 3D sparse range data points. The medical examples are organized as follows: in Fig. 8.4, left to right, the images depict a slice of an MR brain scan in which the shape of interest Â— hippocampus, a gyrus, and a cerebellum Â— has been identified by a neuroscientist via sparsely placed points (only in selected slices of the MR scan) on the shape boundary. The next image shows the collection of these 3D points (in red) and the initialized snake pedal model (in green) followed by images depicting an intermediate stage of fitting and the final fitted model respectively. As evident, the model achieves a visually accurate fit to the data. In the case of a hippocampus, model fitting was performed for 30 left and right brains and computed volumes were validated against manual measurements. The Table 8.1 depicts the ratio of the left to right hippocampal volumes. The degree of match between the manual and computed volume ratios is very encouraging and sets the stage for further research in automatic methods for epilepsy studies.
PAGE 121
109 % (a) (b) (c) (d) (e) (f) (g) (h) Â• (i) 0) (k) (1) Figure 8.4. Fitting results of snake pedal surfaces to a hippocampus (d), a gyrus (h), and a cerebellum (1) with initializations shown in (b), (f), and (j), the intermediate stages are shown in (c), (g), and (k), respectively, (a), (e) and (i) are slices of MRI scans for the hippocampus, gyrus, and cerebellum, respectively.
PAGE 122
110 Table 8.1: Comparison between the manual and computed ratios of the left to right hippocampal volumes Manual Automatic Index DilllU Volume L/R Volume L/R 4897 R 19 3.3345 3.4142 1 4897 L lb 3.3310 0.9895 3.6625 0.9322 4915 R 20 3.8114 3.7209 2 4915 L 19 3.7668 0.9882 3.6412 0.9785 4957 R 19 3.8114 3.7209 3 4957 L 20 3.7668 0.9882 3.6412 0.9785 4816 R 11 3.5521 3.5430 4 4815 L OA 20 3.4093 0.9598 3.4598 0.9765 1020 R 20 3.9881 3.7116 5 1020 L 18 4.1070 0.9710 3.9344 1.0600 4910 R 18 4.2202 3.9125 D A f\i n T 4910 L 16 4.1125 1.0262 3.9549 1.0108 A ^Ci A T> 4784 R 19 4.0915 4.0754 7 4784 L 22 3.8446 0.9494 3.8525 0.9453 4834 R 19 4.4162 3.8423 8 4834 L 18 4.3163 0.9773 4.0691 1.0590 4oUl sx 1 n 19 4.8368 4.3127 9 4501 L 19 3.8990 0.8061 4.0135 0.9306 4078 R 20 3.7788 3.6167 10 4078 L 20 3.3580 0.8886 3.2051 0.8862
PAGE 123
] 111 11 4216 R 18 3.6067 0.9208 3.2892 0.9655 4216 L 17 3.3206 3.1758 12 1103 R 19 4.4856 0.7932 4.5758 0.8103 1103 L 18 3.5579 3.7078 13 1214 R 21 4.3109 1.0069 4.1128 1.0786 1214 L 22 4.3407 4.4360 14 4944 R 18 3.9551 0.9934 3.6416 1.0264 4944 L 18 3.9814 3.7379 15 4940 R 18 3.3589 0.9349 3.0344 0.9732 4940 L 28 3.1401 2.9530 16 4921 R 17 3.6023 0.8939 3.5350 0.9356 4921 L 17 3.2201 3.3074 17 4848 R 22 4.8625 0.8861 4.2088 0.9325 4848 L 19 4.3088 3.9247 18 4455 R 22 4.4133 1.0001 4.2519 1.0265 4455 L 20 4.4178 4.3647 19 4764 R 19 4.2601 0.6153 3.5871 0.6239 4764 L 18 2.6212 2.2381 20 4695 R 18 3.4442 1.0144 3.0730 1.1071 4695 L 18 3.4939 3.4021 < A
PAGE 124
"J 112 21 4689 R 16 2.9520 0.8593 3.0296 0.8868 4689 L 16 2.5368 2.6868 22 4540 R 18 3.3108 0.9939 3.4186 0.8935 4540 L 19 3.2907 3.0548 23 4001 R 21 5.3441 1.0483 4.4065 1.1726 4001 L 20 5.6024 5.1672 24 4038 R 19 3.9391 0.8899 3.6803 0.9650 4038 L 17 3.5053 3.5518 25 4142 R 19 3.0900 0.9885 3.3657 0.8837 4142 L 18 3.0546 2.9742 26 4187 R 18 3.9353 0.9699 3.8577 0.9505 4187 L 17 3.8170 3.6669 27 4415 R 18 3.6514 0.9896 3.3460 1.097 4415 L 19 3.6896 3.6725 The two synthetic data examples are organized as follows: in Fig. 8.5, left to right, the images depict a 3D synthetic data set (in red) along with the model initialization (in green) superimposed, which is followed by images of intermediate stage of fitting and the final model fit. Note that in these synthetic data examples, in one case, the data contains shapes exhibiting large bending and in another a large amount of twist. What we would like to stress here, is the fact that the model did not require any explicit bending or twisting transformations (unlike the deformable t 4 Â•.
PAGE 125
113 (d) (e) (f) Figure 8.5. Fitting results of snake pedal surfaces to a bend (c) and twist (f) shape. Initializations are shown in (a) and (d) followed by the intermediate stages (b) and (e) respectively. superquadric) to achieve the fits. As mentioned earlier, modeling schemes that require explicit incorporation of bending, twisting transforms are highly nonlinear and are known to introduce numerical instabilities into the fitting algorithms. In contrast, our modeling scheme does not require explicit application of such bending and twisting transformations to achieve the desired fit. The sixth and seventh examples depict the snake pedal fitting in 3D to a "vulcau" (Spock form Star Wars) head and a mechanical part (anvil) from range data. The images in Fig. 8.6 are arranged in the same way as the examples in Fig. 8.2. In all the examples, we use a 20 x 40 mesh to discretize the model except in the bent shape example (Fig. 8.5), in which we use a 40 x 80 mesh.
PAGE 126
114 Figure 8.6. Fitting results of snake pedal model to spock(c) and an anvil (f). Initializations are shown in (a) and (d) followed by the intermediate stages of fitting in (b) and (e) respectively.
PAGE 127
115 (a) (b) (c) (d) Figure 8.7. Hybrid geometric active contour model fit examples: (a) is model initializations (snake pedal in green and generator in red), (b) and (c) are intermediate stages of evolution and final fits are shown in (d). 8.3 Fitting Results with the Hybrid Geometric Active Model In this section, we present the levelset implementation of the hybrid geometric active model with an ellipse generator described in Chapter 6. In Fig. 8.7 and Fig. 8.8, left to right, images show model initialization Â— depicting the snake pedal in green and the global shape in red, intermediate stages, and final fits, respectively. The human head image is obtained from a slice of an MRI scan, and the spanner image is obtained from the "Computer Vision Test Images" web site [1]. As is evident, the snake pedal was successful in extracting the shape of interest in these images via the model fitting process wherein an imagebased speed function was used to deter the model evolution. Note that the "global" shape depicted in these examples are scaled versions of the generator. The "global" shape provides us with a "global size" and "global orientation" of the shape of interest. More medical data fitting examples are given in Fig. 8.9. In each row of the figure, from left to right, images show the model initializations, intermediate stages and the final model fits. The first and second examples contain the 2D slices from
PAGE 128
(a) (b) (c) (d) Figure 8.8. Hybrid geometric active contour model fit examples: (a) is model initializations (snake pedal in green and generator in red), (b) and (c) are intermediate stages of evolution and final fits are shown in (d). an abdominal CT scan. In the first row, the model is initialized in the stomach and in the second row, the model is initialized inside the liver metastasis. The third and fourth example contain 2D slices from an MR scan of the human heart. The model is initialized to capture the endocardium and the myocardium structure respectively. The snake pedal is shown in green and the global shape in red. We use imagebased speed function to deter the model evolution in all the four examples. Figure 8.10 depicts four topological change examples. In each row, from left to right, images depict model initialization, intermediate stages of evolution and final fit respectively. In the first example, the concept of a "global" shape may be associated with a "global" orientation of a cluster of shapes (in this case, two). In the second and fourth example, the snake pedals are initialized to include all the objects in the images, then the models shrink and split to fit to all the object boundaries. While in the third example, the snake pedal is initialized as a single small ellipse, as the fitting proceeds, the model expands and splits, and finally fits to all the object contours in
PAGE 129
117 Figure 8.9. Hybrid geometric active contour model fit examples: Left to right, model initialization (snake pedal in green and generator in red), intermediate stages of evolution and final fit.
PAGE 130
118 the whole image. The global shapes of the generator are not shown here since the meaning of the "global" shape in these examples is not very useful. In all these examples, we implemented the levelset form of the equation: = g{x,y) F{k,) JbV
PAGE 131
(a) (b) (c) (d) Â•IÂ» '{SB* D^D D^D D^D D^D (e) (f) (g) (h) (i) (j) (k) (1) y>^ yÂ»> y.> (m) (n) (o) (P) Figure 8.10. Topological change examples with hybrid geometric active contour models: Left to right, model initialization (snake pedal in yellow or green and generator in red), intermediate stages of evolution and final fit.
PAGE 132
CHAPTER 9 CONCLUSIONS AND FUTURE DIRECTIONS In this thesis, we have presented a powerful geometric shape modeling scheme which allows for the representation of global shapes with local detail and permits model shaping as well as topological changes via physicsbased control. Our modeling scheme consists of representing shapes by pedal curves and surfaces. Physicsbased control was introduced for shaping these geometric models by using a dynamic spline to represent the position of the pedal point in the regular pedal curve/surface. The model, dubbed as a "snake pedal," allows for interactive manipulation by applying forces to the snake. A levelset implementation of the geometric active contour facilitates topological changes in the snake pedal model. The hybrid geometric active contour model in essence introduces a vehicle for incorporating a global parameterized shape descriptor into the now popular framework of geometric active contours / surfaces . 9.1 Salient Features of the Model The snake pedal models are bestowed with several unique features which make them very powerful and attractive for use in shape synthesis and recovery applications. A list of these features are: 120
PAGE 133
121 9.1.1 Geometric Model with Physicsbased Control The deformable superquadric model is created by superimposing a dynamic spline on a geometric primitive, namely, a conventional superellipsoid. The snake pedal model is created by applying the pedaling operation to a geometric primitive with respect to a dynamic spline. Both models allow for the representation of global and local shape characteristics of an object, but the snake pedal model is inherently a geometric model, while the deformable superquadric model is a physicsbased model. 9.1.2 Compact Representation Unlike the deformable superquadric model, the snake pedal model does not need to explicitly include the bending, tapering and twisting parameters in its representation to realize the global deformation for the same. Since there are no additional nonlinearities introduced into the model, numerical efficiency and stability can be largely improved in the model fitting to the image data with the snake pedal model. 9.1.3 Automatic Handling of Topological Changes The snake pedal model, embedded in a levelset framework, naturally handles the topological changes while recovering the shapes from image data. We introduce the novel concept of a global model into the now popular PDEbased curve/surface evolution framework, which is very useful in the shape recognition and image indexing tasks.
PAGE 134
122 9.2 Future Directions Our experience with the snake pedal model has revealed several issues that are relevant to the continued development and success of the model. This section summarizes these issues and indicates some potential further research directions. 9.2.1 Invariant Representation A rigid motion invariant representation of the global parameters as well as of the coarse level coefficients of local parameters is required to collect the statistics using a Bayesian framework. If we solve the global parameters in accordance with a certain predefined order while fitting the model, the invariance requirement for the global parameters can be satisfied. Invariant representation of the local representation is currently under investigation and the preliminary results are encouraging. 9.2.2 Introduce Dvnamics in the Representation Dynamics is necessary in the motion tracking and motion synthesis tasks. Our physicsbased control framework can be extended to shape and nonrigid motion estimators by incorporating the time factor into the governing equation. The equation of motion will become a secondorder PDE and efficient numerical techniques need to be developed to solve the problem. 9.2.3 Extend the hvbrid geometric active model to 3D Currently, we only implemented the 2D hybrid geometric active contours. We need to extend the snake pedal model to cope with topological changes in shape in 3D
PAGE 135
123 by introducing the 3D hybrid geometric active surfaces. We also need to experiment with more image data in 3D.
PAGE 136
APPENDIX A DERIVATION OF STIFFNESS MATRIX K In Chapter 2, we obtained the internal deformation energy in an active/deformable model. To obtain the expression for the stiffness matrix K, we employ the finite element method (FEM) to discretize the model in material coordinates u = {u,v). As shown in Fig. 5.2, we use bilinear quadrilateral elements in conjunction with linear triangular elements. The later i sused for representing the regions around the north and south poles. We will first discuss the bilinear quadrilateral elements, and the linear triangular elements can be handled in the similar way. A.l Bilinear Quadrilateral Elements Bilinear quadrilateral elements are quite common in FEM literature [87]. In this section, we present the basic transformation equations between material and element coordinates. Denote the finite element nodal shape functions by , where i = 1, 2, no, and no is the number of nodes associated with each element Ej. For the bilinear quadrilateral elements, we have Ni{u,^) = ^{l+cocj,){l+4^iPi), (A.l) where {uji,ipi) are the reference coordinates of node i as shown in Fig. (A.l). 124
PAGE 137
125 (4) (1,1) Â®(l,l) 0) (3) (1,1) r2) (1,1) (a) Figure A.l. Bilinear quadrilateral element, where {u, v) are material coordinates and {lo, ip) are the reference coordinates. The reference coordinates (a;, ^p) are related to the material coordinates u Â— {u, v) by 2 2 u = {uUc) iJ; = {vVc), (A.2) where [uc, Vc) are the coordinates of the element center, a and b are the length and width of the element. In our implementation, we use a = b. Differentiating Eq. (A.l) with respect to u and v gives (A.3) dNi dNi dto dNi dip 1 ^ d^ = a7 + di^d^ = 2^^^ + (A.4)
PAGE 138
. ^ 126 When integrating function f{u, v) over element Ej, we perform the integration in the reference system: j f{u,v)dudv = j J f{uj,ip)^dudip (A.5) A.2 Derivation of the Local Stiffness Matrix As described in section 5.2.4, we use the membrane plus thin plate energy expressed in Eq.(5.17) as the internal deformation energy imposed on the snake, By using the procedure described by Terzopoulos and Metaxas [78], we can obtain the local stiffness matrix. For the membrane energy part, if we denote the stiffness matrix for the element Ej by G sjjnoxno^ ^j^g entries of can be represented as {kDim = J + ^^) du dv, (A.7) where l,m = 1,2, ...,no. Since we use a quadrilateral element, no = 4 in this case. Using the formula in Eq. (A.3), (A. 4) and (A.5), we obtain the following local stiffness
PAGE 139
127 matrix for the membrane energy D a 2 1 1 2 1 1 1 1 1 1 2 + () 1 1 2 1 1 2 (A.8) Adopting the similar scheme, we obtain the following stiffness matrix Kj for the thin plate energy part in element Ej * ab 11 11 1 1 11 11 1 11 1 (A.9) The summation of and gives the full stiffness matrix for element Ej K^ = K4 + K^ (A.IO) Note that both and K( are symmetric positive definite (SPD) matrixes, therefore their summation, K', is also a SPD matrix. A. 3 Kronecker Tensor Product Representation of the Assembled Stiffness Matrix The matrix discussed earlier is only the local stiffness matrix for element Ej, which describes the stressstrain relationship between all the nodes within element
PAGE 140
128 u 1 2 3 4 El 5 E2 6 E3 7 1 E4 9 Es 10 E6 11 1 V (a) Figure A. 2. A small 3x4 mesh in {u,v) space. Ej. For example, ^3 = 4 means applying one unit displayment to node 3 results in 4 unit force on node 2. We need to assemble the local matrix over all the elements to obtain the entire stiffness matrix K. To achieve this, we introduce a novel method via the Kronecker tensor product operation. A small 3x4 mesh in {u, v) space as shown in Fig. A. 2 is used to illustrate the idea. We will only consider the first part of the (the first term in the summation of Eq. (A. 8)) for explanation, the second part of K4 and can be handled in the same way. In Fig. (A. 2), we assign a number to each node, ranging from 1 to 12, and denote each element by Ei, E2, E^. We choose node 6 for our illustration, and investigate the stressstrain realtionship between node 6 and every other node in the entire mesh. From Eq. (A. 8), observe that applying one unit displacement to node 6 will result in ^(^) x 1 unit, (^f ) J x (1) unit, ^( J) x (2) unit force on node 2, node 1, and node 5, respectively. This concept can be visualized by the nodal molecule
PAGE 141
129 WI b , 6a ' 00 2) (+21 (a) wi b . 6a Â®0 (b) 2 ) ( +2) wib I 1^ wib X X (c) 6a I I 6a Â©Â® (d) (a) Figure A. 3. Membrane molecules representation, as shown in Fig. (A. 3) (a). We can use a similar representation for elements E2, E4 and Â£'5, resulting in the nodal molecules as shown in Figs. (A. 3) (b), (c), and (d), respectively. Note that there is no interaction between node 6 and nodes in elements E3 and Eq. From the nodal molecular representation, when applying one unit displacement to node 6, we can obtain the force resulting on every other node in the entire 3x4 mesh by summing up all the molecular values at the same molecular location. The result is shown in Fig. (A. 4). The assembled molecular representation can be expressed as a Kronecker product fomula wi b 1 2 1 1 4 1 (A.ll)
PAGE 142
130 W b 6a 5)@Â© (a) Figure A. 4. Assembled membrane molecular representation. The assembled stiffness matrix that corresponds to the first part of the membrane energy (first term in the summation of Eq. (A. 8)) is Kr"' = (y^)A.Â®B where (A.12) 2 1 0 1 4 1 0 1 2 and BÂ„ = 2 1 0 1 1 2 1 0 0 1 1 0 2 1 (A.13) We now prove the above conclusion by expanding the Kronecker tensor product expression in (A.12). According to the Kronecker product definition in Eq. (5.25),
PAGE 143
(parti) 2 1 0 1 4 1 0 1 2 2By By 2 1 0 1 By iBy By 0 By '^By 1 2 1 0 0 1 1 0 2 1 1 2 2 42 0 0 2 1 0 1 1 0 01 21 41 01 1 01 84 0 21 01 1 1 0 4 8 4 0 1 2 1 0 01 21 04 8 0 1 04 81 0 21 01 42 0 1 1 0 2 0 1 1 02 4 2 1 01 22 02 4 (A.14)
PAGE 144
132 In this matrix, each entry represents the stressstrain relation between node / and node m, where l,m = 1,2, ...,12. For example, the entries in the 6th column (A;^"'""^)(6 represent the force caused in node / by the unit displacement in node 6. We can see from column 6 that one unit displacement in node 6 contributes (^) X 8 unit force to itself, and (^J) x (4) unit force to node 5 and 7, (^^) x 2 unit force to node 2 and 10, (^^) x (1) unit force to node 1, 3, 5 and 7, and zero force to node 4, 8 and 12. This is exactly what Fig. (A. 4) conveys. For the second part of the membrane energy, performing a similar procedure, we have Kr"^==(^^).B.Â®AÂ„ (A.15) where 1 1 0 1 2 1 0 1 1 and By = 4 10 1 14 10 0 14 1 10 14 (A.16) Since a = 6 in our implementation, the entire stiffness matrix corresponding to the membrane energy is , Â• (A.17)
PAGE 145
133 Similarly, we have the stiffness matrix corresponding to the thin plate energy K, = (^)A,Â®A,. (A.18) In general, for an m x n mesh, denote 2 1 1 4 1 1 4 1 1 2 (A.19) Ay Â— 4 1 1 4 1 Bx = 1 1 1 1 2 1 1 4 1 1 4 e 3?" 1 2 1 1 1 (A.20) (A.21)
PAGE 146
134 By = 2 1 1 2 1 1 1 1 2 1 1 2 (A.22) then the assembled stiffness matrix corresponding to the energy expressed in Eq.(5.17) is: K = + Kt = wi{A^ (g) By + Â® Ay) + W2{A^ <8> Ay), (A.23) where u;i and W2 are tension and rigidity parameters. Note that the difference between the first column/row and last column/row in and A^ or B^; and B^ is due to the different boundary conditions, namely, the natural boundary and the periodic boundary conditions respectively. When employing the same boundary conditions and let m = n, A = Aj; = Aj,, B = B^ = B^, Eq. (A.23) becomes K = w;i(A (g) B + B O A) + u;2(A (g) A). (A.24) This is the same equation as (5.23).
PAGE 147
APPENDIX B RELATIONSHIP OF ~ i From the differential geometry theory, we have the relationship between the curvature of the snake pedal kg and the corresponding higher dimensional function k. = V IIV02 By using V02 = JbV^i, kg can be expressed as: (B.l) ke = V\JB^(t>l (B.2) We denote each entry of Jb as: Jr = Jbu JBi2 (B.3) where Js^^ = J^^i > Js is the symmetric matrix. Denote V01 ' d^i d(j)i dx dy (B.4) 135
PAGE 148
136 Jb can be written as: JbV0i = Jbh Jbi2 JBi2 Jb22 dx oy T ^ L T 36, (B.5) The magnitude of is: (B.6) Therefore, we can represent A;e in terms of by: JbV0i 1 (iB. Jb. + Jb. JbJ + (JL + JIJ ^11^ ]. 5^ [ WBi: + ^B.j "5^5^ + {'JBu'Jb,2 + Jb.^Jb^^) ^J^y gy + (JbÂ„ Jb. + Jb. JbJ + (JI., + JIJ ]. (B.7)
PAGE 149
REFERENCES [1] http://machine_vision.cse.psu.edu/book/testbed/images/. Index of /book/testbed/images, Penn State University Computer Vision Lab. [2] A. A. Amini, T.E. Weymouth, and R.C. Jain. Using dynamic programing for solving variational problems in vision. IEEE Trans, on Pattern Analysis and Machine Intelligence, 12(9):855867, 1990. [3] E. Bardinet, N. Ayache, and L.D. Cohen. Fitting of isosurface using superquadrics and freeform deformations. In Proc. of IEEE Workshop on Biomedical Image Analysis WBIA'94, pages 184193, Seattle, 1994. [4] A.H. Barr. Global and local deformations of solid primitives. Computer Graphics (ACM), 18(3):2130, 1984. [5] M.O. Berger and R. Mohr. Towards autonomy in active contour models. In Proc. 10th International Conference on Pattern Recognition, pages 847851, Atlantic City, NJ, 1990. IEEE Computer Society Press. [6] A. Blake and A. Zisserman. Visual Reconstruction. MIT Press, Cambridge MA 1987. [7] M.P. Do Carmo. Differential Geometry of Curves and Surfaces. PrenticeHall, Inc., New Jersey, 1976. [8] M.P. Do Carmo. Riemannian Geometry. PrenticeHall, Inc., New Jersey, 1992. [9] V. Caselles, F. Catte, T. Coll, and F. Dibos. A geometric model for active contours in image processing. Numerische Mathematik, 66:131, 1993. [10] V. Caselles, R. Kimmel, and G. Sapiro. Geodesic active contours. In Fifth International Conference on Computer Vision, pages 694699, 1995. [11] I. Cohen and L.D. Cohen. Hyperquadric model for 2d and 3d data fitting. In Proc. oflCPR, pages 403405, Piscataway, NJ, 1994. [12] L.D. Cohen. On active contour models and ballons. CVGIP: Image Understanding, 53(2):211218, 1991. 137
PAGE 150
138 [13] T. Cootes, A. Hill, C. Taylor, and J. Haslam. The use of active shape models for locating structures in medical images. Image and Vision Computing, 12(6):355366, 1994. [14] D. DeCarlo and D. Metaxas. Blended deformable models. IEEE Trans. Pattern Analysis and Machine Intelligence, 18:443448, 1996. [15] C.L. Epstein and M. Gage. The curve shorting flow. In A. Chorin and A. Majda, editors, Wave Motion: Theory, Modeling, and Computation. SpringerVerlag, New York, 1987. [16] C.L. Epstein and M.I. Weinstein. A stable manifold therem for the curve shortening equations. Math. Comp., 34, 1980. [17] M. Gage. An isoperimetric inequality with applications to curve shortening. Duke Math., 50(1225), 1983. [18] M. Gage and R.S. Hamilton. The heat equation shrinking convex plane curves. Journal of Differential Geometry, 23:6996, 1986. [19] I.M. Gelfand and S.V. Fomin. Calculus of Variations. PrenticeHall, Englewood Cliffs, NJ, 1963. [20] S. Geman and D. Geman. Stochastic relaxation, gibbs distribution, and bayesian restoration of images. IEEE Trans, on Pattern Analysis and Machine Intelligence, 6(6):721724, 1984. [21] R. Gilmore, M. Childers, C. Leonard, R. Quisling, S. Roper, S. Eisenschenk, and M. Mahoney. Hippocampal volumetrics differentiates patients with temporal lobe epilepsy and extratemporal lobe epilepsy. Archives of Necrology, 52:819824, 1995. [22] M. Grandall, H. Ishii, and RL. Lions. User's guide to viscosity solutions of secondorder partial differential equations. Bull. Am. Math. Soc, 27:167, 1992. [23] A. Gray. Modern differential geometry of curves and surfaces. CRC Press, Boca Raton, 1993. [24] M. Grayson. The heat equation shrinks embedded plane curves to round points. Journal of Differential Geometry, 26:285314, 1987. [25] R.P. Grzeszczuk and D.N. Levin. Segmentation images with stocastically deformable contours. In R.A. Robb, editor, Proc. Third Conf. on Visualization in Biomedical Computing (VBC'94), volume 2359 of SPIE Proc, pages 7289, Bellingham, WA, 1994. [26] S. Han, D.B. Goldgof, and K.W. Bowyer. Using hyperquadrics for shape recovery from range data. In IEEE Proceedings of the 3rd International Conference on Computer Vision, pages 492496, Berlin, 1993.
PAGE 151
139 [27] M. Hebert, K. Ikeuchi, and H. Delingette. A spherical representation for recognition of freeform surfaces. IEEE Pattern Analysis and Machine Intelligence, 17(7):681690, 1995. [28] A. Hill, A. Thornham, and C.J. Taylor. Modelbased interpretation of 3d medical images. In Proc. 4th British Machine Vision Conf. (BMVC'93), pages 339348, Surrey, UK, September 1993. [29] G.E. Hinton, C.K. Williams, and M.D. Revow. Adaptive elstic models for handprinted character recognition. In J.E. Moody, S.J. Hanson, and R.P. Lippman, editors. Advances in Neural Information Processing Systems. Morgan Kauffman, San Mateo, CA, 1992. [30] W.M. Hsu. Direct manipulation of freeform deformation. Computer Graphics, 26(2):177184, 1992. [31] M. Kass, A. Witkin, and D. Terzopoulos. Snakes : A active contour models. Int. Journal of Computer Vision, 1:321331, 1987. [32] S. Kichenassamy, A. Kumar, P. Olver, A. Tannenbaum, and A. Yezzi. Gradient flows and geometric active contour models. In Fifth International Conference on Computer Vision, pages 810815, 1995. [33] B. B. Kimia, A. Tenenbaum, and S. W. Zucker. Toward a computational theory of shape: An overview. Lecture Notes in Computer Science, 427:402407, 1990. [34] R. Kimmel, A. Amir, and A.M. Bruckstein. Finding shortest paths on surfaces using level sets propogation. IEEE Trans, on Pattern Analysis and Machine Intelligence, 17(6):635640, 1995. [35] S.H. Lai and B.C. Vemuri. An o(n) iterative solution to the poisson equation in lowlevel vision problems. In IEEE Conference on Computer Vision and Pattern Recognition, pages 914, 1994. [36] S.H. Lai and B.C. Vemuri. Physicallybased adaptive preconditioning for early vision. In IEEE Workshop on Physicsbased Modeling in Computer Vision, 1995. [37] S.H. Lai and B.C. Vemuri. Shermanmorrisonwoodbury fornularbased algorithms for the surface smoothing problem. LAA, 1997. [38] J.D. Lawrence. Catalog of Special Plane Curves. Dover Publications, Inc. New York, 1972. [39] C. Leonard, C. Williams, R. Nicholls, O. Agee, K. Voeller, J. Honeyman, and E. Staab. Angelman and praderwilli syndrome. A magnetic resonance imaging study of differences in cerebral structure. American Journal of Medical Genetics, 46:2633, 1993.
PAGE 152
140 [40] R.J. LeVeque. Numerical Methods for Conservation Laws. Birkhauser, Boston, 1992. [41] P.L. Lions. Generalized Solutions of HamiltonJ acobi Equations. Pitman Publishing, Boston, 1982. [42] E.H. Lockwood. A book of Curves. Cambridge University Press, 1963. [43] A. Lu and E.L. Wachspress. Solution of lyapunov equations by alternating direction implicit iteration. Computers Math. Applic, 21(9):4358, 1991. [44] J. Malik and P. Perona. Computational model of texture segmentation. In Proc. IEEE Comput Soc Conf Comput Vision Pattern Recognition, pages 326332, Piscataway, NJ, 1989. [45] R. Malladi, J. A. Sethian, and B. C. Vemuri. A topology independent shape modeling scheme. In SPIE Proc. on Geometric Methods in Computer Vision II, volume 2031, pages 246256,. SPIE, July 1993,. [46] R. Malladi, J. A. Sethian, and B.C. Vemuri. Shape modeling with front propagation : A level set approach. IEEE Trans. Pattern Analysis and Machine Intelligence, 17(2):158175, 1995. [47] P.S. Marc, S. Menet, and G. Medioni. Bsnakes: implementation and application to stereo. In lU Workshop, pages 720726, 1990. [48] T. Mclnerney and D. Terzopoulos. Topologically adaptable snakes. In IEEE ICCV, pages 840845, 1995. [49] F. Morgan. Riemannian Geometry. John and Bartlett Publishers, Boston, 1993. [50] T. O'Donnell, T. Boult, and A. Gupta. Global models with parametric offsets as applied to cardiac. In IEEE Proc. of CVPR, pages 293299, 1996. [51] P. J. Olver, G. Sapiro, and A. Tannenbaum. Affine invariant detection: Edges, active contours and segments. In IEEE Proc. on CVPR, pages 520525, June, 1996. [52] S. Osher and J. A. Sethian. Fronts propagating with curvature dependent speed: Algorithms based on halmitonj acobi formulation. Journal of Computational Physics, 79:1249, 1988. [53] D.V. Ouellette. Schur complements and statistics. Linear Algebra and its Applications, 36:187295, 1981. [54] A. Pentland and J. Williams. Good vibrations: Model dynamics for graphics and animation. Computer Graphics, 23(3):215222, 1989. ^
PAGE 153
141 [55] T. Poggio and T. Toree. Illposed problems and regularization analysis in early vision. In Baumann, editor, Proc. AARPA Image Understanding Workshop, pages 257263, New Orleans, LA, 1984. [56] T. Poggio, V. Toree, and C. Koch. Computer vision and regularization theory. Nature, 317(6035), 1985. [57] C.S. Poon, M. Braun, R. Fahrig, A. Ginige, and A. Dorrell. Segmentation of medical images using an active contour model incorporating regionbased image festures. In R.A. Robb, editor, Proc. Third Conf. on Visualization in Biomedical Computing (VBC'94), volume 2359 of SPIE Proc, pages 9097, Bellingham, WA, 1994. [58] G. Sapiro, R. Kimmel, and V. Caselles. Object detection and measurements in medical images via geodesic deformable contours. In SPIE Proc, volume 2573, pages 366378, Bellingham, WA, 1995. [59] T.W. Sederberg and S.R. Parry. Freeform deformation of solid geometric models. In ACM SIGGRAPH, pages 151 160, Dallas, TX, 1986. [60] J. A. Sethian. Curve shorting makes convex curve circular. Invent. Math, 76:357, 1984. [61] J. A. Sethian. Curvature and the evolution of fronts. Commun. Math. Phys., 101:487499, 1985. [62] J.A. Sethian. Numerical algorithms for propogating interfaces: Halmitonjacobi equations and conservation laws. Journal of Differential Geometry, 31:131161, 1990. . * [63] A. A. Shabana. Dynamics of Multihody Systems. Wiley, New York, 1989. [64] J. Shah. Recovery of shapes by evolution of zerocrossings. Technical report, Dept. of Mathematics, Northeastern University, Boston, MA, 1995. [65] J. Shah. A common framework for curve evolution, segmentation and anisotropic diffusion. In IEEE Conf. on Computer Vision and Pattern Recognition, San Francisco, CA, 1996. [66] J.R. Shewchuk. An introduction to the conjugate gradient method without the agonizing pain. Technical Report Tech. Report CMUCS94125, School of Computer Science, Carnegie Mellon University, 1994. [67] F. Solina and R. Bajcsy. Recovery of parameteric models from range images: The case for superquadrics with global deformation. IEEE Trans, on Pattern Analysis and Machine Intelligence, 12:131146, 1990. [68] S.Osher and C.W. Shu. Highorder essentially nonoscillatory schemes for hamiltonjacobi equations. SIAM J. Numerical Analysis, 28(4):907922, 1991.
PAGE 154
142 [69] L.H. Staib and J.S. Duncan. Boundary finding with parametrically deformable models. IEEE Trans, on Pattern Analysis and Machine Intelligence, 14:10611075, 1992. [70] G. Strang. Introduction to Applied Mathematics. WellesleyCambridge Press, Wellesley,MA, 1986. [71] G. Szekely, A. Kelemen, Ch. Brechbuhler, and G. Gerig. Segmentation of 2d and 3d objects from mri volume data using constrained elastic deformations of flexible fourier surface models. Medical Image Analysis, 1(1), 1996. [72] R. Szeliski and D. Tonnesen. Surface modeling with oriented particle systems. Computer Graphics SIGGRAPH, 26{2):1851M, 1992. . . ; [73] R. Szeliski, D. Tonnesen, and D. Terzopolous. Modeling surfaces of arbitary topology with dynamic particles. In Proceedings IEEE CVPR, pages 8287, 1993. [74] H. Tek and B. B. Kimia. Image segmentation by reactiondiffusion bubbles. In Fifth International Conference on Computer Vision, Boston, MA, 1995. [75] D. Terzopoulos. Regularization of inverse visuak problems involving discontinuities. IEEE Trans, on Pattern Analysis and Machine Intelligence, 8(4):413424, 1986. [76] D. Terzopoulos. On matching deformable models to images. In Topical Meeting on Machine Vision, volume 12, pages 160167, 1987. [77] D. Terzopoulos and K. Fleischer. Deformable models. The Visual Computer, 4(6):306331, 1988. [78] D. Terzopoulos and D. Metaxas. Dynamic 3d models with local and global deformations : Deformable superquadrics. IEEE Pattern Analysis and Machine Intelligence, 13(7):703714, 1991. [79] D. Terzopoulos, A. Witkin, and M. Kass. Constraints on deformable models: Recovering 3d shape and nonrigid motion. Artificial Intelligence, 36(1):91123, 1988. [80] A.N. Tikhonov and V.A. Arsenin. Solutions of Illposed Problmes. Winston and Sons, Washington, DC, 1977. [81] B.C. Vemuri and A. Radisavljevic. Multiresolution stochastic hybrid shape models with fractal priors. ACM Trans, on Graphics, pages 177207, 1994. [82] R. Whitaker. Volumetric deformable models. In R.A. Robb, editor, Proc. Third Conf. on Visualization in Biomedical Computing (VBC'94), volume 2359, Rochester, MN, 1994.
PAGE 155
143 [83] J. Wickert. Nonlinear diffusion in scale space: From the continuous to the discrete setting. In M.O. Berger et ai, editor, Proc. ICAOS'96: Images, Wavelets and PDE's, volume 219, pages 111118, Springer, New York, 1996. [84] P.H. William, T.A. Saul, V.T. William, and B.P. Flannery. Numerical Recipe in C. Cambridge University Press, 2 edition, 1992. [85] D.J. Williams and M. Shah. A fast algorithm for active contours and curvature estimation. CVGIP: Image Understanding, 55:1426, 1992. [86] M. Worring, A.W.M Smeulders, L.H. Staib, and J.S. Duncan. Parametrerized feasible boundaries in gradient vector fields. In A.C.F. Colchester and D.J. Hawkes, editors, Information Processing in Medical Imaging: Proc. 13th Int. Conf. (IPMI'QS), Flagstaff, CA, 1993. [87] O.C. Zienkiewicz and R.L. Taylor. The Finite Element Method. McGrawHill, NewYork, NY, 1989.
PAGE 156
BIOGRAPHICAL SKETCH Yanlin Quo was born in Hubei, P. R. China, in 1971. She received her Bachelor of Engineering degree and Master of Engineering degree in electronics engineering from Tsinghua University in July 1993 and July 1995, respectively. Currently, she is pursuing the Doctor of Philosophy degree in electrical and computer engineering at the University of Florida. Her interests include computer vision, signal and image processing, computer graphics, computer simulation, and medical imaging. 144
PAGE 157
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. Baba C. Vemuri, Chairman Associate Professor of Computer and Information Science and Engineering 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. Thomas E. Bullock Professor of Electrical and Computer Engineering 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. in G. Harris Assistant Professor of Electrical and Computer Engineering 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. Stanley Y.W. Svl Professor of Cornputer and Information Science and Engineering
PAGE 158
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. Christina M. Leonard Professor of Neuroscience 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 degye^wpoctor of Philosophy. August 1998 Winfred M. Phillips Dean, College of Engineering Karen A. Holbrook Dean, Graduate School

