UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 APPLICATIONS OF VARIATIONAL PARTIAL DIFFERENTIAL EQUATION MODELS IN MEDICAL IMAGE PROCESSING By FENG HUANG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2004 PAGE 2 Copyright 2004 by Feng Huang PAGE 3 To my parents and my wife. PAGE 4 ACKNOWLEDGMENTS I wish to express my sincere and deepest appreciation to my advisor, Professor Yunmei Chen, for her guidance, advice and support throughout my graduate study at UF. Her enthusiasm, patience and intelligence have helped me achieve more in my studies. I would like to thank other members of my committee (David Wilson, Scott McCullough, Murali Rao, Yijun Liu and Randy Duensing) for their help, valuable suggestions and kind encouragement. I would also like to thank my colleagues for their support and discussions on various topics. Thanks go to Sheshadri Thiruvenkadam, James Akao, Charley Saylor, Weihong Guo, Andrew Rubin, Jungha An, Xiqiang Zheng. Chapter 2 was the result of collaboration with Yunmei Chen, Murali Rao, David Wilson, Sheshadri Thiruvenkadam, Weihong Guo, Jungha An, Hemant D. Tagare and Edward A. Geiser. Several publications based on the work in this chapter. Chapter 3 is a joint work with Yunmei Chen, Randy Duensing, Charley Saylor, James Akao and Scott McCullough. The work was submitted to IEEE Transaction on Medical Imaging. Chapter 4 was coauthored with Yunmei Chen, Randy Duensing, James Akao and Andrew Rubin. The work was submitted to Magnetic Resonance on Medicine. Chapter 5 was written in collaboration with Yunmei Chen, Yijun Liu, Weihong Guo, Qingguo Zeng, Xiqiang Zheng, Guojun He and Xiaolu Yan. This work was supported in part by Invivo Corporation (Gainesville) and McLaughlin Dissertation Fellowship. iv PAGE 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS.................................................................................................iv LIST OF TABLES.............................................................................................................ix LIST OF FIGURES.............................................................................................................x LIST OF ABBREVIATIONS...........................................................................................xii ABSTRACT.....................................................................................................................xiv CHAPTER 1 INTRODUCTION........................................................................................................1 1.1 Motivating Applications.........................................................................................2 1.1.1 Segmentation of Ultrasound Cardiac Image.................................................2 1.1.2 Noise Tomography.......................................................................................2 1.1.3 Calibration of Sensitivity Map.....................................................................3 1.1.4 Intensity Correction......................................................................................4 1.1.5 Fiber Tracking..............................................................................................5 1.2 Review of Conception............................................................................................6 1.2.1 Description of General Tomography Problem.............................................6 1.2.2 MumfordShah Model..................................................................................6 1.2.3 Total Variation Denoising............................................................................7 1.2.4 Level Set.......................................................................................................8 1.2.5 Distance Function.......................................................................................10 1.2.6 Mutual Information....................................................................................11 1.2.7 Inpainting....................................................................................................12 1.3 Numerical Techniques..........................................................................................13 1.3.1 Choice of Heaviside Function....................................................................13 1.3.2 Computation of uudiv ......................................................................13 1.3.3 Initial Signed Level Set and Reinitialization..............................................14 1.3.4 Fast Sweep Method for Level Set Based Optimization..............................15 1.4 Contributions of Our Study..................................................................................16 1.4.1 Segmentation with Prior Information.........................................................16 1.4.2 Inpainting....................................................................................................16 v PAGE 6 1.4.3 Modified MumfordShah Model for Tomography.....................................17 1.4.4 Fiber Tracking............................................................................................18 2 SEGMENTATION WITH PRIOR INFORMATION................................................19 2.1 Introduction...........................................................................................................19 2.1.1 Definition of Segmentation........................................................................19 2.1.2 Importance of Segmentation in Medical Applications...............................19 2.1.3 Review of Segmentation Methods..............................................................21 2.1.4 Remaining Problems in Existing Methods.................................................24 2.2 Method..................................................................................................................24 2.2.1 Generation of Prior Information.................................................................24 2.2.1.1 Shape analysis..................................................................................25 2.2.1.2 Intensity profile................................................................................28 2.2.2 Framework Description..............................................................................31 2.2.2.1 Gradient term....................................................................................32 2.2.2.2 Shape term........................................................................................34 2.2.2.3 Key points term................................................................................36 2.2.2.4 Intensity profile term........................................................................38 2.2.2.5 Proposed model and EulerLagrange equations...............................40 2.2.3 Parameter Choosing....................................................................................42 2.3 Experiments and Results.......................................................................................43 2.3.1 Shape Analysis...........................................................................................44 2.3.2 Intensity Profile..........................................................................................45 2.3.3 Segmentation with Prior Shape Information..............................................47 2.3.3.1 Synthetic image................................................................................48 2.3.3.2 Cardiac ultrasound images...............................................................49 2.3.4 Segmentation with Prior Shape Information and Key Points.....................50 2.3.4.1 Synthetic image................................................................................50 2.3.4.2 Cardiac image...................................................................................52 2.3.5 Segmentation with Prior Shape and Intensity Profiles...............................57 2.4 Conclusion............................................................................................................65 3 APPLICATION OF PDE BASED INPAINTING ON MR PARALLEL IMAGING...................................................................................................................68 3.1 Introduction...........................................................................................................69 3.1.1 Inpainting....................................................................................................69 3.1.2 Sensitivity Maps.........................................................................................69 3.1.3 Intensity Correction Map............................................................................70 3.2 Methods................................................................................................................71 3.2.1 Review of TV Inpainting Model................................................................72 3.2.2 Modification of the Exponent of u in the TV Model.............................73 3.2.3 Modified TV Model for Sensitivity Maps..................................................74 3.2.4 Numerical Implementation Method...........................................................75 3.2.5 Modified TV Model for Intensity Correction Map....................................78 vi PAGE 7 3.3 Experiments and Results.......................................................................................78 3.3.1 Application on Sensitivity Map..................................................................78 3.3.1.1 Phantom with no intrinsic holes.......................................................80 3.3.1.2 Coronal phantom..............................................................................81 3.3.1.3 Neurovascular images......................................................................84 3.3.1.4 Cardiac images.................................................................................85 3.3.2 Application on Intensity Correction Map...................................................88 3.3.2.1 Uniformity scheme...........................................................................88 3.3.2.2 Phantom............................................................................................90 3.3.2.3 Clinical images.................................................................................91 3.4 Conclusion............................................................................................................93 4 MODIFIED MUMFORDSHAH MODEL FOR MEDICAL TOMOGRAPHY......95 4.1 Background of Tomography.................................................................................95 4.1.1 Existing Methods........................................................................................96 4.1.2 Motivation of the Proposed Framework...................................................100 4.2 Proposed Model..................................................................................................101 4.2.1 Description of the Proposed Model..........................................................101 4.2.2 Level Set Forms in Different Cases..........................................................102 4.2.2.1 Binary field model..........................................................................102 4.2.2.2 Piecewise constant model...............................................................103 4.2.2.3 Piecewise texture model.................................................................106 4.2.2.4 General piecewise smooth model...................................................106 4.2.2.5 Model for tomography with guiding image...................................108 4.3 Application on Noise Tomography....................................................................110 4.3.1 Introduction of Noise Tomography..........................................................110 4.3.2 Electronic Noise.......................................................................................110 4.3.3 Intercoil Noise Correlation.......................................................................112 4.3.4 Relation with Tomography Models..........................................................113 4.3.5 Experiments and Results..........................................................................113 4.3.5.1 Experiment in one dimension.........................................................113 4.3.5.2 Experiment in two dimensions.......................................................114 4.3.5.3 Experiment in two dimensions with guiding image.......................116 4.3.5.4 Experiment in two dimensions with guiding image by fast sweeping...............................................................................................117 4.4 Conclusion..........................................................................................................118 5 FIBER TRACKING BY FRONT PROPAGATION...............................................120 5.1 Introduction.........................................................................................................120 5.1.1 Understanding the Basis of DTI...............................................................123 5.1.2 Existing Fiber Reconstruction Methods...................................................125 5.1.2.1 Line propagation algorithms..........................................................126 5.1.2.2 Front evolution methods.................................................................129 5.1.3 Several Considerations in Fiber reconstruction........................................130 5.1.3.1 Termination criteria........................................................................130 vii PAGE 8 5.1.3.2 Effect of noise................................................................................132 5.1.3.3 Branching.......................................................................................132 5.1.3.4 Algorithms to interpolate a tensor field..........................................132 5.1.3.5 Display of experimental tracks.......................................................133 5.1.4 Limitations of Existing Methods..............................................................134 5.2 Methods..............................................................................................................135 5.2.1 Motivations of the Proposed Method.......................................................135 5.2.2 Data Description.......................................................................................138 5.2.2 Front Evolution and Path Record.............................................................138 5.2.4. Tracking Back.........................................................................................140 5.3 Experiments and Results.....................................................................................141 5.4 Conclusion..........................................................................................................143 LIST OF REFERENCES.................................................................................................145 BIOGRAPHICAL SKETCH...........................................................................................156 viii PAGE 9 LIST OF TABLES Table page 21 Experiment results with different for Epi.............................................................59 22 Experiment results with different for Endo..........................................................63 31 The relative error of inpainted results......................................................................81 32 Ghost ratio of coronal phantom................................................................................84 33 Ghost ratio of neurovascular images........................................................................85 34 Ghost ratio of cardiac images by applying highresolution sensitivity maps of first image.................................................................................................................87 35 Ghost ratio of cardiac images by applying lowresolution sensitivity maps of each image................................................................................................................88 ix PAGE 10 LIST OF FIGURES Figure page 11 Level set.....................................................................................................................9 21 Examples of segmentation for cardiac images.........................................................20 22 Cluster and average result........................................................................................45 23 Example of average shapes of abnormal clusters.....................................................45 24 Twenty training images with the segmented epicardiums.......................................46 25 Average images of the above 20 images with the average contour superimposed..47 26 Segmentation with prior shape for synthetic image.................................................49 27 Segmentation with prior shape for cardiac image....................................................50 28 Segmenation with prior shape and key points for synthetic image..........................52 29 Segmentation with prior shape and key points for cardiac image 1.........................54 210 Segmentation with prior shape and key points for cardiac image 2.........................56 211 Segmentation with prior shape and intensity profile for epicardiums.....................61 212 Intensity profile for endocadiums.............................................................................62 213 Segmentation results for endocardiums...................................................................64 31 Inpainting results for the artifical holes....................................................................81 32 Inpainting results for coronal phantom....................................................................83 33 Inpainting results for neurovascular images.............................................................85 34 The flow chart of the intensity correction method...................................................89 35 Intensity correction for phantom..............................................................................92 36 Intensity correction results for clinical neurovascular images.................................93 x PAGE 11 41 Experiment result in 1 dimension...........................................................................114 42 Experiment result in 2 dimensions.........................................................................115 43 Experiment result in 2 dimensions with ideal guiding image................................117 44 Experiment result in 2 dimensions with real MR guiding image...........................117 45 Reconstructed image by fast sweeping..................................................................118 51 White matter structure and its relationship with the information provided by DTI.........................................................................................................................122 52 Linear linepropagation approach..........................................................................127 53 White matter fiber trajectory as a space curve.......................................................128 54 Flow chart of fiber tracking by fast marching on grids..........................................141 55 Fiber tracking results..............................................................................................143 xi PAGE 12 LIST OF ABBREVIATIONS GPR: ground penetrating radar AAM: active appearance models HARDMRI: high angular resolution diffusion magnetic resonance imaging ART: algebraic reconstruction technique BV: bound variation LV: left ventricular CT: computed tomography MAP: maximum a posteriori DTI: diffusion tensor imaging MIIG: mutual information of image geometry DWI: diffusion weighted imaging ECT: emission computed tomography MMSLS: modified mumfordshad model by level set ED: end diastole MRA: magnetic resonance angiography EIT: electrical impedance tomography MRI: magnetic resonance imaging Endo: endocardial MRNT: magnetic resonance image guided noise tomography EPI: echo planar imaging Epi: epicardial NEMA: national electrical manufacturers association ES: end systole NT: noise tomography FBI: fresh blood imaging OSEM: ordered subsets expectation maximization FBP: filtered back projection FIESTA: fast imaging employing steadystate acquisition PDE: partial differential equation PET: positron emission tomography fMRI: functional magnetic resonance imaging RF: resonance frequency FOV: field of view SENSE: sensitivity encoding GKS: Gaussian kernel smoothing xii PAGE 13 SMASH: simultaneous acquisition of spatial harmonics SNR: signal to noise ratio SOS: sum of squares SPECT: single phantom emission computed tomography TE: echo time TR: transmit time TV: total variation xiii PAGE 14 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy APPLICATIONS OF VARIATIONAL PARTIAL DIFFERENTIAL EQUATION MODELS IN MEDICAL IMAGE PROCESSING By Feng Huang December 2004 Chair: Yunmei Chen Major Department: Mathematics Images of various kinds are increasingly important to medical diagnostic processes. Difficult problems are encountered in selecting both the most appropriate image reconstruction methods and postprocessing methods. Reconstruction methods involve producing optimal quality images from raw data. Postprocessing methods involve acquiring the highest quality information from these images. Medical image processing (i.e., reconstruction and post processing) has created tremendous opportunities for mathematical modeling, analysis, and computation. I examined various applications of the variational partial differential equation (PDE) method. This method is one of the most recent and successful approaches to medical image processing. I worked on the following important medical image processing problems: clustering, denoising, segmentation, registration, tomography, inpainting, and field tracking. Novel PDE based models and corresponding numerical methods are introduced for each of these problems. xiv PAGE 15 We proposed an original framework for segmentation that may incorporate various forms of prior information. This information can include shape, key points, or intensity profiles. We also proposed the level set formulation and a numerical algorithm for the model. We proposed a novel inpainting model and applied it to magnetic resonance imaging (MRI). This novel inpainting model considers the automatic choice of diffusion type. A specific modification of this model for sensitivity maps was produced. We also introduced a new framework (the modified mumfordshad model by level set) for tomography problems with extremely noisy and limited data. This algorithm is flexible, it can be used for any geometry with or without prior information. This model is very easy to modify for various cases. Finally, we proposed a novel tracking algorithm. The proposed method solves multidiffusiondirection and branch problems existing in fiber tracking. We applied those models on segmentation of ultrasound cardiac images, MRI brain images, noise tomography, inpainting of sensitivity maps for MRI surface coils, and fiber tracking. Experimental results and computational costs were compared with conventional methods. xv PAGE 16 CHAPTER 1 INTRODUCTION Variational and partial differential equation (PDE) methods have been widely used in image processing in the past few years. Examples include continuous mathematical morphology, shape from shading, invariant shape analysis, image segmentation, tracking, image restoration, stereo, contrast enhancement, and inpainting. The basic idea is to represent an image as an 2 R function. This function usually satisfies a time dependent PDE that characterizes the giving problem. The solution of the differential equation gives the processed image at the scale A PDE can be obtained through direct derivation of the evolution equations such as the classical snakes model by Kass, Witkin and Terzopolous [1]; the anisotropic diffusion model by Perona and Malik [2]; and the image inpainting by Bertalmio et al. [3]. A PDE can also be obtained from variational problems. This is more common in image processing. The basic idea is to minimize an energy functional. t Using variational and PDE methods in image analysis leads to modeling images in a continuous domain. This simplifies the formalism, which become grid independent and isotropic. Conversely, when the image is represented as a continuous signal, PDEs can be seen as the iteration of local filters with an infinitesimal neighborhood. This interpretation of PDEs allows one to unify and classify a number of the known iterated fillters, as well as to derive new ones. Another important advantage of the variational and PDE approach is the possibility of achieving accuracy, stability and high speed with the help of the extensive results of numerical PDE methods. 1 PAGE 17 2 This dissertation shows some applications of variational PDE methods for medical imaging. Chapter 1 presents the motivating applications for our research, reviews medical imaging concepts, and shows numerical techniques used to PDE models. This chapter also proposes the various conclusions and contributions of our study. 1.1 Motivating Applications Partial differential equation based medical image processing has been applied in many fields. Several specific applications motivated our research. Those applications are segmentation of ultrasound cardiac image, Noise Tomography and calculation of sensitivity map of probe coil and fiber tracking. 1.1.1 Segmentation of Ultrasound Cardiac Image The task of determining the epicardial and endocardial borders of the left ventricular (LV) on echocardiography images is essential for quantification of cardifunction. For ultrasound images, it is very difficult because of poor contrast, low signal to noise ratio, information dropout at the borders, and the inhomogeneity of intensities inside and outside the regions between epicardial and endocardial borders. The aim here is to segment the endocardium and epicardium in a twochamber or four chamber image of the heart. The endocardium and epicardium are not distinguishable in the image. Consequently, our task was to find and complete the boundary. 1.1.2 Noise Tomography Noise Tomography, which aims to determine the distribution of the conductivity in the sample, is a new noninvasive medical imaging technique invented by MRI Devices Cooperation. The NT technique is designed to use the correlations in the detected electronic thermal noise in an RF probe array, and to use the relationship between PAGE 18 3 conductivity and the noise power coupled between the sample and probe to measure the electrical conductivity distribution within the sample. Many applications in diagnostic imaging have been documented: gastrointestinal and oesophgeal function [4, 5], hyperor hypothermic treatment of malignant tumors [6, 7] imaging of the head [8, 9] pulmonary function [10], cancer detection [11, 12] measurements of cardiac output [13], locating the focus of epileptic seizures [14], and so on. Noise Tomography generates a threedimensional map of the conductivity of the sample. Since it is based entirely on detection of the thermal noise generated within the body, NT is a completely noninvasive technique; and introduces no external power, chemicals, or radionuclides into the body. However, because of the limitations of collection time and coil design, the data collected for image reconstruction are extremely noisy and limited. This leads to an illposed, big size inversion problem. The aim here is to create an efficient algorithm for tomograpgy with extremely noisy and limited data. 1.1.3 Calibration of Sensitivity Map A sensitivity map deals with the sensitivity information of RF probe coils in an MRI system it does not refer to the state of the object under examination. Knowledge of spatial receiver sensitivity implies information about the origin of detected MR signals, which can be used for image generation. Therefore, samples of distinct information content can be obtained at one time by using distinct probe coils in parallel implying the possibility of reducing scan time in Fourier imaging without having to travel faster in kspace. Parallel imaging techniques, which are based on the fact, have been widely applied since the 1990s of last century. SiMultaneous Acquisition of Spatial Harmonics (SMASH [15]) and Sensitivity Encoding (SENSE [16]) are among the most famous PAGE 19 4 techniques. All of those techniques require highly accurate knowledge of component coil sensitivity information, which is the sensitivity map. Besides this application, the sensitivity map is also crucial to correct the inhomogeniety of surface coil MRI. The raw sensitivity maps are impaired by noise. Straightforward elimination of noise by lowpass filtering results in errors at object edges. Hence denoising is required for raw maps. Some data sets with large areas (such as pulmonary MRA using fresh blood imaging) contribute little or no signal. In this case, there will be some holes, where we have no information of sensitivity at those places, in the sensitivity maps and therefore interpolation techniques are required to fix those holes. Furthermore, to deal with slightly varying tissue configurations and motion, extrapolation over a limited range is also necessary. The aim here is to create a technique to do denoising, interpolation, and extrapolation at the same time for the sensitivity map. 1.1.4 Intensity Correction In the areas of medical imaging, acquired images may be corrupted by poor signal to noise ratio (SNR) or nonuniformities in spatial intensity. Poor SNR can cause blurred or ambiguous feature edges and nonuniformities in spatial intensity. Consequently, structures, textures, contrasts, and other image features may be difficult to visualize and compare both within single images and between a set of images. As a result, attending physicians or radiologists presented with the images may experience difficulties in interpreting the relevant structures. Nonuniformities can hinder visualization of the entire image at a given time and can also hinder automated image analysis. In MR images, the acquired images generally contain intensity variations resulting from the inhomogeneous sensitivity profiles of the surface coil or coils. In general, tissue next to the surface coil appears much brighter than tissue far from the coil. Spatial intensity variations introduced PAGE 20 5 by surface coil nonuniformity hinder visualization because one cannot find a window/level adjustment to encompass the entire field of view. When such images are filmed, the operator tries to select a setting which covers most of the features of interest. Furthermore, uncorrected image inhomogeneity makes it difficult to perform image segmentation and other aspects of image analysis. Therefore, techniques to enhance SNR and correct nonuniformity are necessary for MR images. 1.1.5 Fiber Tracking Diffusionweighted MRI (DWI) adds to conventional MRI the capability of measuring the random motion of water molecules, referred to as diffusion. Water in tissues containing a large number of fibers (such as cardiac muscle and brain white matter). The diffusion is fastest along the direction that a fiber is pointing to, but slowest in the direction perpendicular to it. Water diffuses isotropically in tissues that contain few fibers. DWI renders noninvasively such complex in vivo information about how water diffuses into an intricate 3D representation of tissues. Conventional MRI can be used to investigate spatial relationships among different anatomical regions, but it is unable to infer the connectivity of these regions. DWI provides profound histological and anatomical information about tissue structure, composition, architecture, and organization. Changes in these tissue properties can often be correlated with processes that occur in development, degeneration, disease, and aging, so this method has become more and more widely applied [1720]. Because DWI provides directional information concerning microscopic tissue fiber orientation at the voxel scale, fiber tracking (or fiber reconstruction) provides connection information between voxels based on the information from DWI. PAGE 21 6 1.2 Review of Conception 1.2.1 Description of General Tomography Problem The goal of tomography is to reconstruct object profiles and extract object features given a set of tomographic measurements. In general, the problem of tomography is the problem of recovering (at least approximately) the values of a function from (possibly approximate and/or possibly weighted) an integral of the function over some subset of its domain. I.e. given the integral nidxxfxRMii,...,2,1, and try to recover function on the area In medical tomography problems, is the collected data, andis some fixed function determined by the character of the machine. Integral domain can be different dimensions in different applications. xf iM iR 1.2.2 MumfordShah Model The MumfordShah model was formulated by Mumford and Shah in 1989 [21], originally for an image segmentation problem. The segmentation problem can be defined as follows: Let N R be the domain of a given observed image be the edges of find a decomposition of 0u 0u i and an optimal piecewise smooth approximation uof such that varies smoothly with each 0u u i and rapidly or discontinuously across the boundary of Mumford and Shah proposed the following energy functional to solve it i \220],[dxudxdyuuuE (11) where 0, are fixed parameters, to weight the different terms in the energy. The first term is called the fidelity term. It penalizes the difference between the smooth image and the original image. In the second term, u is the norm of the gradient of u. This term 2L PAGE 22 7 protects the smoothness of u outside the boundary In the last term,  refers to the length of the boundaries in the image. For  ,u a minimizer of the above energy, is an optimal piecewise smooth approximation of the initial, possibly noisy, image will be smooth only outside Theoretical results of the existence and regularity of minimizers of this functional can be found in [21]. Later, Rudin, Osher and Fatemi [22] solved this problem by letting the power of u0 u u u be 1, which makes sure that the smooth is along the level set and therefore preserve and enhance edges. u20 u 0 u0 u dxdy 1.2.3 Total Variation Denoising Total Variation (TV) denoising method is a PDEbased technique that remove noise from images. This method was first introduced by Rubin, Osher and Fatemi [22]. The basic idea is to minimize the following functional: dxdyudxdyu2min (12) The TV functional does not penalize discontinuities in u and thus allows us to recover the edges of the original image. The associated EulerLagrange equation of Equation 12 is 0uuuu (13) Since u is in the denominator, in order to avoid the singularity, it is common to use a slightly perturbed norm 2uu where is a small positive constant, to replace u This is equivalent to minimize the functional udxdy22 (14) PAGE 23 8 Acor and Voger [23] showed that the solutions of the perturbed problems Equation 14 converge to the solution of Equation 12 when 0 As an important improvement of TV functional, Strong and Chan [24] introduced the weighted TV functional dxxuxTV (15) for spatially adaptive (selective) image restoration. The function is chosen so that is larger away from possible edges and smaller near a likely edge. Hence we allow for greater smoothing away from edges and less smoothing at the edges. Certain choices of (x) were given in [24] and [25]. Numerical results are very promising. Chen and Wunderli [26] prove that there exists a unique solution for the weighted TV functional. Later, Chen and Levine improves the model [27]. They change the power of u to be a function which depends on )(xP u This modification makes that one to control the type of diffusion at each location in the image. This model works better for piecewise smooth images and those in which the noise and the edges are extremely difficult to distinguish. 1.2.4 Level Set Level set methods, introduced in 1988 by Osher and Sethian [28] are numerical techniques designed to track the evolution of interfaces. Rather than follow the interface itself, the level set approach takes the original curve (the red one on the left below), and builds it into a surface. That coneshaped surface, which is shown in bluegreen on the right below, has a great property; it intersects the xy plane exactly where the curve sits. The bluegreen surface on the right below is called the level set function, because it PAGE 24 9 accepts any point in the plane as input and hands back its height as output. The red front is called the zero level set, because it is the collection of all points that are at height zero. Figure 11. Level set The idea of moving the front is that instead of moving the red front, which can do all sorts of weird things, one might try and instead move the surface on the right. In other words, the level set function expands, rises, falls, and does all the work. To find out where the front is, get out the saw and cut the surface at zero height. The reason this is a good idea lies in taking a good look at the surface on the right. It will always be "nice"; the red front can get wildly contorted, but our blue level set function is wellbehaved. All the complicated problems of breaking and merging are easily handled. Better yet, everything works for three dimensional surfaces with no change. An equally important advantage is that it is easy to build accurate numerical schemes to approximate the equations of motion. That is because rather than track buoys around, which might end up colliding, one can instead stand on the xy plane and compute the answer. Using terminology from a completely different language, marker particle methods are mantoman coverage, while level set methods are a zone defense. All PAGE 25 10 together, the trick of embedding the front in a higher dimensional function is well worth the added cost (in fact, with some work, that cost may be made the same as that of marker techniques). A typical example of level set function is given by the signed distance function to the curve. In our case, level set is an effective implicit representation for because it allows for automatic change of topology, such as merging and breaking, and the calculations are made on a fixed rectangular grid. Hence, the level set method is applied here for numerical implementation. 1.2.5 Distance Function Given a compact set in space, the distance function for at a point x is defined as the smallest Euclidean distance of x to the points in Fast and accurate computation of the distance function is an important step in many applications. For example, level set methods [28] rely heavily on the distance function to represent an interface, and in particular to keep the level set function well behaved during time evolution. In the surface interpolation model of Zhao et al. [29], a surface interpolating a given points set is constructed by first computing the distance function to In the Island Dynamics model of Chen et al. [30] in materials science, one needs the distance function to help understand the aggregation of atoms into islands. Generally, applications in physics, chemistry, molecular biology, and even urban planning that use point pattern analysis require construction of the distance function. Furthermore, in applications such as ray tracing and surface rendering where one needs information about the normals or other geometric quantities of the surfaces of the given objects, the distance function is essential. PAGE 26 11 There are many different methods for finding distance function. The method chosen for our study is the closest point method [31]. To generate signed distance function, we need to know a point is inside or outside a curve. Tsai [31] introduces one method. One can also use the Jordan curve theorem: Run a semiinfinite ray vertically up from the test point, and count how many edges it crosses. At each crossing, the ray switches between inside and outside. This method is very fast and easy to implement. 1.2.6 Mutual Information Mutual information [3234], or simply information, was introduced by Shannon in his landmark 1948 paper ``A Mathematical Theory of Communication.'' It is a measure of the statistical dependency between two random variables based on Shannon's entropy. Lets first review some probability theorems. Given a sample space and events S E and in the probability of F S E is written EP and the probability of is written The joint probability of F F P E and is the probability that both occur, The F and F E P E and are independent if and only if The F F FPP E and E P FEPFPEor and PF E P The conditional probability of E occurring, given that has occurred, is F FP Fand E P F E P / Bayes Rule: If nE,...,1 E are n mutually exclusive events whose union is the sample space S, and E is any arbitrary event of with nonzero probability, then S EPEPk Ek/ EP E EP k Given random Nvector X with density xp the Shannon entropy is defined as NRdxxpxpXHlog X H is a measure of the uncertainty on the values X may assume. Let's now consider a second random variable Y with probability distribution PAGE 27 12 yp then NNRRdxdyyxpyxpYXH,log,, is the joint entropy of X and Y and is the joint density function. The yx, p YXH, is a measure of the uncertainty on the values X and Y may jointly assume. The YXH, is maximum if X and Y are statistically independent (completely uncoupled): in this case HHYXH Y X The is minimum if X, Y H X is a function of Y (deterministic coupling): in this case Y YHXH XH, X YXH, YXHYHXH, Y, XMI X X Y, YHXH XMI X The definition of Mutual Information between and Y naturally follows from the properties of Mutual Information MI is defined as The MI is clearly a measure of the statistical coupling between and Y The lowest possible value of MI is 0, and this happens if and only if and Y are independent (uncoupled). The highest possible value of MI is and this happens if and only if X is function of Y (deterministic coupling). An important property of MI is that it is invariant to any strictly monotonic transformation of and Y 1.2.7 Inpainting The concept of digital inpainting was first introduced into the image processing community by Prof Sapiro's [3] group at the Department of electrical and computer engineering, University of Minnesota, as inspired by the manual inpainting process of real restoration artists and museum conservators. In nature, inpainting is an image interpolation and restoration problem. Based on the applications, there are numerical methods for inpainting, such as PDE base transportation inpainting [3], jointly continue/interpolate levellines (geometry) and gray values (photometry) in a smooth PAGE 28 13 fashion [35], curvature driven diffusions (CDD), Bayesian framework for image restoration, total variation inpainting [36], elastica inpainting [37], MumfordShahEuler inpainting [38], and so on. 1.3 Numerical Techniques Some numerical techniques are introduced in this section. These numerical techniques are applied to all PDE models in chapters 2,3, and 4. 1.3.1 Choice of Heaviside Function There are several ways to define the Heaviside function According to Chan and Vese [39], the following definition is suggested H zzHarctan2121 ,and 0,'H is defied as 221uu 1.3.2 Computation of uudiv Since u is in denominator, singularity may happen. To avoid singularity, several approaches can be applied. One way is to add a small positive constant in the denominator; however how to choose this constant is a problem. Another approach was introduced in Blomgren et al. [40]. The idea is to let uuw and solve the system wdivuuw to avoid the singularity. PAGE 29 14 Chan and Vese [39] introduced the approach, which was used in our experiments. The following difference schemes are applied ),(,)(2121uuuhuuudivyx huuunijynijxnij11, jijiijyjiijijyjijiijxjiijijxuuuuuuuuuuuu,1,1,,,1,1,,, then ijjijijijiuCCCCuCuCuCuChuudiv43211,41,3,12,1121 where h is the step length. ijaC11 jiaC,121 ijbC13 1,41jibC 2121,1,22huuhuanjinjinijxij 212,1,122huuhubnjinjinijxij 1.3.3 Initial Signed Level Set and Reinitialization Because the energy which is minimized is not convex and also there is no uniqueness for the minimizers, the algorithm may not converge to a global minimizer for a given initial condition. In general, for real images with complicated features, the initial conditions with small initial curves were suggested by Vese and Chan [41]. This is for the following two reasons: first it has the tendency to converge to a global minimizer; and second, the algorithm is much faster. Closest point method by Tsai [31] was used to compute the initial signed level set During the iteration, reinitialization for is necessary at each step, because we work PAGE 30 15 with the regularized function This procedure is standard [39, 42] and prevents the level set function from becoming too flat, or it can be seen as a rescaling. H div 1.3.4 Fast Sweep Method for Level Set Based Optimization Recently, the fast algorithm for level set based optimization [43] was introduced by Bing Song. Applying this method to solve a class of optimization problems with level set representation can improve the computational speed dramatically. Because this method does not solve the EulerLagrange equation and only uses the sign of level set, the implementation is much simpler than the numerical method introduced above. This algorithm is defined for HFmin where H is the Heaviside function; is any functional dependent on F and is the level set function. The outline of the algorithm is as follows: Step 1: Initialize. Construct an initial partition (one part for 0 one part for 0 ) and compute the value of F according to Step 2: Advance. For each point x in image, if the energy F decreases when we change x to x then update this point by xx ; otherwise, x remains unchanged. We sweep the pixels in some prescribed order. For example, in image segmentation, we can sweep the pixels row by row. We can use either GaussSeidel or Jacobi iteration in each sweep. Step 3: Repeat step 2 until the energy F remains unchanged. Because there is no uu term and we do not need to reinitialize the level, computation speed can be dramatically reduced. PAGE 31 16 1.4 Contributions of Our Study 1.4.1 Segmentation with Prior Information Chapter 2 proposes a framework for segmentation that incorporates different prior information. The prior information could be shape, key points or, intensity profile. We also propose the level set formulation and a numerical algorithm for the model. According to the available prior information, this model can easily be modified for the specific priors. Several examples for different prior information are given in the experiment section. The proposed model has five advantages. First, with the prior shape information, this model can reliably segment images in which the complete boundary was either missing, or was low resolution and low contrast. Second, besides segmentation, the algorithm also provides registration estimations of translation, rotation, and scale that maps the active contour to the prior shape. Third, with the help of key points, the problem (which is to determine local shape variations from the average shape) can be improved. Fourth, we used maximizing MIIG rather than MI to match intensity profiles of two images. We did this because that the MIIG takes neighborhood intensity distribution into account and hence gives better a description of the intensity profile than MI does. Fifth, the parameter in the model of active contour with shape can be selected by the model itself, and the segmentation with the optimized parameter arrives at higher image gradients, forms a shape similar to the prior, and captures the prior intensity profile. Experiments on synthetic and clinical images show those advantages. 1.4.2 Inpainting A novel inpainting model was proposed and was applied to MRI. This work made the following contributions: Introduced the inpainting technique into the field of MRI PAGE 32 17 Introduced a novel inpainting model that considers the automatic choice diffusion type Provided a fast and robust technique for calculation of a RF coil sensitivity map Provided a realtime and robust technique for simultaneous correction for nonuniformity and preservation for SNR. All of the major issues involved in the correction of sensitivity maps were effectively and simultaneously dealt with: denoising is accomplished, holes have been fixed, and the original map has been smoothly extended to the entire image. By the experiments with intensity correction maps, it can be concluded that, with the use of inpainted intensity correction maps, the inpainting method is capable of correcting nonuniformity while perfectly preserving SNRs. 1.4.3 Modified MumfordShah Model for Tomography Chapter 4 introduces a novel framework MMSLS for the tomography problem with extremely noisy and limited data. This algorithm is flexible. It can be used for any geometry with or without prior information. It is very easy to change the model for different cases. Furthermore, this approach has no special requirements. Hence it may be easily applied to a variety of tomography problems. The main contributions of this work include providing a solution for medical tomography problem with extremely noisy and limited data; introducing a model for tomography with guiding image; introducing a novel noninvasive technique for conductivity map, which can be used for cancer detection and much more. The MMSLS algorithm has proven itself to provide accurate results. Because the MMSLS algorithm is object based, not pixelbased, this approach, for an equivalent output resolution, is much faster than conventional methods. Moreover, the reconstructed image, using MMSLS, is smooth, and the edge is perfectly preserved. The visualization PAGE 33 18 is therefore much better than when using other approaches. Another impressive advantage of this approach is that it can accomplish segmentation and reconstruction at the same time. 1.4.4 Fiber Tracking A novel tracking algorithm is proposed. The proposed method solves multidiffusiondirection, branches and tracking back problems existing in fiber tracking problems. There are three major differences between this work and other fast marching based methods. First, local polynomial fitting is used to define evolution speed. Second, the closest point method is applied to find the most reasonable path to connect two points. Third, the front propagation method is modified to deal with multidiffusiondirections in some voxels. Hence this method is more efficient way to partially solve multifiber problem. PAGE 34 CHAPTER 2 SEGMENTATION WITH PRIOR INFORMATION In this chapter, we propose a new variational framework for image segmentation that incorporates with prior information. Models are created for segmentation problems with several different kinds of prior information like the following: expected shape, expected shape and a few points on the boundary, expected shape and intensity profiles. This chapter is organized as follows. In the introduction section, the background of segmentation and review of existing segmentation methods are introduced. Then, in the method section, a PDE based segmentation model is proposed with numerical algorithms. In experiment and result section, we show experimental results of the application of the proposed method to synthetic images and cardiac ultrasound images and discussion of the proposed framework ends this chapter. 2.1 Introduction 2.1.1 Definition of Segmentation Segmentation is the process of identifying regions of pixels in an image. The main goal is to divide an image into parts that have a strong correlation with objects or areas in the real world. We can use the gradient of the intensity (edges), or gray level similarities, or other features in the image to perform the segmentation. 2.1.2 Importance of Segmentation in Medical Applications An example of the usefulness of this process is performing measurements of cardiac function (Figure. 21). The object for this image might be to determine the area of blood in the left ventricle. 19 PAGE 35 20 A B Figure 21. Examples of segmentation for cardiac images. A) a 2 dimensional example, B) a 3 dimensional example The region is shown in the left of the above image. This region can be drawn by hand. If images are obtained throughout the cardiac cycle and over the entire left ventricle, then the amount of blood in the left ventricle can be determined at end systole (ES) and end diastole (ED) from which the percentage of blood that is ejected on one heart beat can be determined (ejection fraction). This is an important measurement to cardiologist but would require intensive interaction with a computer by an expert in order to determine the blood volume on a typical data set. These data sets usually contain over 200 images. You might expect that methods to automate this process are both desired and necessary. In addition to obtaining the stroke volume, if the outside wall of the myocardium, epicardium, can also be obtained, then the thickening of the left ventricular wall could also be obtained throughout the cardiac cycle. The differences in thickening over the left ventricular wall can show regions that are not functioning well due to decreased blood flow. Applications also include Multiple Sclerosis (MRI), MR Angiography, Brain Tumor (MRI) and Neuro Applications. PAGE 36 21 2.1.3 Review of Segmentation Methods A large number of segmentation algorithms have been proposed. For the different applications and for the different image types, there are different types of method for segmentation. According to the information used for segmentation, there are three classes of segmentation methods: global knowledge, edgedbased and regionbased. Global knowledge methods rely on knowing something about the image such as the expected pixel intensity. This is possible in CT where the intensity values are in Hounsfield units. Thus, one can look for all pixels of value near 808 for bone. Typically, threshold techniques are used for this process. Most of the methods based on global knowledge of the image are thresholding techniques. Thresholding techniques are effective when the intensity levels of the objects fall squarely outside the range of levels in the background. Because spatial information is ignored, however, blurred region boundaries can create havoc. Edgebased methods find the borders between regions. The edgebased methods [1, 4448] rely on the information of the edges, such as high magnitude of image gradient. Edgebased methods center around contour detection and their weakness in connecting together broken contour lines make them, too, prone to failure in the presence of blurring. Regionbased methods may be more applicable when edges are difficult to determine because of noise. The regionbased methods [21, 39, 4952] make use of homogeneity on the statistics of the regions being segmented. We can use many different properties of the image intensity to define regions such as gray level, color, texture and shape to name just a few. Hybrid techniques using a mix of the methods above are also popular. PAGE 37 22 In numerous medical imaging modalities, the boundaries of anatomical structures of interest may be very difficult to capture. The reasons for this difficulty includes low contrast, low resolution, or even that a portion of the object is not in the field of view. For example, in echocardiography the images are typically degraded by noise and signal drop out. In some cases, an image sequence may even have portions of the myocardium that lie outside of the sector scan so that some segments of an object boundary may not be visible at all. Thus, the segmentation problem for these images is very challenging and cannot be successfully segmented by above methods. In an effort to overcome these difficulties, various techniques have been developed to incorporate prior information of an objects expected shape into the segmentation process. For example, Cootes, et al. [53] and Wang, et al. [54] constructed a statistical model of an objects shape from a set of corresponding points defined on a set of training images and then incorporated this information into a Bayesian formulation to find an objects boundary. Cootes, et al. [55] also used a Gaussian model to fit a training set of corresponding feature points. Later Cootes and Taylor used nonGaussian mixture models to fit data for specific applications [56]. Staib and Duncan [57] specified the shape of the curve by creating statistical priors on the Fourier coefficients of the contour. This prior was incorporated into a segmentation process by a posteriori technique. Yuille, Hallinan, and Cohen [58] and Tagare [59] used a deformable template approach to incorporate shape information into the segmentation process. Recently, statistical shape knowledge has been incorporated into edge based or region based active contours. Leventon, Grimson, and Faugaras [60] and Chen et al. [61] incorporated prior shape information into the geodesic active contour model. Chen et al. [62] used a principal PAGE 38 23 component analysis technique to build a statistical shape model from a training set of curves, where each curve is represented by a signed distance function. A force depending on the shape model was added at each iteration of the levelset surface evolution for the geodesic active contour to force the active contour towards high image gradients and a maximum posteriori estimate of shape and pose. In Tsai et al.[52] the shape model was obtained by averaging the aligned training samples. Their model finds the segmentation transformation which maps the evolving interface to the shape prior by minimizing an energy functional containing the energy of geometric active contour and a measure of the shape disparity between the active contour and the prior. The disparity is measured by using distance function. A number of research groups have incorporated a statistical shape model into a regionbased segmentation scheme including Cremers, Schnorr, and Weickert [51], Tsai et al. [52], and Paragios and Rousson [63]. In the Cremers, Schnorr, and Weickert approach, an energy functional is minimized that includes both a shape energy term corresponding to the Gaussian probability as well as a MumfordShah [21] energy term. Tsai et al. build a parametric shape model similar to the one described by Leventon, Grimson, and Faugaras [60], where parameters are calculated to minimize a region based objective function which provides the segmentation. The shape model of Paragios and Rousson was formed using the level set representation and the variability of the shape. It was designed to maximize the loglikehood function, which is represented by a Gaussian density function where the mean corresponds to the level set representation of the shape and the variance refers to the variation of the aligned samples. Then, the shape model was combined into a geodesic active region model [64] in a variation framework. Leventon, Faugauras, Grimson, and Weikert used a training set to extract PAGE 39 24 prior information of the intensity and curvature profiles of the features. This information was incorporated into the segmentation processing by using an approach similar to the one developed in Leventon et al. [60]. 2.1.4 Remaining Problems in Existing Methods While experimental results have shown the effectiveness of priorbased models in numerous medical applications, many problems remain including the complexity and variability of the images, the accuracy of the measurements obtained, and the rapid computation time required by the user. One practical problem is how to determine the parameter that balances the influences from image information and priors. If the evolution of an active contour is mainly governed by the force depending on the image gradient, it may be sensitive to the initial step or may leak through the boundary where the edge feature is not clearly defined. Conversely, if the force depending on the shape prior is the dominating term, the active contour may not arrive at the boundary of the object of interest even though it has a shape similar to the prior. 2.2 Method In this section, a general varational framework for segmentation with prior information is proposed in 2.2.2. The method for generation of required prior information is explained in 2.2.1. 2.2.1 Generation of Prior Information As we mentioned above, the segmentation problem could be very difficult. Therefore some prior information may be necessary. In this section, the methods for generation of prior information based on training set are discussed. PAGE 40 25 2.2.1.1 Shape analysis To get the prior shape information, an average shape needs to be purified from a set of training shapes. A modified procrustes method is applied here to cluster shapes into groups, and each group generate an average shape. The study of Procrustes shape analysis was provided a solid mathematical foundation in the seminal work of D. Kendall [65], where the contours are considered to be vectors in a high dimensional complex projective space. A comprehensive introduction to the topic can be found in the text by C. Small [66]. When the heart is imaged from the parasternal shortaxis view, it has a simple geometric shape that can be reasonably modeled by a continuous contour written as the union of 4 elliptical arcs [67]. So Procrustes shape analysis can be used in echocardiographic images analysis. The shape of an object is defined in a mathematical context as all the geometrical information that remains after location, scale and rotational effects are filtered out; i.e. an object's shape is invariant under the Euclidean similarity transformations of translation, scaling and rotation. For analysis purposes shape is described by a finite number of points on each object's surface, which are called landmarks. Given two shapes A,B. The first step is to define landmarks of those two shapes. In order to apply procrustes, the center of the parameterized contours are the origin. And those two contours have same number of points, all the angles that between radials connecting origin and the point are same. In standard procrustes, the distance between two contours A,B are defined as Frobenius norm of AB, and when try to minimize the distance by processes of rotation scaling and translation, each step is processed separately. PAGE 41 26 In this modified procrustes, the distance is defined as the area difference and rotation, scaling and translation are processed at the same time. The following definitions and theorems are for 2 dimensional case. Definition 1 Let A, B be two contours, then the area between those two contours are defined as a(A, B)= the area of (A B A B). (21) Definition 2. The shape distance between two contours is defined as the minimum Area between two contours after rotation, transfer and scaling. 21,AA AD ()=min 21,AA )),((21,,ATRAaTR (22) Where is the scalar parameter, R is the rotation matrix, T is the transfer vector. and R can be combined into one matrix Let T, then there are 4 unknowns to get the shape distance. Theorem 1 shows how to solve those 4 unknowns. abba 21tt Theorem 1. Let 2*)(nijaA 2*)(nijbB then achieves at 21,1,122,2,1)()(iiiiibbbbc )) ,*((min,BTARareaTR abbaR , When 21ttT 211111211211111212221121112221*00)(00)(ttbacacaccacacacacaacacacaacniiniiiniiiniiniiiniiiniiiniiiniiiiniiiniiiniiii niiiniiiniiiiiiniiiiiibcbcababcababc12111211212211)()( (23) PAGE 42 27 In the proof of this theorem, ),*(BTARarea was approximated by niiiiiibabtbaaac1121121()(( iibtaa2222)) Given a set of parameterized shapes, a basis is randomly chosen from those shapes. Applying Equation 23 to minimize the distances between the basis and all other shapes by rotation, translation and scaling. Then the average shape is the equally weighted summation of those calibrated shapes. This definition of average curve is simple but not accurate. A better way is to apply distance function here to find the curve, which minimize the summation of distances between the curve and all other curves. The clustering method is named as near basis method. Given a set of shapes, randomly choose one as the basis of first group, minimize the shape distance between this basis and all other shapes by applying Equation 23. All the shapes that have distance less than a certain threshold will be set in the same group as the basis. The threshold hold can be set as a value (<1) multiply the mean of all distances. Then randomly choose one shape from the given shapes that are not in the first group as the basis for the second group. Repeat until all shapes are in some group. Now there are several groups. Compute the average of all of those groups. Then we try to reset the groups according to the averages. Minimize the distance between each contour and each average. A matrix, which has size number of contours multiply by number of groups, is generated. Each curve is reset in the group that has nearest average to it. Repeat group resetting until convergence is approached. Here is no proof for this convergence, but experimentally it converges very quickly. An alternative method we used to create the prior shapes was the selforganizing maps algorithm. If we want to group the training contours C into n nii,...,1, k PAGE 43 28 clusters (say ), we first take three arbitrary contours as the initial contours, denoted by At t iteration, randomly select a contour denoted by 3k3,2,1 0jmj 1 1 tX from the training set, and compare the disparity in shape between 1 tX and each of To do this comparison we first align 2,1jtj 3, m 1 tX to each and denote the aligned by tmj 1tX jjjTtXR11 jtX ~ Then we compute jjtXa jA m,1 ~ defined in (2.1). Suppose ( is the smallest number in 1A 3,2 ,1, jAj We keep and unchanged, and update t m2 t m 3 t1 mby tm11 tX1 t t1 mt m1 ~ 1 where t is a smooth function of ,and decreases to zero as t.After a large number, say of iterations, three average shapes m are generated. Then three clusters are formed by the curves that are closest to the average shapes. The closeness is again measured by the measurement in (2.1). t m N 1jip, 3,,..., 21 NjCimi,...,1 Ii, p i ,pp 1,0* I *C iI* I mi,...,1 Ii, I 2.2.1.2 Intensity profile In this section we introduce our method for generating a intensity model ( average intensity profile across the average shape ) from a set of training images. Let be a set of training segmentations in one cluster, and be the set of images associated with C. Let also Cbe the average shape in the cluster. Our task is to generate an average intensity profile across from the training images If us created just by averaging all the sample images can be blurred, due to the large intensity variations in the training images. To reduce the influence from the outliners we select a subgroup from PAGE 44 29 miIi,...,1, in which the disparity in intensity profiles are relatively small. To do this, we first align each contour mipi,...,1, C to Cby a similarity transformation iiiTR,, that minimizes iCC,* a defined in section (2.1). 0 *,0Cdx 0 *I I x1jC aITxjji1i 11iiiaITxR 1j H YHXHY Next, we examine the disparity in the training intensity profiles across the training segmentation as follows. Let xV (24) be a neighborhood of C. If two images andare related approximated by i j bTRRIjjiii11 (25) for some constants aandb, we may define the closeness/disparity measurement in the intensity profiles near the segmentation Candby 021VjjjidxbTxRI The integral is over Vrather than the entire image domain, since the intensity profiles near the segmentations are more meaningful. If the relation (25) is not valid, maximizing mutual information has been proven to be effective in solving matching problems, in particular in matching multimodality images [33, 34, 68, 69]. One of the advantages of using mutual information is that it does not require an explicit function that relates two images, but only assumes that aligned images explain each other better than when they are not aligned. 0 The mutual information between two random vectors X and Y is defined as YXXMI, (26) where PAGE 45 30 NRdzzpzpZHlog (27) is the Shannon entropy of a random Nvector Z with density zp and NNRRdxdyyxpyxpYXH,log,, (28) is the joint entropy of X and Y and yxp,iI is the joint density function. Consider intensity at each voxel as the realization of a random variable, and all the random variables in an image have the same distribution. By using (26)(28), the common mutual information (MI) of two images and on a region V can be computed by jI 0 2021212121,,log,,,RjiVdidiigpifpiigfpiigfpIIMI (29) where iiiiTxRIxf11 jjjjTxRIxg11 (210) gfp, ,and are computed over V. fp gp 0 Note that if the locations of two points in an image are switched, then the intensity profile is changed, but the density function remains the same. In order to match intensity profiles of two images, we propose to maximize the mutual information of image geometry (MIIG) rather than mutual information (MI). By MIIG of and on V, we mean that iI jI 0 dmdnnxGpmxFpnmxGxFpnmxGxFpIIMIIGRjiV,,log,,,100 (2.11) where hxfhxfxfhxfhxfxF2,,,,2 hxghxgxghxghxgxG2,,,,2 PAGE 46 31 with and f g are defined in (2.10), 54321,,,,iiiiim 109876,,,,iiiiin , and 5421didididm 3di di 109876didididididn Since MIIG uses neighborhood information, MIIG of two images gives better description of the closeness of the intensity profiles of these two images than MI. We now describe the construction of the average intensity profile I across Ideally, we would like to generate *C I by 111*,max0*jjVIIIMIIG (212) where jVIIMIIG,*0 is defined as in (2.11) with xIf* and jjjjTxRIg11 Since this formulation is computationally intensive, we restricted ourselves to the cases where 11111*jjjjjjTxRIaI (213) The weights were determined by using Equation 212. ja 2.2.2 Framework Description To get a desirable segmentation in these cases, prior information of shape, intensity and key points besides the image itself can be applied in this framework. There are 4 terms in the proposed vartional framework. Gradient term shape term, key points term and intensity term. Several notations are defined for the following sections. I : The image which is to be segmented I : The average image generated according to Equation 213 : The image domain of I C : The segmentation in I which is the unknown. PAGE 47 32 *C : The average curve generated according to section (2.2.1.1) u : A Lipschitz function level set such that Cis the zero level set and 0xux is the set inside C zH : the Heaviside function, that is 1 zH if and if 0z 0zH 0 z and zH' (in the sense of distribution) be the Dirac measure. 2.2.2.1 Gradient term Gradient term is general in geometric active contours methods that are widely used for segmentation, that is, for finding object boundaries in images. Geometric active contours are using initialized curves, which dynamically evolve in the image to minimize their energy (or length in some conformal metric). Their energy is designed in such a way that its minimum occurs when the trace of the curve is over points of high gradient in the image. Because such points often define object boundaries, the active contour becomes stationary at the boundary. Let be a parameterized curve with a regular parameter 1,0ppC p To regard the curve as an active contour, we associate energy C CE with it: 10'dppCpCIgCE (214) where Ig can be chosen as 2*11IGIg (215) where is a parameter, and 241xexG so that Ig is strictly positive in homogeneous regions and near zero on the edges [1, 46, 47]. Level set methods are PAGE 48 33 extensively used in active contour evolution because they allow the curve to develop cusps, corners, and topological changes. Notice that the length of the zero level set of uin the conformal metric dppCgds' can be computed by uguuHg Hence the Gradient term in level set form is uIguuE (216) The geometric active contour evolves to minimize this energy and stops when its trace is over points of high gradient. In many images, object boundaries cannot be detected by seeking pixels of high image gradient. Regionbased methods that make use of homogeneity of the statistics of the local features and properties have been developed for such images. Mumford and Shah [21] proposed a variational framework for segmenting an image into homogeneous regions. Zhu and Yuille [49] proposed the region competition algorithm. Their algorithm combines the attractive geometrical features of snake/balloon models and statistical technique of region growing, and was derived by minimizing a generalized Bayes/MDL criterion using the variational principle. The regionbased methods proposed in [39, 50, 51, 70] are an active contour based on minimizing the MumfordShah functional [21]. Recently, the Geodesic Active Region models [64, 71, 72] were proposed. These models combine boundary and regionbased segmentation models. The evolution of the active contour is influenced by a region force that optimizes the segmentation according to the expected intensity of different regions, and a boundary force that contains information regarding the boundaries between the different regions. PAGE 49 34 2.2.2.2 Shape term There is a serious practical problem with active contours they do not become stationary at a boundary if the boundary has segments, which have low gradient. A geometric active contour will often leak through such gaps in the boundary. The problem is that the standard geometric active contours do not have any information about how the gaps are to be bridged. One solution to the problem is to incorporate into the active contour with some prior information about the expected overall shape of the boundary. Then the active contour can compare its shape with the expected shape and bridge the gaps in a meaningful way. A statistical model of shape variation was constructed from a set of corresponding points across the training images [53] [54]. This information was used in a Bayesian formulation to find the object boundary. In [73] a Gaussian model was fit to a training set of corresponding feature points. In [56] mixed models were used to fit to the data for specific applications where the distributions are nonGaussian. In an alternate approach in [57], Staib and Duncan specified the shape of the curve by creating statistical priors on the Fourier coefficients of the contour. This prior was incorporated into segmentation processing by a posteriori technique. In this work, we propose a shape term for incorporating prior shape knowledge in the active contour. We modify the energy function of the contour so that the propagation velocity of the active contour depends on the image gradient as well as the prior shape. The velocity is designed such that the propagation stops when the active contour arrives at high gradients and forms a shape similar to the shape prior. The shape prior is obtained by a modified procrustes method (Section 3.2.2). The modified energy function gives a satisfactory segmentation in the presence of gaps, even when the gaps are a substantial PAGE 50 35 fraction of the overall boundary. We specify the entire prior shape, because a priori, we do not know where the low contrast gaps will occur in the boundary. Our model differs from the model in [60] because we use a variational approach instead of a probability approach. That makes it possible to discuss the existence of the model solution in BV (bounded variation) space. _ To begin the mathematical description, we first specify when two curves have the same shape. Two curves Cand have the same shape, if there exist a scale 1 2C a rotation matrix R (rotation by an angle ), and a translation vector T such that coincides with Let be a curve, called the shape prior, representing the shape we expect of the boundary, and let 1C T RCCnew22 *C Ig be the function defined in Equation 215. To capture the shape prior we would like to find a curve Cand the transformation *C R T ,such that the curve and Care closely aligned. Therefore, to measure the closeness of to we introduce the shape energy term T RC Cnew*C C 102',,,dppCTpRCdTRCE (217) where yxCdyxd,,,* is the distance of the point yx, from In the numerical computation, this distance function is obtained by the fast marching method proposed by Sethian[74]. In level set form, we need measure the similarity of the shapes between the zero level set of uand which can be evaluated by *C *C dxTd2 u Rx Since dpp Cd'2 can be computed by uuHd2 ud2 we have the shape term in level set form uTRxduTRuE2,,, (218) PAGE 51 36 2.2.2.3 Key points term The experimental results have shown that gradient term with shape term provide very promising results in various applications. However, due to the accuracy and efficiency requirement of medical image analysis problems, and the complexity of the medical image, as well as the variability of the anatomic shapes of interest in medical images, how to use shape prior in getting a better segmentation is still a challenging problem. One of the most difficult problems is to determine local shape variations from the average shape. One solution may be using nonrigid registration to assist segmentation [75, 76], but this requires better region or edge information in the image to find the velocity field, and also computationally costs. The aim of key points term is to find a better way to locate boundaries of interest, that has significant signal loss and relatively larger shape variation from the shape prior. In particular, Cardiac ultrasound images exhibit signal drop out so that some segments of the boundary of the myocardium may not be visible in the image. Also ultrasound images often have a limited imaging range and the boundaries are not completely imaged. Moreover, the shapes of the boundaries vary from average shape, that can be constructed from a set of training shape. In particular, in abnormal heart image, the shape distortion of the myocardium from the average shape cant be ignored. However, several points the boundary pass through can easily be given by experts. Therefore we extend the segmentation algorithm developed in Chen et al. [62] by incorporating a few key points on the boundary of interest in addition to shape prior into geometric active contours. We will use the idea of matching nonequivalent shapes by the combination of a rigid transformation and a pointwise local deformation [63, 77] in PAGE 52 37 our modeling. Soatto and Yezzi [77] view a general deformation as the composition of a finite dimensional group action (e.g. rigid or affine transformation) and a local deformation, and introduced a notion of shape average as the entity that separates a group action from deformation. Paragios, Rousson and Ramesh [63] proposed a variational framework for global as well as local shape registration. Their optimization criterion accounts for global (rigid, affine) transformation and local pixelwise deformation. Similar idea was also used in Leventon et al. [78] in defining shape prior models in the level set representations. These ideas can be used to simultaneous approximation and registration of nonequivalent shapes, and tracking moving and deforming objects through time. However, the question of how to determine the local deformation has not been involved in these works. We will employ these ideas in our algorithm by viewing the evolution of active contour as a deformation of the interface, This deformation consists of a rigid transformation and a local deformation. We will use the average shape to determine the rigid transformation that better maps the interface to the prior shape, and use the image gradient and a few key points to determining the local deformation that provides more accurate segmentation. Based on these thought, the proposed variational framework includes a key points term that is able to incorporate a priori knowledge on expected shape and a few points that boundary should pass through into active contours. We modify the energy function of geodesic active contour so that it depends on the image gradient and prior shape, as well as a few prior points. We only need a few points since we have information on expected shape. The modified energy function gives a satisfactory segmentation in the PAGE 53 38 presence of relatively large shape distortion, and gaps that are a substantial fraction of the overall boundary. To combine both prior shape (a global constrain) and prior points (a local constrain) into a single variational framework, we use level set formulation. Suppose, are prior points (given by experts) the boundary passing through. Then the level set u, for mxx,...,1ix 0 .,...,1mi Now let xf be a smooth function defined on the image domain such that 1 0 x f 1 x f for .,...,1,mixxi The function can be obtained by convolution of and xf f where is a function taking value one on the points and zero elsewhere, and f mx x,...,1 is the standard mollifier. Now we can define key points term to be dxxuxfuE221 (219) This term forces the interface passing through the given points since minimizing it with sufficiently small results that umust be close to zero at those points. Note that is nonzero only on the xf neighborhood of the given points, so the Key Points term doesnt affect much the shape of the contour outside the neighborhood of the given points. 2.2.2.4 Intensity profile term Besides using shape prior, Leventon et al. [78]incorporated intensity and curvature priors to segmentation process by an approach similar to the one developed in Leventon et al. [60]. There are two reasons to have intensity profile terms. First, the model with intensity profile terms indeed performs segmentation and registration simultaneously. The registration in the model combines both a rigid transformation and a local PAGE 54 39 deformation. The rigid motion is determined explicitly by shape matching, while the local deformation is determined implicitly by the image gradient and prior intensity profile. Another reason to use prior intensity profiles to assist segmentation is that in some cases information of only the expected shape may not be sufficient to guide the active contour to the boundary of the object of interest. For instance, in some 2chamber cardiac ultrasound images, the image intensities of the myocardium are nonuniform, a significant portion of the border appears at low contrast, and their shapes are not similar to the prior. For this class of images a model incorporating only prior shape in the active contours may not be able to give an accurate segmentation. In this section, we present an alternative method for generating average intensity profiles from a set of training images. To generate an intensity model in a Gaussian model was used to compute the joint distribution of the intensity values and signed distances to the boundary from a set of segmented training images. In ultrasound images signal is partially in the form of speckle. Since the statistics of speckle are nonGaussian modeling the randomness of ultrasound images using a Gaussian model is not appropriate. Our method of generating an intensity model is model free and based on maximizing mutual information of image geometry between the intensity model and aligned training images in section (2.2.1.2). By the reasons discussed in the previous section (2.2.1.2), we propose to use the MIIG of the prior image I and aligned novel image I in a neighborhood of Cas the measurement of the disparity in intensity profiles across the interface and prior shape. Hence we define the intensity profile term to be TxRIxIMIIGEVTR11*,,,0 (220) PAGE 55 40 0V is defined in section (2.4). Notice that MIIG is always positive; we use negative sign to minimize the energy. 2.2.2.5 Proposed model and EulerLagrange equations Consider gradient (Equation 216), prior shape (Equation 218), key points (Equation 219) and intensity profile (Equation 220) terms together, we have the following model for segmentation with prior information in level set form. TxRIxIMIIGdxxuxfuTRxCduuIguVuTR11*32221,,,,,,,2*,2min0321 (221) )3,2,1(0ii are parameters balancing the influences from the image gradient, prior shape, key points and intensity profile. The active contour governed by Equation 221 is forced to arrive at high gradient, form a shape similar to the prior, across the key points and capture the prior intensity profile near the feature. Let xIGxf* TxRIGxg11 dxigifGiigfp2121,1,, where 222221,yxeyx G(According to [69]), then 2211,,diiigfpifpR 1212,,diiigfpifpR dxTxRxgigifGuiigfp2121221.,1,, (2.22) dxTxRxgigifGiigfp1121221.,1,, (2.23) PAGE 56 41 dxRxgigifGTiigfp1121221.,1,, (2.24) Let and notice 0V 221212121,,log,,,RdidiigpifpiigfpiigfpgfMI then 22212221212121,,,,,,log,,,RRdidigpfpgpfpgfpgpfpgfpgfpgpfpgfpdidiigpifpiigfpiigfpgfMI Notice we have 1,212didigfpR 0,,,,,,22212122RRdidigpgpgfpgfpdidigpfpgpfpgfpgpfpgfpgfpgpfpgfp then by Equation 2.22 dxdidiigpifpiigfpTxRgigifGdidiigpifpiigfpiigfpgfMIRR222121212121221212121,,log,1,,log,,, Similarly dxdidiigpifpiigfpTxddRgigifGgfMIR221212111212,,log,1, and 221212111212,,log,1,RdxdidiigpifpiigfpRgigifGTgfMI With the definitions and deduction above, the evolution equations associated with the EulerLagrange equation of Equation 221 are PAGE 57 42 uuudgdivutu2212 (2.25) 0,,0txnu ; xxux,0,0 u (2.26) 0310,0,, tgfMIdxuRxddut (2.27) 0310,0,,tgfMIdxuxddRddut (2.28) 0310,0,,TTtTgfMIdxuddutT (2.29) where R is the rotation matrix in terms of the angle is evaluated at d TRx 2.2.3 Parameter Choosing Generally speaking, the problem of determining the parameters 3,2,10 ii is not trivial. It depends on the specific problem. For a new type of image, it is a time consuming work to find the appropriate parameters. In this section, we present an alternative approach, that is not only able to incorporate both shape and segmentation, but also be able to provide the better estimate for the parameter used in the model. Our model is a coupled optimization problem, which consists of a minimization and a maximization problems. The minimization problem is TRuETRu,,,min,,, with dxxuxfuTRxCduuIguTRuEuTR2221,,,,,,2*,2min,,,321 (230) The maximization problem is PAGE 58 43 TRFTR,,min,, with TxRIxIMIIGTRFV11*,,,0 (231) where MIIG is the mutual information of image geometry defined in section (2.2.1.2), and ,, and Ttogether with uare solutions of Equation 230 corresponding to a fixed R The energy functional in Equation 230 is increasing in Without the joint problem (Equation 231), in Equation 230 takes the smallest value when E 0,0 that reduces to the energy functional for geometric active contour. It will leak through the boundaries with weak gradients. By maximizing the energy functional in Equation 231 over all the possible solutions ,, and Tof Equation 230 corresponding to R we can get an optimal estimate for and hence a better segmentation corresponding to this 2.3 Experiments and Results In this section, the proposed methods are applied on different applications. In real applications, it is possible that there are no all 3 kinds of prior information (shape, key points and intensity profile), any available prior of those 3 can be used for model (221). In the following experiments, we will show the results of shape analysis, generation of intensity profile, segmentation with gradient and prior shape information, segmentation with gradient, prior shape and key points information, segmentation with gradient prior shape and intensity profile information. Most of the experiments are for ultrasound cardiac images. Because the task of determining the epicardial and endocardial borders of the left ventricular (LV) on PAGE 59 44 echocardiography images is essential for quantification of cardifunction. It is very difficult due to poor contrast and information dropout at the borders, the inhomogeneity of intensities inside and outside the regions between epicardial and endocardial borders, and the low signal to noise ratio. 2.3.1 Shape Analysis In this experiment, the shape analysis method introduced in section 2.2.1.1 is applied. Epicardial and endocardial borders segmentations of 167 patients on 4 chamber ultrasound cardiac images, which are given by experts, are used as training set in this experiment. The task of this experiment is to cluster those shapes and find the average shape for each cluster. Among those patients, there are 112 normal and 55 abnormal. For each patient, there are 4 shapes given by expert, namely: ED (end diastole) Endo (endocardial), ED Epi (epicardial), ES (end systole) Endo, ES Epi. Applying Equation 21, 23 and the near basis cluster method. 167 sets of shapes were divided into 17 clusters with the cluster coefficient 0.1. 112 normal sets were in 3 clusters, 55 abnormal sets were in 14 clusters. Figure 22 shows the clustering results and the average shapes of normal patients. Figure 22A is the biggest cluster of normal patients, there are 79 patients in this cluster. The yellow curves are all 79 ED Epi and 79 ES Epi. The red curves are all 79 ED Endo and 79 ES Endo. Those *s show the average contour of ED Epi and ES Epi. Those +s show the average contour of ED Endo and ES Endo. From A, we can see that those contours in the same cluster are very tight, which means that shapes in the same cluster are similar. Figure 22B to C show the average contours of the 3 normal clusters without showing all contours. The green * with yellow curve is ES Epi. The green * with red curve is ED Epi. The blue + with yellow curve is ES Endo. blue + with red curve is ED Endo. B) is the average of 79 PAGE 60 45 patients, C) is the average of 24 patients, D) is the average of 1 patient. We can see the shape in D) is very different than B) and C). Figure 22. Cluster and average result. A) The biggest cluster of normal patients with all 79 sets of curves, B)The average shapes of the first normal cluster, C) The average shapes of the second normal cluster, D) The average shapes of the third normal cluster Figure 23 shows some examples of the average shapes of abnormal clusters. For 55 abnormal patients, 14 clusters were generated with the cluster coefficient 0.1. Again, The green * with red curve is ED Epi. The blue + with yellow curve is ES Endo. blue + with red curve is ED Endo. Figure 23. Example of average shapes of abnormal clusters. A) the average shape of the biggest cluster with 10 patients, B) the average shape of the biggest cluster with 9 patients, C) the average shape of the biggest cluster with 8 patients 2.3.2 Intensity Profile We applied the algorithm in section 2.2.1.2 on 2chamber cardiac ultrasound images. The method was tested against one cluster of a database of 85 apical 2 chamber PAGE 61 46 end diastolic (ED) echo cardiographs. Images acquired retrospectively from 61 normal patients. The images were grouped into three clusters. The 20 images from one of the clusters with expert traced epicardial borders superimposed are displayed in Figure 24. Figure 24. Twenty training images with the segmented epicardiums used as training images and shapes Using Equation 212 described in section 2.2.1.2, the intensity profile for each average shape was computed. The average shape and the associated average intensity *C PAGE 62 47 profile I for this cluster with 20 images is displayed in Figure 25. Figure 25A shows the direct average of those images, i.e. simply set in Equation 213 be 1. Figure 25B shows the average of those images with determined by using Equation 212. From Figure 25, we can see that A is more blurred than B. Moreover, we computed the sum in Equation 212 for the determined by using Equation 212 and The results were 196.73 and 184.08 respectively, confirming our suspicion that the weighted average is better. ja 20//1ja ja ja 20 I C*, d2 ,1 g t x, x0 A B Figure 25. Average images of the above 20 images with the average contour superimposed. A) by equally weighted these aligned images, B) by the proposed method 2.3.3 Segmentation with Prior Shape Information In this section, we assume that only the shape prior information is available. Hence Equation 221 changes to be uTRxuuguuTR2min1,,,,,32 (232) And according to Equation 225 to 229 the evolution equations are uuddivutu212 (233) 0,,0xnu ; uxu0, (234) PAGE 63 48 010,0, tdxuRxddut (235) 010,0,tdxuxddRddut (236) 010,0,TTtdxuddutT (237) We have tested our model on synthetic images and ultrasound images. 2.3.3.1 Synthetic image The aim of the first experiment was to verify that the active contour with the prior shape could fill in the gaps in a boundary. We used the first two terms in Equation 221 for this experiment. Figure 26A shows a curve that we used as the prior shape for this experiment. Figure 26B shows an image of a white semicircular disc on a black background. The boundary of the semicircular disc should be thought of as the boundary of the prior shape after some partial occlusion. The task is to see whether the active contour with prior shape information can utilize this partial boundary while filling in the rest. The active contour was initialized as the ellipse as shown in Figure 26B. Evolving the active contour according to the level set formulation of Equations 234 to 237 with the parameters 1 5.0 (in xg ), 11 02.0 dt 4 .10 0 0 ,T we get the stationary contour C in Figure 26C. The transformation parameters are 0,00 25.1 50 .00 and 13.2 14.2 0 T(pixels). In Figure 26D, the prior shape C, and the contour TRC are shown as the white and the dotted curves. From Figures C and D we can see that the contour Ccaptures the high gradient in the image I and the prior shape even though complete gradient information is not available. *C PAGE 64 49 Figure 26. Segmentation with prior shape for synthetic image. A) The prior shape B) The image with the initialized active contour (the ellipse), C) The final stationary contour, D) The contour (solid curve), and the contour (dotted curve) *C *C TRC* 2.3.3.2 Cardiac ultrasound images The aim of the second experiment was to segment the epicardium (the outer boundary of the myocardium surrounding the left ventricle) in a 4chamber image of the heart (Figure 27B). The epicardium is not completely imaged in the image, and our task is to find and complete the boundary using a prior shape. The prior shape (Figure 27A) is the average epicardial boundary generated in section 2.3.1 of the biggest normal group with 79 patients, The contour shown in Figure 27B was used as an initial contour. This contour was evolved according to the Equations 234 to 237, and it finally stopped at the location of the solid contour in Figure 27C. To validate this result we had the expert manually segment the epicardium. The dotted contour in Figure 26C is the experts epicardium. Our segmentation is close to the experts. The parameters used in this experiment are 1 5.0 (in ), xg 21 , 05.0dt 10 0 0 0,00 T The solutions we obtained are the contour Cin Figure 26C, and the transformation parameters 66.0 10.00 and (pixels). In Figure 26D we present the contour 10.0,88 .00T T RC for inspection. It is similar to the average shape in Figure 27A. *C PAGE 65 50 Figure 27. Segmentation with prior shape for cardiac image. A) A cluster of 79 epicardiums and the average shape B) The initial contour in the ultrasound image, C) The final stationary contour (solid curve) and the experts epicardium (dotted curve), D) The contour *C TRC* 2.3.4 Segmentation with Prior Shape Information and Key Points In this section, we assume that both prior shape and key points are available. Hence the energy functional changes to be dxxuxfuTRxCduuIguuTR2221,,,,,,2*,2min321 (238) And according to Equation 2.25 to 2.29 the evolution equations are uuudgdivutu2212 (239) 0,,0txnu ; xxuxu,0,0 (240) 010,0, tdxuRxddut (241) 010,0,tdxuxddRddut (242) 010,0,TTtdxuddutT (243) 2.3.4.1 Synthetic image The aim of the first experiment is to verify that the active contour with the prior shape and points can fill in the gaps in a boundary in a meaning for way in the sense of passing through the prior points. PAGE 66 51 Figure 28A shows a typical binary image with three points and an ellipse superimposed on it. The ellipse and points are used as the prior shape and points in this experiment, respectively. The object to be segmented partially occluded, and the shape of its boundary is nonequivalent to the prior shape. We want to see whether the active contour with the prior shape and points can utilize the partial boundary while filling in the rest. The active contour was initialized as a solid contour shown in Figure 28C. Evolving the active contour according to Equation 239 to 243 with the parameters 12.01 08.02 5.0 (in xg ), 05.0 dt 10 0 0 ,T, we get the stationary contour (the dotted one) in Figure 28C, and the transformation parameters 0,00 C 91.0 1.00 4 and 3 0,5.0 0 T(pixels). We can see that even though complete gradient information is not available, the contour captures the high gradient in the image C I passes through three prior points, and forms a shape similar but not the same as the prior C. To show the advantage of using prior points we compare the segmentation results obtained by using model (238) with (232). Figure 28B shows the segmentation result by using model (232). In Figure 28B the solid contour represents the shape prior and initial contour, while the dotted contour is the segmentation result. Since the prior points is not incorporated in the model (232) the segmented contour only captures the prior shape and high gradient. It cant accurately capture local shape variation. PAGE 67 52 Figure 28. Segmenation with prior shape and key points for synthetic image. A) An image with the prior shape and three points superimposed, B) The segmentation result using Equation 232 (dotted) and the initial contour (solid), C) The segmentation result using Equation 238 (dotted) and the initial contour (solid). 2.3.4.2 Cardiac image The aim of the second experiment is to segment the endocardium (the inner boundary of the myocardium surrounding the left ventricle) in a twochamber image of the heart (Figure 29A for a typical image). The endocardium is not completely imaged in the image, and its shape is not the same as the average shape (the shape prior). Our task is to determine the endocardium using average shape and five points given by an expert. The prior shape is created by the same way as that in section 2.2.1.1. For this particular problem to create the prior shape, endocardial boundaries were outlined in __ patients by an expert echocardiographer. The boundaries were grouped as explained in section 2.2.1.1. Figure 29B shows the average contour for one of the clusters (the dotted contour), the endocardium outlined by an expert (the solid contour), and five pints on the experts contour, That are in the location where the image gradients are low, and the local shape distortions are larger. Figure 29C presents the image IG* From this PAGE 68 53 image we can see the dropout of image information at several parts of the endocardium.are 1 5.0 (in xg ), 21 05.0 dt 10 0 0 0,00 T. 0024 .10 ,6.24 0 TR, To segment the endocardium in the image shown in Figure 29A, the active contour was initialized as the contour shown in Figure 29A. This contour was evolved according to the Equations 239 to 243 and it finally stopped at the location of the dotted contour in Figure 29D. We also obtained the transformation parameters 1710.00 and T. The solid contour in Figure 29D is the experts endocardium. 3.32 We can see that our segmentation is close to the experts. To see the shape variation between the solution of Equation 238 and average shape we aligned the solution of Equation 238 onto the average shape using the solutions of Equation 238. Figure 29E shows the disparity in shape between these two contours. The dotted contour is the transformed solution of Equation 238, and the solid one is the average shape. We can see that our active contour formed a shape different from the prior one in order to capture the high gradients and given points. Figure 29F provides the segmentation result obtained by using model (232) against the experts endocardium. In this figure the dotted contour is the solution of Equation 232, and the solid contour is the experts endocardium. Comparing Figure 29D with Figure 29F, we can see that the solution of Equation 238 is closer to the experts contour than that of Equation 232. Figure 29G presents the shape comparison between the solution of Equation 232 and the prior. We can see that they are very similar. Since Equation 232 is not able to incorporate prior points, its solution can only capture the high image gradients and the average shape, but it cant provide a desirable segmentation as the experts endocardium. PAGE 69 54 Figure 29. Segmentation with prior shape and key points for cardiac image 1. A) a typical 2chamber ultrasound image with an initial contour; B) Experts endocardium (solid contour), average shape (dotted contour), and five points on the experts contour, C) The image IG* ; (D). The endocardium segmented by using Model 238 (dotted) and the experts contour (solid), (E) The transformed solution of Equation 238 (dotted), and the average shape, F). The endocardium segmented by using Model 232 (dotted) and the experts contour (solid), G) The transformed solution of Equation 232 (dotted), and the average shape. The last experiment is repeating the second experiment in another 2chamber ultrasound image. The aim and the procedure of this experiment are the same as the second one. We list the figures below for the results of this experiment in the same order PAGE 70 55 as above. The segmentation is given in Figure 210D represented by the dotted contour, the transformation parameters are 9883.00 1981.00 and T. From these results and figures we can have the same conclusions as that found in the second experiment. 1.49,7.280 PAGE 71 56 Figure 210. Segmentation with prior shape and key points for cardiac image 2. A) a typical 2chamber ultrasound image with an initial contour, B) Experts endocardium (solid contour), average shape (dotted contour), and five points on the experts contour, C). The image IG* (D). The endocardium segmented by using Model 238 (dotted) and the experts contour (solid), E) The transformed solution of Equation 238 (dotted), and the average shape, F). The endocardium segmented by using Model 232 (dotted) and the experts contour (solid), G) The transformed solution of Equation 232 (dotted), and the average shape. PAGE 72 57 2.3.5 Segmentation with Prior Shape and Intensity Profiles In this algorithm, we assume prior shape and intensity profile is available. Hence the energy functional (224) changes to be TxRIxIMIIGuTRxCduuIguVuTR11*321,,,,,,,*,2min0321 Because the heavy computation duty of MIIG in the evolution and the parameter choice complexity, we used the method introduced in 2.2.3. Then according to Equations 230 and 231 the coupled energy functional pair is uTRxCduuIguTRuEuTR *,2min,,,2,,,, (2.44) TxRIxIMIIGTRFV11*,,,0 (2.45) Notice Equation 244 is same as Equation 232. Hence the associated EulerLagrange equations of Equation 244 are Equations 233 to 237. However, how to solve Equation 245 efficiently still remains an open question. In our experiments we first solved Equation 244 to get a sequence of solutions iiiiTRC ,,, corresponding to a sequence of kii,...,1, Then we computed iiiTRF ,, in Equation 245 for each k i 1. In this computation we partitioned the intensities into 16 bins, and used the discrete form of Shannon entropy to calculate the MIIG. If jjT jRF , is the smallest one among them, we took jjjjTRC ,,, as our model solutions. In this way we could get a better estimate for and hence, a better segmentation, that captures high gradients, shape prior, and intensity profile. However, it may not be the best, since the comparison of the energy values in Equation 245 was only for finitely many s. PAGE 73 58 We applied this algorithm on 2chamber cardiac ultrasound images. The epicardiums and the endocardiums in these images were not completely imaged, and our task was to find and complete the epicardial and endocardial boundaries using prior shape and intensity profile. To create the prior shape, epicardial boundaries were outlined by an expert echocardiographer on 85 images acquired at ED from 62 patients. Using the method described in section 2.2.1.1, the boundaries were grouped into three clusters and the average shape of each cluster was computed. Using the method described in section 2.2.1.2, the intensity profile for each average shape was computed. The average shape and the associated average intensity profile *C I for one cluster is displayed in Figure 25B. To segment the epicardial border in a novel image displayed in Figure 211B, we used the average contour and intensity profile near the contour in Figure 25B as the priors. The active contour was initialized with the ellipse displayed in Figure 211A. Evolving the active contour using Equations 233 to 237 with a fixed we obtained segmentation together with a similarity transformation C TR,, that are the solutions of Equation 244. By varying we generated a sequence of solutions of Equation 244. We chose the optimal value of to be the one maximizing Equation 245. Finally, the solutions for the Model 244 coupled with Equation 245 were chosen as the solutions of Equation 244 corresponding to this optimal The first column of the Table 21 displays 8 different values of By the procedure described above we obtained 8,...,1,,, iTRCjjjj The third and the fourth columns present the TxRI11 xIMIV*,0 and PAGE 74 59 TxRIxIMIIGV11*,004.0 respectively. Since the 4 th column of this table is largest when we selected the solutions of Equation 213 to correspond to this choice of 04.0 The segmentation (solid) corresponding to this is shown in Figure 211B together with the expert traced border (dotted). The distance between the expert and algorithm generated borders are tabulated in column 2 of the table. It is defined as where is the distance function of C, and is our segmentation. The units of the distance are the numbers of the pixels, and the pixel size is 0.62mm0.62mm. From this table we see that the segmentation corresponding to NiiCNpCd1/**Cd1,0,01NppNppC 04.0 is the one having smallest distance from experts contour and largest value of MIIG. This statement is not true for MI. 20 04.0 Table 21. Experiment results with different for Epi dist MI MIIG 0.04 3.1142 1.8894 7.9971 0.20 3.2903 1.9192 7.9794 0.40 3.3552 1.9207 7.9698 2 3.3962 1.9579 7.9623 20 4.7160 1.9207 7.8914 40 6.6697 2.0133 7.7524 0.02 25.4498 0 0 0.004 26.2190 0 0 The second row in Figure 211 displays this experimental result. The segmentations (solid) in Figure 211B and C are the solutions of Equation 244 with 04.0 and respectively. Comparing the segmentation results with the experts borders (dotted) this figure provides visual confirmation of the result presented in the table. To further test the method, we also used the same initial contour and the optimal value determined in the previous test in two additional images. The PAGE 75 60 segmentation results (solid) together with experts borders (dotted) are presented in Figures 211D to G. The segmentations in the left and right columns are the solutions of Equation 244 corresponding to 04.0 and 20 respectively. Comparing the results in the left column of Figure 211 with those in the right column, we observe that 04.0 also provides good segmentations in these two new images indicating that the optimal estimate of from one image can possibly be used for other members of the cluster. Of course, the shapes of the object boundaries and their intensity profiles must be similar. PAGE 76 61 Figure 211. Segmentation with prior shape and intensity profile for epicardiums. A) An ellipse used as the initial contour in our expenments for three images in the following three rows, (B)(G). Each row presents the segmentations (solid, red), experts borders (dotted, black), and average shape (green) in an image. The segmentations in the left column and right column are the solutions of Equation 244 corresponding to 04.0 and 20 respectively. PAGE 77 62 Moreover, the proposed model is less sensitive to the initial step than the edgebased active contours, since the propagation of the contour is influenced by intensity profiles in addition to image gradients. Figure 212. Intensity profile for endocadiums. A) 12 training images with the segmented endocadiums used as training images and shapes, B) Average endocardums with associated intensity profiles generated from the training shapes and images in A) by the proposed algorithms. We also applied our algorithm to segment the endocardiums in ultrasound images. The shape model of endocardium and the prior intensity profiles across the shape model shown in Figure 212B were generated from a set of 12 training shapes and their associated images shown in Figure 212A by using our proposed method explained above. PAGE 78 63 To segment the endocardial border in a novel image displayed in Figure 213B, we used the average contour and intensity profile near the contour in Figure 212B as the priors. Table 22. Experiment results with different for Endo dist MI MIIG 0.1 1.97 1.3711 6.8060 0.2 2.14 1.4390 6.7197 0.01 2.36 1.1653 6.6282 0.001 2.78 1.0783 6.4244 1.0 8.18 1.1719 6.4050 Table 22 shows the choice of .We observed the same phenomina as that in the previous table that the larger the MIIG is, the smaller the distance of the segmentation from the experts contour is. This statement is again not true for MI. From this table we see that the segmentation corresponding to 1.0 0 is the optimal value of and the solution corresponding to this is the model solution of Equation 244. Figure 213 has the same organization as Figure 211. Again, by comparing the left column with the right column, it is evident that the model solutions are closer to the experts contour than the solution of Equation 244 with other value of PAGE 79 64 Figure 213. Segmentation results for endocardiums. A) An ellipse used as the initial contour in our experiments for three images in the following three rows, (B)(G). Each row presents the segmentations (solid, red), experts borders (dotted, black), and average shape (green) in an image. The segmentations in the left column and right column are the solutions of Equation 244 corresponding to 1.0 and 1 respectively. PAGE 80 65 2.4 Conclusion We proposed a framework for segmentation that incorporates different prior information. The prior information could be shape, key points or intensity profile. We also proposed a level set formulation and a numerical algorithm for the active contour. In applications, the framework can be modified to be different model for the priors available. Several examples for different prior information were given in the experiment section. According to experiments for segmentation with prior shape information, the active contour could reliably segment images in which the complete boundary was either missing (the ultrasound images) or was low resolution and low contrast (fMRI). Besides the segmentation, the algorithm also provides estimates of translation, rotation, scale that map the active contour to the prior shape. These estimates are useful in aligning images. In the experiments for segmentation prior shape and key points, we add the prior information on a few points to a variational framework in level set formulation to make up the drawback of those models that use only the gradient and shape information. These prior points are got from expert and are known to be passed by the contour we are looking for. Since the distance of a point on a curve from that curve is zero, by using level set formulation we can incorporate the prior points into active contours with shape by creating an energy term. This energy is designed such that when it is minimized, the surface, where the interface is embedded, takes value close to zero at the prior points. The experiment results showed that the proposed model using both prior shape and points greatly improves the segmentation result comparing with the active contours with shape but without prior points. PAGE 81 66 In the experiment for segmentation that incorporates the shape and intensity priors in active contours, we solve a coupled optimization problem. The energy function in one part of the model is a weighted sum of the energy of the geometric active contour and a shape related energy. The minimizer of this energy function presents the segmentation and the transformation that aligns the segmentation to the prior shape. The second part of the model provides an optimal estimate of the weight used in the first one by maximizing the MIIG of the intensity prior and the aligned novel image near the feature over all the alignments that are the solutions of the first part corresponding to different weights. The improvements of this model over existing active contour algorithms are in two aspects. First, we used maximizing MIIG rather than MI to match intensity profiles of two images. The reason for doing this is that the MIIG takes neighborhood intensity distribution into account, and hence, gives better description of intensity profile than MI. Secondly, the parameter in the model of active contour with shape can be selected by model itself, and the segmentation with the optimized parameter arrives at higher image gradients, forms a shape similar to the prior, and captures the prior intensity profile. Moreover, the coupling idea can be used for any models where the solution is impacted by three forces. We applied our model to the problem of cardiac boundary determination in ultrasound images, for which the methods using edge or region information only can not give a good result. Even the active contours with shape prior struggle with such data. The proposed model was tested against a database of epicardial borders traced by an expert on echocardiographic images acquired from the apical 2chamber view. The preliminary results were encouraging. PAGE 82 67 There are still some drawbacks of this framework. First, because of the shape term this model does not allow big shape difference. However, a significant local shape difference may occur in application. Therefore, we are now working on more nonrigid model. The other drawback is that the initial curve cannot be far from the desired segmentation. The reason is that there may be some regions have high gradient between initial curve and the desired segmentation. In that case, initial curve will be stopped by that high gradient. At last, the computation for MIIG is expensive. Hence, an alternative way for computation of MIIG is necessary. PAGE 83 CHAPTER 3 APPLICATION OF PDE BASED INPAINTING ON MR PARALLEL IMAGING Inpainting is an image interpolation method. PDE based digital inpainting techniques are finding broad applications. This chapter shows the first work in which PDE based inpainting techniques are applied into the field of MR parallel imaging. First ,A novel model and its corresponding numerical method are initially introduced. Then two applications of this model on MR Parallel Imaging are introduced. The first application of this model is for sensitivity maps. Coil sensitivity maps are important in parallel imaging and they often need extrapolation and hole filling (holes being dark regions of low signal in MR images). All of these problems can be solved simultaneously by applying inpainting techniques. Experiments for determining coil sensitivity maps for phantoms, Neurovascular and cardiac MR images demonstrate the speed and accuracy of the proposed model. Images generated using SENSE utilizing inpainted sensitivity maps, raw sensitivity maps and oversmoothed sensitivity maps are compared. The second application is for intensity correction map. The goal of this study is to find a technique which can generate uniform image with optimized Signal to Noise Ration (SNR) for Multichannel RF array coils. Experiments for images from different popular systems show that the inpainted intensity correction map can generate uniform image without losing SNR. 68 PAGE 84 69 3.1 Introduction 3.1.1 Inpainting Image inpainting was originally an artistic phrase referring to an artists restoration of a picture's missing pieces. Computer techniques could significantly reduce the time and effort required for fixing digital images, not only to fill in blank regions but also to correct for noise. Digital inpainting techniques [36, 7982] are finding broad applications such as, image restoration, disocclusion, perceptual image coding, zooming and image superresolution, error concealment in wireless image transmission, etc. Due to its broad range of applications, various methods for inpainting have been developed, ranging from nonlinear filtering methods, wavelets and spectral methods, statistical methods (especially for textures), etc. The most recent approach to nontexture inpainting is based on the PDE method and the Calculus of Variations. According to Chan and Shen [80] PDE based TV models and the MumfordShah model work very well for inpainting problems with a more local nature such as hole filling (holes being regions of low signal in MR images). Hole filling is the most important issue in the correction of MR sensitivity maps (which are generally derived from MR images). 3.1.2 Sensitivity Maps Sensitivity maps are meant to contain the sensitivity information of RF probe coils in MRI systems. Highly accurate knowledge of the spatial receiver sensitivity is required by both SMASH (SiMultaneous Acquisition of Spatial Harmonics [83]) and SENSE (Sensitivity Encoding [16, 84]). Sensitivity maps are also crucial to correcting the inhomogeneities of MRI surface coils. It must be noted that accurate sensitively information can only be obtained where signal is present. Some data sets with large areas contributing little or no signal. Such dark regions (i.e. holes) are common, for example, in PAGE 85 70 pulmonary MRI using Fresh Blood Imaging (FBI). In this case sensitivity map interpolation techniques are required to fix the holes. To deal with slightly varying tissue configurations and motion, extrapolation over a limited range is also necessary. We therefore need a technique, which is capable of simultaneously interpolating and extrapolating. There are some existing techniques for this purpose, such as the polynomial fit procedure [16], thinplate splines [85], wavelets [86],and Gaussian kernel smoothing [87]. Those methods are based on the assumption that the sensitivity map is sufficiently smooth (i.e., containing no sharp variations). By our experiments, this assumption is not always true. If the sensitivity map is piecewise smooth, which is usually true for nonuniform loading, the existing methods are not sufficient. The inpainting model proposed here is a technique to handle interpolation and extrapolation for piecewise smooth maps simultaneously. 3.1.3 Intensity Correction Map In the areas of medical imaging, acquired images may be corrupted by poor signal to noise ratio (SNR) or nonuniformity in spatial intensity. Multichannel RF array coils produces a multitude of individual images. These images may be combined to produce a single composite image with high SNR. The standard technique for this process is the socalled sumofsquares (SOS) algorithm, however, consideration of the correlation of noise between channels also permits a somewhat better algorithm, commonly referred to as the optimal reconstruction [88]. This algorithm produces nearly the highest SNR attainable, at each pixel, given the individual channel images. The limitation of these algorithms is that for coil array elements with very high sensitivities, nonuniformity in the resultant image is likely to occur. To correct nonuniformities, a correction map is usually used to multiply the nonuniform image and generate a uniform one. The PAGE 86 71 correction map could be the inverse of the convolution of nonuniform image with Gaussian filter [87], or the inverse of the surface coil's sensitivity map [89, 90]. Our experiments show the correction map generated by with Gaussian filter could be over smoothed and substantial noise amplification can occur, which hinders the visualization of salient features. Accurate sensitivity map is an appropriate tool for intensity correction. If the sensitivity map is not accurate, then noise may be amplified. However, to generate accurate sensitivity map for intensity correction requires external reference and may take long time. On the other hand, those methods and some other methods provide a low SNR uniform image. Given a low SNR uniform reference and the high SNR sum of squares inpainting technique can be applied to generate a correction map, which can correct the nonuniformity and keep high SNR of sum of squares. 3.2 Methods In this section we first review the conventional Total Variational (TV) inpainting model. Then modified models are proposed for the applications in parallel imaging. The corresponding numerical method is given for implementation. Methods for correction of sensitivity map and intensity correction map are also introduced. It will be helpful to first explain the following terms. The process by which holes are fixed using inpainting is referred to as diffusion. There are two basic types of diffusion, isotropic and anisotropic. In isotropic diffusion, pixels in hole regions are fixed by assigning them a value equal to the isotropically weighted average value of their surrounding neighbors. Isotropically here meaning the weights will be independent of direction from the pixel, but decreasing with distance. Anisotropic diffusion involves assigning these pixels values obtained by extending the intensity contours from nonhole regions into hole regions along the intensity level set. The intensity level set being PAGE 87 72 simply the contour lines of equal intensity in an image, i.e., the set of equiintensity curves. Intensity here refers only to the magnitude of the complex pixel values. 3.2.1 Review of TV Inpainting Model If the inpainting problem is treated as an inputoutput system, then the input is the image u which needs inpainting along with the specification of the hole regions (dark regions), D, in the image. The output is a denoised image without any holes. 0 u The holefree part of the image, where inpainting is not needed, is not supposed to be changed much. Mathematically we wish to keep and uare as close as possible on the parts of the image not contained in D. For this purpose we define, u 0 : the image domain D\ : the area where the image information is contained, which is referred to as known locations. Let us define, DxDxxD\,,10)( (31) Then the desired similarity of and u at known locations can be achieved through minimizing the term, u 0 u0 dxuD2 x This term is called the fidelity term. To reconstruct the missing parts without losing local image texture (i.e., near the holes), it is necessary that the image model smoothly extends the image into the holes and preserves the edges of objects contained in the image. This can be achieved through minimizing [22, 36] dxu (32) PAGE 88 73 This term is used to smooth and extend the image along the intensity level set, hence this term is called the smoothing term. The conventional TV model combines these two terms and is given as, dxudxuuxDuuED20021],[ (33) where is a parameter used to balance the smoothing and fidelity terms. Given and D, the which minimizes the energy functional will be the inpainting result. 0u u 3.2.2 Modification of the Exponent of u in the TV Model The main feature of the conventional TV model is its ability to preserve edges given the fact that the anisotropic diffusion is always perpendicular to [91]. Minimizing u dxu makes sure that the smoothing is along the level set and therefore preserves and enhances edges. However, there are limitations that arise when reconstructing piecewise smooth images or those, in which the noise and the edges are extremely difficult to distinguish, which is often true for sensitivity maps. In smooth regions, isotropic diffusion is appropriate to remove noises. If traditional TV model is applied on to a region where noise and edges are difficult to distinguish, then anisotropic diffusion is applied and those noises are kept. Because the exponent of u determines the type of diffusion (isotropic or anisotropic), this naturally leads to models in which the exponent of u is a function P(X) [27]. More specifically, one might want to consider uPP where, lim 2)(0 ypy and 1)(lim ypy Therefore we may wish to minimize dxuuxP),( (34) PAGE 89 74 where P can either be a constant, 1 P or 2 P or where P can be a function, 0,20,*111)(1),(20DDuGxguxP (35) where are parameters. is used to accentuate the edges between regions of different intensities, larger will better preserve edges. is for the Gaussian function, /4/exp22xG which is used to smooth the image by the standard method of weighted averages. At areas of fairly uniform intensity, u is small, causing P(x) to be large. Hence the diffusion is more isotropic. At edges between regions with very different intensities, u is large so that P(x) will be close to 1; anisotropic diffusion is then applied so that the edges will be preserved. For the edges between D and D\ D is larger than 0, hence P(x) will be set to 2 and image will be isotropically extended. Therefore the goal of choosing P(x) is to control the type of diffusion at specific locations in the image, so that isotropic and anisotropic diffusion are effectively combined (with little threshold sensitivity). Specifically, total variation diffusion is used to preserve edges and create varying combinations of isotropic and anisotropic diffusion in regions that may or may not contain significant features. 3.2.3 Modified TV Model for Sensitivity Maps In order to apply the model to sensitivity maps, the conventional model may be modified based on the particular application. There will be a sensitivity map for each coil j. Let be the image generated by coil j, and let ju jI I be the homogeneous image, then PAGE 90 75 0ju is defined as In order to avoid singularities and to consider every coil, should be replaced by IIj/dx2 uuxD0 xDuDjj0, )((Duxj jjjDdxIIuxj2 (36) If the square root of the sum of squares of the s is used as the homogeneous image I (i.e., sumofsquares image), then the sum of squares of every coils sensitivity map will be unity at each pixel. Therefore jI dxujj221 (37) is to be minimized, where uis the inpainted sensitivity map for the jth probe coil. Then the modified TV model for the sensitivity map is, j dxudxudxIIuuEjjjuxPjjjjjj22),(21221][ (38) where are used to balance the terms in energy functional. Increasing in Eq. 8 will result in a smoother image. Larger makes the square root of the sumofsquares of the sensitivity maps closer to unity. 3.2.4 Numerical Implementation Method The EulerLagrange equation for Equation 38 is 011).()22jjjjPjjjuuuPuIII (39) PAGE 91 76 with boundary condition 0/ nuj for the jth coil. where n is the outward unit normal to the image boundary, To derive the evolution equations corresponding to the EulerLagrange equation, the following difference scheme [92] can be applied in order to discretize Equation 39. Let ube an image, then define lklkklylkklklylklkklxlkklklxuuuuuuuuuuuu,1,1,,,1,1,, and then klyklxkluuu, where correspond to rows and columns of the image respectively. lk, Let be a vector field and let be the step space for this field (generally equal to unity), then according to [92] ),(21vvv h ,)(21hvvvdivyx Applying the difference scheme above, the iteration formula for Equation 39 is found to be, jjDjDjnlkjnlkjnlkjnlknjuCCCCIIIuCuCuCuCujj243212,1,4,1,3,,12,,11111 (310) where klklapC1 lklkapC,1,12 klklbpC3 1,1,4lklkbpC PAGE 92 77 2221,1,22klpnlknlknklxkluuua 222,1,122klpnlknlknklxkluuub with j corresponding to the jth coil, corresponding to rows and columns of the image respectively, and corresponding to the nth iteration. lk, n Pixels in I (the sumofsquares image) with intensity below some particular threshold will be treated as holes, i.e., belonging to D. This threshold can be automatically determined in several ways. One simple approach is finding the maximum pixel intensity in the entire image and setting the threshold to be 0.05 times this value [i.e., 0.05(maximum intensity)]. A good initial guess for ucan reduce the number of iterations involved in applying Equation 310. Raw sensitivity maps for each coil, defined as j DxguessDxxIxIxrawujj,,_ are used as the initial guess for u. In our experiments, the magnitude of above is the average intensity of the phase of is taken to equal the phase of Another issue associated with sensitivity maps is the fact that they contain complex data. Therefore we can either inpaint the vector map directly or inpaint the magnitude and phase part separately. Our experiments showed no explicit improvement through the use of inpainted phase maps. Hence to save inpainting time we only apply inpainting to the magnitude of the sensitivity map, and the phase of is used for j guessjI jI guess jI ju PAGE 93 78 3.2.5 Modified TV Model for Intensity Correction Map To protect the signal to noise ratio, the intensity correction map needs to be isotropically smooth. Hence the inpainting model for intensity correction map is dxudxuuxDuuED220021],[ (311) where is a parameter to balance smooth term and fidelity term. There are many numerical ways to minimize this energy functional. One general way is to solve its Euler Lagrange equation: with boundary condition 0)())((0uIuuxD 0nu By the same difference scheme, it follows the Iteration Scheme 41,1,,1,11,DDnlknlknlknlknlkuuuuu (312) 3.3 Experiments and Results In this section, the proposed inpainting techniques are applied to two applications. One is probe coil sensitivity map, the other one is intensity correction map. Notice that those two kinds of maps are all the division of two images. To avoid singularity, there are some region without values (holes) To fix those holes without losing texture is the main target. And inpainting is good tool to fix holes. Therefore, inpainting technique can be applied to any map, which is division of two images. In all of those experiments, Matlab codes are run on a COMPAQ PC with 1G Hz CPU and 1G RAM. 3.3.1 Application on Sensitivity Map The proposed inpainting model Equation 38 was applied as sensitivity map correction. The corrected sensitivity maps were then used by SENSE in the reconstruction of MR images. The accuracy and efficiency of conventional inpainting technique can be seen in [36, 80, 82]. To show the accuracy of the modified model for PAGE 94 79 sensitivity map, we first apply the proposed the method to a phantom without holes. Some holes were randomly chosen and inpainted. The results were then compared with the original sensitivity map that has no holes. Since the actual sensitivity map is not known in real applications, we were unable to determine directly if an inpainted sensitivity map was good or not. Therefore we applied the inpainted sensitivity maps to reconstruct an MR image and then compare the reconstructed image with the reference image (which is generated by using all of kspace as opposed to SENSE which only utilizes part of the kspace data). Let the phrase intensity difference refers to the difference in magnitudes between the reconstructed and referenceimages at each pixel. We define the ghost ratio as the magnitude of the intensity difference (at each pixel) summed over every pixel in the image divided by the sum of the absolute values of each pixel in the reference image. The ghost ratio is actually the relative error. Notice that the intensity of reconstructed image by SENSE will be lower than reference image because of the lower energy in Kspace. According to Parsevals equation, the intensity of all reconstructed image were calibrated by multiplying reduce factor when calculate the intensity difference in all experiments. An alternative way to fix holes is to apply Gaussian kernel smoothing (12). For the sake of comparison with our method (i.e., the modified TV model), sensitivity maps corrected by Gaussian kernel smoothing (GKS) and raw sensitivity maps are also applied to SENSE. In many real applications, a lowresolution image is used to produce a sensitivity map. According to [93], it is possible to extract the sensitivity map directly from a fully sampled central band of kspace. (The center of the band being the center of kspace, with PAGE 95 80 the band extending along the complete width in the frequency encode direction, but leaving out the top and bottom sections of the kspace.) Since the acquisition of the sensitivity map data and the data to be reconstructed occurs simultaneously, errors due to sensitivity miscalibration are eliminated. An alternative method for acquiring the sensitivity map for dynamic MR (such as a cardiac series) is to apply one highresolution sensitivity map (i.e., utilizing the entire kspace) for a complete time sequence. Both methods are tested in our experiments to show the flexibility of the inpainting technique. The choice of parameters may depend on the system and loading. According to experiments, the model is not very sensitive to the choice of parameters. In our experiments, we set 25.0,1 ,1 0 and 1 3.3.1.1 Phantom with no intrinsic holes The phantom was constructed from a 197mm (7.75) ID by a 180mm long acrylic tube that was filled with a solution of Cu2SO4 (2.0 grams/Liter) and NaCl (4.5 grams/Liter). Fig 1a is a photo of this phantom. Data was then collected by a 1.5 T SIEMENS system (FOV 300 mm, matrix 128128, TR 300 ms, TE 15 ms, flip angle 90, Slice thickness 5 mm, number of averages 1) with an 8channel tuned loop array of coils. The time required for the calculation of 8 sensitivity maps was 4.35 seconds. If we ignore the container there are no holes at all in the original sensitivity map. Figure 31B shows the magnitude of the original sensitivity map of coil 1. 5 holes were then randomly made in the original sensitivity map. Figure 31C shows those holes. Figure 31D shows the inpainted result for coil 1. The local patterns are perfectly preserved in the inpainted result. PAGE 96 81 A B C D Figure 31. Inpainting results for the artifical holes. The maps shown in BD are for coil 1. A) Photo of the phantom, B) magnitude of the reference sensitivity map, C) magnitude of the sensitivity map with 5 randomly chosen holes B) magnitude of the inpainted sensitivity map by the proposed method We define the relative error as the sum over the absolute value of the difference between the inpainted result, and the original sensitivity, for each pixel in the map over the sum of the magnitudes of the original sensitivity map. Table 31 shows the relative error of each channel. The average relative error of all 8 channels using the modified TV method is 0.71%. The average relative error using GKS is 4.45%. Table 31. The relative error of inpainted results coil Relative error by inpainting Relative error by GKS 1 0.49% 3.52% 2 0.40% 4.45% 3 1.84% 7.57% 4 0.70% 4.87% 5 0.46% 3.44% 6 0.29% 2.61% 7 1.08% 4.12% 8 0.44% 4.97% average 0.71% 4.45% 3.3.1.2 Coronal phantom Figure 32 and Table 2 show results for a Coronal Phantom collected by a 1.5 T GE system (FOV 480 mm, matrix 256256, TR 500 ms, TE 13.64 ms, flip angle 90, Slice thickness 3 mm, number of averages 1) with an 8channel Neurovascular Array coil. The time required for the calculation of 8 sensitivity maps was 16.11 seconds. PAGE 97 82 In this experiment, highresolution sensitivity maps are generated by using images reconstructed utilizing the full Kspace for each coil. Lowresolution sensitivity maps are generated by using images reconstructed with the 31 central rows of fully sampled data. Lowresolution sensitivity maps are more likely to be used in real applications. Figure 32A is the Reference Image, which is reconstructed utilizing the full Kspace. Figures 32B to 32D show the magnitude of the sensitivity map of coil 6. Figure 32B is the magnitude of the raw sensitivity map. Figure 32C shows the result of inpainting applied to 32B, while Figures 32D is the result of applying GKS to 32B. From Figures 32C, it can be seen that inpainted sensitivity maps have no holes and thus extend the map to the entire image. Furthermore, in Figures 32C, the local patterns are perfectly preserved. Figures 32E to H show the intensity difference between images reconstructed by SENSE with a Kspace reduction factor of 4 (i.e., using 25% of the full Kspace data) and the reference image, those images use the same gray scale map. Figure 32E shows this intensity difference when the reconstruction is done using highresolution inpainted sensitivity maps. Figure 32F shows this intensity difference when the reconstruction is done using highresolution raw sensitivity maps while 32G shows this intensity difference for GKS corrected sensitivity maps. From Figure 32E to G, it is seen that the result obtained by using inpainted sensitivity maps is much better than that obtained by using either the raw sensitivity maps or the GKS sensitivity maps. Figure 32H shows the intensity difference obtained using lowresolution inpainted sensitivity maps. Figures 32H also shows that even in this case result is still better than that obtained by the use of highresolution GKS sensitivity maps. PAGE 98 83 A B C D E F G H Figure 32. Inpainting results for coronal phantom. A) Reference Image, B) magnitude of the raw sensitivity, C) magnitude of the inpainted sensitivity map, D) magnitude of the sensitivity map by Gaussian kernel smooth, E) The intensity difference by using high resolution inpainted sensitivity map, F) The intensity difference by using raw high resolution inpainted sensitivity map, G) The intensity difference by using high resolution Gaussian Kernel smoothed sensitivity map, H) The intensity difference by using low resolution inpainted sensitivity map, E) to H) use the same gray scale map Table 32 gives the ghost ratios obtained through the use of difference sensitivity maps. The term, R_factor, refers to the reduction factor of Kspace. It can be seen that the results obtained by using the inpainted sensitivity maps have much lower ghost ratios for each reduction factor than the results obtained by using raw sensitivity maps. Notice that the results obtained from the highresolution GKS sensitivity maps are very similar to those obtained from the lowresolution GKS sensitivity maps. The reason for this is that the GKS correction involves an isotropically smooth blurring of the image, which is very similar to the isotropic blurring which results due to an image being low resolution. An image reconstructed through the use of a fully sampled central region of Kspace also results in an isotropically smoothed version of the reference image, since this is also a lowresolution image. (This is true since any point in Kspace is a superposition of signals from each point in the slice.) For a uniform phantom, these two maps(a low PAGE 99 84 resolution GKS map and a highresolution GKS map) may generate similar results. For the same reason that low resolution results in an isotropic smoothing, the images from inpainted lowresolution sensitivity maps are similar to those resulting from GKS sensitivity maps. Table 32. Ghost ratio of coronal phantom R_factor = 2 R_factor = 3 R_factor = 4 Inpainted high resolution 1.06% 3.12% 2.24% Raw high resolution 6.07% 7.41% 6.41% GKS high resolution 2.18% 4.48% 4.76% Inpainted low resolution 2.03% 4% 4.6% Raw low resolution 6.7% 7.96% 7.26% GKS low resolution 2.20% 4.48% 4.72% 3.3.1.3 Neurovascular images Both Figure 33 and Table 33 show the results for Neurovascular Images collected by a 1.5 T GE system (FOV 440 mm, matrix 256192, TR 1000 ms, TE 13.504 ms, flip angle 90, Slice thickness 2 mm, number of averages 2) through the use of a fast spin echo pulse sequence with an 8channel Neurovascular Array coil. The time required for the calculation of 8 sensitivity maps was 15.01 seconds The contents of Figure 33 and Table 33 are similar to that of Figure 32 and Table 32. The reducefactor for Figure 33 E to Figure 33 H is also 4. Table 33 again shows that the results obtained by using inpainted sensitivity maps have much lower ghost ratios for each reduction factor. For Figure 33H, the lowresolution sensitivity maps are generated by using images reconstructed with the 31 central rows of fully sampled kspace data. PAGE 100 85 A B C D E G F G H H Figure 33. Inpainting results for neurovascular images. A) Reference Image, B) magnitude of the raw sensitivity, C) magnitude of the inpainted sensitivity map, D) magnitude of the sensitivity map by Gaussian kernel smooth, E) The intensity difference by using high resolution inpainted sensitivity map, F) The intensity difference by using raw high resolution inpainted sensitivity map, G) The intensity difference by using high resolution Gaussian Kernel smoothed sensitivity map, H) The intensity difference by using low resolution inpainted sensitivity map E) to H) use the same gray scale map Table 33. Ghost ratio of neurovascular images R_factor = 2 R_factor = 3 R_factor = 4 Inpainted high resolution 4.21% 5.14% 6.82% Raw high resolution 13.75% 14.40% 15.62% GKS high resolution 5.91% 7.74% 11.28% Inpainted low resolution 6.61% 8.36% 10.57% Raw low resolution 15.07% 15.49% 17.10% GKS low resolution 8.81% 9.7% 12.76% 3.3.1.4 Cardiac images Table 4 and Table 5 show the results for Cardiac Images collected by a 1.5 T GE system (FOV 240 mm, matrix 224192, TR 4.530 ms, TE 1.704 ms, flip angle 45, Slice thickness 5 mm, number of averages 2) through Fast Imaging Employing SteadyState Acquisition (FIESTA) with a GE 4channel cardiac coil. Breathholds ranged from 10 20 seconds. There are 20 images per heartbeat. Highresolution sensitivity maps (for each coil) are obtained from the first image in this sequence through the use of the full K PAGE 101 86 space. SENSE is then applied to each image in the sequence including the first by using the same sensitivity maps. Table 34 shows the ghost ratios obtained by applying the high resolution sensitivity maps from the first image reconstruction to all 20 images in the same heartbeat. The R_factor of Kspace here is 2. Note the ratios for image 1 are much lower then the others since strictly speaking the sensitivity maps are only correct for this image due to the motion of the heart during the heartbeat. Table 5 shows the results of applying lowresolution sensitivity maps obtained from each image in the series. These are generated from the 51 central rows of fully sampled kspace data [93]. The reduction factor of Kspace is in this case is only 1.6232. The noninteger value is due to rows being skipped only outside the 51 row central band. Table 34 and 35 show that inpainted sensitivity maps are more accurate than both the raw sensitivity maps and GKS sensitivity maps for dynamic image. PAGE 102 87 Table 34. Ghost ratio of cardiac images by applying highresolution sensitivity maps of first image Image Inpainted raw GKS 1 2.19% 9.18% 2.98% 2 11.69% 16.50% 12.04% 3 15.32% 19.43% 16.42% 4 16.10% 20.22% 17.10% 5 15.92% 19.98% 16.68% 6 15.70% 19.53% 16.76% 7 15.59% 19.40% 16.58% 8 15.39% 19.36% 16.20% 9 15.38% 19.26% 16.35% 10 15.57% 19.43% 16.60% 11 15.53% 19.42% 16.28% 12 15.33% 19.23% 16.34% 13 15.23% 19.14% 16.14% 14 15.04% 18.98% 15.80% 15 15.07% 18.92% 16.06% 16 15.08% 18.95% 15.95% 17 15.19% 19.17% 15.98% 18 15.45% 19.46% 16.47% 19 14.93% 19.06% 15.65% 20 11.83% 16.48% 12.47% average 14.38% 18.55% 15.24% PAGE 103 88 Table 35. Ghost ratio of cardiac images by applying lowresolution sensitivity maps of each image Image Inpainted raw GKS 1 10.64% 17.85% 11.96% 2 11.29% 18.57% 12.25% 3 11.08% 18.34% 12.64% 4 11.16% 18.27% 12.58% 5 11.04% 17.99% 11.93% 6 10.62% 17.59% 12.00% 7 10.45% 17.23% 11.66% 8 10.56% 17.18% 11.31% 9 10.28% 17.00% 11.23% 10 10.42% 17.26% 11.62% 11 10.45% 17.25% 11.17% 12 10.37% 17.17% 11.53% 13 10.26% 17.09% 11.37% 14 10.26% 17.13% 10.97% 15 10.19% 17.26% 11.42% 16 10.11% 17.18% 11.00% 17 10.24% 17.34% 11.19% 18 10.63% 17.74% 12.06% 19 10.81% 17.85% 11.59% 20 10.69% 17.92% 11.67% average 10.58% 17.56% 11.66% 3.3.2 Application on Intensity Correction Map 3.3.2.1 Uniformity scheme Figure 34 shows the steps in the correction process described in HeelsBergen et al. [94]. First, the multichannel MRI raw data is duplicated to produce two data files. One data file is used to generate an image with optimized signal to noise ratio but suffering from the resultant nonuniformity. The other data file is used to generate a uniform image but with low SNR. Then each pixel of the uniform image is divided by the corresponding pixel from the high SNR image to generate a raw correction map. Inpainting is then applied to the raw correction map to generate a smooth noiseless and hole free correction map. Lastly, multiply the optimized SNR image with the correction map to generate the final image. Here the calibration of the correction map is done c PAGE 104 89 through inpainting which is capable of denoising, interpolation and extrapolation all at the same time. Figure 34. The flow chart of the intensity correction method. In the following experiments, the sumofsquares is used as the high SNR image with nonuniform intensity. The uniform references were generated by different methods to test the robustness of the proposed method (Equation 36). Dividing each pixel of the uniform image by corresponding pixel from the optimized SNR image can generate the raw correction map. To avoid exaggerated noise, structural and nonstructural regions may be separated based upon the gradient magnitude before the division. A gradient threshold value is set at a desired gradient magnitude level. Pixels having gradient PAGE 105 90 magnitudes at or above the threshold value are considered to meet a first criterion for defining structure in the image, while pixels having gradient magnitudes lower than the threshold value are initially considered nonstructure. The threshold can either be decided by gradient histogram or by operator intervention. Because the raw correction map will be smoothed and inpainted, the structural and nonstructural regions do not to be separated accurately. The result will not be significantly influenced if part of structural regions is wrongly clustered. Because of the separation, some holes could be generated. The proposed model is applied to inpaint the raw correction map. The final image is generated by multiplication of the high SNR image and the correction map. The correction ratio map balances intensity across the image to produce a more uniform result. And the experiments show multiplication with a smooth ratio does not appreciably decrease signal to noise, the result image has also high SNR. In all of those experiments, the whole intensity correction process, which includes every step described above, cost 0.5 to 0.7 seconds in together. 3.3.2.2 Phantom Figure 35 shows the results for an Axis Phantom collected by a 1.5 T GE system (FOV 480 mm, matrix 256256, TR 500 ms, TE 13.64 ms, flip angle 90, Slice thickness 3 mm ) with an 8channel Neurovascular Array coil.. In this experiment, matlab 6.5 was used running on a PC with a 1GHz CPU and with 1G RAM. The whole intensity correction process involving each step described above, cost 0.7 seconds in CPU time. Two sets of data were collected for the same slice to calculate the signal to noise ratio in accord with the National Electrical Manufacturers Association (NEMA) standard. PAGE 106 91 Figure 35A shows the magnitude of the uniform image and Figure 35B shows the corresponding signal to noise ratio. To generate optimized signal to noise ratio image, weighted sum of squared algorithm was used with consideration of the correlation of noise between channels. Figure 35C shows the sum of squared image and Figure 35D shows the corresponding signal to noise ratio. Figure 35E shows the raw correction map and Figure 35F shows the inpainted correction map. Linear interpolation was used to get the initial guess for holes in raw correction map. Figure 35G is the final result image and Figure 35H shows the corresponding SNR. The final result is clearly uniform. And by comparison of Figure 35D and Figure 35H, it can be seen that the signal to noise ratio of the final result has no big difference with the optimized SNR. Hence the uniformity is dramatically improved and the signal to noise ratio is well protected through the proposed technique. 3.3.2.3 Clinical images Figure 36 shows some intensity correction results for clinical Neurovascular Images. Those data collected by a 1.5 T GE system (FOV 440 mm, matrix 256192, TR 1000 ms, TE 13.504 ms, flip angle 90, Slice thickness 2 mm, number of averages 2) through the use of a fast spin echo pulse sequence with an 8channel Neurovascular Array coil. Figure 36A, D and G are the sumof squares that are not uniform but with high SNR. Figure 36B, E and H are the uniform reference images that are uniform but with low SNR. Figure 36C, F and I are the corrected images by the proposed method. The SNRs of the uniform reference at head, neck, chest are 25.8,23.5 and 13.5. The SNRs of the corrected image at head, neck, chest are 44.2, 39.6 and 14.5. PAGE 107 92 Figure 35. Intensity correction for phantom. A) the magnitude of the uniform reference image, B) the corresponding SNR of the reference image, C) the sum of squared image, D) the corresponding SNR of the sumofsquares, E) the raw correction map, F) the inpainted correction map, G) the final intensity corrected image, H) SNR of the intensity corrected image PAGE 108 93 Figure 36. Intensity correction results for clinical Neurovascular images. A,D,G) sum of squares, B,E,H) uniform references, D,F,I) intensity corrected images by the proposed method 3.4 Conclusion A novel inpainting model was proposed and was applied to sensitivity maps. From Figs 1 to 3, it can be concluded that all of the major issues involved in the correction of sensitivity maps are effectively and simultaneously dealt with, denoising is accomplished, holes have been fixed and the original map has been smoothly extended to the entire image. From Tables 1 to 5 it is further concluded that inpainted sensitivity PAGE 109 94 maps produce better results than both raw sensitivity maps and GKS sensitivity maps. The ghost ratios of the results obtained by using inpainted sensitivity maps are only 23% to 61% of the ghost ratios obtained by the use of raw sensitivity maps, and 48% to 94% of the ghost ratios obtained from using GKS sensitivity maps. The authors of [16] suggested that the sumofsquares is applicable only if the object phase is sufficiently smooth so as to ensure that the corrected maps will be largely unaltered. With the application of the proposed model, this restriction becomes unnecessary since isotropic and anisotropic diffusion are effectively combined. The inpainting technique is fast. The average process time for a 256 256 image is less than 2 seconds. The proposed model functioned well for each type of image tested, phantom, neurovascular, and cardiac, and for both the highresolution and lowresolution methods of generating sensitivity maps. The proposed method does not need extra data for uniformity as some other methods do. By the experiments for intensity correction maps, it can be concluded that, with the use of inpainted intensity correction maps, this method is capable of correcting nonuniformity while perfectly preserving SNRs. PAGE 110 CHAPTER 4 MODIFIED MUMFORDSHAH MODEL FOR MEDICAL TOMOGRAPHY In this chapter a novel objectbased variation framework for Medical Tomography problems with extremely noisy and limited data is proposed. MumfordShah Model by using Level Set, one of the most recent and very successful approaches in medical image processing, is introduced into the tomography field. In the method section, simultaneous tomographysegmentation Models for different cases are given. A tomography with guiding image model is also introduced. The application of a Fast Algorithm for Level Set Based Optimization onto the proposed models ends the method section. In the application section, a new noninvasive medical imaging techniqueNoise Tomography is presented. The application of the proposed method onto Noise Tomography is discussed. Experimental results and comparison with conventional methods are also reported in this section. A conclusion of the proposed framework ends this chapter. 4.1 Background of Tomography hT is study is for medical tomography problem with very limited and noisy data. The problem of reconstructing physically inaccessible objects from tomographic data arises in many medical imaging applications, such as Xray Computed Tomography (xray CT), Positron Emission Tomography (PET), Single Photon Emission Computed Tomography (SPECT), Electrical Impedance Tomography (EIT) and so on. The image reconstruction problem arising in those applications is commonly referred to as tomography, which is an inversion problem from a mathematical point of view. The 95 PAGE 111 96 inverse problems that arise in these applications are typically illposed and require solving large dimensional estimation problem. 4.1.1 Existing Methods To solve an illposed problem, a typical approach is to parameterize the region under investigation in terms of pixel values and employ regularization techniques in the reconstruction procedure to combat the illposedness and possible modeling errors. Tikhonov Regularization and Total Variation Method are wellknown examples using this approach, which were applied in medical tomography [95, 96]. Based on the assumption that the image should be smooth, Tikhonov Regularization method adds the smooth term as the Regularization. To preserve the edges better, Total Variation Method changes the power of gradient in the smooth term from 2 to 1. Those methods work well and it can solve the noisy data problem. However, regularization methods often contain a larger number of pixels to reconstruct and result in a large dimensional estimation problem, which is time consuming. Moreover, in the cases when only extremely limited number of data are available, this type of pixelbased approaches often produces poor reconstruction. To solve the large dimensional estimation problem in medical imaging, many successful numerical algorithms were created. Algebraic reconstruction technique(ART [97]) was the original reconstruction method used in medical imaging. The basic idea of this method is that makes an initial guess at the solution, compute projections based on the guess and then refine the guess based on the weighted difference between the actual projections and the desired projections. The method works, but is slow and susceptible to noise. With higher level of noise, the convergence is very slow, and is not guaranteed. When a full set of high signal to noise ratio (SNR) data are available, filtered back PAGE 112 97 projection (FBP) [98100] is the most important algorithm in 2D tomography. This algorithm is essentially a numerical implementation of the Radon inversion formula. And, some simple approaches avoiding singular integrals may be used. The discrete implementation depends on the scanning geometry, i.e. the way the data is sampled. This method is widely used in 2D tomography applications. For CT, FBP is the standard reconstruction algorithm. For PET [101, 102], the reconstruction method most commonly used is also the FBP algorithm. However, this method does not work well for low SNR data. The other class of tomography methods is statistical methods. Ordered subsets expectation maximization (OSEM) [103, 104]and maximum a posteriori (MAP) [101, 105] are the most famous ones. Statistical methods are very successful in PET and SPECT. However, the EMalgorithm has a limitation that is sensitive to data noise. And those algorithms require that the distribution of measured data is known. Besides those methods, Layer Stripping Method is also an interesting method. This method was first described by Cheney et al. [106]. This technique is unusual in that it reconstructs images from the periphery inward. It is claimed that for EIT, this method achieves a more accurate conductivity contrast than other methods. However, no available layerstripping algorithm works well on completemodel data. [107, 108] Data in tomography problems can sometimes be extremely noisy and limited. In there cases the algorithms above cannot accurately solve those tomography problems. Objectbased reconstruction methods provide an alternative approach. The basic idea of objectbased reconstruction method is to significantly reduce the number of unknowns in a reconstruction problem by representing the object field with a low order parametric PAGE 113 98 model and then recast the original reconstruction problem into a parameter estimation problem. These parameters typically contain internal intensity values and important geometric features of target objects, such as boundaries, which are commonly used for object classification and recognition. Miller et al. [109] have proposed a reconstruction approach which uses object boundaries to describe the object shape, and a set of basis functions to represent the texture of the object as well as the background. However, Millers method can only handle single simplyconnected objects, while in practice, the ability to handle multiple objects is often required. Moreover, the position and shape of the initial curve have significant impact on the convergence speed of the solution. Later, Feng [110] used Curve Evolution Approach to overcome those drawbacks and some successful applications on CT and ground penetrating radar (GPR) were achieved. But this algorithm requires the number of different objects, which are not available in many medical imaging applications. Also, the texture of each of object is required, it is not necessary to know the basis of the texture for piecewise smooth image. Yu el al. [111] proposed an nonlocal boundary information into the regularization method for edgepreserving tomographic reconstruction. Yu et al assume the image is smooth except at the boundary of each object, and function which is 0 near the boundary curves and 1 otherwise, is used to control the smoothness. This model is actually same as MumfordShah model. The function works in the same way as the derivative of heaviside function by Chan et al. [39]. Chan et al. [112] proposed level set and total variation regularization for inverse problems with discontinuous coefficients. Authors regularized the recovered image by the total variation norm, which indirectly controls both the jumps in the coefficient and the length of the level sets. Authors noted that in the standard level h h PAGE 114 99 set methodology, only the length part is regularized, which is insufficient for our purpose because we need to control the jumps as well, due to the illposedness of the underlying inverse problem. Another advantage of controlling the TVnorm instead of just the length of the level sets is that triple junctions are allowed to have arbitrary angles. Through the use of a novel variational multiple level set function approach borrowed from image processing, Chan et al. recover coefficients with multiple constant values associated with the discontinuous regions. If the identified image hasconstant values, then we just need level set functions. Moreover, the exact number of constant values need not be known a priori an upper bound suffices. If we use more level set functions than are actually needed, the redundant regions will disappear or merge with the other regions during the iterative process. n n2log Objectbased reconstruction methods provide a direction for tomography problems with noisy and limited data. And fortunately, the reconstructed images are piecewise statistical homogeneous in many applications. This property is true for most of the natural images. And sometimes, there is a reference image which can used as guidance to generate a new image. Those statistical properties and guiding image can be used to dramatically reduce the number of unknowns. Hero et al. [113] proposed a Bspline Model to interpolate the MRI boundary and applied this boundary information to reconstruct emission computed tomography (ECT). This method solves resolution mismatch problem, which can produce severe bias due to blurring of an organs tracer intensity across the organ boundary. This method works only for single organ case. However, we can learn from this paper that to avoid the bias caused by resolution PAGE 115 100 mismatch it is better to incorporate the MRI side information as an integral part of the image reconstruction process. 4.1.2 Motivation of the Proposed Framework Based on the ideas above, we propose a MumfordShah model based framework for tomography problems with noisy and limited data. MumfordShah model is region based and can detect the edges that separate homogeneous regions [111]. Therefore, we propose models based on MumfordShah model and using statistical properties to solve tomography problems with noisy and limited data. To overcome the complicated topology, multiphase level set method is applied here [112]. Hence this framework is named as MMSLS (Modified MumfordShah by Level Set). In framework, we assume there are some prior information we know about the reconstructed image, for example, smoothness of the boundary, piecewise constant/smooth/texture of the image, a guiding highresolution image. When one or more of those prior information are available, the framework can be applied. This work is based on the combination of [110113]. And we extends previous approaches in a number of ways. First, this framework works for piecewise smooth image through MumfordShah model. The Feng et al. [110] model need prior information of the texture basis which is not required in this model. The Chan et al. [112] model is based on total variation and the diffusion is anisotropic. We prefer the MumfordShad model because of the computation efficiency. One may argue that the computation for curve length term in MumfordShad model is time consuming, but some new techniques [43, 114, 115] for segmentation can be applied here to dramatically reduce the computation complexity. One may also argue that the anisotropic diffusion may be better than the isotopic diffusion in MumfordShad model. But for most of the medical images, PAGE 116 101 in the segmentation boundary, the image is isotropically smooth. In case the image is piecewise constant, the proposed framework is same as the models in Chan et al. [112] and Feng et al. [110]. The application of new numerical techniques in Chan [114] for segmentation is the second extension of this work. Because of application of the multiphase level set, which is originally used for segmentation in Chan and Vese [50], the proposed framework requires loglevel set when there are different objects, which is much less than level sets required in [110]. The application of fast sweeping method in [114] for MumfordShah image segmentation dramatically reduces the computation complexity. Third, not like in Hero et al. [113] when we incorporate the MRI side information as an integral part of reconstruction, there could be more than one object in the MR guiding image. The reason is the level set can take care of complex topology. n2 n n 4.2 Proposed Model 4.2.1 Description of the Proposed Model In the objectbased tomography problem of our interest, we want to determine the position and shape of the object according to measured data. A region based level set approach is adopted here and energy functional is developed to attract the zero level set to the object boundary: )()(),(12\2122HdxfMdxxfxRfEiii (41) Where 1 and 2 are parameters to balance those 3 terms, iis for different measurement. The first term is the fidelity term. As a tomography problem, the reconstructed function should satisfy xf ndxxfxRMii,...,2,1, i Because most of the natural images are smooth in each object, the second term is to make sure the PAGE 117 102 reconstructed image is piecewise smooth. The last term is to smooth object edges, since most of the natural objects have smooth edges. xf 1 H is the one dimensional Hausdorff measure, which generalizes the length notion for regular curves. This term penalize the length of the boundary and hence smooth the boundary. This model is actually a Mumfordshad model based implementation of formula (4)(6) in Yu and Fessler [111] Equation 41 gives a simpler numerical method than the Yu and Fessler model [111] wwlwx\ 4.2.2 Level Set Forms in Different Cases In this part, several cases of the image are considered. Specific models and corresponding EulerLagrange Equations are given. Several numerical issues are also discussed. 4.2.2.1 Binary field model Assume that there are only two intensities in an image, i.e. several homogeneous objects with same unknown constant intensity x are scattered in a homogeneous background with intensity y. Prior knowledge of number of objects is not required, and the background does not have to be continuous. Since now the function f is a twovalue function, the smooth term is no longer needed. Level set function is applied here to implement the Hausdorff measure 1 H [41]. Let xxxxx,0,0,0 where w is the set of regions with value x. Then the energy functional is PAGE 118 103 dpHMdpHRydpHRxyxEiiii221**),,( (42) where H is the Heaviside function 0,00,1zzH and the solution of f is. 0,0,)(yxzf Let and iiidpHRF iiidpHR1 G, notice H where zHdzd z in the sense of distributions then the EulerLagrange Equations are iiidivRMyGxFdtd*)( iiiiiMyGxFFdtdx)( iiiiiMyGxFGdtdy)( where 0,'H 4.2.2.2 Piecewise constant model This model is for the case that there are several homogeneous regions (phases) with different constant value of in the image and the maximum possible number n of phases is given. IX Since now the function f is piecewise constant, the smooth term is no longer needed. For n phases, level set functions nmlog Ri : are needed. The union of the zerolevel sets of i will represent the edges in the segmented image. The vector PAGE 119 104 level set function m ,...,1)(mH and the vector Heaviside function ),...,()(1HH whose components are only 1 or 0, are applied to define the segments or phases in the domain in the following way: two pixels (x1 ,y1) and (x2 y2) in will belong to the same phase or class, if and only if In other words, the classes or phases are given by the level sets of the function i.e. one class is formed by the set ))2,2((yx)(H ))1,1((HyxH consyxtan, Hyx, t_ IX ijk ijkX I IX 32 1111 HH H 321 0111 HH H111X 011X lnIIRXm21*iXzf II), Il XE( ) Hvector (one phase or class contains those pixels (x, y) of having the same value ) ,(yxH ). Characteristic function I and constant values can be defined as the followed sample. In case, there are 9 phases, then 3 level sets are needed to describe it. And use and to denote and i is for 1 j is for 2 k is for 3 i, j, k can only be either 0 or 1, 1 means H, 0 means 1H. e.g. Correspondingly is the conductivity of area 1 >0 2 >0 3 >0 is the conductivity of area 1 <0 2 >0 3 >0. Then the energy functional is miildpHMdp122 And the solution of f is iz ( In case, there are 4 phases, the energy functional is PAGE 120 105 ),(IIXE dxdyHHRXdxdyHHRXiii211021111 2221002101111iiiMdxdyHHRXdxdyHHRX (43) 21HH The following is one numerical method for Equation 43. Define iidxdyHHRT2111 dxdyHHRTii21101 dxdyHHRTii21011 dxdyHHRTii210011 jiiiijijMTXDIF, then EulerLagrange Equations are .11111DIFdivt 20010201111HXXHXXRii Similarly, .22122DIFdivt 10001110111HXXHXXRii For 11X 0110000010110101111TTMTXTXTXTX Similar formulas for 000110,,XXX Hence PAGE 121 106 MTMTMTMTXXXXTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT 00011011000110110000010010001100000101011001110100100110101011100011011110111111 4.2.2.3 Piecewise texture model In some cases, the intensity of each object is not constant but follows certain texture, i.e. the intensity function of object I is jjIjIIxBxfjIBI where provide the basis for describing the object I texture [109, 116, 117], and are the corresponding expansion coefficients for the object I. Suppose the texture basis for those objects are known and the upper bound of the number of objects is also known, then the energy functional is xBjI x j miilnIljjIjIIlIjIdpHMdpxBREm12122),( (44) Now, the unknowns are the coefficients of those bases. And the solution of f is IjjIjIzzBzf)( 4.2.2.4 General piecewise smooth model In case no prior information about number of objects and no prior information about intensity, the problem can be solved using only two level set functions. The idea of twolevel set is from Vese and Chan [41]. They successfully applied multiphase level set to solve segmentation problem. Because there is no prior information available, the advantage of reducing unknowns from number of pixels to number of parameters is lost. However, because of the isotropic diffusion, it is still better than total variation method PAGE 122 107 on speed. Furthermore, this model can do segmentation and reconstruction simultaneously. Define function 00,00,00,00,21212121andfandfandfandff Using the Heaviside function, the relation between f and the level set functions can be expressed by 212121211111HHfHHfHHfHHff Then the energy functional is ),,(21 uE llldxdyHHfRdxdyHHfR21211 222121111lllMdxdyHHfRdxdyHHfR (45) dxdyHHudxdyHHu212121211 dxdyHHudxdyHHu21212121111 2212HH And define dxdyHHuRdxdyHHuRDIFiiii21211 EMMdxdyHHuRdxdyHHuRiiii 2121111 We deduce the EulerLagrange equations for u uDIFRiii1 in 0,021 0nu on 0,021 and 0,021 where PAGE 123 108 iiiMdxdyuRDIF Similarly we have the EulerLagrange equations for u uu,, uDIFRiii1 in 0,021 0nu on 0,021 and 0,021 uDIFRiii1 in 0,021 0nu on 0,021 and 0,021 uDIFRiii1 in 0,021 0nu on 0,021 and 0,021 the EulerLagrange equations for 1 and 2 221112111HuDIFuRdivvEtiii 2212211HuDIFuRHuDIFuRiiiiii 2211HuDIFuRiii Similarly 12122222HuDIFuRdivtiii 1211211HuDIFuRHuDIFuRiiiii 1211HuDIFuRiii 4.2.2.5 Model for tomography with guiding image In some medical imaging applications, image of one physical property is given and the image of another physical property is wanted. For example, the MRI guided Noise Tomography (MRNT). In MRNT the MR image, which is the map of proton density, is PAGE 124 109 given, and the map of conductivity is wanted. In this case, geometry information, like shape and position, of each phase can be achieved by segmenting the given guiding image. Then using that information to compute wanted value for each object. If we do segmentation and reconstruction separately, resolution mismatch may occur. This can produce severe bias due to blurring of an organs tracer intensity across the organ boundary. The proposed model below combines segmentation and reconstruction together and it is beneficial to both segmentation the given image and reconstruction for the new image with noisy and limited data. Denote u given guiding image Iu the segmented guiding image for class I Parameters 1 2 trade off between terms Other notations are same as above. By assuming the wanted value is piecewise constant and there are at most n phases, then the energy functional is ilInIImiiinIIiIIIIdxdyuuHMdARXuXEmm,2212112221*),,( (46) And the solution of f is iizXzf ,)( The first term is still the fidelity term, this term is mainly for reconstruction. The second term penalize the length of the boundary. This term works for both segmentation and reconstruction. The third term is for segmentation. This term is borrowed from [39], which is for MumfordShad model based segmentation. Here we assume the guing image is also piecewise constant. This model solves the tomography problem with guiding PAGE 125 110 images. It is also simultaneously capable of not only segmentation for both guiding and reconstructed images but also denoising the guiding image. 4.3 Application on Noise Tomography 4.3.1 Introduction of Noise Tomography Noise Tomography, which aims to determine the distribution of the conductivity in the sample, is a new noninvasive medical imaging technique invented by MRI Devices Cooperation. The NT technique is designed to use the correlations in the detected electronic thermal noise in an RF probe array and the relationship between conductivity and the noise power coupled between the sample and probe to measure the electrical conductivity distribution within the sample. Many applications of conductivity distribution in diagnostic imaging have been documented. For example, gastrointestinal and oesophgeal function [4, 5], hyperor hypothermic treatment of malignant tumors [6, 7], imaging of the head [8, 9, 118], pulmonary function [10], cancer detection [11, 12, 119], measurements of cardiac output and investigation to locate the focus of epileptic seizures [13]. 4.3.2 Electronic Noise Historically, the electronic thermal noise generated within the sample has not been encoded during the Magnetic Resonance Imaging (MRI) process. The pulse sequences and gradient encoding do not encode the noise but the new methods of partially parallel imaging (SENSE, SMASH, etc.,[16] [15]) demonstrate that the probe sensitivity patterns are encoding functions for the MRI signal. The sensitivity patterns also encode the noise to approximately the same extent as they encode the signal. This is true because the electronic thermal noise of the tissue is normally dominant and the physical size and structure of the coil dictate the acquisition site of the 'body' noise. Although the regions PAGE 126 111 of MRI signal acquisition and noise acquisition are similar they are not the same, because the probe's magnetic field is associated with MRI signal acquisition whereas the probe's electric field is associated with noise acquisition. If several probes and receiver channels are used in MRI signal acquisition, the noise in a given channel can be assumed to be significantly localized by the electric field profile of the probe associated with that channel. This is a critical concept because it indicates that noise can have an associated information content. The association exists because the noise originates in tissues of the body and is acquired in a predictable and repeatable way by means of a local probe. The NT technique refines existing technology by measuring only the covariance of noise within a multiple receiver system. This is important because, in principle, the amount of nonsample noise is largely irrelevant since correlating the two noise outputs eliminates the uncorrelated noise by averaging. To determine the distribution of the conductivity in the sample, the NT technique uses the nonresonant electronic thermal noise of the tissue sample, detected by an array of RF probes, and the known, or calculated overlap of the electromagnet fields of those probes. Electronic thermal noise is a basic property of all matter and has been well characterized in both onedimensional (Johnson/Nyqusit) and threedimensional (Black Body Radiation) cases. This NT technique is based on a combination of the fluctuation dissipation theorem and the principle of reciprocity. The objective of the NT technique is to use the noise correlation between the different channels in an array of RF probes to plot the conductivity distribution within the sample. When an oscillating RF probe is placed near a conducting sample, energy in the probe is dissipated in the sample and the impedance of the probe is changed. This change PAGE 127 112 in impedance is a function of the internal impedance of the probe and of the coupling between the sample and the probe. Thevenin's theorem shows that the thermal noise power detected across a circuit is proportional to the equivalent resistance of the circuit. Therefore, if an RF probe is moved with respect to the sample and the resulting change in noise power is measured, the change in coupling between the sample and the coil may be determined. A similar effect on the impedance of the circuit and the detected noise power will occur if the probe is kept in a fixed location and the conductivity of the sample changed. Utilizing these principals, a onedimensional conductivity map of the sample may be constructed by simply moving the probe at a fixed distance from the surface of the sample and measuring the changes in the sample noise. Adding an array of coils, which are measured in parallel, will increase the sampling rate and therefore the imaging rate. This arrangement, however, will determine only one piece of information per coil position. In order to determine n independent parameters about the sample, a minimum of n independent coil positions must measured. 4.3.3 Intercoil Noise Correlation In order to increase the sample information obtained from any given probe position and to reduce the effects of nonsample noise sources, the noise correlation between different coils in an array of coils will be measured. Using the normalized noise covariance will provide up to n(n1)/2 different pieces of sample information from each position of the coil array and will reduce the effects of nonsample related detector noise sources. The noise correlation between channels is proportional to: dVrErErkj PAGE 128 113 Where E I E k are the electric field resulting from the two probes at the point r in space and (r) is the conductivity at the point r in space. 4.3.4 Relation with Tomography Models For this Noise Tomography problem, this distribution of conductivity is wanted. Hence is the unknown function. The noise correlation dVrErErMkjjk rErErRkjjk generated by coil i and coil j is the collected data, is the fixed function determined by the character of the machine. The goal here is to find conductivity map with given and jkM jkR 4.3.5 Experiments and Results All programs were produced by using Matlab, and ran on a Compaq PC with 1GHz CPU speed. 4.3.5.1 Experiment in one dimension The phantom sample is a disc with a diameter of 16 cm and a width 5 cm. Figure 41 shows the result of the experiment using MMSLS without using any prior information. The red curve shows where the phantom disc actually is, and the blue curve shows the result by using the MMSLS to binary field model. The unit of xaxis is the centimeter. By using the real position and disc width, the ideal norm of RXM is calculated as 28.2831. The resulting norm, using MMSLS, is 33.6771. The CPU time for this reconstruction was 6.875 seconds. PAGE 129 114 Figure 41. Experiment result in 1 dimension 4.3.5.2 Experiment in two dimensions In this experiment, there are only 231 equations with which to reconstruct the image. But to achieve resolution of 2mm2mm, 10000 unknowns need to be solved. Furthermore, the measured data itself is very noise because the data is noise correlation. Therefore this inverse problem is extremely ill posed. Figure 42 shows the experiment result. Figure 42A is the MR image of the phantom. The phantom is constructed of a 197mm (7.75) ID by 180mm long acrylic tube that is filled with a solution of Cu 2 SO 4 (2.0 grams/Liter) and NaCl (4.5 grams/Liter). The conductivity of this solution is on the order of 0.8S/m. There are three smaller acrylic tubes, each of which is 180mm long, contained within the 197mm OD tube. Figure 42B shows the photo of the phantom. If the container is ignored, the digitized real image should look like the image in Figure 42C. Figure 42D is the result by using MMSLS to binary field model. The reconstructed image has a resolution of 1.97mm1.97mm. The CPU time is 25 seconds. The solved relative conductivity of saline solution and water is 1.0 and 0.2, respectively. The black dots (around the periphery of the blue area) are the result of automatic segmentation. If it is assumed that the relative conductivity of distilled water is 0, and the relative conductivity of saline solution is 1, then the norm of RXM is 1.9712e006. To compare the result of MMSLS and traditional algorithms, Figure 42E and 42F show the PAGE 130 115 result by using pixelwise method Total Variation and ART. The reconstructed image (Figure 36) by Total Variation Method has a resolution of 3.97mm3.97mm. The CPU time is 189 seconds. The reconstructed image (Figure 42F) by ART has a resolution of 14.1mm*14.1mm. The CPU time is 32 seconds. A B C D E F Figure 42. Experiment result in 2 dimensions. A) MRI Image of Phantom, B) NT Phantom Photo, C) Digitized NT Phantom, D) Reconstructed image by MMSLS, E) Reconstructed image by Total Variation, F) Reconstructed image by ART In this experiment, the model for Total Variation algorithm is dxxfMdxxfxRfElll)()()(22 subject to ub lb xf PAGE 131 116 where ub and lb are the upper bound and least bound of xf and are given as prior information. And the method introduced in [97] is applied for the ART algorithm. 4.3.5.3 Experiment in two dimensions with guiding image In this experiment MR image is used as the guiding image and tried to generate more accurate conductivity map with the same tomography data in the previous experiment. The model for tomography with guiding image is applied in this experiment. Because there are 3 kinds of objects: Cu, and containers, 2 level sets are used to describe the geometry properties of those 3 objects. 42SO NaCl Figures 43 display the results based on the ideal guiding image. Figure 43A is the ideal guiding Image. Figure 43B through D is the segmentation results of the ideal guiding Image according to vector level set function defined in 2.3.2. Figure 43E is the reconstructed conductivity distribution map. The solved relative conductivity of saline solution, water and plastic is 1.0000, 0.2805 and 0.0399. Figures 44 display the results based on the real MR image of the phantom. Figure 44A is the real MR Image. Figure 44B to 44D is segmentation result based on the real MR Image. Figure 44E is the reconstructed conductivity distribution map. The solved relative conductivity of saline solution, water and plastic is 1.0000, 0.1873 and 0.1230. With the same tomography and an extra guiding image, the result 44 E is much better than 42 D. Furthermore, even though the MR image used here is very nonuniform, but the MMSLS to Tomography with guiding image model still can segment it accurately. It can be seen that this model can not only make the reconstruction better with the guiding image but also can segment the guiding image better with the tomography data. PAGE 132 117 A B C D E Figure 43. Experiment result in 2 dimensions with ideal guiding image. A) Ideal MR Image, B) Segmented Object 1, C) Segmented Object 2, D) Segmented Object 3, E) Reconstructed conductivity distribution map A B C D E Figure 44. Experiment result in 2 dimensions with real MR guiding image. A) real MR Image, B) Segmented Object 1, C) Segmented Object 2, D) Segmented Object 3, E) Reconstructed conductivity distribution map 4.3.5.4 Experiment in two dimensions with guiding image by fast sweeping In this experiment, the fast algorithm for Level Set Based Optimization [43] is applied to MMSLS binary field model. The computation time is significantly reduced. And the reconstruction result is accurate. PAGE 133 118 Figure 45 shows the result by using ideal guiding image. The CPU time is 2.047 seconds. The solved relative conductivity of saline solution, water & plastic (water and plastic can not distinguished by binary field model) is 0.918, 0.289. Figure 45. Reconstructed image by fast sweeping 4.4 Conclusion A novel algorithm MMSLS for Tomography problem with extremely noisy and limited data is introduced. This algorithm is flexible. It can be used for any geometry with or without prior information. It is very easy to change the model for different cases. Furthermore, this approach has no special requirements. Hence it may be easily applied to a variety of Tomography problems. The MMSLS algorithm has proven itself to provide accurate results. The reconstructed image by MMSLS has higher resolution and better visualization then conventional method. And the object property (relative conductivity in our experiments) is also accurately calculated. Because the MMSLS algorithm is object based, not pixelbased, this approach, for an equivalent output resolution, uses no more than 5% of the time used by traditional pixelbased methods in the 2D experiment. If fast sweeping is applied, the reconstruction PAGE 134 119 time can be no more than 0.3% of the time used by traditional method. Actually, it is also the key idea of this work that change a pixelbased problem with numerous unknowns to be a object based problem with only few unknown parameters. Moreover, the reconstructed image, using MMSLS, is smooth, and the edge is perfectly preserved. The visualization is therefore much better than when using other approaches. And mesh generation is no longer crucial, because a finer mesh can be generated with only a small addition to the calculation. Another impressive advantage of this approach is that it can accomplish segmentation and reconstruction at the same time. Actually, the zero level set is the desired segmentation. This is a very useful feature in clinical application. PAGE 135 CHAPTER 5 FIBER TRACKING BY FRONT PROPAGATION In this chapter, front propagation method is applied to fiber tracking for DWMRI. The proposed method solves multidiffusiondirection, branches and tracking back problems existing in fiber tracking problems. 5.1 Introduction Diffusionweighted MRI (DWMRI) is a technique used for relating image intensities to the relative mobility of endogenous tissue water molecules. Two equal and opposite large magnetic field gradients are applied before the data are acquired. In the absence of any motion of the water molecules, the magnetic spins will be dephased by the first gradient, and then completely rephased by the second gradient. However, because the actual motion of the water molecules is that of a random walk, the second gradient will not completely rephase the spins. The signal intensity will therefore be exponentially attenuated proportional to the mean diffusion length. Areas with relatively high mean diffusion length will appear dark on the MRI images relative to areas with low mean diffusion length. Diffusion Tensor Imaging (DTI) is an application of diffusion imaging where several sets of DWMRI are acquired with the diffusion gradients applied in different directions. This technique enables the detection of diffusion anisotropy in various mediums such as brain white matter. The diffusive properties of an anisotropic medium can be described with a 3 X 3 symmetric tensor. The eigenvalues of the diffusion tensor are the diffusion coefficients in the three principal directions of diffusivity, and the 120 PAGE 136 121 eigenvector corresponding to the largest eigenvalue is the main diffusivity direction in the medium. It is then possible to construct image maps from the main diffusivity direction, and for various anisotropy indices calculated from the eigenvalues. The anisotropy indices range from 0, in the case of completely isotropic diffusion, to 1, in the case of completely anisotropic diffusion. And the index is called fractional anisotropy (). Let fa 1 2 3 be the three eigenvalues and 1 2 3 be the three eigenvectors; fractional anisotropy may be defined as 21 22 2 23222312321fa Experimental evidence has shown that water diffusion is anisotropic in organized tissues such as muscles or brain white matter. In the last decade, the quantitative description of this anisotropy with DTI has become well established in the research environment and its first applications in the clinic are now being reported. These applications employ both the directional anisotropy that can be measured by DTI as well as removal of this anisotropy through the use of the tensor trace. For example, DTI is presently being explored as a research tool to study brain development, multiple sclerosis, amyotrophic lateral sclerosis (ALS), stroke, schizophrenia and reading disability, while trace imaging has become an essential part of clinical acute stroke assessment. As first shown by Basser et al. [120], the diffusion ellipsoid obtained from DTI can not only provide a quantitative orientationindependent measure of diffusion anisotropy, but also the predominant direction of water diffusion in image voxels. However, it was not until the end of the decade and the beginning of the new millenium that the first successful in vivo fiber tracking results were published [121123]. The reason for this time lag between the availability of the tensor diagonalization techniques that provide the vector PAGE 137 122 information and the actual mapping of the fibers is due to the inherent complexity of connecting these macroscopic voxelbased vectors in a reproducible threedimensional manner. At present, new tracking algorithms are being developed rapidly. In order to better utilize this promising technology, it is important to understand the basis of the anisotropy contrast in DTI and the limitations imposed by using a macroscopic technique to visualize microscopic restrictions. These basics are briefly reviewed, while a more comprehensive overview is presented by Beaulieu et al. [124] Figure 51. White matter structure and its relationship with the information provided by DTI. A) anisotropy map, B) colorcoded map, C) vector map, D) a pixel contains bundles of axons and neuroglial cells, E) axon filled with neuronal filaments PAGE 138 123 5.1.1 Understanding the Basis of DTI DTI essentially provides two types of information about the property of water diffusion; the extent of diffusion anisotropy and its orientation. By assuming that the largest principal axis of the diffusion tensor aligns with the predominant fiber orientation in an MRI voxel, we can obtain 2D or 3D vector fields that represent the fiber orientation at each voxel. The 3D reconstruction of tract trajectories, or tractography, is a natural extension of such vector fields. Before further describing tractography, it is important to discuss what exactly DTI measures and how the data relate to the tract trajectories we are trying to derive from the measurement. Figure 51 is the Figure 1 in [125]. This figure explains the basis of DTI. As we mentioned above, in DWI, areas with relatively high mean diffusion length will appear darker on the MR images compared to areas with low mean diffusion length. Figure 51 A) is a slice of anisotropy image based on DWI. DTI is a set of DWI acquired with the diffusion gradients applied in different directions. In typical DTI measurements, the voxel dimensions are in the order of 1 mm and DTI measures the averaged diffusion properties of water molecules inside it. Figure 51B and C show the tracts that are running along various directions that are large enough to discern visually. Figure 51B uses color to show the directions. Please notice the colored arrows beside B. The red indicates fibers running along the rightleft direction, green inferiorsuperior, and blue anteriorposterior. Figure 51C is the vector field map. Very often, image resolution is sufficiently high for the white matter tracts to contain several voxels. The white matter tracts, in turn, consist of densely packed axons (neuronal projections) in addition to various types of neuroglia and other small populations of cells. Figure 51D shows a cell. Inside the voxel, water molecules are distributed between these cell types and the PAGE 139 124 extracellular space (805% are intracellular). Thus, even a voxel within a single white matter tract consists of very inhomogeneous environment, and water molecules are likely to experience high anisotropy as judged from the cytoarchitecture of the axon. Inside an axon, water molecules are surrounded by high concentration of neuronal filaments, which are polymers of protein molecules. Figure 51E shows the neuronal filaments in an axon. Each monomer protein has a molecular weight of 50 kDa. The neuronal filament is far larger than that, and multiple filament bundles are densely packed in the axon. The axonal membrane as well as the wellaligned protein fibers within an axon restricts water diffusion perpendicular to the fiber orientation, leading to anisotropic diffusion. Myelin sheaths that surround the axons may also contribute to the anisotropy for both intraand extracellular water. These contributions have been studied in detail by Beaulieu [124]using nonmyelinated axons, showing that the contribution of myelin may be significant, but that the axonal contribution dominates. In addition to the inhomogeneity in terms of the cellular environment of water molecules, inhomogeneity of axonal orientation within one voxel is also an important factor to be considered. Despite this multidirectional environment, previous DTI studies have shown that water diffusion in many regions of the white matter is highly anisotropic (Figure 51A) and that the orientation of the largest principal axis aligns to the predominant axonal orientation. However, when studying axonal architecture using DTI, it is very important to understand the limitations arising from the inhomogeneity of the water environment. First, the conventional DTI data acquisition and processing methods may not be able to properly handle a voxel containing more than one population of axonal tracts with different orientations. This issue will be discussed in more detail below. Second, DTI PAGE 140 125 cannot provide information on cellularlevel axonal connectivity. Multiple axons from individual cells may merge into or branch out from one voxel. Within a voxel, cellularlevel axonal information of multiple compartments is averaged. The third important limitation is that afferent and efferent pathways of axonal tracts cannot be judged from the direction of water diffusion. To resolve the problem of multiple fiber orientations within a single voxel, high angular resolution diffusion imaging (HARD) was proposed by Tuch et al. in [126, 127]. In HARD, a geodesic, high value diffusion gradient sampling scheme is applied. In regions of fiber crossing the diffusion signal exhibited multiple local maxima/minima as a function of diffusion gradient orientation, indicating the presence of multiple intravoxel fiber orientations. The multimodality of the observed diffusion signal precluded the standard tensor reconstruction, so instead the diffusion signal was modeled as arising from a discrete mixture of Gaussian diffusion processes in slow exchange, and the underlying mixture of tensors was solved for using a gradient descent scheme. The multitensor reconstruction resolved multiple intravoxel fiber populations corresponding to known fiber anatomy.[127] b We will first discuss recent 3D tract reconstruction techniques that employ 3D vector fields obtained from DTI measurements. Because of the advantages of HARDMRI, we propose our fiber reconstruction method for HARDMRI. 5.1.2 Existing Fiber Reconstruction Methods Currently, there are several different approaches to reconstruct white matter tracts, which can be roughly divided into three types. Techniques classified in the first category are based on line propagation algorithms that use local tensor information for each step of the propagation. The main differences among techniques in this class stem from the way PAGE 141 126 information from neighboring pixels is incorporated to define smooth trajectories or to minimize noise contributions. [128130]The second type of approach is front evolution methods. In this kind of method, all possible directions are tried, the result obtained is first the front, and then the track corresponding to that front [131, 132]. The basic ideas of those two kinds of methods are introduced here. 5.1.2.1 Line propagation algorithms The most intuitive way to reconstruct a 3D trajectory from a 3D vector field is to propagate a line from a seed point by following the local vector orientation. Notice, if a line is propagated simply by connecting pixels, which are discrete entities, the vector information contained at each pixel may not be fully reflected in the propagation. In the simple example illustrated in Figure 52A, axonal tracts are running along 30from the vertical line. When applying the discrete pixel connection approach, a judgment has to be made about which pixel should be connected (for instance, is the 30vector angle pointing at pixel (1, 2) or (2, 2)?). No matter what the judgment is, it should be clear that this simple pixel connection scheme cannot represent the real tract even in such a simple case. Because in this simplest method, only the node (discrete integers in coordinates) has direction and the fiber only passes through the nodes, hence the field is discrete number field. The simplest way to convert the discrete voxel information into a continuous tracking line is to linearly propagate a line, in a continuous number field. It is called continuous field because it is assumed that the whole voxel has the same direction and any point (continuous in coordinates) has direction. This conversion from the discrete to continuous number field is shown in Figure 52B. PAGE 142 127 Figure 52. Linear linepropagation approach. Doubleheaded arrows indicate fiber orientations at each pixel. Tracking is initiated from the center pixel. A) discrete number field The shaded pixels are connected, B) continuous number field. A line instead of a series of pixels is propagated. C) interpolation approach to perform nonlinear line propagation. Large arrows indicate the vector of the largest principal axis. For every step size, a distanceweighted average of nearby vectors is calculated. The simple linear approach demonstrated in Figure.52A, B can be modified to create a smooth (curved) path, which should be more accurate when the curvature of the constructed line is steep with respect to imaging resolution. A schematic diagram of a simple interpolation approach that can achieve such a smooth path is shown in Figure. 52C. In this example, the line is propagated with a small, predefined step size. Whenever it moves to a new coordinate, a distanceaveraged vector orientation is calculated. In this most simple example, the average distance between two pixels closest to the new coordinate was seed to draw a smooth line between pixels. In more rigorous approaches, diffusion tensors, rather than the vector of the largest principle axis, are interpolated at each coordinate as the line is being propagated. Notice Figure 52C used a different vector field than A and B; the reason is that the result of interpolation approach is same as continuous propagation method for the vector field in B. PAGE 143 128 With interpolated diffusion tensors, we have the following propagation method. Let white matter fiber tract trajectory be represented as a 3D space curve, i.e., a vector, r(s), parameterized by the arc length, s, of the trajectory. The Frenet equation describing the evolution of r(s) is stdssdr (51) where tis the unit tangent vector to s sr at s. These vectors are depicted in Figure. 53. Figure 53. White matter fiber trajectory as a space curve s r. The local tangent vector, is identified with the eigenvector, 1st 1 1sr associated with the largest eigenvalue of the diffusion tensor, D, at position 1s r Since the normalized eigenvector, 1 associated with the largest eigenvalue of the diffusion tensor, D, lies parallel to the local fiber tract, the direction in coherently organized white matter. Within acceptable experimental error, several groups have confirmed this to be true in the human heart. A key idea in this fiber tract following algorithm is to equate the tangent vector, t(s), and the unit eigenvector, 1 calculated at position, r(s): srst1 (52) PAGE 144 129 Combine (1) and (2) we have srdssdr1 with subject to 00rr (53) Where specifies a starting point on the fiber tract. This procedure can now be repeated starting at the new point, and can be iterated to predict the location of discrete points along the fiber trajectory, 0r 1sr sr 5.1.2.2 Front evolution methods Front evolution methods are based on the stepwise evolution of a front from an initial seed point. The speed and direction of evolution of the front are calculated from the local diffusion information, such that evolution is fastest along the direction of fastest diffusion. Parker, WheelerKingshott and Barker first introduced the front evolution method in [131]. The method in [131] is based on the fast marching method [133]. There are 4 steps in this method. First step is the evolution of a front from a seed point using a variant of the fast marching method at a rate governed by a speed function. rFrrrTr'' T Where is what we want to find, which is the propagation front.' rT r is chosen to be the neighbor for which is most closely aligned with the direction of normal to the front, is the speed function, which is defined as 'rr rn rF ',',min1111111rrrnrrnrrF (54) or rnrrFrF ','min122 (55) PAGE 145 130 Where r1 is the dominant tensor direction at pixel r The second step is the generation of paths from all points in a given data set to the seed point. dTrT rAmin where is the path. The third step is the creation of a connectivity metric, which assigns a measure of connectivity for a given path based on the worst case along the length of the path of the property to which they are sensitive. Define TTF1minmin1 (56) or 12minw (57) The last step is the selection of a subset of the paths as being reasonable pathways of connection, which are points with a high connectivity to the seed point. Tournier et al. proposed another front evolution method in [132]. Although this technique holds certain similarities to the Fast Marching Algorithm proposed by Parker et al., here the front is not evolved in a voxelbyvoxel manner, and there is no gradient descent to establish the actual tracks. Hence this method is faster than Parkers. 5.1.3 Several Considerations in Fiber reconstruction For any fiber reconstruction method, there are several general considerations. 5.1.3.1 Termination criteria Line propagation or front evolution must be terminated at some points. The most intuitive termination criterion is the extent of anisotropy. In a low anisotropy region, such as gray matter, there may not be a coherent tract orientation within a pixel and the orientation of the largest principal axis is more sensitive to noise errors (for isotropic PAGE 146 131 diffusion, the anisotropy information is dominated by noise and becomes purely random). The fractional anisotropy of the gray matter is typically in the range 0.1.2. Therefore, one simple termination approach is to set the threshold for tracking termination at 0.2. Another important criterion is the angle change between pixels. For the linear line propagation model, large errors occur if the angle transition is large. Even for the interpolation approach, it should be noted that the diffusion tensor calculation assumes that there is no consistent curvature of axonal tracts within a voxel. The presence of curvature violates the assumption that the diffusion process along any arbitrary axis is Gaussian, thereby invalidating the routine tensor calculation. Therefore, it is preferable to set a threshold that prohibits a sharp turn during line propagation. The significance of this angletransition threshold depends on the particular trajectories of tracts of interest and the image resolution. An image resolution of 1 mm, for example, is high enough to smoothly reconstruct the curvature of trajectories of major tracts in the brainstem and the corticocortical association fibers connecting the functional regions of the brain. Under such favorable conditions, the angles between connected vectors are small and the termination criteria are dominated by the magnitude of the fractional anisotropy. However, for smaller tracts in environments that are structurally highly convoluted, such as subcortical Ufibers, the same resolution may be too coarse to smoothly represent the trajectories and, thus, analysis for the errors in tensor calculations and fiber tracking due to the curvature becomes more important. While at this moment there are no comprehensive simulation studies in this regard, some preliminary data on the effect of the turning radius on the tracking results in particular examples can be found in recent reports. PAGE 147 132 5.1.3.2 Effect of noise 3D vector field obtained from DTI contains noise and, as a consequence, the calculated vector direction may deviate from the real fiber orientation. One of the drawbacks of line propagation methods is that noise errors accumulate as the propagation becomes longer. The extent of these errors as a function of signaltonoise ratio (SNR) has been evaluated for linear and interpolated models. The results show a dependence on the shape of the trajectory, anisotropy, resolution and the particular interpolation method used. The noise effect can be reduced using smoothing or interpolation techniques, which averages vector or tensor information among neighboring pixels with a cost in the reduction in effective resolution. More information can be found in [134]. 5.1.3.3 Branching White matter tracts often have extensive branching, which renders tracking computationally complex. For example, bifurcation of a line during propagation is already a mathematically involved issue. From a programming point of view, merging two lines rather than splitting a line into two can much more easily handle this problem. In this approach, tracking is initiated from all pixels within the brain and tracking results that penetrate the pixel of interest are kept. In other words, instead of using the pixel of interest as a seed pixel, all pixels in the brain are used as seed pixels. 5.1.3.4 Algorithms to interpolate a tensor field In case the resolution of MRI image is low, it is necessary to interpolate the tensor field to get the smooth fiber. There are several ways to generate the interpolated tensor field. Since the tensor field is from MR images, the first possible object to interpolate is the raw diffusion weighted MR images [123], and then use the interpolated high PAGE 148 133 resolution MR images to calculate the tensor field. This is the most defensible method; however there are considerable computational expenses. The second choice is matrix interpolation, wherein the tensor matrix, D, is calculated once per sample point, and then interpolated componentwise to produce a tensor at intermediate locations. Were the components of the tensor matrix simply linear combinations of the measured channels, then matrix interpolation would be equivalent to raw image interpolation. However, in reality, the components of the tensor matrix are not linear combinations of the measured channels. The last choice is to interpolate eigenvalues or eigenvectors. Kindlmann, Weinstein and Hart [135] called this kind of interpolation eigensystem interpolation. This has immediate relevance for the barycentric map methods of assigning color and shading, and the littensor shading model. This style of interpolation has the tremendous benefit of doing all the eigensystem calculations once as a preprocess, and storing the results at each dataset point. Despite its obvious benefits, there is a subtle issue of correspondence, which complicates eigensystem interpolation. 5.1.3.5 Display of experimental tracks Track display is important for understanding the 3D shape, anatomical location, relative spatial relationships, and internal fiber structure of pathways. The precise anatomical location of tracks can be visualized using 2D anatomical overlays of tracks onto anatomical background images. If any bead in a pathway is found to lie inside a given voxel, that voxel is assigned the designated color of that pathway. Background images were either I0 (from DTMRI) or coregistered 3D T1weighted MPRAGE images. Visualization of wholebrain track data in 3D can be difficult because of overlapping tracks and limitations in computer graphics. Thus tracks representing PAGE 149 134 pathways of interest were selectively extracted and displayed. From the wholebrain data, tracks are extracted if they pass into a spatial selection volume (SSV). This selection procedure is passive and does not alter the quality of the tracks. The SSVs can be combined with logical operations to obtain a pathway or its component subsets. 5.1.4 Limitations of Existing Methods Linear propagation methods suffer from three major disadvantages. Firstly, there is no (or at best limited) natural mechanism to allow for the branching of tracts from a single start point (an anatomically reasonable occurrence at the spatial resolution of DTI data, is seen, for example, in the corona radiata); (which means) meaning that connectivity is restricted to a representation as a onetoone mapping between points in different regions. It is possible to produce a limited impression of diverging pathways with streamlinelike methods by utilizing multiple interpolated start points within the start voxel. However, the density of traces after the point of branching or fanning will reduce as the traces diverge when using these techniques, thus undersampling (in general in a nonuniform manner) subsequent connected regions. This change in density will not in general represent any reduction in packing density of the underlying axons. Therefore, this kind of method cannot take advantage of high angular resolution data. The second major restriction is that there is no attempt to determine how reasonable, or likely, any path is in representing a true pathway of connections as demonstrated by Pierpaoli et al. [136]. Furthermore, this approach suffers from a high susceptibility to noise; the error caused by noise could be accumulated during the propagation. Front propagation methods overcome those shortcomings. First, it allows for the branching of tracts from a single start point, because the front can propagate to any number of directions. Second, the fiber is along the fastest direction of the front and the PAGE 150 135 speeds in other directions are also known, hence it is easy to determine how reasonable or likely, any path is in representing a true pathway of connection. Third, because the evolving element is the front instead of one point on a line in the frontpropagation methods, the error caused by local noise may be remitted by propagation from other points. However, in Parker et al. [131], the front evolution result is a time map, T, from the seed. To find the paths of connection, a gradient descent has to be established. This computation for determination of connection paths is not simple. Tournier et al. [132] introduced another front propagation method using a fiber orientation probability density function. This method is not the voxelbyvoxel manner, hence the front can evolve to different directions from one seed simultaneously, but the speeds in different directions are different. The speeds are decided by fiber orientation probability density functions. It should be more accurate because it is not voxelbyvoxel manner, however the computation complexity could be dramatically increased. Another drawback of the method of Tournier et al. [132] is that the result is only the front, no true path of connection is available. 5.2 Methods 5.2.1 Motivations of the Proposed Method In this work, we are trying to (partially) solve the problems in existing fiber reconstruction methods. Since front propagation method can solve the main problems in linear propagation method (see 5.1.4), we also apply front propagation method to track the fibers. To avoid computation complexity [132], we still use voxelbyvoxel method. That is, we assume the fiber always passes through the nodes. If resolution is low and kinks in PAGE 151 136 the reconstructed path are evident, interpolation techniques should provide an advantage in accuracy. On the other hand, if the resolution is sufficiently high, the noninterpolating techniques offer faster calculation. According to Mori and Ziji [125], when using DTI images with 1 1 3 mm resolution, the stretch of the corticospinal tract between the pons and the motor cortex has a rather small average angle of 7.8 5.8(SD) between connecting pixels, for which the effect of interpolation should be minimal. In case the resolution is not high enough, we need to do interpolation according to 5.1.3.4 before fiber reconstruction. To find the true connection path, the closest point method [31] is applied in our algorithm. In this method, the closest path from the seed to any points in propagated region of the front is recorded during the propagation; then there is no more effort to find the connection path after front propagation. Here we assume the reconstructed fiber tends to minimize its length. Therefore, the true connection path can be found in no effort in our algorithm and the calculation for gradient descent [131] is avoided. Also, a novel definition for propagation speed is applied in our algorithm. According to recent literature [125, 131, 137], the reconstructed fiber should have the following 3 properties: first, the tangential direction of the reconstructed fiber should be same as the direction of water diffusion which is the vector field given by DWI; second, the reconstructed fiber should have no sharp turns that tend to minimize its length; and third, fibers do not pass through gray matter. To achieve the first two properties, Parker et al. [131] define the propagation speed as Equation 54 or 55 in 5.1.2.2. Formula 54 ensures that front evolution will occur most rapidly if both r1 and '1r are close to being colinear with each other. The incorporation of the term 'r 2F in Equation 55 PAGE 152 137 provides a memory of previous iterations of the front position. This can be interpreted as embedding connectivity information already established between the start point and in the value of Those definitions are very reasonable. But we notice, first, only one tensor direction 2F r1 is considered in formula (54) and (55). In case there are more than one fiber in one voxel, this approach has no solution. Therefore, in our algorithm all directions from HARD MRI data are considered. Second, we notice that formula (54) and (55) use the worst case speed. People may argue that it should be better if we consider all those terms together instead of simply choosing the smallest number among them. To address this objection, we propose to use a polynomial to describe all terms together. Given a point r and a tensor direction r1 at that point, we want to find the next point where the fiber should go. Since each neighbor of r may be a candidate, for each neighbor i r we construct local fiber connecting r and the neighbor by a polynomial and force it to pass through those two points with the tangential direction r1 at r and ir1 at i r then the local fiber follows the first property. To follow the second property, we calculate the length of those curves defined by those polynomials and find the shortest one. And then we choose the corresponding i r as the next point of r Because the length of the curve is decided by curvature and torso, the shortest curve has most likely no sharp turn and hence follows the second property of the reconstructed fiber. To smoothly connect those piecewise polynomials, we need to avoid sharp turns at the connection nodes, hence the angle between the direction in previous point of r and the direction of the next point of r is calculated, and a threshold is set to avoid sharp turn. The third property of the reconstructed fiber is used as a termination condition in our algorithm. Our algorithm is based on the idea of fast marching, and the termination PAGE 153 138 condition of fast marching method is when the set of neighborhoods becomes empty [74]. We delete an entry from the set of neighborhoods when the next propagating point of this entry is isotropic or when a sharp turn happens. To deal with noise problem and interpolation before fiber reconstruction, we proposed a method in [134]. Please refer to this paper for details. To display the tracking result, 1.5 T MR image is used as background. 5.2.2 Data Description There are four sets of data used as input for the proposed method. First, the threedimensional diffusion direction field generated with DWI data. In [134], we reported the details. For the tracking algorithm, we assume that the input data has low noise, high resolution, and more than 1 directions in some voxels. Second, the segmentation result of gray matter was generated with 1.5 T MRI data. The segmentation result is used as termination result in our algorithm. Third, seeds, the start region of the fibers, were given by experts; fourth, the targets (i.e. the interested regions the fibers may pass through) were given by experts. 5.2.2 Front Evolution and Path Record The first step toward determining paths of connection to a given start, or seed, involves the growth of a volume from this point. This is achieved using a modified version of the fast marching algorithm, a rapid implementation of boundaryvalue level sets methods [138]. We are primarily interested in the behavior of the surface of the volume (the front). The speed, at which the front propagates from the start point(s), is linked to the information contained in the field. Each iteration of the front position evolution involves the determination of the speed at candidates, which are the neighbors of the front. Let PAGE 154 139 r : the current point which is on the front r : the previous point of r on the fiber '' r : the next possible point on the fiber neighboring r ird, : the i th fiber direction at r krd,' : the k th fiber direction at' r which is used to connect r and r jirrC,,'', : the cubic polynomial that pass through r and '' r the tangential direction of is ji,,'' rrC, ird, at r and jrd,'' at '' r Then the speed from jirrS,,'', r to '' r with ird, and jrd,'' direction is defined as ''',,,',,'',,,'',rrLengthikrrCLengthjirrCLengthjirrS (58) and the speed from '',rrS r to '' r is defined as constjirrSrrSji,,'',min1'',, (59) Cubic polynomial C is used here to approximate the local fiber path, so the first property of the reconstructed fiber is satisfied. The definition of Equation 51 is a measure of the length of the local fiber. The definition of Equation 52 is the inverse of the minimum length of local fibers which connects jirr,,'', r and '' r So the second property of the reconstructed fiber is satisfied.is a constant number used to exaggerate the difference of speeds. For 3D, there are 26 neighbors for each voxel. The neighbors which have already in the same fiber or which are not gray matter are not considered as candidates. Here the third property of the reconstructed fiber is satisfied. Notice the speed const PAGE 155 140 is defined for each direction of r the missing branch problem of [131] mentioned in previous section can be avoided. In each iteration, we calculate the speeds for all neighbors of the current point r and then find the greatest speed among all neighbors. Notice all neighbors are considered, not only the neighbors of r but also neighbors of each front point]; hence the branching problem and noise accumulation problem in linear propagation method are solved. If the greatest speed is greater than a threshold, then add the neighbor ' r which has the greatest speed, into the fiber. Save r and '' r and the directions used to generate the greatest speed for the next iteration and tracking back in the future. r is called the closest point of '' r So r is the closest point of r We save all the speed information for the neighbor speed contest in the future. If the greatest speed is less than the threshold, then change the current point r to be another possible point on the front and mark r as an impossible point. We continue this process, until no possible points available or the fiber arrives at a point in the target region. We repeat the front evolution process for each seed. 5.2.4. Tracking Back After the process of front evolution, the fiber starts from one seed may have many branches, see Figure 51E. We need find the shortest path to connect the target and seed. Based on the closest point idea (or greedy algorithm), we start from the point r that the fiber touches in the target. We then read the saved record of the fiber to find the closest point r of r save it, and then find the closest point of r and save it too. We repeat this process until we track the fiber back to the seed. This is the shortest path connecting the seed and the target point. Figure 54 shows the flow chart of the proposed algorithm. PAGE 156 141 Figure 54. Flow chart of fiber tracking by fast marching on grids 5.3 Experiments and Results The proposed tracking method was applied to human MRI data. Given one region of seeds and two regions of targets, fibers starting from the seeds will touch one of those two regions or touch none of them. The purpose of this experiment is to find the ratio of touching and the path of those fibers. PAGE 157 142 DWI brain data were acquired using a GE 3.0tesla scanner using a singleshot spin echo EPI sequence. The scanning parameters for the DWI acquisition are: repetition time(TR) = 1000ms, echo time(TE) = 85 ms, the field of view (FOV) = 220 mm 220 mm. 24 axial sections covering the entire brain with the slice thickness = 3.8mm and the intersection gap = 1.2mm. Because it is required that the fiber has to pass through grids, this algorithm is reasonable when high resolution DWI is available. In general case, it is difficult to get the high resolution DWI, hence the DWI data for this experiment is interpolated. The interpolated and denoised tensor direction field is generated according to [134]. Dr Yijun Liu provided the seeds and target regions. There are 149 start points. The proposed tracking method was applied. The total tracking time is 240 seconds. Figure 55 shows the data and results. Out of 149 seeds, 92 touched target 1, 11 touched target 2, with a ratio 8.4:1. Figure 55A is the seeds region, the background is MR Image of that slice. Figure 55B is the mask for seeds slice, fiber can not go out of the mask because there is no fiber in white matter, this mask is used as a termination condition. Figure 55C is the map of target regions. The yellow region is target 1. The pink region is target 2. The background is the MR Image of this slice. Figure 55D is the mask on targets slice. Figure 55E is one example of the reconstructed fiber. We can see that the reconstructed fiber has several branches. It shows that the proposed method solves the branching problem. Figure 55F shows the correspondence between seeds and targets. The fibers start from yellow region arrive target 1, the fibers start from red region arrive target 2. The fibers that start from the blue region miss both targets. Figure 55G and H show all reconstructed fibers in 3D anatomic images. The lower slice is the seed slice. PAGE 158 143 The upper slice is the fa map of targets slice. The two targets are shown by different colors in the map. Fibers arrive different targets are shown by different colors. G and H are different views for the same image. af Figure 55. Fiber tracking results. A) the seeds region, B) the mask for seeds slice, C)the map of target regions, D)the mask on targets slice, E) one example of reconstructed fiber, F) the correspondence between seeds and targets, G) and H) all reconstructed fibers in 3D anatomic images 5.4 Conclusion A novel tracking algorithm is proposed. To solve the problems in existing front evolution methods, there are three major differences between this work and other front evolution methods. First, local polynomial fitting is used to define evolution speed. Second, the closest point method [31] is applied to find the most reasonable path to connect two points. Third, fast marching method is modified to deal with multidiffusiondirections in some voxels. The proposed method solves multidiffusiondirection, branches and tracking back problems existing in fiber tracking problems. There is one parameter (the threshold for PAGE 159 144 speed) to choose, hence it is easy to use. The method is fast: it only spent several minutes in all of our experiments. Experiments on simulated and human brain data show the effectiveness of this method. The drawback of this algorithm is that it requires that the fiber to pass through grids which is not reasonable when the image resolution is low. A more general algorithm will be shown in a separate work. PAGE 160 LIST OF REFERENCES [1] M. Kass, A. Witkin, and D. Terzopoulos, "Snakes: Active contour models," International Journal of Computer Vision, vol. 1, pp. 321331, 1987. [2] P. Perona and J. Malik, "Scalespace and edge detection using anisotropic diusion," IEEE Trans. Pattern Anal. and Mach. Intell., vol. 12, pp. 629639, 1990. [3] M. Bertalmio, G. Sapiro, V. Caselles, and C. Ballester, "Image inpainting," presented at SIGGRAPH Computer Graphics Proceedings, 2000. [4] A. J. Baxter, Y. F. Mangnall, E. J. Loj, B. H. Brown, D. C. Barber, A. G. Johnson, and N. W. Read, "Evaluation of applied potential tomography as a new noninvasive gastric secretion test," Gut, vol. 30, pp. 17301735, 1988. [5] S. J. Watson, R. H. Smallwood, B. H. Brown, P. Cherian, and K. D. Bardhan, "Determination of the Relationship Between the pH and Conductivity of Gastric Juice.," Physiol. Meas., vol. 17, pp. 21 27, 1996. [6] F. A. Duck, Physical Properties of Tissue A Comprehensive Reference Book. London: Academic Press, 1990. [7] P. H. Moller, K. G. Tranberg, B. Blad, P. H. Henriksson, L. Lindberg, L. Weber, and B. R. R. Persson, "EIT for measurement of temperature distribution in laser thermotherapy.," in Clinical and Physiological Applications of Electrical Impedance Tomography., vol. 24, D. Holder, Ed. London: UCL Press, 1993, pp. 217226. [8] F. J. McArdle, B. H. Brown, and A. Angel, "Imaging cardiosynchronous impedance changes in the adult head.," in Clinical and Physiological Applications of Electrical Impedance Tomography, vol. 20, D. Holeder, Ed. London: UCL Press, 1993, pp. 177184. [9] D. S. Holder, "Physiological constraints to imaging brain function using scalp electrodes," in Clinical and Physiological Applications of Electrical Impedance Tomography, vol. 21, D. Holder, Ed. London: UCL Press, 1993, pp. 185200. [10] N. D. Harris, A. J. Suggett, D. C. Barber, and B. H. Brown, "Applications of applied potential tomography (APT) in respiratory medicine," Clin. Phys. Physiol. Meas., vol. 8, pp. 155165, 1987. 145 PAGE 161 146 [11] B. Blad and B. Baldetorp, "Impedance spectra of tumor tissue in comparison with normal tissue: A possible clinical application for electrical impedance tomography," Physiological Measurement, vol. 17 Suppl 4A, pp. 105115, 1996. [12] N. Chauveau, "Ex vivo discrimination between normal and pathological tissues in human breast surgical biopsies using bioimpedance spectroscopy," Annals of the New York Academy of Sciences, vol. 873, pp. 4250, 1999. [13] J. C. Buell, "A practical, costeffective, noninvasive system for cardiac output and hemodynamic analysis," American Heart Journal, vol. 116, pp. 657664, 1988. [14] D. Weinstein, L. Zhukov, and G. Potts, "Localization of Multiple Deep Epileptic Sources in a Realistic Head Model via Independent Component Analysis," School of Computing ,University of Utah, Salt Lake City UUCS2000004, 2000. [15] D. Sodickson and W. Manning, "Simultaneous acquisition of spatial harmonics (SMASH): Ultrafast imaging with radiofrequency coil arrays," Magn Reson Med, vol. 38, pp. 591 603, 1997. [16] K. P. Pruessmann, M. Weiger, M. B. Scheidegger, and P. Boesiger, "SENSE: Sensitivity Encoding for Fast MRI," Magn Reson Med, vol. 42, pp. 952 962, 1999. [17] D. LeBihan, "Diffusion, perfusion and functional magnetic resonance imaging," J Mal Vasc, vol. 20, pp. 20314, 1995. [18] M. E. Moseley, Y. Cohen, J. Mintorovitch, L. Chileuitt, J. Suruda, D. Norman, and P. Weinstein., "Evidence of anisotropic selfdiffusion in cat brain," presented at 8th Annual Meeting of SMRM, Amsterdam, 1989. [19] M. E. Moseley, J. Kucharczyk, H. S. Asgari, and D. Norman, "Anisotropy in diffusion weighted MRI," Magn. Reson. Med, vol. 19, pp. 321326, 1991. [20] P. J. Basser and C. Pierpaoli, "Microstructual and physiological features of tissues elucidated by quantitative diffusion tensor MRI," Magn. Reson. Med, vol. 111(B), pp. 209219, 1996. [21] D. Mumford and J. Shah, "Optimal approximation by piecewise smooth functions and associated variational problems," Communications in Pure and Applied Mathematics, vol. 42, pp. 557685, 1989. [22] L. Rudin, S. J. Osher, and E. Fatemi, "Nonlinear total variation based noise removal algorithms," Physica D, vol. 60, pp. 259268, 1992. [23] R. Acar and C. R. Vogel, "Analysis of total variation penalty methods for illposed problems," Inverse Problems, vol. 10, pp. 12171229, 1994. PAGE 162 147 [24] D. Strong, "Adaptive total variation minimization image restoration," UCLA, 1997. [25] T. Chan and D. Strong, "Relation of Regularization Parameter and Scale in Total Variation Based Image Denosing," Dept. of Math, UCLA, Los Angeles CAM 967, 1996. [26] Y. Chen and W.Wunderli, "Adaptive total variation for image restoration in BV space," J. Math. Anal. Appl, vol. 272, pp. 117137, 2002. [27] Y. Chen, S. Levine, M. Rao, and J. Stanich, "Functions of p(x)Growth in Image Restoration," submitted to SIAM J. Appl Math, 2003. [28] S. Osher and J. A. Sethian, "Fronts propagating with curvature dependent speed: algorithms based on HamiltonJacobi formulations," J. Comput. Phys., vol. 79, pp. 1249, 1988. [29] H.K. Zhao, S. Osher, and R. Fedkiw, Fast surface reconstruction using the level set method," presented at IEEE Workshop on Variational and Level Set Methods (VLSM'01), Vancouver, Canada, 2001. [30] S. Chen, B. Merriman, M. Kang, R. E. Caflisch, C. Ratsch, L. Cheng, M. Gyure, R. P. Fedkiw, and S. Osher, "Level set method for thin film epitaxial growth," J. Comput. Phys., vol. 167, pp. 475500, 2001. [31] Y. R. Tsai, "Rapid and Accurate Computation of the Distance Function Using Grids," UCLA CAM, Los Angeles 0036, 2000. [32] B. Pompe, P. Blidh, D. Hoyer, and M. Eiselt, "Using mutual information to measure coupling in the cardiorespiratory system," IEEE Eng Med Biol Mag, vol. 17, pp. 3239, 1998. [33] A. Collignon, F. Maes, D.Vandermeulen, P. Suetens, and G. Marchal, "Automated multimodality image registration using information theory," presented at Processing in medical Images, Brest, France, 1995. [34] G. Hermosillo, C. ChefdHotel, and O. Faugeras, "A variational methods for multimodal image matching," Intl. J. Comp. Visi, vol. 50, pp. 329345, 2002. [35] C. Ballester, M. Bertalmio, V. Caselles, G. Sapiro, and J. Verdera, "Fillingin by joint interpolation of vector fields and grey levels," IEEE Trans. Image Process, vol. 10, pp. 12001211, 2001. [36] T. F. Chan and J. Shen, "Mathematical models for local nontexture inpaintings," SIAM J. Appl Math, vol. 62, pp. 10191043, 2001. [37] T. F. Chan, S. H. Kang, and J. Shen, "Euler's elastica and curvature based inpaintings," SIAM J. Appl. Math, vol. 63, pp. 564592, 2002. PAGE 163 148 [38] S. Esedoglu and J. H. Shen, "Digital inpainting based on the MumfordShahEuler image model," European J Appl Math, vol. 13, pp. 353370, 2002. [39] T. F. Chan and L. A. Vese, "Active contours without edges," IEEE Trans. Med. Imag., vol. 10, pp. 266277, 2001. [40] P. Blomgren, T. F. Chan, P. Mulet, and C.K.Wong, "Total Variation Image Restoration: Numerical Methods and Extensions," presented at IEEE International Conference On Image Processing, 1997. [41] L. A. Vese and T. F. Chan, "A multiphase level set framework for image segmentation using the Mumford and Shah model," UCLA CAM, Los Angles 0125, 2001. [42] M. Sussman, P. Smereka, and S. Osher, "A level set approach for computing solutions to incompressible twophase flow," J. Comput. Phys., vol. 5, pp. 146159, 1994. [43] B. Song and T. F. Chan, "A Fast Algorithm for Level Set Based Optimization," SIAM J. on Scientific Computing, 2003. [44] V. Caselles, F. Catte, T. Coll, and F. Dibos, "A geometric model for active contours in image processing," Numerische Mathematik, vol. 166, pp. 131, 1993. [45] R. Malladi, J. Sethian, and V. B., "Shape modeling with front propagation:A level set approach," IEEE Trans Pattern Anal machine Intell, vol. 17, pp. 158175, 1995. [46] V. Caselles, R. Kimmel, and G. Sapiro, "Geodesic active contours," International Journal of Computer Vision, vol. 22, pp. 6179, 1997. [47] S. Kichenassamy, A. Kumar, P. Olver, A. Tannenbaum, and A. Yezzi., "Gradient flows and geometric active contour models," presented at IEEE ICCV, Boston, 1995. [48] A. Yezzi, S. Kichenassamy, A. Kumar, P. J. Olver, and A. Tannenbaum, "A geometric snake model for segmentation of medical imagery," IEEE Trans. Med. Imag., vol. 16, pp. 199209, 1997. [49] S. C. Zhu and A. Yuille, "Region competition:Unifying snakes, region growing, and Bayes/MDL for multiband image segmentation," IEEE PAMI, vol. 18, pp. 884900, 1996. [50] T. F. Chan and L. Vese, "A level set algorithm for minimizing the MumfordShah functional in image processing," presented at IEEE Computer Society Proceedings of the 1st IEEE Workshop on Variational and Level Set Methods in Computer Vision, Vancouver, Canada, 2001. PAGE 164 149 [51] D. Cremers, C. Schnrr, and J. Weickert, "Diffusionsnakes combining statistical shape knowledge and image information in a variational framework," presented at IEEE Intl. Workshop on Variational and Levelset Methods, Vancouver, 2001. [52] A. Tsai, A. Yezzi, W. Wells, C. Tempany, D. Tucker, A. Fan, E. Grimson, and A. Willsky, "Modelbased curve evolution technique for image segmentation," presented at CVPR, Hawaii, 2001. [53] T. Cootes, C. Beeston, G. Edwards, and C. Taylor, "Unified framework for atlas matching using active appearance models," presented at Int'l Conf. Inf. Proc. in Med. Imaging, 1999. [54] Y. Wang and L. H. Staib, "Boundary Finding with Correspondence Using Statistical Shape Models," presented at CVPR, Santa Barbara, CA, 1998. [55] T. Cootes, C. Taylor, D. Cooper, and J. Graham, "Active shape modeltheir training and application," Computer Vision and Image Understanding, vol. 61, pp. 3859, 1995. [56] T. Cootes and C. Taylor, "Mixture model for representing shape variation," Image and Vision Computing, vol. 17, pp. 567574, 1999. [57] L. Staib and J. Duncan, "Boundary finding with parametrically deformable contour methods," IEEE TransPatt Analysis and Mach Intell, vol. 14, pp. 10611075, 1992. [58] A. Yuille, P. W. Hallinan, and D. S. Cohen, "Feature extraction from faces using deformable templates," Int J Computer Vision, vol. 8, pp. 99111, 1992. [59] H. D. Tagare, "Deformable 2D Template Matching Using implementation of the MumfordShah functional for image segmentation, denoising, interpolation, and magnification," IEEE Trans. on Image Proc., vol. 10, pp. 11691186, 2001. [60] M. E. Leventon, E. Grimson, and O. Faugeras, "Statistical shape influence in geodesic active contours," presented at IEEE Conf. CVPR (2000), 2000. [61] Y. Chen, S. Thiruvenkadam, F.Huang, H. Tagare, D. Wilson, and A. Geiser, "On the Incorporation of shape priors into geometric active contours," presented at IEEE Workshop in Variational and Level Set Methods in Computer Vision, Vancouver, Canada, 2001. [62] Y. Chen, H. Tagare, S. R. Thiruvenkadam, F.Huang, D. Wilson, K. S. Gopinath, R. W. Briggs, and A. Geiser, "Using prior shapes in geometric active contours in a variational framework," International Journal of Computer Vision, vol. 50, pp. 315328, 2002. PAGE 165 150 [63] N. Paragios, M. Rousson, and V. Ramesh, "Matching Distance Functions: A ShapetoArea Variational Approach for GlobaltoLocal Registration," presented at ECCV, Copenhagen, Denmark, 2002. [64] N. Paragios and R. Deriche, "Geodesic Active Contours for Supervised Texture Segmentation," presented at CVPR, Cofu, Greece, 1999. [65] D. G. Kendall, "Shapemanifolds, Procrustean metrics, and complex projective spaces," Bulletin London Mathematics Society, vol. 16, pp. 81121, 1984. [66] C. G. Small, The Statistical Theory of Shape Springer Series in Statistics. New York: SpringerVerlag, 1996. [67] E. A. Geiser, D. C. Wilson, D. Wang, D. A. Conetta, J. D. Murphy, and A. D. Hutson, "Autonomous epicardial and endocardial boundary detection inechocardiographic shortaxis images," Journal of the American Society of Echocardiography, vol. 11, pp. 338348, 1998. [68] P. Thevenaz and M. Unser, "Optimization of mutual information for multiresolution image registration," IEEE Trans. on Image Proc., vol. 9, pp. 20832099, 2000. [69] P. Viola and W. M. Wells, "Alignement by maximization of mutual information," IJCV, vol. 24, pp. 137154, 1997. [70] A. Tsai, J. A. Yezzi, and A. S. Willsky, "Curve evolution implementation of the MumfordShah functional for image segmentation, denoising, interpolation, and magnificati.," IEEE Trans. on Image Proc., vol. 10, pp. 11691186, 2001. [71] N. Paragios, "Geodesic Active Regions and Level Set Methods: Contributions and Applications in Artificial Vision," in School of Computer Engineering. Nice, France: University of Nice/Sophia Antipolis, 2000. [72] N. Paragios, "A varitional approach for segmentation of the left ventricle in MR cardiac images.," presented at 1st IEEE Workshop on Variational and Level Set Methods in Computer Vision, Vancouver, B C, 2001. [73] T. Cootes, C. Taylor, D. Cooper, and J. Graham, "Active shape model their training and application," Computer Vision and Image Understanding, vol. 61, pp. 3859, 1995. [74] J. A. Sethian, "Fast marching methods," SIAM Review, vol. 41, pp. 199235, 1999. [75] B. C. Vemuri, Y. Chen, and Z. Wang, "Registration Assisted Image Smoothing and Segmentation," presented at ECCV, Copenhagen, Denmark, 2002. PAGE 166 151 [76] F. Richard and L. Cohen, "A New Image Registration Technique with Free Boundary Constraints: Application to Mammography," presented at ECCV, Copenhagen, Denmark, 2002. [77] S. Soatto and A. J. Yezzi, "DEFORMOTION: Deforming Motion, Shape Average and the Joint Registration and Segmentation of Images," presented at ECCV, Copenhagen, Denmark, 2002. [78] M. E. Leventon, O. Faugeras, E. Grimson, and W. Wells, "Level set based segmentation with intensity and curvature priors," presented at Mathematical Methods in Biomedical Image Analysis, 2000. [79] C. Ballester, M. Bertalmio, V. Caselles, G. Sapiro, and J. Verdera, "Fillingin by joint interpolation of vector fields and grey levels," IEEE Trans. on Image Proc., vol. 10, pp. 12001211, 2001. [80] T. F. Chan, J. Shen, and L. Vese, "Variational PDE models in image processing," Amer. Math. Soc. Notice, vol. 50, pp. 1426, 2003. [81] T. F. C. a. J. Shen, "Nontexture inpainting by curvature driven diffusions," J. Visual Comm. Image Rep, vol. 12, pp. 436449, 2001. [82] V. Caselles, J. M. Morel, and C. Sbert, "An axiomatic approach to image interpolation," IEEE Trans. Med. Imag., vol. 7, pp. 376386, 1998. [83] D. Sodickson, Spatial Encoding Using Multiple RF Coils: SMASH Imaging and Parallel MRI: John Wiley & Sons Ltd, 2000. [84] K. P. Pruessmann, M. Weiger, M. B. Scheidegger, and P. Boesiger, "Coil sensitivity encoding for fast MRI," presented at 6th Intl. Soc. Mag. Reson. Med, Sydney, Australia, 1998. [85] M. Angel, G. Ballester, Y. Machida, Y. Kassai, Y. Hamamura, and H. Sugimoto, "Robust Estimation of Coil Sensitivities for RF Subencoding Acquisition Techniques," presented at 9th Intl. Soc. Mag. Reson. Med, 2001. [86] F. H. Lin, Y. J. Chen, J. W. Belliveau, and L. L.Wald, "Estimation of coil sensitivity map and correction of surface coil magnetic resonance images using wavelet decomposition," presented at 9th Intl. Soc. Mag. Reson. Med, Brighton,UK, 2001. [87] M. Cohen, R. DuBois, and M. Zeineh, "Rapid and effective correction of RF inhomogeneity for high field magnetic resonance imaging," Hum Brain Mapp, vol. 10, pp. 204211, 2000. [88] P. B. Roemer, W. A. Edelstein, C. E. Hayes, S. P. Souza, and O. M. Mueller, "The NMR phased array," Magn Reson Med, vol. 16, pp. 192225, 1990. PAGE 167 152 [89] J. W. Murakami, C. E. Hayes, and E. Weinberger, "Intensity correction of phased array surface coil images," Magn Reson Med, vol. 35, pp. 585590, 1996. [90] K. P. Pruessmann, M. Weiger, M. B. Scheidegger, and P. Boesiger, "Coil sensitivity maps for sensitivity encoding and intensity correction," presented at 6th Intl. Soc. Mag. Reson. Med, Sydney, 1998. [91] S. O. L. Rudin, and E. Fatemi., "Nonlinear total variation based noise removal algorithms.," Physica D, vol. 60, pp. 259268, 1992. [92] G. Aubert and L. Vese, "A variational method in image recovery," SIAM J. Numer. Anal, vol. 34, pp. 19481979, 1997. [93] C. A. McKenzie, E. N. Yeh, M. A. Ohliger, M. D. Price, and D. K. Sodickson, "SelfCalibrating Parallel Imaging With Automatic Coil Sensitivity Extraction," Magn Reson Med, vol. 47, pp. 529538, 2002. [94] V. HeelsBergen, R. Teunis, Fuderer, and Miha, "Magnetic resonance imaging method and apparatus." USA: US Philips Corporation, 1996. [95] K. Paulson, W. Breckon, and M. Pidcock, "Optimal measurements in electrical impedance tomography," IEEE EMBS, vol. 14, pp. 1730, 1992. [96] M. Vauhkonen, D. V. asz, P. A. Karjalainen, E. Somersalo, and J. P. Kaipio, "Tikhonov Regularization and Prior Information in Electrical Impedance Tomography," IEEE Trans. Med. Imag., vol. 17, pp. 285293, 1998. [97] G. T. Herman and L. Meyer, "Algebraic reconstruction techniques can be made computationally efficient," IEEE Trans. Med. Imag., vol. 12, pp. 600609, 1993. [98] A. M. Cormack, "Representation of a function by its line integrals, with some radiological applications I," J. Appl. Physics, vol. 34, pp. 27222727, 1963. [99] A. M. Cormack, "Representation of a function by its line integrals, with some radiological applications II," J. Appl. Physics, vol. 35, pp. 195207, 1964. [100] S. R. Deans, The Radon Transform and some of its Applications. New York: Wiley, 1983. [101] J. Qi and R. Leahy, "Resolution and noise properties of MAP reconstructions in fully 3D PET," IEEE Trans. Med. Imag., vol. 19, pp. 493506, 2000. [102] T. Hebert and R. Leahy, "A generalized EM algorithm for 3D Bayesian reconstruction from Poisson data using Gibbs distributions," IEEE Trans. Med. Imag., vol. 8, pp. 194202, 1989. [103] T. K. Moon, "The ExpectationMaximization algorithm," IEEE Trans Signal Processing, vol. 13, pp. 4760, 1996. PAGE 168 153 [104] J. A. Bilmes, "A gentle tutorial of the EM algorithm and its application to parameter estimation for Gaussian mixture and hidden Markov model," University of California, Berkeley Technical Report ICSITR97021, 1997. [105] R. Leahy and J. Qi, "Statistical approaches in quantitative positron emission tomography," Statistics and Computing, vol. 10, pp. 147165, 2000. [106] M. Cheney and D. Isaacson, "A layerstripping algorithm for impedance imaging," IEEE EMBS, vol. 13, pp. 34, 1991. [107] J. Sylvester, "A convergent layer stripping algorithm for radially symmetric impedance tomography problem," Comm. Partial Differential Equations, vol. 17, pp. 19551994, 1992. [108] E. Somersalo, M. Cheney, D. Isaacson, and E. L. Isaacson, "Layerstripping: A direct numerical method for impedance imaging," Inverse Problems, vol. 7, pp. 899926, 1991. [109] E. L. Miller, M. Kilmer, and C. Rappaport, "A new shapebased method for object localization and characterization from scattered field data," IEEE Trans. Geoscience and Remote Sensing, vol. 38, pp. 16821696, 2000. [110] H. Feng, W. C. Karl, and D. A. Castanon, "A curve evolution approach to objectbased tomographic reconstruction," IEEE Trans. Med. Imag., vol. 12, pp. 4457, 2003. [111] D. F. Yu and J. A. Fessler, "Edgepreserving tomographic construction with nonlocal regularization," IEEE Trans. Med. Imag., vol. 21, pp. 159173, 2002. [112] T. F. Chan and X. C. Tai, "Level set and total variation regularization for elliptic inverse problems with discontinuous coefficients," Journal of Computational Physics, vol. 193, pp. 4066, 2003. [113] A. O. Hero, R. Piramuthu, J. A. Fessler, and S. R. Titus, "Minimax emission computed tomography using high resolution anatomical side information and Bspline models," IEEE Transactions on Information Theory, vol. 45, pp. 920938, 1999. [114] T. F. Chan and S. Esedoglu, "A Multiscale Algorithm for MumfordShah Image Segmentation," UCLA, Los Angeles cam0377, Dec 17, 2003 2003. [115] F. Gibou and R. Fedkiw, "Fast hybrid kmeans level set algorithm for segmentation," Preprint, in review, pp. http://math.stanford.edu/~fgibou 2003. [116] A. Kaup and T. Aach, "Coding of segmented images using shape independent basis functions," IEEE Trans. Med. Imag., vol. 7, pp. 937947, 1998. PAGE 169 154 [117] M. Eden, M. Unser, and R. Leonardi, "Polynomial representation of pictures," Signal Process, vol. 10, pp. 385393, 1986. [118] G. Cusick, D. S. Holder, A. Birkett, and K. G. Boone, "A system for impedance imaging of epilepsy in ambulatory human subjects," Innov. Tech. Biol. Med., vol. 15(Suppl. 1), pp. 3439, 1994. [119] M. Radai, S. Abboud, and M. Rosenfeld, "Evaluation of impedance technique for detecting breast carcinoma using a 2D numerical model of the torso," Annals of the New York Academy of Sciences, vol. 873, pp. 360 369, 1999. [120] P. Basser, J. Mattiello, and D. LeBihan, "Estimation of the effective selfdiffusion tensor from the NMR spin echo," J. Magn. Reson. B, vol. 103, pp. 247, 1994. [121] S. Mori, B. Crain, V. Chacko, and P. C. M. van Zijl, "Three dimensional tracking of axonal projections in the brain by magnetic resonance imaging," Ann. Neurol., vol. 45, pp. 265269, 1999. [122] R. Xue, P. C. M. van Zijl, B. J. Crain, M. Solaiyappan, and S. Mori, "In vivo threedimensional reconstruction of rat brain axonal projections by diffusion tensor imaging," Magn. Reson. Med, vol. 42, pp. 11231127, 1999. [123] T. E. Conturo, N. F. Lori, T. S. Cull, E. Akbudak, A. Z. Snyder, J. S. Shimony, R. C. McKinstry, H. Burton, and M. E. Raichle, "Tracking neuronal fiber pathways in the living human brain.," presented at Natl Acad. Sci. USA, 1999. [124] C. Beaulieu, "The basis of anisotropic water diffusion in the nervous system," NMR Biomed, vol. 15, pp. 435., 2002. [125] S. Mori and P. C. M. v. Ziji, "Fiber tracking: principles and strategies A technical review," NMR in Biomedicine, vol. 15, pp. 468480, 2002. [126] D. S. Tuch, R. M. Weisskoff, J. W. Belliveau, and V. J. Wedeen, "High angular resolution diffusion imaging of the human brain," presented at 7th Intl. Soc. Mag. Reson. Med, Philadelphia, 1999. [127] D. S. Tuch, T. G. Reese, M. R. Wiegell, N. Makris, J. W. Belliveau, and V. J. Wedeen, "High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity," Magn Reson Med, vol. 48, pp. 577582, 2002. [128] P. J. Basser, S. Pajevic, C. Pierpaoli, J. Duda, and A. Aldroubi, "In vivo fiber tractography using DTMRI data," Magn Reson Med, vol. 44, pp. 625632, 2000. [129] N. F. Lori, E. Akbudak, J. Shimony, T. Cull, A. Snyder, R. Guillory, and T. Conturo, "Diffusion tensor fiber tracking: Methods, experimental results, and error simulations," presented at Diffusion NMR and MRI: From the Single Molecule to the Entire Human Brain,, 2001. PAGE 170 155 [130] N. F. Lori, E. Akbudak, J. S. Shimony, T. S. Cull, A. Z. Snyder, R. K. Guillory, and T. E. Conturo, "Diffusion tensor fiber tracking of human brain connectivity: aquisition methods, reliability analysis and biological results," NMR in Biomedicine, vol. 15, pp. 493515, 2002. [131] G. J. M. Parker, C. A. M. WheelerKingshott, and G. J. Barker, "Estimating distributed anatomical connectivity using fast marching methods and diffusion tensor imaging," IEEE Trans. Med. Imag., vol. 21, pp. 505512, 2002. [132] J. D. Tournier, F. Calamante, D. G. Gadian, and A. Connely, "A novel fibretracking technique: front evolution using a fibre orientation probability density function," presented at ISMRM 10th Scientific Meeting & Exhibition, Hawaii, 2002. [133] J. A. Sethian, "A fast marching level set method for monotonically advancing fronts," Proc. Nat. Acad. Sci., vol. 93, pp. 15911595, 1996. [134] Y. Chen, Y. J. Liu, G. J. He, W. Guo, Q. G. Zeng, F. Huang, and X. L. Yan, "Estimation, Smoothing, and Charaterization of Appranet Diffusion Coefficient Profiles from High Angular Resolution DWI," presented at 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, Washington ,DC, 2004. [135] G. Kindlmann, D. Weinstein, and D. Hart, "Strategies for direct volume rendering of diffusion tensor fields," IEEE Transactions on Visualization and Computer Graphics, vol. 6, pp. 124138, 2000. [136] C. Pierpaoli, A. Barnett, S. Pajevic, R. Chen, L. Penix, A. Virta, and P. Basser, "Water diffusion changes in Wallerian degeneration and their dependence on white matter architecture," NeuroImage, vol. 13, pp. 11741185, 2001. [137] D. S. Tuch, "Mapping cortical connectivity with diffusion MRI," presented at IEEE ISBI, Washington, D.C, 2002. [138] J.A.Sethian, "A fast marching level set method for monotonically advancing fronts," Proc. Nat. Acad. Sci., vol. 93, pp. 15911595, 1996. PAGE 171 BIOGRAPHICAL SKETCH I was born in Anhui province, China, on March 10, 1973. At the age of seventeen, I finished high school and entered the mathematics Department of Anhui University in the fall of 1989. I received my bachelors degree in science in July 1993 (as I turned 20 years of age.) In 1996, I obtained a masters degree of engineering in the field of computer aided geometric design (CAGD) from Hefei University of Technology. In the fall of 1996, I took a job at Beijing Forestry University as a mathematics lecturer, where I worked until the 1999. I joined the Ph.D. program of the Department of Mathematics at the University of Florida in August 1999. The next year, I began work with Dr. Chen on medical image processing. It became apparent to me that this field was the most exciting thing I had encountered in my whole life. Since the spring of 2001, I have worked in Diagnostic imaging, Invivo Corporation, Gainesville on a number of joint projects. 156 