Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UFE1001162/00001
## Material Information- Title:
- EFFICIENT AND ROBUST ALGORITHMS FOR MULTIMODAL IMAGE REGISTRATION
- Copyright Date:
- 2008
## Subjects- Subjects / Keywords:
- Algorithms ( jstor )
Computerized axial tomography ( jstor ) Coordinate transformations ( jstor ) Datasets ( jstor ) Estimators ( jstor ) Image rotation ( jstor ) Images ( jstor ) Images of transformations ( jstor ) Mathematical independent variables ( jstor ) Mathematical robustness ( jstor )
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright the author. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Embargo Date:
- 8/8/2012
- Resource Identifier:
- 81266610 ( OCLC )
## UFDC Membership |

Downloads |

## This item is only available as the following downloads: |

Full Text |

PAGE 1 EFFICIENT AND R OBUST ALGORITHMS F OR MUL TIMOD AL IMA GE REGISTRA TION By JUNDONG LIU A DISSER T A TION PRESENTED TO THE GRADUA TE SCHOOL OF THE UNIVERSITY OF FLORID A IN P AR TIAL FULFILLMENT OF THE REQUIREMENTS F OR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORID A 2002 PAGE 2 Cop yrigh t 2002 b y Jundong Liu PAGE 3 This is dedicated to m y paren ts, brother and sisters. PAGE 4 A CKNO WLEDGMENTS I ha v e man y p eople to thank. First I express m y gratitude to Dr. Baba C. V em uri, for his excellen t academic guidance and unconditional supp ort throughout the y ears. I thank him for the n umerous long discussions and constructiv e commen ts on m y dissertation. It w as a great pleasure for me to conduct this dissertation under his sup ervision. I w ould also lik e to thank Dr. Sahni, Dr. Da vis, Dr. Rangara jan and Dr. Rao for their willingness to serv e on m y committee. I am indebted to Dr. Jose Marro quin for pro viding v aluable ideas and commen ts during the course of this w ork. My appreciation and sp ecial thanks also go to Dr. F rank Bo v a and Dr. Lionel Bouc het for pro viding the medical image data for giving useful commen ts regarding to the soft w are system presen ted in this dissertation. Needless to sa y , I am grateful for the supp ort of m y colleagues and friends at the Computer and Information Science and Engineering Departmen t at the Univ ersit y of Florida: Jun Y e, Zhizhou W ang, Eric Sp ellman, Neeti V ohra and F ei W ang. I w ould also lik e to tak e this opp ortunit y to express m y sp ecial thanks to m y paren ts, brother and sisters, who nev er hesitate to deliv er their lo v e, encouragemen t and supp ort to me. This researc h w as supp orted in part b y the gran ts NSF I IS9811042, NIH R O1-RR13197 and NIH R01-NS42075. iii PAGE 5 T ABLE OF CONTENTS page A CKNO WLEDGMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii ABSTRA CT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi CHAPTER1 INTR ODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Problem Description . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Image Registration . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1 Denition . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.2 T ransformations . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Thesis Con tribution . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 RELA TED W ORK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1 F eature-Based Metho ds . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.1 P oin t-Based F eatures . . . . . . . . . . . . . . . . . . . . . 10 2.1.2 Con tour and Surface Based F eatures . . . . . . . . . . . . . 10 2.1.3 Edge/Ridge Based Metho ds . . . . . . . . . . . . . . . . . 12 2.1.4 Reductiv e Metho ds . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Direct Metho ds . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.1 Image Similarit y Measures . . . . . . . . . . . . . . . . . . 14 2.2.2 Dra wbac ks of MI-based Metho ds and the Related Impro v emen ts . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2.3 Computational Eciency . . . . . . . . . . . . . . . . . . . 17 2.2.4 Robustness . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2.5 In terp olation Artifacts . . . . . . . . . . . . . . . . . . . . 18 2.2.6 Extensions to Cop e with Nonrigid Motion . . . . . . . . . . 18 2.3 In tro duction of Our Metho d . . . . . . . . . . . . . . . . . . . . . 19 3 INV ARIANT IMA GE REPRESENT A TION . . . . . . . . . . . . . . . . 22 3.1 Lo cal F requency Represen tation . . . . . . . . . . . . . . . . . . . 23 3.1.1 Gab or Filters . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.1.2 Computing Lo cal-frequency Represen tations . . . . . . . . 26 4 COMPUT A TION EFFICIENCY ISSUE . . . . . . . . . . . . . . . . . . 29 4.1 Matc hing Lo cal F requency Represen tations with Eciency . . . . 30 iv PAGE 6 4.1.1 Spline Based Represen tation of the Registration T ransformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.1.2 3D Rigid T ransformation . . . . . . . . . . . . . . . . . . . 31 4.2 Numerical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.2.1 3D Ane Motion . . . . . . . . . . . . . . . . . . . . . . . 35 5 R OBUSTNESS ISSUE . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 5.1 Robust Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . 37 5.2 Related W ork . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5.3 Robustly Matc hing Lo cal F requency Represen tations . . . . . . . 39 5.3.1 V ector-v alued Lo cal F requency with Ane Motion . . . . . 40 5.3.2 Lo cal F requency Magnitude with Rigid Motion . . . . . . . 42 5.3.3 Numerical Implemen tation . . . . . . . . . . . . . . . . . . 43 6 NONRIGID MOTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 6.1 Numerical Implemen tation . . . . . . . . . . . . . . . . . . . . . . 47 6.2 Adaptiv e Time Step t . . . . . . . . . . . . . . . . . . . . . . . . 48 6.3 Incremen tal Nonrigid Up dating . . . . . . . . . . . . . . . . . . . 49 7 EXPERIMENT RESUL TS . . . . . . . . . . . . . . . . . . . . . . . . . . 50 7.1 Exp erimen ts with ESD-based Measure . . . . . . . . . . . . . . . 50 7.1.1 Syn thesized Motion: Estimation Results . . . . . . . . . . . 50 7.1.2 Real 3D Rigid Motion: Estimation Results . . . . . . . . . 52 7.1.3 Nonrigid Motion: Syn thetic and Real Data Example . . . . 53 7.2 Exp erimen ts with the L 2 E Measure . . . . . . . . . . . . . . . . . 56 7.2.1 Robust Prop ert y of L 2 E Measure . . . . . . . . . . . . . . . 56 7.2.2 Comparison of Magnitude and V ector-v alued Lo cal F requency Sc hemes . . . . . . . . . . . . . . . . . . . . . . . 59 7.2.3 Comparison with NMI using Syn thetic MRI Data . . . . . 60 7.2.4 Comparison with 3D CT-MR Data . . . . . . . . . . . . . . 62 8 CONCLUSIONS AND FUTURE W ORK . . . . . . . . . . . . . . . . . . 69 8.1 Ma jor Con tributions . . . . . . . . . . . . . . . . . . . . . . . . . 69 8.2 F uture W ork . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 BIOGRAPHICAL SKETCH . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 v PAGE 7 Abstract of Dissertation Presen ted to the Graduate Sc ho ol of the Univ ersit y of Florida in P artial F ulllmen t of the Requiremen ts for the Degree of Do ctor of Philosoph y EFFICIENT AND R OBUST ALGORITHMS F OR MUL TIMOD AL IMA GE REGISTRA TION By Jundong Liu August 2002 Chair: Baba C. V em uri Ma jor Departmen t: Computer and Information Science and Engineering The goal of image registration is the alignmen t of t w o or more images of the same scene or ob ject. It is one of the most widely encoun tered problems in a v ariet y of elds (including medical image analysis, remote sensing, satellite imaging, optical imaging, etc.). Numerous algorithms for registering m ultimo dal image data ha v e b een rep orted in these areas. Among them, m utual information (MI) based metho ds are widely accepted as the most eectiv e w a y of handling rigid m ultimo dal image registration problems. But MI metho ds ha v e their dra wbac ks. They are computationally in tensiv e, they lac k robustness, and they do not handle non-rigid motion w ell. The suite of algorithms presen ted in this dissertation is an attempt to pro vide an alternativ e to the MI based algorithms. The primary con tributions of this dissertation are (1) a v ector-v alued lo cal frequency image represen tation that is in v arian t to brigh tness and con trast v ariations in an image; (2) a cost function | for computing the lo cal frequency represen tations of misaligned images | based vi PAGE 8 on the exp ectation of the squared dierence and a signican t mo dication of the Newton metho d to n umerically optimize the cost function in a computationally ecien t manner; and (3) a statistically robust cost function whose optimization yields an estimate of the registration parameters required to align the giv en image pair. The k ey adv an tage of this tec hnique is its abilit y to cop e with alignmen t of image pairs that ha v e a signican t nono v erlap. Exp erimen ts w ere carried out using real clinical CT and MR brain scans to ev aluate the p erformance of our algorithms. They clearly demonstrated the higher time-eciency of our ESD metho d and the sup erior robustness of our L 2 E metho d o v er the m utual information based sc hemes. vii PAGE 9 CHAPTER 1 INTR ODUCTION Medical image analysis tec hnology , with image segmen tation, image matc hing /registration, motion trac king and the measuremen t of anatomical and ph ysiological parameters as the main researc h areas, has seen a tremendous amoun t of gro wth o v er the past decade. The w ork describ ed in this thesis is concerned with the problem of automatically aligning 3D medical images. In particular it deals with the task of aligning image acquired using dieren t imaging devices, con taining b oth corresp onding and complemen tary structure, commonly referred to as the problem of m ulti-mo dalit y image registration. 1.1 Problem Description Image registration is one of the most widely encoun tered problems in a v ariet y of elds including but not limited to medical image analysis, remote sensing, satellite imaging, optical imaging, etc. A p ossible denition of this problem is: Determine the co ordinate transformation, or mapping, relating dieren t views of the same or similar ob jects. These views ma y arise from: The same ob ject \imaged" with dieren t sensors. The same ob ject \imaged" rep eatedly with the same sensor. Multiple similar ob jects all \imaged" with the same sensor. A single ob ject \imaged" with a sensor, and a mo del (matc hing) Medical imaging pla ys a v ery imp ortan t role in mo dern clinical practice. Often, a single mo dalit y alone can not pro vide adequate information ab out a patien t's condition, so they are imaged b y a second and sometimes more mo dalities. Dieren t mo dalities, as sho wn in Figure 1.1 (from the homepage of the 1 PAGE 10 2 Figure 1.1: Image from dieren t mo dalities, from left to righ t: CT, MR, PET, SPECTImage Science Institute, Departmen t of medical Imaging. Av ailable online at h ttp://www.isi.uu.nl/Researc h/Registration /registration-frame.h tml), can usually pro vide dieren t, complemen tary or partially o v erlapping asp ects of the anatom y under examination. So if the information from them can b e com bined, it will greatly facilitate the clinicians to p erform surgical planning, diagnosis and treatmen t. Among the widely used mo dalities, X-ra y , CT (computed tomograph y), MRI (magnetic resonance imaging) and US (ultrasound) depict anatomical information, while PET (p ositron emission tomograph y), SPECT (single-photon emission computed tomograph y), fMRI (functional MRI) pro vide information on the metab olism of the underlying anatom y . Figure 1.1 sho ws images from 4 dieren t mo dalities: CT, MR, PET, SPECT, illustrating the dieren t t yp es of information they con tain. As w e can see from the gure, CT is sensitiv e to b one densit y , MR to tissue densit y whereas PET and SPECT depict the ph ysiology . Figure 1.2 (from the homepage of the Image Science Institute Departmen t of medical Imaging. Av ailable online at h ttp://www.isi.uu.nl/Researc h/Registration/ registration-frame.h tml) also sho ws ho w m ulti-mo dal data can b e used to pro vide a b etter understanding of ph ysiology of the h uman brain aided b y the presence of precise anatomical information. It sho ws the righ t and left hemispheres of a brain, segmen ted from an MR image. Information from a SPECT scan is o v erla y ed on the PAGE 11 3 Figure 1.2: An example of a m ultimo dal visualization cortex. Color enco ding is used to indicate the amoun t of cerebral blo o d p erfusion: from ligh t gra y for lo w p erfusion, via y ello w, to red for high p erfusion. As w e can see, the picture sho ws an area with increase blo o d p erfusion in the righ t hemisphere and there is indicativ e of a pathology on the righ t hemisphere. Ho w ev er, the m ulti mo dalit y images are usually acquired with dieren t devices, at dieren t times, so inevitably , there will b e some motion b et w een them. This mak es accurate geometrical registration a prerequisite step for the eectiv e fusion of the information from m ulti mo dalit y images. 1.2 Image Registration 1.2.1 Denition Let I 1 ( x; y ; z ) and I 2 ( x; y ; z ) denoted t w o images. The relationship b et w een these t w o images can b e expressed as I 2 ( x; y ; z ) = g ( I 1 ( T ( x; y ; z ))) (1.1) where T is a 3D spatial-co ordinate transformation and g is an in tensit y transformation. The registration problem is to nd the optimal spatial and in tensit y transformation so that the images are matc hed w ell. Finding the parameters of the optimal geometric co ordinate transformation is generally the k ey to an y registration PAGE 12 4 problem while the in tensit y transformation is not alw a ys the task in terest. In this thesis, w e will primarily fo cus on the spatial co ordinate registration. 1.2.2 T ransformations The transformations can b e classied in to global and lo cal transformations. A global transformation is giv en b y a single equation whic h maps the en tire image. Lo cal transformations map the image dieren tly dep ending on the spatial lo cation and are th us more dicult to express succinctly . The most common global transformations are rigid, ane and pro jectiv e transformations. A transformation is called rigid if the distance b et w een p oin ts in the image b eing transformed is preserv ed. A rigid transformation can b e expressed as: 8><>: u ( x; y ) = ( cos ( ) x sin ( ) y + d x ) x v ( x; y ) = ( sin ( ) x + cos ( ) y + d y ) y (1.2) where u ( x; y ) and v ( x; y ) denote the displacemen t at p oin t ( x; y ) along the X and Y directions; is the rotation angle, and ( d x ; d y ) the translation v ector. A transformation is called ane when an y straigh t line in the rst image is mapp ed on to a straigh t line in the second image with parallelism b eing preserv ed. In 2D, the ane transformation can b e expressed as: 8><>: u ( x; y ) = ( a 11 x + a 12 y + d x ) x v ( x; y ) = ( a 21 x + a 22 y + d y ) y (1.3) where 264 a 11 a 12 a 21 a 22 375 (1.4) denotes an arbitrary real-v alued matrix. Scaling transformation, whic h has a transformation matrix of s 1 0 0 ^ s 2 and shearing transformation, whic h has a matrix PAGE 13 5 of ( 1 s 3 0 1 ) are t w o examples of ane transformation, where s 1 ; s 2 and s 3 are p ositiv e real n um b ers. A more in teresting case, in general, is that of a planar surface in motion view ed through a pinhole camera. This motion can b e describ ed as a 2D pro jectiv e transformation of the plane 8><>: u ( x; y ) = m 0 x + m 1 y + m 2 m 6 x + m 7 y +1 x v ( x; y ) = m 3 x + m 4 y + m 5 m 6 x + m 7 y +1 y (1.5) where m 0 ... m 7 are the global parameters. In clinical practice, the most commonly used global registration transformations are rigid and ane. F or the brain images tak en from the same patien t using dieren t devices, e.g, CT/MR, or the same devices but at dieren t times, usually a rigid transformation is adequate to explain the v ariation b et w een them. When t w o images depicting the same scene are tak en from the same viewing angle but from dieren t p ositions, i.e, the camera w as zo omed in/out, or rotated around its optical, an ane transformation is required to matc h these t w o images. In this thesis, w e will b e mainly dealing with these t yp e of global transformations. When a global transformation do es not adequately explain the relationship of a pair of input images, a lo cal transformation ma y b e necessary . T o register an image pair tak en at dieren t times with some p ortion of the b o dy exp erienced gro wth, or to register t w o images from dieren t patien ts falls in to this lo cal transformation registration category . A motion (v ector) eld is usually used to describ e the c hange/displacemen t in lo cal transformation problem. 1.3 Thesis Con tribution Numerous algorithms for registering m ulti-mo dal image data ha v e b een rep orted in literature [11, 2, 3, 4, 5 , 6]. Among them, m utual information (MI) based metho ds are widely accepted as the most eectiv e w a y of handling rigid PAGE 14 6 m ulti-mo dal image registration problems. But MI metho ds ha v e their inheren t dra wbac ks, including b eing computationally in tensiv e, lac king in robustness, and lac k of ease in handling nonrigid motion. The suite of algorithms presen ted in this dissertation is an attempt to pro vide an alternativ e to the MI based metho ds. The ma jor con tribution of this thesis can b e summarized as: A v ector-v alued lo cal frequency image represen tation that is in v arian t to brigh tness and con trast v ariations in an image. A cost function { for computing the lo cal frequency represen tations of misaligned images { based on the exp ectation of the squared dierence and a signican t mo dication of the Newton metho d to n umerically optimize the cost function in a computationally ecien t manner. A statistically robust cost function whose optimization yields an estimate of the registration parameters required to align the giv en image pair. The k ey adv an tage of this tec hnique is its abilit y to cop e with alignmen t of image pairs whic h ha v e a signican t non-o v erlap. 1.4 Thesis Outline The rest of the thesis is organized as follo ws: Chapter 2 reviews the related w ork in handling the m ulti-mo dal image registration problem. This eld has receiv ed considerable atten tion o v er the past decade and no attempt is made here to b e exhaustiv e. Rather, the predominan t approac hes to the registration problem, mostly with medical applications are describ ed, fo cusing on the assumptions that are made and the practical adv an tages and disadv an tages of eac h approac h. Chapter 3 describ es the lo cal frequency represen tation whic h w e are using to retriev e the common information residing in the m ulti-mo dal image pair. The adv an tage/justication as w ell as the pro cedure to compute the lo cal frequency represen tation will b e describ ed in this c hapter. PAGE 15 7 Chapter 4 describ es the approac h used to attac k the computational eciency problem. A similarit y measure Exp ectation of Squared Dierence (ESD) together with a mo died Newton metho d will b e presen ted to sho w the impro v emen t made in ac hieving a signican t sp eed-up in the registration pro cedure. The exp erimen tal comparison with resp ect to a v ery p opular metho d Mutual information based sc heme, will b e giv en in Chapter 7. Chapter 5 presen ts our approac h to handle dieren t eld of view (F O V) problem. A robust estimator | in tegral squared error (ISE), also kno wn as the L 2 error ( L 2 E ) is in tro duced for ac hieving the registration. Chapter 6 describ es the preliminary results of applying our lo cal frequency in a lev el set form ulation [7] of the nonrigid motion alignmen t problem. Chapter 7 presen ts the exp erimen tal results in supp ort of the algorithms describ ed in the previous c hapters. Most sp ecically , the results of comparison b et w een ESD and m utual information, L 2 E and normalized m utual information in aligning misaligned m ulti-mo dal image pairs. Finally the lev el set ev olution metho d for nonrigid motion estimation are all presen ted in this c hapter. Chapter 8 summarizes the main con tributions of this thesis and discusses the directions for future researc h. PAGE 16 CHAPTER 2 RELA TED W ORK T o solv e an y image registration problem, v e asp ects [8] m ust b e considered: 1) image v ariations, 2) feature space, 3) similarit y measure, 4) searc h space and 5) searc h strategy . V ariations are the dierences in lo cations and v alues of pixels b et w een t w o images. Generally image v ariations can b e classied in to three t yp es. The rst one is the dieren t co ordinate system used in the pro cess of acquisition whic h causes images to b e misaligned. F or this case, a geometric transformation accoun ts for the discrepancy . Brigh tness Constancy equation, whic h assumes the apparen t brigh tness of mo ving ob jects remains constan t, holds and n umerous optical ro w tec hniques can b e applied to solv e the problem. The second t yp e is due to not only the co ordinate mis-alignmen t, but the in tensit y v ariations. The third t yp e is caused b y lo cal c hanges, suc h as an ev olving pathology , surgery to implan t or resect an organ from a b o dy , etc. The second and third t yp e of v ariations mak e the brigh tness constancy assumption in v alid and are more dicult to handle than the rst. Kno wledge of the t yp e of image v ariations eects the c hoice of feature space, the set of features extracted from the image and used to carry out the matc hing task. The feature ma y b e ra w image in tensit y v alue; geometric features: salien t p oin ts, edges, con tours, surfaces; statistical features suc h as momen t in v arian t and cen troids; or high lev el structural features whic h con tain meaningful in v arian t information. Similarit y measures are the criteria used to determine the co ordinate and p ossibly the in tensit y transformation b et w een the image pair b eing registered. The similarit y measures can b e dened either to the extracted features or the 8 PAGE 17 9 in tensit y directly . T ypically they include geometric correlations, sum of squared dierence, sum of absolute dierence, phase correlation; statistical measures, suc h as corresp onding in tensit y v ariance, momen ts of join t probabilit y distribution, and m utual information. Searc h space is the class of p ossible transformations the image pair can p ossibly dier b y . Both global and lo cal transformation are included in the class. A global mapping implies that the whole image undergo es the same mapping and its t ypical forms are ane, rigid and p olynomial. A lo cal mapping, on the other hand, is a function of p osition and eac h pixel undergo es dieren t motions. Lo cal motion can b e observ ed in cases suc h as, heart motion, tumor gro wth etc. Searc h strategies are the w a ys to carry out the searc hing for the optimal solution. T ypically they include hierarc hical or m ulti-resolution tec hniques, relaxation, generalize Hough transformation, linear programming and nonlinear programming. In the follo wing, w e will review the literature on existing computational algorithms that ha v e b een rep orted for ac hieving m ulti-mo dal registration. W e will p oin t out their limitations and hence motiv ate the need for a new and ecien t computational algorithm for ac hieving our goal. Broadly sp eaking, image registration metho ds can b e classied in to t w o classes [9] namely , feature-based and direct metho ds. In the former, features are extracted from the t w o images to b e registered and matc hed to estimate the transformation b et w een the t w o data sets. In the latter, the co ordinate transformations b et w een the images are determined directly from the image data or a deriv ed image-lik e represen tation of the same. 2.1 F eature-Based Metho ds This section reviews approac hes based on the extraction and alignmen t of corresp onding features. PAGE 18 10 2.1.1 P oin t-Based F eatures One of the most straigh t forw ard metho ds is to nd the corresp onding landmarks in the t w o images, and then obtain the registration b y forcing these landmark p oin ts in to alignmen t. An easy w a y to nd the transformation b et w een the iden tied p oin t sets is to minimize the mean square distance giv en b y , D is ( T ) = K X i =1 ( x i T ( y i )) 2 : (2.1) where X = f x i g and Y = f y i g , i = 1,..., K are the p oin ts selected from the t w o input images. A w ell-established tec hnique to solv e the ab o v e equation is the use of singular v alue decomp osition (SVD) [10]. T o iden tify the corresp onding landmarks is the rst step to carry out the p oin t based registration. One of the most rexible approac hes to image registration has b een based on man ual p oin t landmark iden tication; salien t and accurately lo catable p oin ts of morphology of the visible anatom y , usually iden tied in teractiv ely b y the user [1, 11 ]. Landmarks can also b e iden tied automatically from the image based on some geometric information [11], i.e., p oin ts at the lo cus of the optim um of some geometric prop ert y , e.g. lo cal curv ature extrema, corners, etc., generally lo calized in an automatic fashion [12]. Landmark-based registration is widely applicable in the sense that, at least in theory , it can b e applied to an y image, regardless what the ob ject or sub ject is, it do es not dep end on an y sp ecic imaging proto col, th us can b e used on an y mo dalit y . 2.1.2 Con tour and Surface Based F eatures P ellizari et al. [13] w ere the rst one using surface matc hing to register dieren t mo dalit y images. The surfaces w ere formed b y delineating the external PAGE 19 11 skin surface in eac h slice of the MR and PET scans. The skin surfaces w ere then used to estimate the transformation parameters b y minimizing a least square criterion. The minimization is carried out n umerically using P o w ell's metho d. In order to nd the corresp ondence in eac h iteration, they made an assumption that the t w o surfaces are predominan tly spherical. The task of extracting corresp onding shap es in their metho d is quiet easy , but the accuracy of the rigid alignmen t is not guaran teed due to the skin surface deformation and so is the robustness of the matc hing due to the symmetries and lac k of the ne structure in the skin surface. A n um b er of p o w erful, alternativ e surface/feature-based strategies w ere no w b eing dev elop ed as w ell. F or sets of delineated p oin ts, man y sa w the Besl and Mac k a y iterativ e closest p oin t (ICP) app oarc h as an attractiv e idea [14]. Ho w ev er, this strategy could b e quiet sensitiv e to c hanges in feature lo cations and false p ositiv es in the detection of closest p oin ts. This problem w as partially xed in the adaptiv e-threshold w ork of F eldmar and Ay ac he [15]. In addition, Rangara jan et al. [16] noted the problems with robustness to outliers and prop osed an alternativ e sc heme based on assignmen t matrices and robust statistics kno wn as robust p oin t matc hing (RPM). The ICP and RPM strategies implicitly sim ultaneously solv e the corresp ondence and matc hing problem. Deform-able mo dels later w ere used to facilitate the mo deling of the surface shap e in an image. Li et al. [17] presen ted t w o con tour-based metho ds whic h use region b oundaries and other \strong" edges as matc hing primitiv es. The matc hing algorithm is based on the c hain-co de correlation and other shap e similarit y criteria suc h as in v arian t momen ts. Ev ans et al. [18 ] dev elop ed a registration sc heme based on appro ximating the 3D w arp b et w een the mo del and target image b y a 3D thin plate spline tted to landmarks. Szeliski et al. [19 ] dev elop ed a fast matc hing algorithm for registration of 3D anatomical surfaces found in 3D medical image data. They use a distance metric minimization to nd the optimal transformation PAGE 20 12 b et w een the surfaces. The k ey dierence in their sc heme from earlier metho ds is the use of tensor pro duct splines to represen t dieren t registration transformations, namely rigid, ane, trilinear, quadratic etc. In addition, they also in tro duced the no v el concept of o ctree splines for v ery fast computation of distance b et w een surfaces. Da v atzik os et al. [20] in tro duced a t w o stage brain image registration algorithm. The rst stage in v olv ed using activ e con tours to establish a one-toone mapping b et w een cortical and v en tricular b oundaries in the brain and in the second stage, an elastic deformation transformation that determines the b est corresp ondence b et w een the iden tied con tours in the t w o data sets is determined. F eldmar et al., [15] dev elop ed a no v el surface to surface nonrigid registration sc heme, using lo cal ane transformations. 2.1.3 Edge/Ridge Based Metho ds One of the dra wbac ks of the surface based metho ds is they essen tially lose v ery imp ortan t information b y extracting only the information from surface, with most of the in ternal tissue information b eing discarded. Main tz et al. [21 ] compared edge [22] and ridge-based [23, 24, 25 ] registration of CT and MR brain images. They describ e a no v el metho d that mak es use of fuzzy edgeness and ridgeness images in ac hieving the registration b et w een MR and CT data sets. F or more on feature-based metho ds, w e refer the reader to the surv ey b y Main tz et al. [26]. All of these approac hes ha v e one commonalit y , i.e., they need to detect features/surface/con tours in the images and hence the accuracy of registration is dictated b y the accuracy of the feature detector. Also, additional computational time is required for detecting these features prior to application of the actual registration sc heme. PAGE 21 13 2.1.4 Reductiv e Metho ds Principal-axis and momen ts-based metho ds are the ma jor examples of the reductiv e metho ds, within whic h the image cen ter of gra vit y and its principal orien tation (principal axes) are used as the matc hing en tit y to carry out the registration [28]. Sometimes, higher-order momen ts are also computed and used in this pro cess. Hill et al. [11] examined the eects of misregistration on the feature space and prop osed the use of the third order momen t as a measure for MR/CT alignmen t, whic h has b een sho wn mathematically to b e related to information theory measure. The results based on the momen ts-based metho ds are not v ery accurate, but this do es not prev en t p eople from widely using those metho ds, esp ecially in the registration problems that do not require high accuracy , b ecause of the automatic and v ery fast nature of its use, and easy implemen tation. 2.2 Direct Metho ds Direct metho ds subsume the approac hes op erating directly on the image grey v alues, without prior feature extraction. Pixel/v o xel based metho ds using the full image grey v alues are the most activ e direction of curren t researc h. They do not start b y reducing the gra y-lev el image to relativ ely sparse extracted information, but use all a v ailable information throughout the whole pro cedure. One straigh tforw ard in tensit y v alue based approac h is the lo cally adaptiv e correlation windo w sc heme [29 ], whic h can ac hiev e high accuracy . But it is computationally demanding [30 ]. A more general sc heme than windo w-based correlation approac h is the optical ro w form ulation, in whic h the problem of registering t w o images is treated as equiv alen t to computing the ro w b et w een the data sets. There are n umerous tec hniques for computing the optical ro w from a pair of images [31 , 32 , 33, 34, 30, 35, 36 ]. Zhao et al. [27 ] presen t an optical ro w based registration sc heme for aligning/registering coronal sections to form a 3D stac k PAGE 22 14 of registered auto-radiographic images. In their algorithm, they assume the ro w is parameterized b y a general ane parameters and use linear constrain ts. The registration is p erformed in 2D yielding a stac k of registered 2D coronal slices. 2.2.1 Image Similarit y Measures Another class of in tensit y v alue based approac hes is based on statistical similarit y measures. W e can express an o v erall image similarit y D for a transformation T b et w een the t w o images, as sum of individual measuremen t similarities S ( : ), at K dieren t p oin ts, D ( T ) = K X i =1 S ( I t ( x i ) ; I s ( T ( x i ))) (2.2) where I t ( x i ) and I s ( T ( x i )) are corresp onding v alues. Cross-Correlation Based Metho ds. In the simplest case, the gra y v alues of input image pair con tain a strong linear relationship. Misalignmen t then can b e remo v ed b y maximizing the correlation of the t w o images. More sp ecially , for a 3D target image I t and a source image I s , the cross-correlation function measures the similarit y o v er all transformations and is giv en b y: C ( T ) = P X I t ( X ) I s ( T ( X )) p ([ P X I s ( T ( X ))] 2 ) (2.3) If the source image matc hes the reference image exactly , barring an in tensit y scale factor, at a transformation T , the cross-correlation C ( T ) will ha v e its p eak at C ( T ). V ariance of In tensit y Ratio. This is the rst and simplest statistical measure prop osed b y W o o ds et al. [2 ] for registering PET and MRI images. The algorithm is based on an assumption that there will b e a predominan tly one-to-one mapping of v alues in the images at registration, then for all p oin ts in an image with a giv en v alue, the v ariance of all corresp onding PAGE 23 15 v alues in the other image will b e minimized at registration. More sp ecically , the algorithm is giv en as follo ws: When the pair of PET-MRI images are accurately registered with a transformation T , then for eac h p osition i , the ratio of the corresp onding v alues is ev aluated: r i ( T ) = I t ( X i ) I s ( T ( X i )) (2.4) Let r b e the standard deviation of this ratio r i , ev aluated o v er all v o xel pairs. Based on the author's assumption, this r is minimized at correct registration. Join t In tensit y Distribution. Lev en ton et al. [3] prop osed a metho d based on matc hing the join t in tensit y distribution of curren t input image with the prior join t in tensit y distribution obtained from training data sets. Tw o metho ds are used to mo del the join t in tensit y distribution of the training data, mixture of Gaussians and P arzen windo wing. The tted Gaussians corresp onds to dieren t anatomical structures and pro vide a rough segmen tation for the registered image pair. T o matc h unregistered input image pair, the algorithm computes the log lik eliho o d of the t w o images and the b est registration is reac hed at its maximal. Results of applying of this algorithm to SPGR and MR scans demonstrate an impressiv e p erformance in term of sp eed and region of con v ergence, with a sub-v o xel registration accuracy . The dra wbac k of this approac h lies in the fact that registered training data are required to get the prior distribution of the join t in tensit y , whic h is not alw a ys a v ailable. Mutual Information. Curren tly the most p opular approac h is based on the concept of maximizing m utual information rep orted in Viola and W ells [4] and W ells et al. [5], Collignon et al. [37 ] and Studholme et al. [6]. Mutual information is a measure from the eld of information theory . F or t w o images, m utual information is a measure of ho w w ell one image PAGE 24 16 explains the other, or vice v ersa. In other w ords, giv en one image, ho w m uc h information do es that giv e ab out the other image? The primary computational task in the m utual information based tec hniques is the estimation of probabilit y densit y functions and their deriv ativ es from the giv en data sets. Tw o original w orks [4, 37] used dieren t approac hes to compute the join t probabilit y densit y . While P arzen Windo w approac h [38 ] has long b een established as a standard metho d for densit y estimation, as used in Viola et al.'s implemen tation [4], using histogram to appro ximate the true probabilit y of densit y function (p df ), as used in [37] gained the p opularit y due to its computation eciency and easy implemen tation. Mutual information measure is dened b y: M I ( I s ; I t ) = H ( I s ) + H ( I t ) H ( I s ; I t ) ; (2.5) with M I ( I s ; I t ) the m utual information of t w o images I s and I t . H stands for en trop y , whic h is a measure of disp ersion of a probabilit y distribution. T o register the images the m utual information is to b e maximized. Lo oking at the denition, this means that the join t en trop y H ( I s ; I t ) is minimized. The adv ances of m utual information based metho ds reside not only in the impressiv e accuracy in the rep orted registration exp erimen ts, but also the generalit y the MI metho ds can pro vide. V ery few assumptions ev er made regarding the relationship that exists b et w een the image in tensities, so m utual information is esp ecially suited for m ulti-mo dalit y matc hing and that it is completely automatic. 2.2.2 Dra wbac ks of MI-based Metho ds and the Related Impro v emen ts The dra wbac ks of original MI-based metho ds [4, 37 ] can b e summarized as follo ws: PAGE 25 17 Lac k of computational eciency . This is partially due to the complexit y of MI's ob jectiv e function. Lac k of robustness, the accuracy dep ends on the initial guess and is sensitiv e to the c hange of eld of view [39] and image in terp olation eect. Since m utual information (MI) w as prop osed to handle registration problem, a considerable amoun t of researc h has b een done to impro v e the p erformance of MIbased metho ds. In the follo wing, w e will fo cus on the researc h w ork rep orted on, computational eciency , robustness to outliers and eect of c hoice of in terp olation used. 2.2.3 Computational Eciency All relev an t results in published literature rep ort the sp eed-up factors with resp ect to the P o w ell's metho d used for optimizing the MI-based matc hing criterion. Maes et al. [ ? ] claimed that a factor of 3 on the a v erage can b e ac hiev ed o v er P o w ell's metho d with the do wnhill simplex, conjugate-gradien t and Lev en b ergMarquardt metho ds using a t w o-lev el m ulti-resolution sc heme. Thev enaz et al. [53] used a spline-based optimizer and claimed a sp eed-up factor of 6 with resp ect to P o w ell's metho d ([37 ]) and Netsc h et al. [54] ac hiev ed a sp eed-up factor of 7 o v er the w ork rep orted in Studholme et al. [6], who used the P o w ell's metho d as w ell. The sp eed-up factor of the algorithm dev elop ed in this thesis (see c hapter 4) with resp ect to P o w ell's sc heme [37] on the a v erage is 13, whic h is a m uc h larger factor than what the rest of the aforemen tioned metho ds ac hiev ed. Note that this comparison of sp eed-up factor is mac hine indep enden t. 2.2.4 Robustness Studholme et al. [6] p oin ted out the standard m utual information is sensitiv e to the eld of view of the scanner used in the image acquisition, namely , with dieren t c hoices of eld of view, the maximization pro cess ma y lead to an incorrect PAGE 26 18 alignmen t. The authors then extended the m utual information to a normalized m utual information to alleviate this problem. Instead of only computing the in tersection of region I s and I t , the whole area of I s + I t is used to a v oid the eect of dieren t R OI: N M I ( I s ; I t ) = ( H ( I s ) + H ( I t )) =H ( I s ; I t ) (2.6) The exp erimen t results [6 ] demonstrated the impro v emen t made b y applying NMI o v er standard MI in term of handling image pair acquired b y scanners with some non-o v erlap in the F O Vs. 2.2.5 In terp olation Artifacts An ideal matc hing metric is a smo oth, unimo dal function of the registration parameter T , whose global maxim um p oin t corresp onds to the correct registration T . Ho w ev er, for the real computer-based implemen tation, w e ha v e to deal with the discrete v alues of data size. In terp olation step is required to equalize the v o xel sizes of the t w o images. F or en trop y based registration measures suc h as m utual information, the registration function of these equalized images can con tain in terp olation artifacts [40 ]: patterns of lo cal extrema. In terp olation-induced artifacts imp ede the registration optimization pro cess and more imp ortan tly , they rule out sub-v o xel accuracy . T o o v ercome this problem, sp ecial in terp olators or c hange of the MI energy form ulation are needed [41]. 2.2.6 Extensions to Cop e with Nonrigid Motion Most m utual information based algorithms for image registration in literature ha v e b een form ulated for global parameterized motion with the exception of the recen t w ork rep orted in Mey er et al. [42] wherein ane transformations as w ell as thin-plate spline w arps are handled. The rep orted registration results are quite impressiv e but the CPU execution times are quite high { of the order of sev eral hours for estimating thin-plate w arps. The problem of b eing able to PAGE 27 19 handle lo cal deformations in a m utual information framew ork is a v ery activ e area of researc h and some recen tly rep orted results on this problem ma y b e found in [3, 43 , 44 , 45, 46, 47]. 2.3 In tro duction of Our Metho d This thesis is an attempt to pro vide an alternativ e to MI metho ds to o v ercome the inheren t dra wbac ks in MI. The lo cal frequency image represen tation is obtained b y ltering the image with a Gab or lter tuned to a certain frequency/orien tation and then computing the gradien t of the phase of the ltered image. Note that the Gab or lter is a complex v alued function in the spatial domain and is obtained b y mo dulating a Gaussian function where the mo dulating factor is a complex exp onen tial. When dealing with complex v alued functions, phase is dened as the ar ctan ( r eal par t=imag inar y par t ) and lies b et w een 0 and 2 . Th us, the phase of a complex v alued function will con tain a wrap around whenev er it exceeds 2 . In order to a v oid this wrap around, one can use the gradien t of phase that (in practise) can b e computed directly from the function without ha ving to compute the phase itself. In this thesis, instead of using global phase of a function/image { whic h is obtained b y taking the F ourier transform of a function/image and discarding the amplitude { w e use its lo cal phase and its gradien t that b ears the adv an tage of b eing tunable to an y frequency/orien tation. Suc h tuning is particularly useful in capturing an y predominan t frequencies and/or orien tations presen t in an image. In addition, lo cal phase/frequency represen tation is relativ ely in v arian t to c hanges in illumination conditions whic h cause gra y v alue c hanges. What do es this mean in the con text of m ulti-mo dal registration? Because, m ultimo dal images of the same underlying scene/ob ject ha v e completely dieren t in tensities, a lo cal phase/frequency based represen tation seems apt to capture the salien t common features b et w een the t w o mo dalities. it This is the premise underlying the dev elopmen t of our m ulti-mo dal registration algorithm describ ed in PAGE 28 20 this thesis (see Chapter 3). Once, w e compute this lo cal frequency represen tation for eac h of the t w o (source and target) images to b e registered, w e are ready to nd the registration transformation whic h will b est matc h these represen tations. Sev eral matc hing criteria ma y b e dened. The rst matc hing measure w e used is based on exp ectation of the squared dierences (ESD) b et w een the transformed source and the target (lo cal frequency) represen tations. This matc hing criteria is minimized o v er all rigid/ane transformations using a v ery fast n umerical tec hniques dev elop ed b y V em uri et al. [9]. The n umerical algorithm is a mo dication of the standard quasi-Newton sc heme and the mo dication in v olv es pre-computation of the Hessian at the optim um prior to kno wing the optim um. This n umerical algorithm has a quadratic con v ergence rate and has t wice the con v ergence range of the standard Newton sc heme. This means that w e can start with an initial guess of 0 (zero v ector) for the estimation of the transformation parameter v ector in most cases ev en for large rotations and translations. In addition, there is a big computational adv an tage in pre-computing the Hessian of the ob jectiv e function b eing minimized. On the con trary , in the standard Newton metho d or in the La v en b erg-Marquardt metho d [60], one needs the Hessian computation at ev ery iteration whic h pro v es to b e exp ensiv e as w ell as unstable. The ESD approac h ga v e v ery impressiv e results in b oth timing and accuracy with real image. But with real CT-MR data that ha v e non-signican t o v erlap, ESD's p erformance degraded signican tly unless considerable prepro cessing in v olving elimination of some of the non-o v erlapping regions w as done. As an alternativ e matc hing criteria, w e designed a robust measure based on the in tegral of the squared error (ISE or L 2 E) b et w een a Gaussian mo del of the residual and its true densit y function [61]. The residual here refers to the dierence b et w een the lo cal frequency represen tations of the transformed (b y an unkno wn transformation) source and target data. This robust estimation framew ork aords PAGE 29 21 the adv an tage of not ha ving to deal with additional tuning parameters that w ould b e required when using inruence functions to ac hiev e robustness. W e minimize this L 2 E matc hing criteria o v er all rigid/ane transformations using a sequen tial quadratic programming tec hnique. W e presen t results of applying this algorithm to MR T1/T2 as w ell as MR/CT images that are misaligned due to real motion of the patien t during imaging. One of the k ey strengths of our metho d is that { due to the form ulation b eing set in a robust framew ork { it can handle data sets that do not ha v e signican t o v erlap as is the case in most practical situations. T o ac hiev e nonrigid registration, w e can use the same lo cal-frequency represen tations of the t w o images to b e registered and apply the tec hnique of lev el-set based ev olution (describ ed in the c hapter 6) b et w een the t w o lo cal-frequency represen tations to ac hiev e the nonrigid registration. W e presen t syn thetic and real data exp erimen ts to demonstrate the p erformance of our nonrigid registration algorithm. The k ey idea used w ere here in ac hieving v ery large deformations is the use of a Lagrangian (mo ving) grid implemen tation of the go v erning equation of motion. An example of nonrigid registration for in ter-mo dalit y (regular MR and con trast agen t injected MR breast image pair) data are presen ted in the exp erimen tal results c hapter. PAGE 30 CHAPTER 3 INV ARIANT IMA GE REPRESENT A TION T o design a registration algorithm for m ulti-mo dal images, t w o fundamen tal questions ha v e to b e answ ered: 1) what is a go o d image represen tation to w ork with? Namely , what represen tation will bring out the common information b et w een the t w o m ulti-mo dal images, while suppressing the irrelev an t information. 2) what is an appropriate similarit y measure to matc h the images with the selected represen tation. Numerous researc h articles ha v e b een published on retrieving information shared b y the m ulti-mo dal data sets while reducing or eliminating the (imaging) sensor dep enden t bac kground information. In [48], the authors presen ted an energy-image represen tation based on directional-deriv ativ e lters. F our directional deriv ativ e lters, orien ted in the horizon tal, v ertical, and the t w o diagonal directions, are applied to the ra w image and the deriv ativ e image is squared yielding an \energy" image . Th us the directional information is preserv ed in this energy represen tation. This form ulation also requires explicit directional lters and explicit ltering with Gaussian functions to create a p yramid. Considering the inheren t c haracteristic of the registration problem as w ell as the computational asp ects, it is reasonable to list follo wing p oin ts as the desired prop erties for the represen tations used for m ulti-mo dalit y registration problem: T o capture the common information b et w een the image pairs. Main tain as m uc h information as p ossible; in ligh t of this, con tin uous/dense represen tations are more desirable than discrete/sparse forms. 22 PAGE 31 23 Should facilitate m ulti-resolution computation, b ecause this can not only mak e it p ossible to cop e with large motions, but also b o ost the con v ergence sp eed as w ell. Must main tain directional information, since, without it, false corresp ondences ma y result for orien ted patterns. 3.1 Lo cal F requency Represen tation Based on these guidelines, w e select the Lo cal F requency , whic h is obtained b y Gab or ltering the input images and then computing the gradien t of the phase, as the in v arian t image represen tation to tac kle the registration problem. This lo cal frequency represen tation has all the adv an tages of an \energy" represen tation and additionally can b e tuned to an y desired frequency/orien tation; hence it facilitates con trol of the alignmen t pro cess. In this thesis, w e estimate the lo cal frequency to b e the one corresp onding to the maximal lo cal energy computed from a suite of frequency tuned Gab or lters and w e call it the dominan t lo cal frequency . In order to dene the lo cal frequency represen tation, w e will rst in tro duce the concept of an analytic signal/function in 1D whic h can b e generalized to higher dimensions [49]. Giv en a signal f in 1D, its analytic signal is dened as f A = f if H i , where f H i is the Hilb ert transform of f , giv en b y f H i = 1 = R + 1 1 ( f ( ) = ( x )) d : The transformation from the real signal to its corresp onding analytic signal can b e regarded as the result of con v olving the real signal with a complex lter, whic h is called a quadrature lter. The argumen t of f A is referred to as the lo cal phase of f , whic h is dened in the spatial domain. The spatial deriv ativ e of lo cal phase is called lo cal or instan taneous frequency [49]. Consulting the desirable c haracters list in 3, the follo wing prop erties of lo cal frequency mak e it a candidate for an in v arian t image represen tation [49 ] in matc hing m ulti-mo dal image data sets: Lo cal frequency is in v arian t to con trast rev ersal. PAGE 32 24 Lo cal frequency estimates are relativ ely in v arian t to signal energy . Lo cal phase v aries in generally the same manner regardless of whether the signal v ariations are small or big, so lo cal phase and frequency are relativ ely insensitiv e to c hanges in illumination conditions. In our dominan t v ector v alued Lo cal F requency represen tation, directional information is main tained. The standard deviation of the Gab or lters, whic h will b e describ ed later, is a built-in scale parameter to facilitate m ulti-resolution computation. Lo cal phase/frequency of an image can b e used to discriminate b et w een a v ariet y of line and edge structures whic h giv e rise to particular orien tation estimates. Its use for the analysis of lo cal image structure has b een explored extensiv ely [50, 49 ]. 3.1.1 Gab or Filters The Gab or Filter is a w ell-kno wn quadrature lter. It ac hiev es the theoretical lo w er b ound of the uncertain t y principle [49] i.e., Gab or functions are optimal in terms of their space-frequency lo calization. The complex 2-D Gab or functions ha v e the follo wing general form (see Bo vik et al. [51] for details on follo wing denitions): h ( x; y ) = g ( x 0 ; y 0 ) exp (2 j ( U x + V y )) ; (3.1) where ( j = p 1), ( x 0 ; y 0 ) T = R ( x; y ) with R b eing the 2D rotation matrix and g ( x; y ) = ( 1 2 2 ) exp ( ( x ) 2 + y 2 2 2 ) : (3.2) Therefore, h ( x; y ) is a complex sin usoid mo dulated b y a 2-D Gaussian function whose asp ect is giv en b y , scale parameter b y , and the ma jor axis making an angle with the x-axis. F or = 1, g ( x; y ) is circularly symmetric and hence need not b e sp ecied. The transfer function i.e., the frequency resp onse of the Gab or PAGE 33 25 10 20 30 40 50 60 10 20 30 40 50 60 10 20 30 40 50 60 10 20 30 40 50 60 0 10 20 30 40 50 60 70 0 20 40 60 80 -1.5 -1 -0.5 0 0.5 1 1.5 x 10 3 0 10 20 30 40 50 60 70 0 20 40 60 80 1 0.5 0 0.5 1 1.5 2 x 10 3 Figure 3.1: In tensit y(top) and p ersp ectiv e (b ottom) plot of the real (left) and imaginary (righ t) parts of a v ertically orien ted 2D Gab or lter function in eqn. 3.1 is giv en b y H ( u; v ) = exp f 2 2 2 [( u 0 U 0 ) 2 2 + ( v 0 V 0 ) 2 ] g (3.3) where ( u 0 ; v 0 ) = ( ucos + v sin; usin + v cos ) and ( U 0 ; V 0 ) is a similar rotation of the cen tral frequency ( U; V ). Hence, H ( u; v ) happ ens to b e a Gaussian bandpass lter with its minor axis at an angle from the u-axis, an asp ect parameter of 1 = , a radial cen tral frequency ! = p ( U 2 + V 2 ) and orien tation = tan 1 ( V =U ) with resp ect to the u -axis. Figure 3.1 is an in tensit y plot of the real and imaginary parts of a 2D Gab or lter. Similarly , the complex 3-D Gab or lters, whic h are used in our implemen tation, ha v e the follo wing general form: h ( x; y ; z ) = g ( x 0 ; y 0 ; z 0 ) exp (2 j ( U x + V y + W z )) ; (3.4) where ( j = p 1), ( x 0 ; y 0 ; z 0 ) T = R ( x; y ; z ) T with R b eing a 3D rotation matrix and g ( x; y ; z ) = 1 (2 ) 3 = 2 3 exp x 2 + y 2 + z 2 2 2 : (3.5) PAGE 34 26 Therefore, h ( x; y ; z ) is a complex sin usoid mo dulated b y a 3D Gaussian that is assumed to b e circularly symmetric but in general need not b e. Also, w e assumed that x = y = z = . 3.1.2 Computing Lo cal-frequency Represen tations The pro cess of computing the lo cal frequency represen tation can b e ac hiev ed in three steps: generate a set of Gab or lters G k ( x; y ; z ) tuned to spatial frequencies ! k , for k = 1 ; :::N ; let q ( k ) + ( x; y ; z ) and q ( k ) ( x; y ; z ) b e the result of con v olution of the image I with the real and imaginary parts of G k resp ectiv ely . q ( k ) + ( x; y ; z ) = I n r eal ( G k ) and q ( k ) ( x; y ; z ) = I n imag ( G k ). Compute the output magnitude M k ( x; y ; z ) = q ( q ( k ) + ( x; y ; z )) 2 + ( q ( k ) ( x; y ; z )) 2 . Select the lter k max with maxim um output magnitude for eac h pixel and compute the lo cal phase gradien t (lo cal frequency estimator) for eac h of these lters using the follo wing equation: r k max ( x; y ; z ) = q ( k max ) + ( x; y ; z ) r q ( k max ) ( x; y ; z ) q ( k max ) ( x; y ; z ) r q ( k max ) + ( x; y ; z ) ( q ( k max ) + ( x; y ; z )) 2 + ( q ( k max ) ( x; y ; z )) 2 : (3.6) Where, k max ( x; y ; z ) = ar ctan ( q ( k max ) ( x; y ; z ) =q ( k max ) + ( x; y ; z )). Note that k max need not b e computed explicitly . Set the v ector-v alued lo cal frequency estimator F ( x; y ; z ) = r k max ( x; y ; z ) The lo cal frequency magnitude is giv en b y F = jj F ( x; y ; z ) jj . Ideally all lters whic h constitute a basis of the frequency space should b e used, but in practice, w e usually select a subset whic h co v ers a part of space near the origin. Figure 3.2 depicts a pair of slices from a set of MR T1 and T2 w eigh ted images and the asso ciated v ector v alued (second ro w) as w ell as scalar v alued (third ro w) lo cal frequency represen tations. Line segmen ts are used to indicate b oth the direction and magnitude of the lo cal frequency v ector in second ro w of 3.2. Note PAGE 35 27 that lo cal frequency is the same for corresp onding structures in the 2 images, ev en when edge strength and con trast ha v e signican t dierences. Another notable prop ert y of a lo cal frequency image represen tation is its scalabilit y: the Gaussian scale parameter in equation 3.5 can b e v aried to directly generate an image scale space represen tation. Scale space prop erties ha v e long b een exploited in the computer vision comm unit y b oth from a computational and lo calization of features p oin t of view. Using scale-space i.e., m ulti-resolution pro cessing leads to computationally ecien t algorithms and alleviates the problem of lo cal-minima en trapmen t when used in the con text of non-con v ex function optimization e.g., the registration problem describ ed in the next section. PAGE 36 28 50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 local freq mag 50 100 150 200 250 50 100 150 200 250 local freq mag 50 100 150 200 250 50 100 150 200 250 Figure 3.2: Images and their corresp onding lo cal frequency represen tations. Left column: T1 w eigh ted, Righ t column: T2 w eigh ted MR data. In tensit y image (rst ro w), corresp onding lo cal frequency v ector images (second ro w) and magnitude images (third ro w). PAGE 37 CHAPTER 4 COMPUT A TIONAL EFFICIENCY ISSUE F or an y registration algorithm, the second comp onen t need to b e considered is a similarit y measure to matc h the presen tation selected. In this thesis, w e dev elop ed t w o similarit y measures to handle dieren t practical demand. W e will presen t the measure to handle the time eciency issue in this c hapter, and the robustness issue in next c hapter. Along with the increasing resolution of mo dern medical imaging devices, computational eciency has b een a more and more imp ortan t issue in 3D medical image pro cessing. Mutual information based metho ds, as w e in tro duced in c hapter 2, ha v e b een kno wn as computationally slo w. The reason is partially due to the complexit y of their ob jectiv e functions, whic h mak es a lot of high eciency n umerical solutions not suited to b e used to estimate the registration. Our implemen tation of one of the original MI w orks, Collignon sc heme [37], has a CPU time ranges from 20 70 min utes for data of size (512 512 128) on SGI on yx mac hine. Curren tly the fast implemen tation of 3D m ulti-mo dal algorithm, other than ours, is rep orted b y Netsc h et al. [54], whic h claims it can run 7 times than Collignon metho d [37]. Our w ork, presen ted in this c hapter, used the Exp ectation of the Squared Dierence as the matc hing measure together with a mo died Newton metho d can run 13 times faster than Collignon metho d [37] with comparable accuracy . 29 PAGE 38 30 4.1 Matc hing Lo cal F requency Represen tations with Eciency T o matc h the lo cal-frequency represen tations of the image pair, w e adopt the hierarc hical spline-based metho d dev elop ed in [9]. Let the lo cal frequency representations of the m ulti-mo dal image pairs dier b y a lo cal displacemen t whic h can b e parameterized as a global rigid transformation, then, the follo wing equation holds for the lo cal-frequency magnitude represen tations: F 1 ( X + T ) = F 2 ( X ) where, F 1 ( ; ) and F 2 ( ; ) are the 3D lo cal frequency magnitude image represen tations computed from the m ulti-mo dal data sets resp ectiv ely and T = ( u; v ; w ) is the 3D displacemen t eld. Our goal is to minimize the exp ectation of the squared dierences [9] giv en b y E sd ( T ) = E f [ F 1 ( X + T ) F 2 ( X )] 2 . This is a nonlinear function of T and its minimization can only b e ac hiev ed using n umerical tec hniques. W e prop ose to represen t the displacemen t eld in a B-spline basis whic h will implicitly incorp orate smo othness assumption on the deformation eld and accrues computational adv antages b ecause the displacemen t need only b e estimated at a few con trol v ertices and in terp olated elsewhere. 4.1.1 Spline Based Represen tation of the Registration T ransformation Let the displacemen t eld b e T = ( u ( x; y ; z ) ; v ( x; y ; z ) ; w ( x; y ; z )) T and the computed eld b e ( ^ u; ^ v ; ^ w ), whic h is represen ted b y B-splines with a small n um b er of con trol p oin ts ^ T = ( ^ u j ; ^ v j ; ^ w j ), where subscript j denotes the con trol p oin t on a B-spline con trol mesh. Then, the displacemen t at a v o xel lo cation i is giv en b y u ( x i ; y i ; z i ) = X j ^ u j q ij ; v ( x i ; y i ; z i ) = X j ^ v j q ij and w ( x i ; y i ; z i ) = X j ^ w j q ij : (4.1) Where, q ij = B j ( x j ; y j ; z j ) = B ( x ^ x j ; y ^ y j ; z ^ z j ). In our implemen tation, w e let B ( x; y ; z ) = (1 j x j )(1 j y j )(1 j z j ) on [ 1 ; 1] 3 . Let F m denote an in termediate lo cal-frequency represen tation obtained b y transforming F 1 using the B-spline represen tation of the transformation and giv en b y F m ( X i ; ^ T ) = F 1 ( x i + P j ^ u j q ij ; y i + PAGE 39 31 P j ^ v j q ij ; z i + P j ^ w j q ij ). Where, X i = ( x i ; y i ; z i ), ^ T = ( ^ u 1 ; ^ v 1 ; ^ w 1 ; :::; ^ u n ; ^ v n ; ^ w n ) T is the estimated displacemen t eld at the n con trol p oin ts and q ij are the basis functions as b efore. Rewriting the error criterion E sd ( ^ T ) w e ha v e, E sd ( ^ T ) = E f ( F m ( X i ; ^ T ) F 2 ( X i )) 2 g : (4.2) F or a fast solution to this non-linear optimization problem, a mo died quasiNewton sc heme dev elop ed in [9] is used. 4.1.2 3D Rigid T ransformation In man y applications, the whole image undergo es the same t yp e of motion, in whic h case, it is p ossible to parameterize the displacemen t eld b y a small set of parameters describing the global motion, e.g., rigid transformations in this case. A 3D rigid transformation ma yb e parametrically expressed with six unkno wn parameters to express the rotation and translation (three Euler angles and three translation comp onen ts). The displacemen t of v o xel ( x; y ; z ) is then dened as follo ws: 266664 u ( x; y ) v ( x; y ) w ( x; y ) 377775 = 266664 r 11 r 12 r 13 r 21 r 22 r 23 r 31 r 32 r 33 377775 266664 x yz 377775 + 266664 d 1 d 2 d 3 377775 266664 x y z 377775 (4.3) Where the (3 ; 3) matrix with en tries ( r 11 ; :::; r 33 ) is the 3D rotation matrix determined b y the three rotation angles 1 ; 2 ; 3 along the x, y and z axis resp ectiv ely , and d 1 ; d 2 and d 3 are the translation comp onen ts in the three axis directions resp ectiv ely . Let m = ( r 11 ; r 12 ; r 13 ; d 1 ; r 21 ; r 22 ; r 23 ; d 2 ; r 31 ; r 32 ; r 33 ; d 3 ) T . T o compute an estimate of the rigid motion, w e rst dene the displacemen ts ^ u j = ( ^ u j ; ^ v j ; ^ w j ) at PAGE 40 32 the con trol p oin ts in terms of the global motion parameters, i.e., ^ u j = 266664 ^ x j ^ y j ^ z j 1 0 0 0 0 0 0 0 0 0 0 0 0 ^ x j ^ y j ^ z j 1 0 0 0 0 0 0 0 0 0 0 0 0 ^ x j ^ y j ^ z j 1 377775 m 266664 ^ x j ^ y j ^ z j 377775 = j m ^ X j : (4.4) Where j represen ts the (3 ; 12) matrix in equation 4.4. Substituting in to equation (4.2), the error criterion E sd ( ^ T ) for the rigid motion case b ecomes E sd ( ^ T ) = E f ( F 1 ( X i + X j w ij ( j m ^ X j )) F 2 ( X i ) 2 ) g : (4.5) 4.2 Numerical Solution W e no w briery describ e the mo died Newton sc heme dev elop ed in [9, 55 ] and emplo y it to nd the n umerical solution to the equation 4.2. The primary structure of the algorithm is giv en b y the follo wing iteration form ula ^ T k +1 = ^ T k H 1 ( ^ T = ^ T ) g ( ^ T k ) ; (4.6) where H is the Hessian matrix and g is the gradien t v ector. Unlik e in a t ypical Newton iteration, in the ab o v e equation, observ e that the Hessian is alw a ys computed at the optim um ^ T = ^ T instead of the iteration p oin t ^ T k . So, one of the k ey problems is ho w to compute the Hessian at the optim um prior to b eginning the iteration. In general, the Hessian of the ob jectiv e function E sd ( ^ T ) with resp ect to the parameters ^ T is not indep enden t of the parameters of the transformation with the exception of the case in v olving only translation. Ho w ev er, in the general case, a clev er tec hnique in v olving the use of inno v ations (incremen tal transformations) with a mo ving co ordinate system f X k = ( x; y ; z ) k g and an in termediate motion v ector ~ T w as in tro duced in [55] and adopted in [9] to dev elop the form ulas for pre-computing the Hessian. This in termediate transformation v ector giv es the motion from f X k g PAGE 41 33 to f X k + 1 g . X k = h ( X 0 ; ^ T k ); (4.7) X k +1 = h ( X k ; ~ T k +1 ) = h ( X 0 ; ^ T k +1 ); (4.8) F m ( X 0 ; ^ T k +1 ) = F m ( X 0 ; ^ T k ~ T k +1 ) (4.9) Instead of dieren tiating the error criterion with resp ect to ^ T , the error criterion no w has to b e dieren tiated with resp ect to ~ T . The inno v ation ~ T k +1 is dened as ~ T k +1 = ~ H 1 ~ g ( ^ T k ) ; (4.10) where, ~ H is the Hessian matrix at the optim um, ~ g is the gradien t v ector at the iteration p oin t ^ T k . The Hessian and the gradien t v ector for the rigid motion case can b e written do wn as follo ws: ~ H = 2 E ( A @ I 1 ( X ) @ X @ I 1 ( X ) @ X T A T ) : (4.11) where, A = @ h ( X ; T ) @ T T = 2666666666666664 P j w j ^ y j P j w j ^ z j 0 P j w j ^ x j 0 P j w j ^ z j 0 P j w j ^ x j P j w j ^ y j P j w j 0 0 0 P j w j 0 0 0 P j w j 3777777777777775 : (4.12) and the gradien t v ector at the optim um is ~ g ( ^ T k ) = 2 E ( F m F 2 ) NM T @ F 2 ( X ) @ X ; (4.13) PAGE 42 34 where, the matrices M and N are as follo ws M = @ X k @ X = 266664 1 + P j ( w j ) x 0 ^ u kj P j ( w j ) y 0 ^ u kj P j ( w j ) z 0 ^ u kj P j ( w j ) x 0 ^ v k j 1 + P j ( w j ) y 0 ^ v k j P j ( w j ) z 0 ^ v k j P j ( w j ) x 0 ^ q k j P j ( w j ) y 0 ^ q k j 1 + P j ( w j ) z 0 ^ q k j 377775 (4.14) and N = @ X k +1 @ ~ T k +1 T ~ T k +1 =0 = 2666666666666664 P j w j ^ y k j P j w j ^ z k j 0 P j w j ^ x kj 0 P j w j ^ z k j 0 P j w j ^ x kj P j w j ^ y k j P j w k j 0 0 0 P j w k j 0 0 0 P j w k j 3777777777777775 : (4.15) Note that the gradien t v ector, ~ g ( ^ T k ), can not b e pre-computed. The estimate of the transformation at step k+1 is no w up dated b y the follo wing up date equation: ^ T k +1 = ~ T k +1 ^ T k = f ( ^ T k ; ~ T k +1 ) : (4.16) Where f is a function that dep ends on the t yp e of motion (rigid, ane, etc.) mo del used. Since the Hessian matrix can b e pre-computed, the n umerical cost p er iteration step of the mo died Newton sc heme is determined b y the calculation of the gradien t v ector. A signican t reduction of the calculation can b e ac hiev ed if at eac h step of the iteration, instead of expressing the gradien t v ector, the partial deriv ativ e of the transformed mo del image F m ( X ; ^ T k ; ~ T k +1 ) with resp ect to the motion v ector ~ T k +1 , w e no w express it in terms of the corresp onding deriv ativ es of the image F 2 ( X ) with resp ect to co ordinates f X g . Therefore, the partial deriv ativ es can b e determined in adv ance and within eac h iteration they only need PAGE 43 35 to b e m ultiplied with the dierence of the image F m ( X ; ^ T k ) and F 2 ( X ). This is precisely what w as done in the ab o v e gradien t equation 4.13. Our algorithm for computing the global 3D rigid transformation consists of the follo wing steps, 1. Precompute the Hessian at the optim um ~ H using equation (4.11). 2. A t iteration k , compute the gradien t v ector using equation (4.13). 3. Compute the transformation ~ T k +1 using equation (4.10). 4. Up date the motion parameter ^ T k +1 using the follo wing up date equations, 2666666666666664 ^ 1 ^ 2 ^ 3 ^ d 1 ^ d 2 ^ d 3 3777777777777775 k +1 = 2666666666666664 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 ~ r 11 ~ r 12 ~ r 13 0 0 0 ~ r 21 ~ r 22 ~ r 23 0 0 0 ~ r 31 ~ r 23 ~ r 33 3777777777777775 k 2666666666666664 ^ 1 ^ 2 ^ 3 ^ d 1 ^ d 2 ^ d 3 3777777777777775 k + 2666666666666664 ~ 1 ~ 2 ~ 3 ~ d 1 ~ d 2 ~ d 3 3777777777777775 k +1 : (4.17) where the (3,3) submatrix with en tries ( ~ r 11 ; :::; ~ r 33 ) is the 3D rotation matrix determined b y the three incremen tal angles ~ 1 ; ~ 2 ; ~ 3 . The mo died Newton metho d discussed ab o v e has man y adv an tages o v er the standard Newton sc heme for solving nonlinear equations. The ma jor computational adv an tage stems from the fact that the Hessian is pre-computed. In addition, the mo died Newton metho d has a broader range of con v ergence than the standard Newton sc heme, as w as sho wn in [55 ]. This allo ws the algorithm to allo w for large motions ev en with an initial guess that is far from the optim um. Finally , the mo died Newton sc hemes is quadratically con v ergen t [55 ]. 4.2.1 3D Ane Motion If the motion b et w een the t w o m ultimo dal images is an ane transformation, then the cost function has to b e c hanged to E sd ( T ) = E f [ F 1 ( X ) T AA T F 1 ( X + T ) PAGE 44 36 F 2 ( X ) T F 2 ( X )] 2 g . In the app endix, w e sho w that the Hessian of the cost function can b e written in a manner that allo ws a large part of it to b e pre-computed and th us accruing signican t sa vings in computational costs. The parameter up date algorithm is as in the rigid motion case except that the Hessian at the optim um can not b e fully pre-computed in this case. PAGE 45 CHAPTER 5 R OBUSTNESS ISSUE 5.1 Robust Estimation When the elds of view (F O Vs) for the image pairs b eing registered are signican tly dieren t, as illustrated in Figure 5.1, the ESD based metho d presen ted in previous c hapter will fail to giv e an accurate estimation. Similarly , the regular m utual information based metho ds also suer from the same problem. Curren tly the most p opular w a y to handle this problem is due to Studhlome et al. [6], who extended the regular MI to a F O V in v arian t similarit y measure { normalized m utual information (NMI). In this c hapter, w e prop ose another statistically robust measure In tegrated Square Error (ISE) or L 2 E to handle this problem. The comparison with NMI as w ell as other robust estimators is giv en in the exp erimen tal result c hapter. Figure 5.1: MR and CT image pair with dieren t elds of view. 37 PAGE 46 38 5.2 Related W ork One of the oldest robust metho ds used in computer vision eld is the Hough transform [56]. The basic idea is to map data in to the parameter space and seek the clusters to explain the data. Hough transform has b een successfully used in some feature detection tasks, including line tting, ellipse tting etc. Another old robust metho d is called regression diagnostics [57]. The idea is to iterativ ely detect p ossible outliers and reject them through the analysis of globally tted mo del. The third group, also the curren tly p opular estimators are called M-estimators. These estimators minimize the sum of a symmetric, p ositiv e-denite function ( r i ) of the residuals r i , with a unique minim um at r i = 0 (where residual is dened as the dierence b et w een the data and the t). Sev eral functions ha v e b een prop osed whic h reduce the inruence of large residual v alues of the estimated t. Hub er [58] emplo y ed the squared error for small residuals and the absolute error for large residuals. Andrews [59] used a squared sine function for small and a constan t for large residuals. In computer vision applications, one of the widely used M-estimators is Loren tzian function. Blac k et al. [35] dev elop ed a framew ork based on the Loren tzian function for robust estimation of optical ro w. The estimators using Sum of Squared Dierence(SSD) and the Loren tzian function ha v e a common form: min n X i =1 ( r i ) (5.1) where r i is the residual of the i th datum. In the registration con text, r i = I 1 ( x i ) I 2 ( t ( x i )). is usually a symmetric, p ositiv e-denite function and it is equal to ( r i ) 2 in the SSD case. When using the Loren tzian function as the estimator, it tak es the form of c 2 2 l og (1 + ( r i c ) 2 ), where c is a constan t. PAGE 47 39 5.3 Robustly Matc hing Lo cal F requency Represen tations T o matc h the lo cal frequency represen tations of the image pair, w e dev elop a statistically robust matc hing criteria based on minimization of the in tegral squared error (ISE) also kno w as the L 2 error or simply L 2 E b et w een a Gaussian mo del of the residual and the true densit y function of the residual. It w as sho wn in Scott [61] that minim um distance estimators, including the L 2 E , are inheren tly robust without requiring the need to sp ecify an y tuning parameters found in other robust metho ds. Figure 5.2 (from [61]) is an example sho wing the abilit y of L 2 E in rejecting outliers. The underlying histogram of the data is mixture of 2 Gaussian comp onen ts (0 : 8 N (0 ; 1) + 0 : 2 N (5 ; 1)). A practical prop ert y of L 2 E is that, it tends to select the largest comp onen t presen t in the data if the comp onen ts are assumed to b e w ell separated and if the n um b er of comp onen ts tted is smaller than required. The unique minimizer of the L 2 E mo del is sho wn in Figure 5.2 along with the MLE t for a comparison. If w e assume that the smaller of the t w o comp onen ts in the data w as attributed to the outliers in the data then, the L 2 E t has rejected the outliers successfully . There are t w o p ossible approac hes in emplo ying the lo cal frequency represen tations to handle the image registration problem: one is use the v ector-v alued phase gradien t represen tation and p erform the v ector matc hing; the other is use lo cal frequency magnitude as the underlying basis for scalar matc hing. In the follo wing section, w e will giv e the cost/ob jectiv e functions for these t w o dieren t sc hemes, and empirically sho w that the computation using the lo cal frequency magnitude tends to giv e more robust results under the L 2 E framew ork. PAGE 48 40 Figure 5.2: The parametric results of L 2 E and MLE in tting a single t w o parameter Normal distribution to data from a mixture of Gaussians. 5.3.1 V ector-v alued Lo cal F requency with Ane Motion Let the m ulti-mo dal image pairs to b e registered, dier b y an ane co ordinate transformation, then, the follo wing equation holds for the corresp onding v ectorv alued lo cal-frequency represen tations: A T F 1 ( AX + b ) = F 2 ( X ) + ( X ) (5.2) where F 1 ; F 2 are v ector elds corresp onding to the dominan t lo cal frequency estimates computed from the m ulti-mo dal (in our case, MR T1 and T2 w eigh ted) data sets resp ectiv ely; X = ( x; y ; z ); matrix A and v ector b corresp ond to the parameters of the unkno wn ane transformation that maps image 2 in to image 1, and ( X ) is the residual v ector at lo cation X . Note that an ane transformation distorts the lo cal frequency , e.g., an image that lo oks lik e cos( ! T X + c ), under an ane transformation of co ordinates will b e mapp ed to cos ( ! T AX + ! T b + c ), whic h PAGE 49 41 corresp onds to a lo cal frequency equal to A T ! . F or this reason, F 1 is prem ultiplied b y A T in Eq. (5.2). The residual v ector eld is assumed to b e comp osed of indep enden t, iden tically distributed, v ector{v alued random v ariables, whose distribution will b e mo deled b y a m ultiv ariate Gaussian with mean and co v ariance matrix 2 = 3 I ( I is the iden tit y matrix). Our goal is to minimize the L 2 E measure giv en b y min E ( ) = Z Z Z f g ( = ) h ( ) g 2 d 1 d 2 d 3 (5.3) where g ( : ) is the m ultiv ariate Gaussian function mo deling the densit y of the residual error, = [ ; ; A ; b ] b eing the v ector describing the Gaussian densit y parameters and , and the ane transformation parameters A ; b , and h is the true unkno wn densit y of the residual error term. Expansion of the in tegrand leads to t w o terms that are dep enden t on and a third term h 2 ( : ) indep enden t of , whic h can b e ignored from the minimization. The rst term in the expansion is R g 2 ( : ) and the second term is 2 E h g ( := ) i.e., the exp ectation of g ( : ) with resp ect to h , the true densit y of the residual. The rst term b eing a Gaussian can b e ev aluated in closed form, and w e can use the follo wing un biased estimator for the second term, 2 N N X i =1 g ( A T F 1 ( AX i + b ) F 2 ( X i ) = ) (5.4) f or i = 1 ; ::; N l attice points: Th us, the minimization using the L 2 E criterion is giv en b y min E ( ) = 1 8 p 3 3 2 N N X i =1 exp f j A T F 1 ( AX + b ) F 2 ( X ) j 2 2 2 g (5.5) T o estimate , w e solv e the minimization problem in equation (5.5) n umerically using a v arian t of the sequen tial quadratic programming approac h (see [62 ]). Eq. (5.5) ma y b e in terpreted as a robust estimation criterion for mo del (5.2). It PAGE 50 42 diers from the standard SSD approac h, in that the quadratic error terms are replaced b y robust p oten tials (in this case, in v erted Gaussians), so that large errors |due, for example, to dierences in F O V's | are not unduly o v erw eigh ted, but rather are treated as outliers and giv en a small w eigh t. The rst term of this error function p enalizes small v alues for the parameter , whic h denes the minimal error magnitude that is to b e considered as an outlier. Th us, this parameter need not b e set b y the user (unlik e in standard robust estimation), but rather it is automatically adjusted in a dynamic w a y b y the optimization metho d. This allo ws the metho d to a v oid b eing trapp ed in lo cal minima, ev en for large transformations, while retaining a high lev el of precision in the determination of the transformation parameters. 5.3.2 Lo cal F requency Magnitude with Rigid Motion If the lo cal frequency represen tations of the m ulti-mo dal image pairs (MR-MR pair acquired under dieren t imaging proto cols) dier b y a displacemen t whic h can b e parameterized as a rigid motion, then, the follo wing equation holds for the v ector-v alued lo cal-frequency represen tations: R T F 1 ( RX + b ) = F 2 ( X ) + ( X ) (5.6) Where R is a rotation matrix. Corresp ondingly the follo wing relation holds for the lo cal frequency magnitude: jj F 1 ( RX + b ) jj = jj F 2 ( X ) jj + 2 ( X ) (5.7) where the residual error eld 2 is assumed to ha v e a distribution mo deled b y a 1D Gaussian with mean and standard deviation of . jj F 1 ( ) jj and jj F 2 ( ) jj are the 3D lo cal frequency magnitude image represen tations computed from the PAGE 51 43 m ultimo dal data sets resp ectiv ely , The L 2 E measure for this sc heme is giv en b y: min T ;; E ( T ; ) = 1 2 p 2 N N X i =1 exp f ( jj F 1 ( RX + b ) jj jj F 2 ( X ) jj ) 2 2 2 g (5.8) T o estimate the parameterized transformation, w e solv e the minimization problem in equation (5.8) n umerically . It is p ossible to use a v ariet y of n umerical metho ds to ac hiev e this n umerical optimization starting with simple gradien t descen t and mo ving to w ard more sophisticated metho ds suc h as the non-linear conjugate gradien t, the Quasi-Newton metho ds or the sequen tial quadratic programming { SQP { (if the p enalt y function formed from the Lagrangian (see b elo w) is assumed to b e non-dieren tiable) (see [62] for details). Of course, the gradien t descen t and the non-linear conjugate gradien t algorithms m ust b e preconditioned for impro v ed p erformance/con v ergence. In this pap er, w e presen t exp erimen tal results that w ere obtained using a v arian t of the sequen tial quadratic programming sc heme. 5.3.3 Numerical Implemen tation The Numerical implemen tation w as ac hiev ed using a constrained minimizer of the equation 5.8 and 5.5. The constrained used w as that 0 < < L , with L b eing a p ositiv e constan t. W e will briery describ e the n umerical sc heme used and refer the in terested reader to [62] for details. Let g j ( ) ; j = 1 ; ::; r b e the set of inequalit y constrain ts in the minimization of E ( ). The exact p enalt y function E + cP is rst formed where, c > 0, P ( ) = max f g 0 ( ) ; g 1 ( ) ; :::; g r ( ) g . Similar to the gradien t pro jection metho ds, in the SQP , the descen t direction is computed b y solving a quadratic programming subproblem. The iterativ e descen t is giv en b y: k +1 = k + k d k (5.9) PAGE 52 44 Where, k is a non-negativ e scalar step size and d k is the descen t direction obtained b y solving a quadratic program in ( d ; ), min ( r E ( k )) t d + (1 = 2) d t H k d + c (5.10) sub ject to g j ( k ) + ( r g j ( k )) t d < ; j = 0 ; ::; r : (5.11) Where, the Hessian matrix H k is required to b e a symmetric p ositiv e denite matrix. The step size can b e c hosen using a v ariet y of sc hemes of v arying sophistication ranging from the basic minimization rule, limited minimization rule of the Armijo rule. See [62 ] for details. One can use the MA TLAB optimization to olb o x whic h con tains v arious t yp es of optimizers including the ab o v e describ ed metho d. Curren tly , w e are w orking on a public domain C language implemen tation of a feasible sequen tial quadratic programming algorithm from Univ ersit y of Maryland (The FSQP home page in the Institute for Systems Researc h in the Univ ersit y of Maryland. Av ailable online at h ttp://www.isr.umd.edu/Labs/CA CSE/FSQP/) and exp ecting to get similar results with a shorter execution time. PAGE 53 CHAPTER 6 NONRIGID MOTION T o handle large nonrigid deformations (those that can not b e handled b y the metho ds in v olving minimization of equation 4.2), w e use the lev el-set based registration metho d describ ed in [7]. W e assume that the nonrigid motion is appro ximated b y piecewise rigid motions. Th us, since the lo cal frequency magnitude represen tation of the images are in v arian t to rigid motions, w e can use these represen tations in computing the nonrigid displacemen t elds. If the deformation b et w een the data sets is large, one ma y ha v e to reinitialize the lo cal-frequency represen tations after ev ery few iterations. The main idea of the lev el-set based nonrigid registration metho d is to ev olv e the lev el-sets of one represen tation in to the other. The lo cal-frequency magnitude represen tations of the m ulti-mo dal data set forms the input to the lev el-set based image registration metho d and the go v erning partial dieren tial equation (PDE) for the ev olution of the lev el-sets is giv en b y , F t ( X ; t ) = ( F 2 ( X ) F ( X ; t )) kr F ( X ; t ) k w ith F ( X ; 0) = F 1 ( X ) (6.1) Where, F 1 and F 2 are the lo cal-frequency image represen tations of the m ulti-mo dal image pair resp ectiv ely and X = ( x; y ; z ). The v ariable t here is used to indicate the time parameter in the ev olution. A pro of of existence and uniqueness of the solution to a PDE of this form where F 1 and F 2 are assumed Lipsc hitz con tin uous functions is giv en in [7]. Note that the equation 6.1 basically morphs the source represen tation in to the target but do es not compute an y underlying geometric (co ordinate) transformation that migh t b e prev alen t b et w een the t w o data sets. W e determine this geometric transformation b et w een the t w o data sets via the use of 45 PAGE 54 46 the follo wing PDE: ~ V t = ( F 2 ( X ) F ( ~ V ( X ))) r F ( ~ V ( X )) kr F ( ~ V ( X )) k w ith ~ V ( X ; 0) = ~ 0 (6.2) where ~ V = ( u; v ; w ) T is the displacemen t v ector at X and the op eration ~ V ( X ) = ( x + u; y + v ; z + w ). Equation (6.2) ma y b e used to ev olv e feature maps and ma y also b e applied to w arp a giv en segmen tation on to an unkno wn data set. Since the gradien t is highly sensitiv e to the noise, I ( X , t ) is con v olv ed with the Gaussian k ernel b efore taking the gradien t. This leads to the follo wing mo dication in the previous equation ! V t = [ I 2 ( X ) I 1 ( ! V ( X ))] r ( G I 1 ( ! V ( X )))) k r ( G I 1 ( ! V ( X ))) k (6.3) w ith ! V ( X ; 0) = ! 0 T o mak e the computation of the lev el-set motion stable in nature when, ( G I 1 ( ! V ( X ))) is near to zero, a stabilizing factor is added to previous equation ! V t = [ I 2 ( X ) I 1 ( ! V ( X ))] r ( G I 1 ( ! V ( X ))) k r ( G I 1 ( ! V ( X ))) k + (6.4) w ith ! V ( X ; 0) = ! 0 where is a small p ositiv e n um b er. F urther, a smo othness constrain t is added to the lev el-set motion mo del whic h mo dies the equation (6 : 4) as follo ws ! V t = [ I 2 ( X ) I 1 ( ! V ( X ))] r ( G I 1 ( ! V ( X ))) k r ( G I 1 ( ! V ( X ))) k + + 0BBBB@ u v w 1CCCCA (6.5) w ith ! V ( X ; 0) = ! 0 where denotes the Laplacian op erator. In the implemen tation, the smo othness constrain t is added externally b y con v olving the motion eld with a Gaussian lter after eac h iteration. PAGE 55 47 6.1 Numerical Implemen tation T o solv e (6.5), V em uri et al. [7] ha v e suggested the use of the minmo d nite dierence sc heme b y Kimmel et al. [63] to ac hiev e the lev el-set based optical ro w. The minmo d function is dened as follo ws m ( x; y ) = 8><>: sig n ( x )min ( j x j ; j y j ) if xy > 0 0 if xy 0 (6.6) Dening C ij k = G ( I 1 ) i u;j v ;k w , w e get ( C x ) ij k = m ( D + x C ij k ; D x C ij k ) ( C y ) ij k = m ( D + y C ij k ; D y C ij k ) (6.7) ( C z ) ij k = m ( D + z C ij k ; D z C ij k ) where D + x , D x , D + y , D y , D + z and D z are dened as follo ws D + x E ij k = E i +1 ;j;k E i;j;k x ; D x E ij k = E i;j;k E i 1 ;j;k x D + y E ij k = E i;j +1 ;k E i;j;k y ; D y E ij k = E i;j;k E i;j 1 ;k y (6.8) D + z E ij k = E i;j;k +1 E i;j;k z ; D z E ij k = E i;j;k E i;j;k 1 z PAGE 56 48 Using these denitions and a forw ard dierence in time, (6.5) can b e rewritten in a discretized form: 0BBBB@ u s +1 ij k v s +1 ij k w s +1 ij k 1CCCCA = 0BBBB@ u sij k v s ij k w s ij k 1CCCCA + 0BBBB@ u sij k v s ij k w s ij k 1CCCCA + (6.9) t ( I 2 ) ij k ( I 1 ) i u sij k ;j v s ij k ;k w s ij k q m 2 ( D + x C ij k ; D x C ij k ) + m 2 ( D + y C ij k ; D y C ij k ) + m 2 ( D + z C ij k ; D z C ij k ) + 0BBBB@ m ( D x C ij k ; D x C ij k ) m ( D + y C ij k ; D y C ij k ) m ( D + z C ij k ; D z C ij k ) 1CCCCA w ith 0BBBB@ u 0ij k v 0 ij k w 0 ij k 1CCCCA = ! 0 This nite dierence sc heme is kno wn to preserv e lo cal min/max for the equation, ensuring stable solutions for the lev el set motion mo del. 6.2 Adaptiv e Time Step t T o stabilize the n umerical in tegration pro cess, Y e in [64] prop ose the use of adaptiv e time step sc heme. Let I 2 ( X ) b e the target image, I ( X ; t ) b e the source image at time t and E ( X ) = G I ( X ). The time step t can b e calculated as follo ws t = 8><>: 1 max fj H u j + j H v j + j H w jg if max( j H u j + j H v j + j H w j ) 6 = 0 k otherwise (6.10) PAGE 57 49 where u = E x ; v = E y ; w = E z H u = u p u 2 + v 2 + w 2 + ( I 2 ( X ) I ( X ; t )) H v = v p u 2 + v 2 + w 2 + ( I 2 ( X ) I ( X ; t )) (6.11) H w = w p u 2 + v 2 + w 2 + ( I 2 ( X ) I ( X ; t )) k = small p ositiv e n um b er (6.12) and max( j H u j + j H v j + j H w j ) is computed in eac h iteration using the maxim um absolute v alue of H u ; H v and H w dened on the 3D lattice ( H u ; H v and H w is computed for eac h grid p oin t). 6.3 Incremen tal Nonrigid Up dating In a Lagrangian grid implemen tation, in eac h iteration, the grid is transformed b y the curren t transformation and then resampled in the transformed space and the displacemen t eld is reinitialized to zero. Strictly sp eaking, this need not b e done at ev ery iteration but in practise, it suces to p erform this only after ev ery few iterations. F or more details on the up wind dierencing sc heme used, w e refer the reader to [7 ]. PAGE 58 CHAPTER 7 EXPERIMENT RESUL TS In this c hapter, w e presen t the the results of applying our ESD as w ell as L 2 E measure to a large n um b er of data sets and demonstrate the sup eriorit y of our approac hes o v er the p opular m utual information based metho ds. The exp erimen ts in section 7.1 sho w the impro v emen t made b y using our ESD measure in shortening CPU execution time without losing accuracy . The sp eed-up factor, as eviden t from the results, is around 13. Exp erimen ts in section 7.2 sho w a b etter p erformance of our L 2 E sc heme o v er normalized MI metho d and the adv an tages reside not only in the robustness but also the computational eciency . 7.1 Exp erimen ts with ESD-based Measure In this section, w e presen t four sets of exp erimen ts on in ter-mo dalit y image registration based on our ESD measure. The rst t w o sets demonstrate the results for syn thesized and real rigid registrations and the third and fourth set demonstrate nonrigid registrations on syn thesized 2D images and real 3D v olumes resp ectiv ely . In the rigid registration exp erimen ts, for syn thesized and the real misalignmen ts, w e compare the computed registrations with the ground truth registrations whic h in the real mis-alignmen t case are obtained from a man ual alignmen t pro cess b y an "exp ert". In the nonrigid registration exp erimen t, v alidation is done only qualitativ ely via visual insp ection. 7.1.1 Syn thesized Motion: Estimation Results W e no w presen t results of using the lo cal-frequency magnitude based image represen tation of the t w o input m ulti-mo dal image data sets, to our parameterized 50 PAGE 59 51 motion estimator for the rigid registration case and to the lev el-set based motion estimator for the nonrigid case. In ter mo dalit y syn thetic rigid motion b et w een MR and CT: T able 7.1 sho ws the estimated rigid (rotation + translation) motion b et w een an MR and CT brain scan acquired from a single sub ject. The left column in the table sho ws the kno wn rigid transformation applied to the source images. The rst three terms are rotation angles ab out the x,y and the z axes resp ectiv ely . Other three are translations along the (x,y ,z) directions resp ectiv ely . The middle column sho ws the corresp onding estimates of these parameters. RMS errors for the rotation and translation parameters are depicted in the righ t column. The accuracy of the ac hiev ed results as eviden t is quite high in all these exp erimen ts. The MR and CT image sizes in this exp erimen t w ere b oth (512 ; 512 ; 128). The n um b er of Gab or lters used for computing the lo cal frequency represen tation w as 7 and the a v erage CPU time for this pro cess w as 106 sec on an R10000 single CPU of an SGI On yx. In the estimation of the global rigid motion, the n um b er of con trol p oin ts w e used w as 8 and the n um b er of lev els in the p yramid represen tations { for the m ulti-resolution framew ork based implemen tation { w as set to 3. The a v erage CPU time for the motion estimation step w as 97sec on the same mac hine. As seen from the table, the estimated rigid motion is v ery accurate. T able 7.1: 3D rigid motion estimates for MR-CT data with kno wn motion b et w een the data. The n um b ers in the paren thesis corresp ond to ( x , y , z , t x , t y , t z ) T rue motion(degree, v o xel) Estimated motion (degree, v o xel) RMSE (0 ; 0 ; 0 ; 10 ; 8 ; 3) (0 ; 0 ; 0 ; 10 : 004 ; 8 : 045 ; 2 : 94 ) (0 : 00 ; 0 : 434) (5 : 73 ; 5 : 73 ; 5 : 73 ; 4 ; 5 ; 6) (5 : 67 ; 5 : 65 ; 5 : 67 ; 3 : 90 ; 5 : 06 ; 5 : 95) (0 : 067 ; 0 : 073) (5 : 73 ; 11 : 46 ; 17 : 20 ; 9 ; 6 ; 3) (5 : 61 ; 11 : 92 ; 17 : 13 ; 8 : 89 ; 5 : 972 ; 3 : 127) (0 : 277 ; 0 : 098) (11 : 46 ; 11 : 46 ; 11 : 46 ; 7 ; 8 ; 9) (11 : 28 ; 11 : 28 ; 11 : 46 ; 6 : 80 ; 8 : 11 ; 9 : 04) (0 : 147 ; 0 : 114) (22 : 92 ; 0 : 0 ; 0 : 0 ; 0 ; 0 ; 0) (22 : 46 ; 0 : 00 ; 0 : 00 ; 0 : 34 ; 0 : 5 4 ; 0 : 10 ) (0 : 266 ; 0 : 387) (0 : 0 ; 0 : 0 ; 17 : 19 ; 0 ; 0 ; 0) ( 0 : 00 ; 0 : 00 ; 17 : 01 ; 0 : 05 ; 0 : 13 ; 0 : 209) (0 : 104 ; 0 : 373) ( 11 : 46 ; 5 : 73 ; 11 : 46 ; 14 ; 4 ; 3) ( 11 : 34 ; 5 : 67 ; 11 : 40 ; 13 : 82 ; 3 : 71 ; 5 : 31) (0 : 085 ; 1 : 348) PAGE 60 52 7.1.2 Real 3D Rigid Motion: Estimation Results In this section, w e presen t the p erformance of our Exp ectation of squared dierence based metho d (ESD) on data con taining real rigid mis-alignmen ts. F or the purp ose of comparison, w e implemen ted the w ell kno wn Mutual Information based metho d (MI) presen ted in Collignon et al. [37 ]. The results from eac h algorithm are compared to ground truth registrations obtained man ually b y an \exp ert", whose man ual registrations are curren tly used in clinical practise at the Univ ersit y of Florida Brain Institute. As will b e seen from the results describ ed b elo w, the k ey adv an tage of our metho d o v er the widely used MI-based metho ds is that, w e are able to obtain comparable accuracy within m uc h shorter CPU time. W e tested the our algorithm and the MI-based metho d on MR-CT data from ten dieren t sub jects. The MR-CT pairs w ere miss-aligned due to motion of the sub ject. The CT image w as of size (512,512,120) while the MR image size w as (512,512,142) and the v o xel dimensions w ere (0 : 46 ; 0 : 46 ; 1 : 5) and (0 : 68 ; 0 : 68 ; 1 : 05) for CT and MR resp ectiv ely . T able 7.2 summarizes the results of the comparison on the ten data sets. The ro w lab eled \T rue" in the table sho ws the ground truth rotation and translation parameters (as assessed b y the lo cal exp ert) of eac h data set. W e should men tion that the lo cal-exp ert's registrations are curren tly in clinical use on a daily basis for neuro-surgery . Both rotation and translation parameters are in (x,y ,z) order. The ro w lab eled, \ESD Err." depicts the absolute error b et w een the estimated parameters and the ground truth using our metho d, while \MI Err." indicates the corresp onding errors for the MI metho d. As eviden t, b oth the metho ds ha v e ac hiev ed a rather high accuracy in the registrations. F or most of the cases, the a v erage error of rotation angle is less than 0.6 degree and the translation is less than 0.6 millimeter. F rom the table, w e can see that the translation errors from our metho d are sligh tly w orse than those from the MI metho d. But, b y p erforming a PAGE 61 53 studen t-t test on the dierence b et w een the motion estimates from our metho d and the MI metho d, w e found the probabilit y of the dierence in error estimates from the t w o metho ds b eing upp er b ounded b y , (0 : 32 mm; 0 : 36 mm; 0 : 31 mm ) is 0 : 9 and (0 : 134 ; 0 : 188 ; 0 : 216) is 0 : 9 for the translation and rotation parameters resp ectiv ely . Th us, in most cases, the dierence in motion estimates obtained from the t w o metho ds is not signican t. The righ tmost column of the table sho ws the CPU time for b oth metho ds to register these large data sets on a single R10000 pro cessor of the SGI-On yx. As eviden t, our metho d is computationally more ecien t than the normalized MI based metho ds. 7.1.3 Nonrigid Motion: Syn thetic and Real Data Example Figure 7.1 is an example sho wing the abilit y of our metho d to handle large nonrigid deformation. Figures 7.1(a) and 7.1(b) depict the source and target 2D images resp ectiv ely used in this exp erimen t. As is eviden t, the in tensit y proles of the t w o images are v ery distinct from one another. The lo cal frequency magnitude is used as the underlying represen tation to compute the nonrigid deformation eld b et w een the t w o images, whic h is depicted in 7.1(c). The result of applying this ro w eld to the source image is sho wn in 7.1(d) As is visually eviden t, the deformed source matc hes the geometry of the target quite accurately . Note that w e do not deal with transforming the in tensit y of the source to that of the target. As men tioned earlier, w e used the Lagrangian framew ork for implemen tation of the go v erning PDE for the nonrigid motion. Strictly sp eaking, the grid should b e reinitialized after ev ery iteration but this pro v es to b e computationally exp ensiv e and not necessary since the deformation b et w een t w o iterations is quite small and in practise, one can aord to reinitialize after sev eral iterations that w ould sucien tly deform the grid. Re-initialization here refers to a zeroing of the deformation eld after resampling of the transformed initial grid. PAGE 62 54 (a) (b) (c) (d) Figure 7.1: An example of nonrigid motion. (a) Source image, (b) target image, (c) deformed mesh and (d) transformed source image. W e no w presen t a real data example in v olving the a pair of T1 and T2 w eigh ted MR images from dieren t sub jects. The estimated v ector eld is applied to the one of the t w o giv en images lab eled as the source image to obtain a transformed image whic h is \closer" to the target image than the original source image. The v ector eld c haracterizing the transformation is computed from equation (6.2. W e sup erimp ose the edge map of the target on to this transformed image to visually depict the accuracy of our registration algorithm. In this exp erimen t, w e PAGE 63 55 (a) (b) (c) (d) (e) (f ) (g) (h) Figure 7.2: Non-Rigid Registration of an MR-MR brain scan, example 1. The rst ro w: (a) source image, (b) target image, (c) globally transformed source image with sup erimp osed target edge map, (d) nonrigidly transformed source with sup erimp osed target edge map. The second ro w is the zo om-in v ersion of the (c) and (d): (e) left upp er and (f ) cen tral parts of (c); (g) left upp er (h) cen tral parts of (d). rst estimated the global ane transform b et w een the pair of input images b eing registered using the rigid/ane motion estimation algorithm describ ed in section 4.1. The lo cal frequency maps of these globally transformed source and target images are then used as input to the nonrigid lev el-set based registration algorithm describ ed in the c hapter 6. Figure 7.2(a) and 7.2(b) depict the source and target images to b e registered. Figure 7.2(c) depicts the rigidly transformed source and the nonrigidly transformed source image is sho wn in Fig. 7.2(d). 7.2(c) and 7.2(d) ha v e b een sup erimp osed with the edge map of target image to visually demonstrate the accuracy of the registration. Figures 7.2(e) and 7.2(f ) depict the residual mis-alignmen t after the application of the global ane transform to the source using a close up view at t w o dieren t lo cations in the source. Figures 7.2(g) and 7.2(h) depict the close up view of the alignmen t ac hiev ed after the application of a PAGE 64 56 nonrigid deformation to the rigidly transformed source image at the same lo cations as in 7.2(e) and 7.2(f ). As is eviden t, the mis-alignmen t has b een dramatically reduced after the application of the nonrigid deformation. 7.2 Exp erimen ts with the L 2 E Measure In this section, w e presen t four sets of exp erimen ts based on our L 2 E measure. The rst set constitutes of a 2D example to depict the adv an tage of L 2 E o v er, the Sum of Squared Dierence (SSD), the Loren tzian estimator and the normalized MI. The second set con tains exp erimen ts with syn thetic 2D data to assess the adv an tages of using magnitude vs. v ector-v alued frequency represen tations as a basis in the L 2 E based matc hing tec hnique describ ed earlier. The third set of exp erimen ts w as conducted with 3D MR T1 and T2 w eigh ted data sets obtained from the Mon treal Neurological Institute (MNI) database. The data sets w ere articially misaligned b y kno wn rigid transformations and our algorithm as w ell as the normalized MI w ere used to estimate the transformation and their p erformances w ere compared. The fourth set con tains exp erimen ts using 3D CTMR data sets obtained from the RadioSurgery Lab at the Univ ersit y of Florida Brain Institute. The comparison b et w een our algorithm and the normalized MI w as conducted in the same w a y as in the third exp erimen t men tioned earlier. 7.2.1 Robust Prop ert y of L 2 E Measure In this section, w e demonstrate the robustness prop ert y of L 2 E and hence justifying the use of the L 2 E measure in the registration problem. This is ac hiev ed b y sho wing the sup eriorit y of L 2 E o v er the SSD and a widely used M-estimator: the Loren tzian function, and then the normalized MI measure. As w e in tro duced in c hapter 5, the estimators using SSD and the Loren tzian function ha v e a common form: min n X i =1 ( r i ) (7.1) PAGE 65 57 (a) (b) Figure 7.3: (a) The source image, a MR slice with half image b eing cut. (b) The target image, obtained b y rotating the original source image. where r i is the residual of the i th datum. In the registration con text, r i = I 1 ( x i ) I 2 ( t ( x i )). is usually a symmetric, p ositiv e-denite function and it is equal to ( r i ) 2 in SSD case. When using the Loren tzian function as the estimator, it tak es the form of c 2 2 l og (1 + ( r i c ) 2 ), where c is a constan t. In order to compare the robustness prop ert y of L 2 E v ersus the estimators based on SSD and the Loren tzian function, w e designed a series of exp erimen ts as follo ws: with a 2D MR slice as the source image, the target image is obtained b y applying an ane transformation to the source image. Outliers are created b y setting the in tensit y v alues of half of the source image to zero. The source and target image pair are sho wn in Figure 7.3. With 15 randomly generated ane transformations, w e applied L 2 E, together with SSD and the Loren tzian function to estimate those motion parameters. T able 1 sho ws the statistics of errors resulting from the 3 dieren t metho ds. In eac h cell, the left 2 b y 2 matrix is the ane matrix, while the righ t v ector sho ws the translations in x and y directions. As exp ected, SSD, b eing sensitiv e to outliers, pro duced v ery unreliable results. The second ro w and the third ro w of table 7.3 sho w the mean and standard deviation of the ane errors estimated with Loren tzian estimator and L 2 E, resp ectiv ely . It is apparen t from this exp erimen t that L 2 E is b etter than the Loren tzian; further, out of the 15 trials, the Loren tzian PAGE 66 58 estimator failed 5 times while L 2 E failed only once (\failed" here means that the results had unacceptably large errors). If w e only coun t the cases whic h ga v e reasonable results, as sho wn in the 4th and 5th ro ws, L 2 E and the Loren tzian ha v e comparable p erformances, b oth b eing v ery accurate. So w e can dra w the follo wing conclusion: L 2 E and the Loren tzian estimator are comparable in terms of accuracy , but L 2 E is less prone to lo cal minima en trapmen t. The reason for the increased robustness of L 2 E is b ecause it has built in an adaptiv e tuning of the parameter: at the b eginning of the optimization, is large, so that the global searc h is more eectiv e; as the iterations pro ceed, is automatically adjusted, so that at the end one has a small whic h p ermits a precise appro ximation. T o compare L 2 E and the normalized MI measure, w e used the same image data as in the ab o v e exp erimen ts. Due to the clipping of the source image, w e migh t in tro duce a large 'syn thetic ap erture' in to the image. This padding of large regions of pixels with iden tical v alues (zero es) can pro duce v ery abnormal join t distributions, and therefore, aect the accuracy of the results for the normalized MI metho d. W e explored the eects of the clipping on L 2 E and normalized MI b y carrying out exp erimen ts as follo ws: randomly generate 15 rigid transformations and apply them on the target image, and then estimate the motion parameters based on t w o dieren t sc hemes: one is to carry out the computation based only on the o v erlapping region, and the other is to base it on the complete, zero-padded image. T able 7.4 depicts the mean error obtained from L 2 E and normalized MI from these exp erimen ts. In eac h cell, the rst v alue sho ws the rotation error and the latter t w o sho w the translation errors in x and y directions. As can b e observ ed, normalized MI gets go o d results only when the o v erlapping-only sc heme is used, while L 2 E yields v ery accurate results with b oth the sc hemes. This not only illustrates the robustness of L 2 E, but implies a computational adv an tage as w ell; PAGE 67 59 using L 2 E, the task of nding the o v erlapping region can b e omitted, and this computation migh t b e computationally in tensiv e, particularly if the clipping shap e is v ery complicated, e.g., a p olygon, or a p olyhedron. 7.2.2 Comparison of Magnitude and V ector-v alued Lo cal F requency Sc hemes As men tioned in the earlier c hapters, there are t w o p ossible represen tations that can b e used when dealing with lo cal frequency to tac kle the registration problem: one is to directly matc h the lo cal frequency v ectors; the other is to use the lo cal frequency magnitude as the underlying represen tation. W e designed 2 exp erimen ts to explore the prop erties of b oth represen tation sc hemes. This comparison can p ossibly b e generalized to other v ector/scalar represen tations. As eviden t from these t w o exp erimen ts, eac h represen tation has its o wn adv an tage o v er the other under sp ecic circumstances. The rst exp erimen t is designed to sho w the sup eriorit y of v ector-v alued lo cal frequency o v er frequency magnitude. Consider the t w o ob jects sho wn in gure 4 The circular region of gure 4.a has a frequency of 0.1 along y direction, and 0 along x direction. Figure 4.b is a rotated v ersion of gure 4.a where the rotation angle is 30 an ti-clo c kwise. As in tro duced in ab o v e section, the circular region of gure 4.b has a lo cal frequency of cos (30 ) sin (30 ) sin (30 ) cos (30 ) ( 0 : 0 0 : 1 ) = ( 0 : 05 0 : 0866 ) except at the b oundary . Based on the v ector-v alued lo cal frequency , our estimation of the rotation angle is (29 : 997 ), whic h is v ery close to the true v alue. Ho w ev er, b ecause g 4.b is just the rotated v ersion of g 4.a , they ha v e the same lo cal frequency magnitude, and hence, it is imp ossible to estimate the rotation information from the transformed images. In the second exp erimen t, w e examine the relativ e accuracy when the in tensit y of one of the images is distorted b y a smo oth m ultiplicativ e bias: giv en the target image I 1 , the source image I 2 is generated b y applying a kno wn co ordinate transformation and then the target image is altered using I 1 ( y ; x ) = ( x; y ) I 1 ( x; y ). In this exp erimen t, the function ( x; y ) is assumed to b e a linear function ( x; y ) = PAGE 68 60 (a) (b) Figure 7.4: (a) A round region with a lo cal frequency of (0.1, 1). (b) obtained b y rotating (a) 30 an ti-clo c kwise. (a) and (b) ha v e the same lo cal frequency magnitudes. ( x + y ) = ( w + h ), while w and h are the width and heigh t of the image resp ectiv ely . T able 7.5 sho ws the estimation results obtained using t w o metho ds with dieren t v alues of . As it is eviden t, the ane parameter estimation using v ector-v alued lo cal frequency is quite sensitiv e to the in tensit y c hanges, while lo cal frequency magnitude yields more accurate estimates ev en when is quite large. F rom the ab o v e exp erimen ts and discussion, w e can see that there is a tradeo b et w een the usage of v ector-v alued and lo cal frequency magnitude. In less noisy en vironmen ts, including the follo wing real 3D m ulti-mo dal data sets, w e used the v ector-v alued lo cal frequency as the underlying represen tation. 7.2.3 Comparison with NMI using Syn thetic MRI Data In this section, w e demonstrate the algorithm p erformance for in ter-mo dalit y rigid registrations. All the examples con tain syn thesized misalignmen ts applied to MR data sets from the brain w eb site at the Mon treal Neurological Institute [65]. for comparison purp oses, w e ha v e implemen ted the SSD algorithm for matc hing the lo cal frequency represen tations as w ell as the MI and normalized MI for matc hing the original image pairs. W e will ho w ev er presen t only the results of comparison with the normalized MI whic h p erforms b etter than the SSD and the ordinary MI sc hemes. F or the normalized MI as w ell as for our metho d, w e compare the PAGE 69 61 computed registrations with the ground truth registration whic h is kno wn for these data sets as they w ere articially misaligned b y kno wn transformations. As will b e seen from the results describ ed b elo w, the k ey adv an tage of our metho d o v er the state-of-the-art in mi-based metho ds namely , the normalized MI-based metho d is that w e are able to handle m uc h larger non-o v erlapping areas b et w een the t w o data sets b eing matc hed. Figure 5 depicts an example of the t w o MR data sets obtained using the T1 and T2 w eigh ting resp ectiv ely . The rst ro w of the gure depicts corresp onding slices from the T1 w eigh ted MR brain scan in the sagittal, coronal and axial views. Second ro w sho ws a similar arrangemen t of images for the T2 w eigh ted MR scan. Note that the eld of view (F O V) for the data sets are signican tly non-o v erlapping. The non o v erlap w as sim ulated b y simply cutting appro ximately 50% of the T2 w eigh ted image (target image). Third ro w depicts the transformed source (T1 w eigh ted) image along with an edge map of the target (T2 w eigh ted) sup erimp osed on the transformed source (T1 w eigh ted) image. As eviden t, the registration is visually quite accurate. The exp erimen t design w as as follo ws: in ev ery exp erimen t, w e applied a kno wn transformation to the T2 MR image, then w e use normalized MI and our metho d to estimate the transformation. 3 sets of 15 randomized transformations ha v e b een used. They are normally distributed around the v alues of (10 ; 5 mm ), (15 ; 10 mm ) and (30 ; 15 mm ) resp ectiv ely , with standard deviations of 5 and 2 mm for rotation and translation resp ectiv ely . The MR images are of size (256, 256, 200). F or a quan titativ e assessmen t of the accuracy of the registration, w e computed the mean and standard deviation of the errors for the 3 rotation angles and the 3 translations. It should b e noted that in all the three sets of exp erimen ts, our L 2 E-based metho d has yielded sup erior p erformance. PAGE 70 62 The n umerical sc heme w e used to implemen t normalized MI is based the P o w ell's metho d, whic h is used in Collignon's w ork [37]. T o compute the discrete histogram, w e used 256 in tensit y bins for eac h mo dalit y . When the image is transformed, tri-linear in terp olation is used to a v erage neigh b oring v o xels. In order to sp eed up the con v ergence and a v oid b eing trapp ed in lo cal minima, a m ulti-resolution sc heme is used in our implemen tation, and for these particular exp erimen ts, w e used a 3 lev el p yramid. Regarding the L 2 E implemen tation, one needs to pa y some atten tion on the direction (sign) of lo cal frequency . A v ector frequency ( v x ; v y ; v z ) lo oks, lo cally , exactly the same as another v ector ( v x ; v y ; v z ), and this am biguit y ma y cause problems in the registration algorithm. T o eliminate it, w e mapp ed the source lo cal frequency F 1 and the transformed target lo cal frequency A T F 2 ( AX + b ) in to the same half space, b y setting F = F , if F x < 0, where F x is the x-comp onen t of the lo cal frequency . The n um b er of Gab or lters used in these exp erimen ts w as 7. T ables 7.6, 7.7 and 7.8 summarize the results of applying our L 2 E algorithm and the normalized MI algorithm to the 3 sets of misaligned MR T1 and T2 w eigh ted image pairs. The tables depict, the maximal v alue, mean and standard deviation of the error v ectors from L 2 E estimation(left column) and normalized MI metho d (righ t column). Within ev ery \cell", the upp er ro w sho ws the 3 error v alues for rotation angles, and the lo w er ro w, the errors for the translation parameters. As eviden t, for most of the estimated parameters, w e obtained smaller maximal errors with L 2 E in comparison to normalized MI. In addition, w e obtained smaller mean in rotation errors along with comparable translation errors. 7.2.4 Comparison with 3D CT-MR Data In this section, w e demonstrate the algorithm p erformance for in ter-mo dalit y rigid registrations using CT-MR data sets. The original aligned CT-MR image PAGE 71 63 pair w ere obtained from the RadioSurgery Lab at the Univ ersit y of Florida Brain Institute. The alignmen t w as ac hiev ed man ually b y an exp ert whose alignmen ts are used in daily clinical practise. The exp erimen ts w ere conducted exactly the w a y as sho wn in T1/T2 examples: w e articially misaligned the CT-MR images b y kno wn rigid transformations, whic h are normally distributed around the v alues of (10 ; 5 mm ), (15 ; 10 mm ) and (30 ; 15 mm ) resp ectiv ely , with standard deviations of 5 and 2 mm for rotation and translation resp ectiv ely . F rom table 7.9, 7.10 and 7.11, one can see that o wn conclusion still holds, i.e, for most of the estimated parameters, our metho d pro duced smaller maximal error, smaller mean in rotation errors and comparable translation errors. PAGE 72 64 T able 7.2: 3D rigid motion estimates for ten MR-CT data sets. Set Item Rotation (degree) T ranslation (mm) CPU(mins.) 1 T rue (5 : 455 1 : 146 14 : 003) (3 : 822 9 : 254 3 : 1094) ESD Err. (0 : 524 0 : 027 0 : 509) (0 : 229 0 : 729 0 : 525) 3 : 44 MI Err. (0 : 714 0 : 755 0 : 004) (0 : 013 0 : 223 0 : 441) 46 : 38 2 T rue (1 : 765 2 : 023 12 : 284) (6 : 661 2 : 340 6 : 280) ESD Err. (0 : 249 0 : 146 0 : 1965) (0 : 335 0 : 259 0 : 133) 3 : 04 MI Err. (0 : 453 0 : 025 0 : 109) (0 : 401 0 : 002 0 : 130) 52 : 10 3 T rue ( 5 : 099 7 : 357 18 : 472) (9 : 790 0 : 901 0 : 228) ESD Err. (0 : 222 0 : 606 0 : 825) (0 : 715 0 : 183 0 : 689) 4 : 49 MI Err. (0 : 137 0 : 533 0 : 218) (0 : 127 0 : 297 0 : 012) 53 : 32 4 T rue ( 5 : 581 2 : 865 21 : 675) (10 : 561 4 : 306 18 : 874) ESD Err. (0 : 023 0 : 175 0 : 145) (0 : 759 0 : 175 0 : 017) 3 : 74 MI Err. (0 : 539 0 : 254 1 : 025) (0 : 121 0 : 253 0 : 415) 56 : 23 5 T rue ( 2 : 498 7 : 540 23 : 737) (3 : 249 2 : 425 3 : 734) ESD Err. (0 : 055 1 : 173 2 : 471) (0 : 757 0 : 706 0 : 744) 4 : 05 MI Err. (0 : 795 1 : 313 1 : 479) (0 : 334 0 : 115 0 : 584) 71 : 25 6 T rue ( 1 : 381 3 : 386 30 : 028) (0 : 6022 7 : 4366 7 : 125) ESD Err. (0 : 328 0 : 303 0 : 003) (0 : 823 0 : 233 0 : 111) 4 : 07 MI Err. (0 : 031 0 : 498 0 : 195) (0 : 173 0 : 191 0 : 313) 29 : 43 7 T rue (13 : 922 3 : 965 18 : 719) (8 : 700 7 : 328 22 : 421) ESD Err. (1 : 320 2 : 721 1 : 611) (1 : 296 1 : 268 0 : 752) 3 : 84 MI Err. (0 : 487 1 : 461 1 : 833) (1 : 060 0 : 818 0 : 179) 60 : 96 8 T rue (14 : 024 3 : 569 21 : 778) (0 : 120 12 : 970 9 : 870) ESD Err. (0 : 419 0 : 946 0 : 6356) (0 : 382 2 : 079 1 : 081) 3 : 97 MI Err. (0 : 115 1 : 177 2 : 681) (0 : 500 1 : 525 0 : 770) 70 : 19 9 T rue (5 : 357 2 : 796 20 : 174) (5 : 217 2 : 611 1 : 156) ESD Err. (0 : 588 0 : 008 0 : 481) (0 : 721 1 : 191 0 : 593) 4 : 01 MI Err. (0 : 713 0 : 193 0 : 943) (0 : 237 1 : 039 0 : 409) 42 : 23 1 0 T rue ( 5 : 477 1 : 364 20 : 730) (0 : 081 3 : 989 0 : 658) ESD Err. (0 : 023 0 : 091 0 : 398) (0 : 707 0 : 547 0 : 355) 3 : 99 MI Err. (0 : 452 0 : 333 0 : 372) (0 : 251 0 : 461 0 : 053) 38 : 12 M ean ESD Err. (0 : 375 0 : 619 1 : 070) (0 : 477 0 : 737 0 : 505) 3 : 86 MI Err. (0 : 443 0 : 654 0 : 886) (0 : 322 0 : 492 0 : 331) 52 : 02 S T D ESD Err. (0 : 388 0 : 840 1 : 435) (0 : 357 0 : 617 0 : 341) 0 : 39 MI Err. (0 : 269 0 : 503 0 : 880) (0 : 296 0 : 485 0 : 241) 13 : 45 PAGE 73 65 T able 7.3: Comparison of estimation errors b et w een SSD, L 2 E and Loren tzian function based estimator. mean standard deviation SSD (0 : 0693 0 : 1320 3 : 0775) (0 : 1268 0 : 3626 38 : 7502) (0 : 0430 0 : 0836 1 : 860) (0 : 1128 0 : 0399 1 : 9884) Loren tzian (0 : 0197 0 : 0832 0 : 9194) (0 : 1167 0 : 0840 9 : 5005) (0 : 0428 0 : 1142 1 : 2636) (0 : 1624 0 : 1180 13 : 3713) L 2 E (0 : 0085 0 : 0188 0 : 1337) (0 : 0293 0 : 0224 1 : 1088) (0 : 0190 0 : 0410 0 : 3068) (0 : 0053 0 : 0485 2 : 3473) Lor entian ( w or k ing cases onl y ) (0 : 0001 0 : 0002 0 : 0001) (0 : 0020 0 : 0010 0 : 0057) (0 : 0020 0 : 0031 0 : 0023) (0 : 0033 0 : 0013 0 : 0089) L 2 E ( w or k ing cases onl y ) (0 : 0001 0 : 0000 0 : 0001) (0 : 0002 0 : 0002 0 : 0059) (0 : 0012 0 : 0000 0 : 0021) (0 : 0004 0 : 0023 0 : 0120) T able 7.4: Comparison of estimation errors b et w een L 2 E and normalized MI Ov erlapping region only Whole image L 2 E (0 : 014 0 : 005 mm 0 : 010 mm ) (0 : 017 0 : 002 mm 0 : 034 mm ) normalized MI (0 : 165 0 : 119 mm 0 : 115 mm ) (0 : 935 0 : 634 mm 0 : 303 mm ) T able 7.5: Estimation accuracy comparison b et w een magnitude and v ector form of lo cal frequency represen tations (see text for details). estimation using magnitude lo cal frequency estimation using v ector form ideal 0 : 0040 0 : 0010 0 : 0310 0 : 0020 0 : 0070 0 : 0450 0 : 0000 0 : 0000 0 : 0010 0 : 0010 0 : 0000 0 : 0010 1 : 0 0 : 0150 0 : 0040 0 : 0308 0 : 0190 0 : 0100 0 : 1668 0 : 2764 0 : 069 0 : 0176 0 : 1937 0 : 2314 0 : 000 2 : 0 0 : 0150 0 : 0041 0 : 0312 0 : 0240 0 : 0207 0 : 1723 0 : 3605 0 : 0351 0 : 0148 0 : 1937 0 : 1791 0 : 0016 4 : 0 0 : 0160 0 : 0050 0 : 0430 0 : 0134 0 : 0002 0 : 0472 0 : 3588 0 : 0402 0 : 0138 0 : 1987 0 : 2001 0 : 0300 9 : 0 0 : 0782 0 : 0866 0 : 0857 0 : 0071 0 : 0102 0 : 0317 0 : 1783 0 : 1221 0 : 082 0 : 3002 0 : 4321 0 : 0600 T able 7.6: Comparison of estimation errors b et w een L 2 E and normalized MI for T1/T2 data set 1 (10 ; 5 mm ) L 2 E normalized MI maximal error (0 : 2082 0 : 1832 0 : 0930) (0 : 2430 0 : 2380 0 : 2415) (0 : 3729 0 : 2135 0 : 5543) (0 : 4440 0 : 3328 0 : 1639) error mean (0 : 1566 0 : 1468 0 : 0510) (0 : 2430 0 : 2380 0 : 2415) (0 : 2416 0 : 1535 0 : 2306) (0 : 2193 0 : 2048 0 : 1099) error STD (0 : 0359 0 : 0257 0 : 0283) (0 : 0955 0 : 0707 0 : 0463) (0 : 0921 0 : 0555 0 : 1909) (0 : 0360 0 : 0294 0 : 0100) PAGE 74 66 Figure 7.5: Registration of a T1-T2 w eigh ted brain MRI with large non-o v erlap. First ro w: T1 w eigh ted image view ed from (a) sagittal, (b) coronal and (c) axial planes. Second ro w: T2 w eigh ted image. The third ro w is the transformed T1 image with sup erimp osed edge maps from T2 image. PAGE 75 67 T able 7.7: Comparison of estimation errors b et w een L 2 E and normalized MI for T1/T2 data set 2 (20 ; 10 mm ) L 2 E normalized MI maximal error (0 : 3531 0 : 1780 0 : 0993) (0 : 2056 0 : 2846 0 : 4830) (0 : 5757 0 : 4560 0 : 7410) (0 : 3960 0 : 6220 0 : 2238) error mean (0 : 1166 0 : 1362 0 : 0743) (0 : 0725 0 : 1101 0 : 3428) (0 : 3480 0 : 2302 0 : 3749) (0 : 1990 0 : 1353 0 : 1429) error STD (0 : 0459 0 : 0280 0 : 0229) (0 : 0309 0 : 0670 0 : 0992) (0 : 1573 0 : 1158 0 : 1090) (0 : 1290 0 : 2240 0 : 0660) T able 7.8: Comparison of estimation errors b et w een L 2 E and normalized MI for T1/T2 data set 3 (30 ; 15 mm ) L 2 E normalized MI maximal error (0 : 2080 0 : 1834 0 : 0855) (0 : 2340 0 : 4260 0 : 4600) (0 : 5022 0 : 3771 0 : 3830) (0 : 2480 0 : 1860 0 : 2245) error mean (0 : 1011 0 : 1760 0 : 0486) (0 : 1386 0 : 2315 0 : 2918) (0 : 3117 0 : 2390 0 : 2142) (0 : 1431 0 : 1495 0 : 1883) error STD (0 : 0726 0 : 0466 0 : 0227) (0 : 0795 0 : 1465 0 : 1200) (0 : 0649 0 : 0845 0 : 1369) (0 : 0888 0 : 0468 0 : 0601) T able 7.9: Comparison of estimation errors b et w een L 2 E and normalized MI for CT-MR data set 1 (10 ; 5 mm ) L 2 E normalized MI maximal error (0 : 5072 0 : 8040 0 : 2624) (0 : 5535 0 : 4500 0 : 3850) (0 : 698 0 : 4753 1 : 1776) (0 : 6090 0 : 4280 0 : 3200) error mean (0 : 1721 0 : 1449 0 : 1574) (0 : 3111 0 : 2192 0 : 2036) (0 : 5213 0 : 1868 0 : 6425) (0 : 3229 0 : 2071 0 : 2244) error STD (0 : 1920 0 : 2084 0 : 0878) (0 : 0335 0 : 1426 0 : 1143) (0 : 2132 0 : 1949 0 : 3547) (0 : 1197 0 : 0788 0 : 3810) T able 7.10: Comparison of estimation errors b et w een L 2 E and normalized MI for CT-MR data set 2 (20 ; 10 mm ) L 2 E normalized MI maximal error (0 : 2354 0 : 4566 0 : 2292) (0 : 4190 0 : 4720 0 : 3800) (0 : 5700 0 : 6988 1 : 2426) (0 : 5932 0 : 9820 0 : 3938) error mean (0 : 1891 0 : 2832 0 : 1223) (0 : 3293 0 : 1220 0 : 2945) (0 : 3078 0 : 3231 0 : 4695) (0 : 1712 0 : 3514 0 : 2638) error STD (0 : 0713 0 : 2056 0 : 1008) (0 : 1029 0 : 0654 0 : 0741) (0 : 2261 0 : 1764 0 : 3711) (0 : 2096 0 : 3211 0 : 1024) PAGE 76 68 T able 7.11: Comparison of estimation errors b et w een L 2 E and normalized MI for CT-MR data set 3 (30 ; 15 mm ) L 2 E normalized MI maximal error (0 : 1462 0 : 4396 0 : 2072) (0 : 4982 0 : 752 0 : 4566) (0 : 7410 0 : 9690 1 : 3680) (0 : 8049 0 : 5900 0 : 6329) error mean (0 : 0740 0 : 3720 0 : 1372) (0 : 2774 0 : 4090 0 : 2763) (0 : 3648 0 : 3884 0 : 7565) (0 : 2014 0 : 2263 0 : 3342) error STD (0 : 054 0 : 0677 0 : 0808) (0 : 1404 0 : 207 0 : 1387) (0 : 2800 0 : 2744 0 : 5541) (0 : 2830 0 : 4024 0 : 2043) PAGE 77 CHAPTER 8 CONCLUSIONS AND FUTURE W ORK Tw o fundamen tal questions for m ulti mo dalit y image registration problem are: 1) what is a go o d image represen tation to w ork with? 2) what is an appropriate similarit y measure for matc hing the t w o images based on the selected represen tation? Mutual information is curren tly the most p opular metho d b eing used in this area. Ho w ev er MI metho ds ha v e their inheren t dra wbac ks, whic h can b e summarized as lac k of computational eciency , lac k of robustness and the lac k of ease in coping with nonrigid. Additionally , MI is not a metric and hence do es not satisfy the standard prop erties of a metric. 8.1 Ma jor Con tributions This dissertation is an attempt to nd an alternativ e in place of MI metho ds. Lo cal frequency is a desirable represen tation used for m ulti-mo dal image registration due to the fact that: it is (1) insensitiv e to the c hange of brigh tness and con trast v ariations; (2) main tains directional information and (3) naturally allo ws for pro cessing the data at dieren t scales/resolutions. Computational eciency is a v ery imp ortan t problem is 3D medical image applications. In this dissertation, w e dev elop ed a fast algorithm that in v olv es minimizing the exp ectation of the squared dierence b et w een the lo cal-frequency represen tations of the source and target images. Exp erimen ts w ere carried out using real clinical CT and MR brain scans. Our algorithm's p erformance is comparable to the results obtained from a Mutual Information-based algorithm in the con text of accuracy of estimated rigid transforms but is m uc h faster in computational sp eed. 69 PAGE 78 70 Most existing metho ds in the literature are unable to cop e with registration of image pairs with large non-o v erlapping eld of view (F O V). The robust algorithm w e presen ted to handle this problem in v olv es minimizing the in tegral of the squared error (ISE or L 2 E) b et w een a Gaussian mo del of the residual and its true densit y function. W e exp erimen tly sho w that our L 2 E based sc heme yields b etter accuracy o v er the normalized MI. In summary , the ma jor con tribution of this thesis can b e summarized as: 1 A lo cal frequency image represen tation used to capture common information. 2 An ESD measure together with a mo died Newton metho d to ac hiev e a v ery fast 3D m ulti-mo dal image registration implemen tation. 3 An no v el L 2 E measure handle the dieren t elds of view problem. 8.2 F uture W ork 1. F urther researc h on m ulti-mo dal nonrigid motion estimation. Up to no w, w e obtained some preliminary results applying our lo cal frequency for nonrigid motion estimation, but the v alidation is based only on visual insp ection. F urther study can b e carried out in follo wing directions: Nonrigid deformation estimation based on prop erties of lo cal frequency . Lo cal frequency has some w ell-established mathematical prop erties whic h can b e used to describ e the nonrigid motion in the realit y . In tegrate nonrigid deformation estimation in to L 2 E framew ork. Curren tly the lo cal nonrigid deformation estimation pro cess is separated from global estimation. If w e can in tegrate it in to the robust L 2 E framew ork, it will pro vide the users with a highly-coupled system and will mak e it p ossible to handle lo cal and global deformations in a single robust framew ork. 2. F urther exploration of L 2 E measure. PAGE 79 71 W e are the rst one to in tro duce the L 2 E measure in to the computer vision area. It is w orth while to apply it on to more applications to exploit its p o w er. The planned examination include: Apply L 2 E to some other general computer vision problems, suc h as curv e tting, to see if an y impro v emen t to the state of art can b e made. Under registration con text, carry out broader comparison with more robust estimators, including W elsc h estimators, least median of square [66 ] etc. PAGE 80 APPENDIX Let I 1 and I 2 b e a pair of images b eing registered. Let us assume that the co ordinate transformation b et w een them is mo deled b y an ane transformation T ( X ) = AX + b , and I 1 ( T ( X )) = I 2 ( X ). Then, the corresp onding v ector-v alued lo cal frequency magnitude represen tations are related b y: A T ~ F 1 ( AX + b ) = ~ F 2 ( X ); (1) W e in tro duce an in termediate image I m suc h that I m ( X ; ^ T ) = I 1 ( ^ A X + ^ b ) (2) where the sym b ol ^ is used to indicate computed v alues whic h are appro ximations to the true v alues. V ector-v alued lo cal frequency magnitude represen tation ~ F m of I m satises: ~ F m ( X ; ^ T ) = ^ A T ~ F 1 ( ^ A X + ^ b ); (3) Our energy function for ane registration is then giv en b y: E sd ( ^ T ) = E jj ~ F m ( X ; ^ T ) jj 2 jj ~ F 2 ( X ) jj 2 2 (4) where the jj ~ F 1 jj and jj ~ F 2 jj are just the scalar-v alued lo cal frequency magnitude represen tations used in this pap er. The gradien t of the energy function with resp ect to the parameter v ector ^ T is: @ E sd @ ^ T = 2 E jj ~ F m ( X ; ^ T ) jj 2 jj ~ F 2 ( X ) jj 2 @ jj ~ F m ( X ; ^ T ) jj 2 @ ^ T !! (5) 72 PAGE 81 73 the Hessian of E sd is then giv en b y: H ( ^ T ) = @ 2 E sd ( ^ T ) @ ^ T 2 (6) = 2 E 0@ @ jj ~ F m ( X ; ^ T ) jj 2 @ ^ T @ jj ~ F m ( X ; ^ T ) jj 2 @ ^ T ! T + jj ~ F m ( X ; ^ T ) jj 2 jj ~ F 2 jj 2 @ 2 jj ~ F m ( X ; ^ T ) jj 2 @ ^ T 2 ! Note that at the optim um ~ F m ( X ; ^ T ) = ~ F m ( X ; T ) = ~ F 2 ( X ), hence, H ( T ) = @ 2 E sd ( T ) @ ^ T 2 = 2 E @ jj ~ F m ( X ; T ) jj 2 @ ^ T ( @ jj ~ F m ( X ; T ) jj 2 @ ^ T ) T ! (7) F or simplicit y , w e only sho w the deriv ation for 2D case, where ~ F m = ( F m 1 ; F m 2 ) T , but ev erything carries o v er to the 3D case. Con tin uing with equation (7) w e ha v e: H ( T ) = 2 E @ ( F 2 m 1 ( X ; T ) + F 2 m 2 ( X ; T )) @ ^ T @ ( F 2 m 1 ( X ; T ) + F 2 m 2 ( X ; T )) @ ^ T T ! (8) = 8 E @ F m 1 ( X ; T ) @ ^ T ; @ F m 2 ( X ; T ) @ ^ T ( F m 1 ( X ; T ) ; F m 2 ( X ; T )) T ( F m 1 ( X ; T ) ; F m 2 ( X ; T )) @ F m 1 ( X ; T ) @ ^ T ; @ F m 2 ( X ; T ) @ ^ T T ! = 8 E @ F m 1 ( X ; T ) @ ^ T ; @ F m 2 ( X ; T ) @ ^ T ( F 21 ( X ) ; F 22 ( X )) T ( F 21 ( X ) ; F 22 ( X )) @ F m 1 ( X ; T ) @ ^ T ; @ F m 2 ( X ; T ) @ ^ T T ! = 8 E F 2 21 ( X ) @ F m 1 ( X ; T ) @ ^ T @ F m 1 ( X ; T ) @ ^ T T + F 21 ( X ) F 22 ( X ) @ F m 1 ( X ; T ) @ ^ T @ F m 2 ( X ; T ) @ ^ T T + F 22 ( X ) F 21 ( X ) @ F m 2 ( X ; T ) @ ^ T @ F m 1 ( X ; T ) @ ^ T T (9) + F 2 22 ( X ) @ F m 2 ( X ; T ) @ ^ T + @ F m 2 ( X ; T ) @ ^ T T ! PAGE 82 74 This expression in general is not indep enden t of the transformation on parameters. W e will therefore consider computation of the Hessian with resp ect to the incremen tal transformations ~ T as discussed in section 4.1. Before getting in to an y details, let's list some kno wn facts rst: Since ~ F m ( X ; ^ T ) = ^ A T ~ F 1 ( ^ A X + ^ b ); (10) Where ^ T = ( ^ t 0 ; ^ t 1 ; ^ t 2 ; ^ t 3 ; ^ t 4 ; ^ t 5 ) T , ^ A = ^ t 0 ^ t 1 ^ t 3 ^ t 4 and ^ b = ( ^ t 2 ^ t 5 ) T . A t iteration N , where the co ordinate system transforms from X N to X N + 1 , the follo wing relation holds: 8><>: F N +1 m 1 = ~ t N +1 0 F N m 1 + ~ t N +1 3 F N m 2 ; F N +1 m 2 = ~ t N +1 1 F N m 1 + ~ t N +1 4 F N m 2 ; (11) Supp ose w e reac h the optim um at iteration N , then the transformation ^ T N will b e equal to ^ T N +1 and the incremen tal transformation v ector ~ T N +1 is giv en b y: ~ T N +1 = ( ~ t N +1 0 ; ~ t N +1 1 ; ~ t N +1 2 ; ~ t N +1 3 ; ~ t N +1 4 ; ~ t N +1 5 ) = (1 ; 0 ; 0 ; 0 ; 1 ; 0 ) No w let's calculate the partial deriv ativ e of F m 1 ( X ; ^ T ) with resp ect to eac h comp onen t of ~ T at the optim um. First, let us examine the deriv ativ e with resp ect to ~ t 0 at step N + 1: @ F N +1 m 1 @ ~ t 0 N +1 = @ ( ~ t 0 N +1 F N m 1 + ~ t 3 N +1 F N m 2 ) @ ~ t 0 N +1 (12) = F N m 1 + ~ t 0 N +1 @ F N m 1 @ ~ t 0 N +1 + ~ t 3 N +1 @ F N m 2 @ ~ t 0 N +1 (13) = F N m 1 + ~ t 0 N +1 @ F N m 1 @ ~ t 0 N +1 Let 8><>: x N +1 = ~ t 0 N +1 x N + ~ t 1 N +1 y N + ~ t 2 N +1 ; y N +1 = ~ t 3 N +1 x N + ~ t 4 N +1 y N + ~ t 5 N +1 (14) PAGE 83 75 then 8><>: @ x N +1 @ ~ t 0 N +1 = x N ; @ x N +1 @ ~ t 1 N +1 = y N ; @ x N +1 @ ~ t 2 N +1 = 1 ; @ x N +1 @ ~ t 3 N +1 = 0 ; @ x N +1 @ ~ t 4 N +1 = 0 ; @ x N +1 @ ~ t 5 N +1 = 0 @ y N +1 @ ~ t 0 N +1 = 0 ; (15) 8><>: @ y N +1 @ ~ t 1 N +1 = 0 ; @ y N +1 @ ~ t 2 N +1 = 0 ; @ y N +1 @ ~ t 3 N +1 = x N ; @ y N +1 @ ~ t 4 N +1 = y N ; @ y N +1 @ ~ t 5 N +1 = 1 ; @ x N @ x N +1 = 1 ; @ y N @ y N +1 = 1 (16) Th us, @ F N +1 m 1 @ ~ t 0 N +1 = F N m 1 + ~ t 0 N +1 @ F N m 1 @ x N @ x N @ x N +1 @ x N +1 @ ~ t 0 N +1 (17) = F N m 1 + 1 @ F N m 1 @ x N 1 x N Since, 8><>: F m 1 ( X N ; ^ T N ) = ^ t N0 F 11 ( X N ) + ^ t N3 F 12 ( X N ); F m 2 ( X N ; ^ T N ) = ^ t N1 F 11 ( X N ) + ^ t N4 F 12 ( X N ) (18) @ F N +1 m 1 @ ~ t 0 N +1 = ( ^ t N0 ; ^ t N3 )( F N 11 ; F N 12 ) T + ( ^ t N0 ; ^ t N3 )( @ F N 11 @ x N x N ; @ F N 12 @ x N x N ) T (19) Similarly , @ F N +1 m 2 @ ~ t 0 N +1 = @ ( ~ t 1 N +1 F N m 1 + ~ t 4 N +1 F N m 2 ) @ ~ t 0 N +1 (20) = ~ t 1 N +1 @ F N m 1 @ ~ t 0 N +1 + ~ t 4 N +1 F N m 2 @ ~ t 0 N +1 (21) = ~ t 4 N +1 @ F N m 2 @ ~ t 0 N +1 (22) = @ F N m 2 @ x N @ x N @ x N +1 @ x N +1 @ ~ t 0 N +1 (23) = @ F N m 2 @ x N x N = ( ^ t N1 ; ^ t N4 )( @ F N 11 @ x N x N ; @ F N 12 @ x N x N ) T PAGE 84 76 Using the same pro cedure, w e can deriv e the expressions for the deriv ativ e with resp ect to the other parameters of ~ T N +1 . These are simply summarized as follo ws: @ F N +1 m 1 @ ~ t 1 N +1 = ( ^ t N0 ; ^ t N3 )( @ F N 11 @ x N y N ; @ F N 12 @ x N y N ) T ; (24) @ F N +1 m 2 @ ~ t 1 N +1 = ( ^ t N0 ; ^ t N3 )( F N 11 ; F N 12 ) T + ( ^ t N1 ; ^ t N4 )( @ F N 11 @ x N y N ; @ F N 12 @ x N y N ) T (25) @ F N +1 m 1 @ ~ t 2 N +1 = ( ^ t N0 ; ^ t N3 )( @ F N 11 @ x N ; @ F N 12 @ x N ) T ; (26) @ F N +1 m 2 @ ~ t 2 N +1 = ( ^ t N1 ; ^ t N4 )( @ F N 11 @ x N ; @ F N 12 @ x N ) T (27) @ F N +1 m 1 @ ~ t 3 N +1 = ( ^ t N1 ; ^ t N4 )( F N 11 ; F N 12 ) T + ( ^ t N0 ; ^ t N3 )( @ F N 11 @ y N x N ; @ F N 12 @ y N x N ) T (28) @ F N +1 m 2 @ ~ t 3 N +1 = ( ^ t N1 ; ^ t N4 )( @ F N 11 @ y N x N ; @ F N 12 @ y N x N ) T ; (29) @ F N +1 m 1 @ ~ t 4 N +1 = ( ^ t N0 ; ^ t N3 )( @ F N 11 @ y N y N ; @ F N 12 @ y N y N ) T ; (30) @ F N +1 m 2 @ ~ t 4 N +1 = ( ^ t N1 ; ^ t N4 )( F N 11 ; F N 12 ) T + ( ^ t N1 ; ^ t N4 )( @ F N 11 @ y N y N ; @ F N 12 @ y N y N ) T (31) @ F N +1 m 1 @ ~ t 5 N +1 = ( ^ t N0 ; ^ t N3 )( @ F N 11 @ y N ; @ F N 12 @ y N ) T ; (32) @ F N +1 m 2 @ ~ t 5 N +1 = ( ^ t N1 ; ^ t N4 )( @ F N 11 @ y N ; @ F N 12 @ y N ) T (33) F rom the ab o v e deriv ations, w e observ e that @ F N +1 m 1 @ ~ t N +1 k ( k = 0 ; 1 ; :::; 5) all tak e the form of ~ L 1 ~ G 1 ( X N ) + ~ L 2 ~ G 2 ( X N ), where ~ L 1 and ~ L 2 are t w o-comp onen t PAGE 85 77 ro w v ectors con taining an y t w o of ^ t Ni ( i = 0 ; 1 ; :::; 5) ; ~ G 1 and ~ G 2 are t w o-comp onen t column v ectors con taining comp onen ts prop ortional to r x F 11 or r y F 12 , where r x and r y represen t the x and y comp onen ts of the gradien t v ector. Similarly w e will use ~ P 1 ~ Q 1 ( X N ) + ~ P 2 ~ Q 2 ( X N ) to denote the form of @ F N +1 m 2 @ ~ t N +1 k ( k = 0 ; 1 ; :::; 5). Since the Hessian with resp ect to incremen tal transformation ~ T at the optim um is giv en b y: @ 2 E sd ( T ) @ ~ T 2 = 2 E @ jj ~ F m ( X ; T ) jj 2 @ ~ T ( @ jj ~ F m ( X ; T ) jj 2 @ ~ T ) T ! = 8 E F 2 21 ( X ) @ F m 1 ( X ; T ) @ ~ T @ F m 1 ( X ; T ) @ ~ T T + F 21 ( X ) F 22 ( X ) @ F m 1 ( X ; T ) @ ~ T @ F m 2 ( X ; T ) @ ~ T T + F 22 ( X ) F 21 ( X ) @ F m 2 ( X ; T ) @ ~ T @ F m 1 ( X ; T ) @ ~ T T + F 2 22 ( X ) @ F m 2 ( X ; T ) @ ~ T @ F m 2 ( X ; T ) @ ~ T T ! An y particular elemen t ( i; j ) of the Hessian matrix tak es the form ~ H ij = E 8 F 2 21 ( X N ) @ F N +1 m 1 @ ~ t i @ F N +1 m 1 @ ~ t j + E 8 F 21 ( X N ) F 22 ( X N ) @ F N +1 m 1 @ ~ t i @ F N +1 m 2 @ ~ t j + E 8 F 22 ( X N ) F 21 ( X N ) @ F N +1 m 2 @ ~ t i @ F N +1 m 1 @ ~ t j + E 8 F 2 22 ( X N ) @ F N +1 m 2 @ ~ t i @ F N +1 m 2 @ ~ t j (34) PAGE 86 78 After a tedious y et a straigh tforw ard deriv ation, w e get: ~ H ij = E 8 F 2 21 ( X N )( ~ L 1 i ~ G 1 i ( X N ) + ~ L 2 i ~ G 2 i ( X N )) (35) ( ~ G 1 j ( X N ) T ~ L T1 j + ~ G 2 j ( X N ) T ~ L T2 j ) + E 8 F 21 ( X N ) F 22 ( X N )( ~ L 1 i ~ G 1 i ( X N ) + ~ L 2 i ~ G 2 i ( X N )) ( ~ Q 1 j ( X N ) T ( ~ P 1 j ) T + ~ Q 2 j ( X N ) T ( ~ P 2 j ) T ) + E 8 F 22 ( X N ) F 21 ( X N )( ~ P 1 i ~ Q 1 i ( X N ) + ~ P 2 i ~ Q 2 i ( X N )) ( ~ G 1 j ( X N ) T ( ~ L 1 j ) T + ~ G 2 j ( X N ) T ( ~ L 2 j ) T ) + E 8 F 2 22 ( X N )( ~ P 1 i ~ Q 1 i ( X N ) + ~ P 2 i ~ Q 2 i ( X N )) ( ~ Q 1 j ( X N ) T ( ~ P 1 j ) T + ~ Q 2 j ( X N ) T ( ~ P 2 j ) T ) Expanding the righ t hand side of the equation (36) will yield 16 terms. T aking the rst of these terms w e get: E 8 F 2 21 ( X N ) ~ L 1 i ~ G 1 i ( X N ) ( ~ G 1 j ( X N )) T ~ L T1 j = ~ L 1 i E 8 F 2 21 ( X N ) ~ G 1 i ( X N ) ( ~ G 1 j ( X N )) T ~ L T1 j W e assume that the exp ectation v alue is indep enden t of the co ordinate system and therefore, E 8 F 2 21 ( X N ) ~ G 1 i ( X N ) ( ~ G 1 j ( X N )) T = E 8 F 2 21 ( X ) ~ G 1 i ( X ) ( ~ G 1 j ( X )) T This term is only dep enden t on the lo cal frequency represen tations, F 11 ( X ), F 12 ( X ), from source image I 1 ( X ), and F 21 ( X ), F 22 ( X ) from target image I 2 ( X ) , so it is a function of I 1 ( X ) and I 2 ( X ) and therefore can b e pre-computed. W e use s 1 ( I 1 ; I 2 ) to denote it. The situation is the same for the other 15 pro duct terms, so o v erall, the Hessian of E sd with resp ect to the incremen tal transformation ~ T can b e computed using: ~ H ij = k =16 X k =1 a k ( ^ T N ) s k ( I 1 ; I 2 ) c k ( ^ T N ) T (36) PAGE 87 79 where a k and c k are functions whic h tak e ^ T N as argumen ts and generate ~ L 1 i , ~ L 1 j , ~ L 2 i , ~ L 2 j , ~ G 1 i ; ::: etc. F rom ab o v e deriv ation,w e can see that the Hessian with resp ect to ~ T is not totally pre-computable. But most of computational exp ense in computing ~ H can b e reduced b y pre-computing s k ( I 1 ; I 2 ). The Hessian with resp ect to ~ T at eac h iteration p oin t ^ T k is giv en b y: ~ H k + 1 ij = k =16 X k =1 a k ( ^ T k ) s k ( I 1 ; I 2 ) c k ( ^ T k ) T (37) The transformation v ector can then b e up dated using the follo wing equations: ~ T k +1 = ( ~ H k +1 ) 1 ~ g ( ^ T k ); ^ T k +1 = ~ T k +1 ^ T k (38) The comp osition dep ends on the t yp e of motion mo del used, whic h is an ane mo del in this app endix. PAGE 88 REFERENCES [1] D. L. G. Hill, D. J. Ha wk es, J. E. Crossman, M. J. Gleeson, T. C. S. Co x, E. C. Bracey , A. J. Strong, and P . Gra v es, \Registration of MR and CT images for skull base surgery using p oin tlik e anatomical features", British Journal of R adiolo gy , 64(767):1030-1035, 1991. [2] R. W o o ds, J. Maziotta and S. Cherry \MRI-PET registration with automated algorithm," J. Comput. Assis. T omo gr. , 17, pp. 536-546, 1993. [3] M. E. Lev en ton and W. E. L. Grimson, \Multimo dal v olume registration using join t in tensit y distributions," First Intl. Conf. on MICCAI , Cam bridge, MA, pp. 1057-1066, 1998. [4] P . A. Viola and W. M. W ells, \Alignmen t b y maximization of m utual information," in Fifth ICCV , MIT, Cam bridge, MA, pp. 16-23, 1995. [5] W. M. W ells I I I, P . Viola, and H. A tsumi, \Multi-mo dal V olume Registration b y Maximization of Mutual Information," Me dic al Image A nalysis , 1(1), pp. 35-51, 1997. [6] C. Studholme, D. L. G. Hill and D. J. Ha wk es, \Automated 3D registration of MR and CT images in the head," Me dic al Image A nalysis , 1(2), pp. 163-175, 1996. [7] B.C. V em uri, J. Y e, Y. Chen, and C. M. Leonard, \A Lev el-set based approac h to image registration," IEEE Workshop on MMBIA , June 10-12, 2000, Hilton Head, SC, 2000. [8] L. G. Bro wn, \A Surv ey of Image Registration T ec hniques," A CM Computing Survey , V ol. 24, No.4, Dec, 1992, pp. 325-376, 1992. [9] B. C. V em uri, S. Huang, S. Sahni, and C. M. Leonard, \An ecien t motion estimator with application to medical image registration," Me dic al Image A nalysis , 2 (1), pp. 79-98, 1998. [10] K. Arun, T. Huang, and S. Blostein: \Least Squares Fitting of t w o 3-D Sets," IEEE T r ansactions on Pattern A nalysis and Machine Intel ligenc e 9(5):698-700, 1987. [11] D. L. G. Hill, D. J. Ha wk es, N. Harrison and C. Ru, \A Strategy for Automated Multimo dalit y Image Registration Incorp orating Anatomical Kno wledge and Imagery Characteristics", Pr o c. Information Pr o c essing in Me dic al Imaging , pp. 182-196, 1993. 80 PAGE 89 81 [12] Z. He, J. Maublan t, J. Cauvin and A. V eyre, \Reorien tation of the Left V en tricular Long-axis on My o cardial T ransaxial T omograms b y a Linear Fitting Metho d," J. Nucl. Me d. , 32, pp. 1794-1800, 1991. [13] C. A. P elizarri, G. T. Y. Chen, and D. R. Sp elbring , \Accurate ThreeDimensional Registration of CT, PET and MR Images of the Brain," J. Comput. Assist. T omo gr. , 13(1), pp. 20-26, 1989. [14] P . J. Besl and N. D. Mac k a y , \A Metho d for Registration of 3D Shap es," IEEE T r ansactions Pattern A nalysis and Machine Intel ligenc e , v ol. 14, no. 2, pp. 239-254, F eb, 1992. [15] J. F eldmar and N. Ay ac he, \Lo cally Ane Registration of F ree-form Surfaces," Pr o c. of IEEE CVPR , pp. 496-501, Seattle, W A, 1994. [16] A. Rangara jan, H. Ch ui, E. Mjolsness, S. P appu, L. Da v ac hi, P . GoldmanRakie, and J. Duncan, \A Robust P oin t Matc hing Algorithm for Autoradiograph Alignmen t," Me dic al Image A nalysis , v ol. 4, no. 1, pp. 379-398, 1997. [17] H. Li, B. S. Manjunath, S. K. Mitra, \A Con tour-Based Approac h to Multisensor Image Registration," IEEE T r ansactions On Image Pr o c essing , V ol. 4, No.3, pp. 320-334, Marc h 1995. [18] A. C. Ev ans, \W arping of Computerized 3D A tlas to Matc h Brain Image V olumes for Quan titativ e Neuroanatomical and F unctional Analysis," Pr o c. SPIE Me dic al Imaging V , v ol. 1445, pp. 236-246, 1991. [19] S. La v allee, R. Szeliski and L. Brunie, \matc hing 3D smo oth surfaces with their 2D pro jections using 3D distance maps," in Pr o c. of SPIE , V ol. 1570, pp. 322-336, 1991. [20] C. Da v atzik os and J. L. Prince, \Brain Image Registration Based on Curv e Mapping," IEEE WBIA , pp. 245-254, Seattle, W A, 1994. [21] J. B. A. Main tz, P . A. v an den Elsen and M. A. Viergev er, \Comparison of Edge-based and Ridge-based Registration of CT and MR brain images," Me dic al Image A nalysis , 1(2), pp. 151-161, 1996. [22] G. Borgefors, \Hierarc hical Chamfer Matc hing: A P arametric Edge Matc hing Algorithm," IEEE T r ansactions on Pattern A nalysis and Machine Intel ligenc e , 10, pp. 849-865, 1988. [23] A. Gueziec and N. Ay ac he, \Smo othing and Matc hing of 3D Space Curv es," Robb. R. A. (ed.), Visualization in Biome dic al Computing , Pr o c. SPIE , 1808, pp. 259-273. SPIE Press, Bellingham, W A, 1992. [24] O. Monga, S. Bena y oun and O. D. F augeras, \F rom P artial Deriv ativ es of 3D Densit y Images to Ridge Lines," Robb. R. A. (ed.), Visualization PAGE 90 82 in Biome dic al Computing , Pr o c. SPIE , 1808, pp. 118-127. SPIE Press, Bellingham, W A, 1992. [25] P . A. V an Den Elsen, \Multimo dalit y Matc hing of Brain Images," Ph.D. Thesis , Utrec h t Univ ersit y , The Netherlands, 1993. [26] J.B. Main tz and M. A. Viergev er, \A Surv ey of Medical Image Registration," Me dic al Image A nalysis V ol. 2, pp. 1-36, 1998. [27] W. Zhao, T. Y. Y oung, and M. D. Ginsb erg, \Registration of threedimensional autoradiographic images b y the disparit y analysis metho d," IEEE T r ansactions on Me dic al Imaging , V ol. 12, No. 4, pp. 782-791, 1993. [28] N. Alp ert, J. Bradsha w, D. Kennedy and J. Correia, \The Principal Axis T ransformation A Metho d for Image Registration," J. Nucl. Me d. , 31, pp. 1717-1722, 1990. [29] M. Okuomi and T. Kanade, \A Lo cally Adaptiv e Windo w for Signal Matc hing," International Journal of Computer Vision , 7(2), pp. 143-162, 1992. [30] R. Szeliski and J. Coughlan, \Spline-Based Image Registration," T ec hnical Rep ort 94/1, Digital Equipmen t Corp oration, Cam bridge Researc h Lab, April 1994. [31] B. P . K. Horn and B. G. Sc h unk, \Determining Optical Flo w". A rticial Intel ligenc e , V ol. 17, pp. 185-203, 1981. [32] J. L. Barron, D. J. Fleet, and S. S. Beauc hemin, \P erformance of Optical Flo w T ec hniques," International Journal of Computer Vision , 1(12):43-77, 1994. [33] J. H. Duncan and T. C. Chou \On the Detection of Motion and the Computation of Optical Flo w," IEEE T r ansactions on Pattern A nalysis and Machine Intel ligenc e , 14(3):346-352, 1992. [34] S. H. Lai and B. C. V em uri, \Robust and Ecien t Computation of Optical Flo w," IEEE Symp osium on CVS , pp. 455{460, Miami, Florida, 1995. [35] M. J. Blac k and P . Anandan, \A F ramew ork for Robust Estimation of Optical Flo w," Pr o c. International Confer enc e on Computer Vision (ICCV) , pp. 231-236, Berlin, German y , 1993. [36] S. Gupta and J. Prince, \On V ariable Brigh tness Optical Flo w for T agged MRI," Pr o c. of IPMI , pp. 323-334, Dordredc h t, Klu w er, 1995. [37] A. Collignon, F. Maes, D. V andermeulen, P . Suetens and G. Marc hal, \Automated Multimo dalit y Image Registration Using Information Theory", In Bizais, Y., Barillot, C., and Di paoloa, R., (eds.), Pr o c. IPMI , pp. 263-274. Klu w er, Dordrec h t, 1995. PAGE 91 83 [38] R. O. Duda and P . E. Hart, Pattern Classic ation and Sc ene A nalysis , John Wiley and Sons, 1973. [39] C. Studholme, D. L. G. Hill and D. J. Ha wk es, \An o v erlap in v arian t en trop y measure of 3D medical image alignmen t," Pattern R e c o gnition , V ol. 32, pp. 71-86, 1999. [40] J. P . W. Pluim, J. B. A. Main tz and M. A. Viergev er \In terp olation artifacts in m utual information based image registration," Computer Vision and Image Understanding , 77(2), pp. 211-232, 2000. [41] E. Haac k e, Z. Liang, \Challenges of Imaging Structure and F unction with MRI," IEEE Engine ering in Me dicine and Biolo gy , Sept., pp. 55 62, 2000. [42] C. R. Mey er, J. L. Bo es, P . H. Bland, K. R. Zasadn y , P .V. Kison, K. Koral, K.A. F rey , and R. L. W ahl, \Demonstrating the accuracy and clinical v ersatilit y of m utual information for automatic m ultimo dalit y image fusion using ane and thin-plate spline w arp ed geometric deformations," Me dic al Image A nalysis , V ol. 1., No. 3, pp. 195-207, 1997. [43] N. Hata, W. M. W ells I I I, M. Halle, S. Nak a jima, P . Viola, R. Kikinis, F. A. Jolesz, \Image Guided Microscopic Surgery System using Mutual-Information based Registration," in Visualization in Biome dic al Computing , Ham burg, German y , pp. 317-326, 1996. [44] R. Bansal and L. H. Staib, \A No v el Approac h for Registration of 2D P ortal and 3D CT Images for T reatmen t Setup V erication in Radiotherap y ," First Intl. Conf. on MICCAI , Cam bridge, MA, pp. 1075-1086, 1998. [45] T. Gaens, F. Maes, D. V andermeulen and P . Suetens, \ Nonrigid Multimo dal Image Registration using Mutual Information," First Intl. Conf. on MICCAI , Cam bridge, MA, 1998. [46] A. Ro c he, G. Malandain, X. P ennec and N. Ay ac he, \The correlation ratio as a new similarit y measure for m ulti-mo dal image registration," First Intl. Conf. on MICCAI , Cam bridge, MA, pp. 1115-1124, 1998. [47] G. P . P enney , J. W eese, J. Little, P . Desmedt, D. L. G. Hill and D. Ha wk es, \A Comparison of similarit y measures for use in 2D-3D medical images," First Intl. c onfer enc e on MICCAI , Cam bridge, MA, pp. 1153-1161, 1998. [48] M.Irani and P .Anandan, \Robust Multi-sensor Image Alignmen t," ICCV 1998 , Bom ba y , India, pp. 959-965, 1998. [49] G.H.Granlund and H.Kn utsson, Signal Pr o c essing for Computer Vision . Klu w er Academic Publishers, AH Dordrec h t, Netherlands, 1995. PAGE 92 84 [50] J. L. Marro quin, M. Servin and R. Ro driguez-V era, 1997, \Adaptiv e quadrature lters and the reco v ery of phase from fringe pattern images," JOSA , 14(8), pp. 1742-1752, 1997. [51] A.C. Bo vik, M. Clark and W.S.Geisler, "Multic hannel T exture Analysis Using Lo calized Spatial Filters". IEEE T r ansactions Pattern A nalysis and Machine Intel ligenc e ,V ol. 12, No. 1, pp. 55-71, 1990. [52] F. Maes, D. V andermeulen and P . Suetens, \Comparativ e Ev aluation of Multiresolution Optimization Strategies for Multimo dalit y Image Registration b y Maximization of Mutual Information," Me dic al Image A nalysis , v ol. 3, no. 4, pp. 373-386, 1999. [53] P . Thev enaz, M. Unser, \Optimization of Mutual Information for Multiresolution Image Registration", IEEE T r ansactions on Image Pr o c essing , v ol. 9, no. 12, pp. 2083-2099, 2000. [54] T. Netsc h, P . Rosc h, A. V. Muiswink el and J. W eese, \T o w ards Real-Time Multi-Mo dalit y 3D Medical Image Registration", ICCV 2001 , v ol. 1, pp. 718-725, July , 2001. [55] N. Diehl and H. Burkhardt, \Motion Estimation in Image Sequences," High Pr e cision Navigation , pp. 297-312, Springer-V erlag, New Y ork, 1989. [56] J. Illingw orth and J. Kittler, \A Surv ey of the Hough T ransform," Comput. Vision Gr aph. Image Pr o c ess. , 44:87-116, 1988. [57] D. A. Belsley , E. Kuh and R. E. W elsc h, \Regression Diagnostics", John Wiley and Sons, USA, 1980. [58] P . J. Hub er, \Robust Statistics", Wiely: New Y ork, 1981. [59] D. F. Andrews, \A Robust Metho d for Multiple Linear Regression", T e chnometrics , 16:523-531, 1974. [60] P . E. Gill, W. W righ t and M. H, Pr actic al Optimization . Academic Press, London, 1981. [61] D. W. Scott, \P arametric mo deling b y minim um L 2 error," T ec hnical Rep ort 98-3, Dept. of Statistics, Rice Univ ersit y , 1998. [62] D. P . Bertsek as, Nonline ar Pr o gr amming , A thena Scien tic Publishers, 1999. [63] R. Kimmel, A. Amir and A. M. Bruc kstein, \Finding Shortest P aths on Surfaces using Lev el Sets Propagation", IEEE T r ansactions on Pattern A nalysis and Machine Intel ligenc e , 17(6):635-640, 1995. [64] J. Y e, \A Lev el-set based approac h to image registration", Master Thesis , Univ ersit y of Florida, 2000. PAGE 93 85 [65] Sim ulated brain database [Online]. Av ailable : h ttp://www.bic.mni.mcgill.ca/brain w eb/, 2002. [66] Z. Zhang, \P arameter Estimation T ec hniques: A T utorial with Application to Conic Fitting", Image and Vision Computing Journal , V ol.15, No.1, pp. 59-76, 1997. PAGE 94 BIOGRAPHICAL SKETCH Jundong Liu w as b orn in Xi Ning, P . R. China. He receiv ed his Bac helor of Science degree from the Computer Science Departmen t at W uhan Univ ersit y , P . R. China in 1993. He earned his Master of Science degree from the Computer Science Departmen t at P eking Univ ersit y , P . R. China, in 1996. He receiv ed his Do ctor of Philosoph y degree in computer and information science and engineering from Univ ersit y of Florida, Gainesville, in August 2002. His researc h in terests include medical imaging, image pro cessing, computer vision, computer graphics, pattern recognition and computational geometry . PAGE 95 I certify that I ha v e read this study and that in m y opinion it conforms to acceptable standards of sc holarly presen tation and is fully adequate, in scop e and qualit y , as a dissertation for the degree of Do ctor of Philosoph y. Baba C. V em uri, Chair Professor of Computer and Information Science and Engineering I certify that I ha v e read this study and that in m y opinion it conforms to acceptable standards of sc holarly presen tation and is fully adequate, in scop e and qualit y , as a dissertation for the degree of Do ctor of Philosoph y. Sarta j K. Sahni Distinguished Professor of Computer and Information Science and Engineering I certify that I ha v e read this study and that in m y opinion it conforms to acceptable standards of sc holarly presen tation and is fully adequate, in scop e and qualit y , as a dissertation for the degree of Do ctor of Philosoph y. Anand Rangara jan Asso ciate Professor of Computer and Information Science and Engineering I certify that I ha v e read this study and that in m y opinion it conforms to acceptable standards of sc holarly presen tation and is fully adequate, in scop e and qualit y , as a dissertation for the degree of Do ctor of Philosoph y. Tim Da vis Asso ciate Professor of Computer and Information Science and Engineering I certify that I ha v e read this study and that in m y opinion it conforms to acceptable standards of sc holarly presen tation and is fully adequate, in scop e and qualit y , as a dissertation for the degree of Do ctor of Philosoph y. Murali Rao Professor of Mathematics PAGE 96 This dissertation w as submitted to the Graduate F acult y of the College of Engineering and to the Graduate Sc ho ol and w as accepted as partial fulllmen t of the requiremen ts for the degree of Do ctor of Philosoph y. August 2002 Dean, College of Engineering Dean, Graduate Sc ho ol |