<%BANNER%>

Neuronal fiber tracking in Dt-Mri

University of Florida Institutional Repository

PAGE 1

NEUR ONAL FIBER TRA CKING IN DT-MRI By TIM E. MCGRA W A THESIS PRESENTED TO THE GRADUA TE SCHOOL OF THE UNIVERSITY OF FLORID A IN P AR TIAL FULFILLMENT OF THE REQUIREMENTS F OR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORID A 2002

PAGE 2

Cop yrigh t 2002 b y Tim E. McGra w

PAGE 3

F or m y wife, Jo, and m y mother, P atti.

PAGE 4

A CKNO WLEDGMENTS I w ould lik e to express thanks to ev ery one who made this researc h p ossible. First, thanks go to Dr. Baba C. V em uri for serving as the c hairman of m y thesis committee, and b eing a v ailable for consultation and advice. Thanks go, also, to Dr. J org P eters for serving on m y committee and imparting the graphics kno wledge required for the visualization asp ect of this pro ject. Also, man y thanks go to the mem b ers of the Computer Vision, Graphics and Medical Imaging (CV GMI) group, esp ecially Zhizhou W ang, Jundong Liu and F ei W ang, for their assistance in preparing this w ork. I am grateful to the facult y mem b ers of the UF Mathematics Departmen t, Dr. Y unmei Chen, for serving on m y committee, and Dr. Murali Rao, for advisemen t and encouragemen t. I w ould also lik e to ac kno wledge the McKnigh t Brain Institute for pro viding data, analysis and v alidation of results. Thanks go to Dr. T om Mareci, Dr. Stev e Blac kband, Dr. P aul Reier, and Ev eren Ozarslan. This researc h w as funded in part b y the NIH gran t R O1-NS42075. iv

PAGE 5

T ABLE OF CONTENTS page A CKNO WLEDGMENTS . . . . . . . . . . . . . . iv LIST OF FIGURES . . . . . . . . . . . . . . . . vii KEY TO ABBREVIA TIONS . . . . . . . . . . . . . viii ABSTRA CT . . . . . . . . . . . . . . . . . . ix CHAPTERS 1 INTR ODUCTION . . . . . . . . . . . . . 1 1.1 Ov erview . . . . . . . . . . . . . . 1 1.2 Con tributions . . . . . . . . . . . . . 2 1.3 Outline of thesis . . . . . . . . . . . . 3 2 BA CK GR OUND O VER VIEW . . . . . . . . . . 5 2.1 DT-MRI Data Acquisition . . . . . . . . . 5 2.1.1 Ov erview of Diusion . . . . . . . . . 8 2.1.2 Ov erview of MR Imaging . . . . . . . . 9 2.1.3 DT-MRI Acquisition . . . . . . . . . 10 2.1.4 Stejsk al-T anner Equation . . . . . . . . 11 2.2 PDE Based Scalar-V alued Image Restoration . . . . 12 2.2.1 Linear Filtering . . . . . . . . . . 13 2.2.2 Nonlinear Filtering . . . . . . . . . . 14 2.2.2.1 P erona-Malik . . . . . . . . . 15 2.2.2.2 T ensor anisotropic diusion . . . . . 16 2.2.2.3 ALM . . . . . . . . . . . 17 2.2.3 V ariational F orm ulation . . . . . . . . 17 2.2.3.1 Mem brane spline . . . . . . . . 18 2.2.3.2 Thin-Plate spline . . . . . . . 19 2.2.4 T otal V ariation Image Restoration . . . . . 19 2.3 PDE Based V ector-V alued Image Restoration . . . . 21 2.3.1 V ector-V alued Diusion . . . . . . . . 22 2.3.2 Riemannian Metric Based Anisotropic Diusion . . 22 2.3.3 Beltrami Flo w . . . . . . . . . . . 23 2.3.4 Color T otal V ariation . . . . . . . . . 23 v

PAGE 6

3 DT-MRI IMA GE RESTORA TION . . . . . . . . . 25 3.1 Noise Mo del . . . . . . . . . . . . . 25 3.2 Restoration F orm ulation . . . . . . . . . . 26 4 NEUR ONAL FIBER TRA CKING . . . . . . . . . 28 4.1 Ov erview of Neuronal Fib er T rac king . . . . . . 28 4.2 F orm ulation . . . . . . . . . . . . . 30 5 NUMERICAL METHODS . . . . . . . . . . . 33 5.1 Data Denoising . . . . . . . . . . . . 33 5.1.1 Fixed-P oin t Lagged-Diusivit y . . . . . . 33 5.1.2 Discretized Equations . . . . . . . . . 34 5.2 Fib er Regularization . . . . . . . . . . . 35 5.2.1 Lagged-Diusivit y . . . . . . . . . . 36 5.2.2 Crank-Nic holson Metho d . . . . . . . . 36 6 NEUR ONAL FIBER VISUALIZA TION . . . . . . . 38 6.1 Rendered Ellipsoids . . . . . . . . . . . 38 6.2 Stream tub es . . . . . . . . . . . . . 38 6.3 Line In tegral Con v olution . . . . . . . . . . 38 6.4 P articles . . . . . . . . . . . . . . 40 7 EXPERIMENT AL RESUL TS . . . . . . . . . . 42 7.1 Data Denoising . . . . . . . . . . . . 42 7.2 Fib er T rac king . . . . . . . . . . . . . 43 7.2.1 LIC . . . . . . . . . . . . . . 43 7.2.2 Stream tub es . . . . . . . . . . . 44 7.2.3 P articles . . . . . . . . . . . . 44 8 CONCLUSIONS AND FUTURE W ORK . . . . . . . 48 8.1 Conclusions . . . . . . . . . . . . . 48 8.2 F uture W ork . . . . . . . . . . . . . 48 8.2.1 Noise Mo del F or Measured Data . . . . . . 48 8.2.2 Quan titativ e V alidation . . . . . . . . 49 8.2.3 Robust Regression . . . . . . . . . . 49 8.2.4 Non-T ensor Mo del of Diusion . . . . . . 49 8.2.5 Automatic Region-Of-In terest Extraction . . . 49 APPENDIX DERIV A TION OF EULER-LA GRANGE CONDITIONS . 50 REFERENCES . . . . . . . . . . . . . . . . . 54 BIOGRAPHICAL SKETCH . . . . . . . . . . . . . . 59 vi

PAGE 7

LIST OF FIGURES Figure page 2.1 Diusion Ellipsoid . . . . . . . . . . . . . . 8 2.2 TV( f 1) > TV( f 2) = TV( f 3) . . . . . . . . . . . 20 2.3 Noisy image (left) and restored (righ t) b y TV norm minimization. . 21 6.1 LIC visualization of syn thetic feld. . . . . . . . . . 39 7.1 F A image of coronal slice of ra w rat brain data. . . . . . . 42 7.2 F A results for smo othed data. . . . . . . . . . . . 42 7.3 LIC fb er tracts in coronal slice of smo othed rat brain data. . . . 43 7.4 LIC fb er tracts in axial slice of smo othed rat brain data. . . . 43 7.5 Comparison of ruoroscopic image (left) with extracted stream tub es (righ t). . . . . . . . . . . . . . . . . . 44 7.6 Axial view of stream tub es in corpus callosum. . . . . . . 45 7.7 Details of corpus callosum. . . . . . . . . . . . . 46 7.8 Sequence from Fib er visualization (2 seconds b et w een images). . . 47 vii

PAGE 8

KEY TO ABBREVIA TIONS CNS: cen tral nerv ous system DT: diusion tensor DTI: diusion tensor imaging F A: fractional anisotrop y LIC: line in tegral con v olution MRI: magnetic resonance imaging PDE: partial dieren tial equation RF: radio frequency R OI: region of in terest TV: total v ariation viii

PAGE 9

Abstract of Thesis Presen ted to the Graduate Sc ho ol of the Univ ersit y of Florida in P artial F ulllmen t of the Requiremen ts for the Degree of Master of Science NEUR ONAL FIBER TRA CKING IN DT-MRI By Tim E. McGra w Decem b er 2002 Chair: Baba C. V em uri Ma jor Departmen t: Computer and Information Sciences and Engineering Diusion tensor imaging (DTI) can pro vide the fundamen tal information required for viewing structural connectivit y Ho w ev er, robust and accurate acquisition and pro cessing algorithms are needed to accurately map the nerv e connectivit y In this thesis, w e presen t a no v el algorithm for extracting and visualizing the b er tracts in the CNS, sp ecically the spinal cord. The automatic b er tract mapping problem will b e solv ed in t w o phases, namely a data smo othing phase and a b er tract mapping phase. In the former, smo othing is ac hiev ed via a w eigh ted TV-norm minimization, whic h striv es to smo oth while retaining all relev an t detail. F or the b er tract mapping, a smo oth 3D v ector eld indicating the dominan t anisotropic direction at eac h spatial lo cation is computed from the smo othed data. Neuronal b ers are traced b y calculating the in tegral curv es of this v ector eld. ix

PAGE 10

Results are expressed using three mo des of visualization. Line in tegral con v olution (LIC) pro duces an orien ted texture whic h sho ws b er path w a ys in a planar slice of the data. A stream tub e map is generated to presen t a three-dimensional view of b er tracts. Additional information, suc h as degree of anisotrop y can b e enco ded in the tub e radius, or b y using color. A particle system form of visualization is also presen ted. This mo de of displa y allo ws for in teractiv e exploration of b er connectivit y with no additional prepro cessing. x

PAGE 11

CHAPTER 1 INTR ODUCTION 1.1 Ov erview The understanding of neurological function, and medical diagnosis of disease and injury require kno wledge of the structural organization of the brain. Magnetic Resonance Imaging pro vides a non-in v asiv e means of studying anatomical connectivit y This connectivit y in the form of axonal nerv e b er bundles, can only b e measured indirectly The presence of these b ers m ust b e inferred from the b eha vior of w ater molecules near the tissue. Diusion tensor magnetic resonance imaging (DT-MRI) is a metho d of measuring the rate of w ater diusion in biological structures. A diusion tensor at eac h lo cation on a regular lattice describ es a v olumetric a v erage of the directional properties of diusion within eac h v o xel. The observ ation that diusion is anisotropic in areas of white-matter b er bundles allo ws trac king of b ers through the lattice. Ho w ev er, acquisition noise corrupts the data measuremen ts. V oltage v ariations in the receiving coil of the MRI mac hine due to thermal noise are a ma jor source of signal degradation. The noise is accepted to b e zero-mean, and additiv e in nature. In tegrating a b er path o v er a noisy eld is an ill-p osed problem. The goals of the w ork presen ted in this thesis are to trace b ers through an ideal, unkno wn eld giv en only noisy observ ations, and to visualize the b ers in a useful w a y 1

PAGE 12

2 The b er tracts are estimated in a t w o stage pro cess, (a) smo othing of the ra w directional images acquired for v arying magnetic eld strengths and estimating the diusion tensor, D, eld o v er the 3D image lattice, (b) computing the dominan t eigen v ector eld from this regularized D and estimating regularized streamlines/in tegral curv es as the desired b er tracts. Computed b ers m ust b e visualized in a w a y that is useful for anatomical study or medical diagnosis. Existing v ector eld visualization tec hniques, suc h as streamlines/stream tub es can b e used to visualize b ers. These tec hniques can also b e mo died to con v ey more information ab out diusion b y incorp orating quan tities deriv ed from the tensor. Here, w e adapt line in tegral con v olution (LIC), stream tub es, and particles to the task of tensor eld visualization. W e will review the literature on the ra w data smo othing, streamline generation, and visualization tec hniques where appropriate. 1.2 Con tributions Although is w ell kno wn that the DT-MRI data requires denoising prior to b er trac king, previous w ork has concen trated on smo othing the v ector eld of dominan t diusion directions, or ev en the diusion tensor itself. Emplo ying these metho ds requires that the diusion tensor, and p ossibly eigensystems of the tensor b e computed from noisy data. Results of these calculations ma y b e meaningless due to prop erties of the noise, and the diusion tensor calculation. W e prop ose to smo oth the data b efore these calculation to pro vide for more ph ysically meaningful results.

PAGE 13

3 A new application of particle-based visualization is presen ted. This tec hnique has previously b een used in ruid mec hanics for visualization of v elo cit y elds. Here w e will adapt the tec hnique to tensor eld visualization. This tec hnique has no asso ciated prepro cessing time, and allo ws real-time in teractiv e visualization of the reconstructed data. Man y other tec hniques require extensiv e computing time for streamline in tegration prior to visualization. 1.3 Outline of thesis The rest of this thesis is organized as follo ws. In Chapter 2 w e will discuss the diusion pro cess and giv e a short in tro duction to the DT-MRI data acquisition pro cess. An o v erview of PDE based image denoising will follo w. Some linear and nonlinear tec hniques for scalar and v ector-v alued images will b e review ed. V ariational tec hniques for form ulating image denoising, and the TV norm will b e presen ted. Chapter 3 will describ e our v ector-v alued image restoration tec hnique, a w eigh ted TV norm minimization. The noise mo del used to form ulate this tec hnique is discussed, as w ell as the motiv ation for mo difying the TV norm usually used for v ector-v alued images. In Chapter 4 w e presen t the pro cess of trac king neuronal b er bundles through the restored DT-MRI data. Previous metho ds of b er trac king will b e discussed. A quan tication of diusion tensor anisotrop y will b e giv en. In this c hapter streamline regularization will b e presen ted as a w a y of enforcing a smo othness constrain t on the b ers.

PAGE 14

4 Chapter 5 will sho w the linearization and discretization of the PDEs in v olv ed in data denoising and b er regularization. The n umerical metho ds used to solv e linear equations will b e describ ed. In Chapter 6 w e will presen t visualization tec hniques for DT-MRI data. These include rendered ellipsoids, stream tub es, LIC and particles. In Chapter 7 exp erimen tal results will b e rep orted. The data denoising results for spinal cord and brain data sets will b e presen ted. Fib er trac king results using sev eral visualization tec hniques will b e presen ted. Conclusions and p ossible directions for future w ork will b e discussed in Chapter 8.

PAGE 15

CHAPTER 2 BA CK GR OUND O VER VIEW 2.1 DT-MRI Data Acquisition F undamen tal adv ances in understanding living biological systems require detailed kno wledge of structural and functional organization. This is particularly imp ortan t in the nerv ous system where anatomical connections determine the information path w a ys and ho w this information is pro cessed. Our curren t understanding of the nerv ous system is incomplete b ecause of a lac k of fundamen tal structural information [21 ] necessary to understand function. In addition, understanding fundamen tal structural relationships is essen tial to the dev elopmen t and application of therapies to treat pathological conditions (e.g. disease or injury). Ho w ev er, most imaging metho ds giv e only an anatomically isolated representation of living tissue b ecause the images do not con tain connectivit y information. Suc h information w ould allo w the iden tication and correlation of system elemen ts resp onding during function. F or example in brain trauma, the relationship b et w een anatom y and b eha vior will only b ecome apparen t when w e are able to discriminate the aeren t nerv e b er path w a ys transmitting the sensation from a stim ulus to the brain or the eeren t path w a ys transmitting impulses from the brain area con trolling b eha vior. 5

PAGE 16

6 Researc h during the last few decades has sho wn that the cen tral nerv ous system is able to adapt to c hallenges and reco v er some function. Ho w ev er, the structural basis for this adaptiv e abilit y is not w ell understo o d. F or the en tire cen tral nerv ous system, understanding and treating ev olving pathology suc h as brain trauma, dep ends on a detailed understanding of the anatomical connectivit y c hanges and ho w they relate to function. Recen tly MR measuremen ts ha v e b een dev elop ed to measure the tensor of diusion. This pro vides a complete c haracterization of the restricted motion of w ater through the tissue that can b e used to infer tissue structure and hence b er tracts. In a series of pap ers, Basser and colleagues [2, 3, 4, 5, 6 36 ] ha v e discussed in detail general metho ds of acquiring and pro cessing the complete apparen t-diusion-tensor of MR measured translational self-diusion. They sho w ed that directly measured diusion tensors could b e recast in a rotationally in v arian t form and reduced to parametric images that represen t the a v erage rate of diusion (tensor trace), diusion anisotrop y (relationship of eigen v alues), and ho w the diusion ellipsoid (eigen v alues and eigen v ectors) can b e related to the lab oratory reference frame. The parametric images [5, 35] of v olume ratio, fractional anisotrop y and lattice anisotrop y index represen t scalar measures of diusion that are indep enden t of the lab reference frame and sub ject orien tation. Therefore, these measures can b e used to c haracterize the tissue pathology e.g., isc hemia, indep enden t of the sp ecic frame of reference used to acquire the images. The dev elopmen t of diusion tensor acquisition, pro cessing, and analysis metho ds pro vides the framew ork for creating b er tract maps based on this complete

PAGE 17

7 diusion tensor analysis [19, 24 29 30]. This has b een used to pro duce b er tract maps in rat brains [30, 47 ] and to map b er tracts in the h uman brain [24] then the rst steps w ere tak en to relate this structural connectivit y to function [19]. The directional prop erties of diusion can b e c haracterized b y a diusion tensor, a 3 3 symmetric matrix of real v alues. In order to calculate the 6 indep enden t comp onen ts of the tensor, the sub ject is imaged in sev en dieren t directions with sev eral magnetic eld strengths. The relationship S = S 0 exp ( P ij b ij D ij ) allo ws the diusion tensor, D and the T-2 w eigh ted image S 0 to b e calculated giv en the samples S and eld gradien t strength, b Previous w ork has concen trated on smo othing the eld of eigen v ectors of D More recen t w ork [14 ] has form ulated regularization tec hniques for the tensor eld itself, ev en constraining the resulting tensors to b e p ositiv e-denite. Here w e will tak e the approac h of smo othing the observ ed v ector-v alued image S prior to calculating D This decision to smo oth in this manner will b e justied in Chapter 3. In summary the anisotrop y of w ater translational diusion can b e used to visualize structure in the brain and pro vides the basis for a new metho d of visualizing nerv e b er tracts. Initial results ha v e b een v ery encouraging and suggest that this approac h to b er mapping ma y b e applied to a wide range of studies in living sub jects. Ho w ev er, it is essen tial to optimize the acquisition and pro cessing algorithms for b er tract mapping and v alidate the results relativ e to kno wn measures of b er tracts.

PAGE 18

8 n rFigure 2.1: Diusion Ellipsoid 2.1.1 Ov erview of Diusion Random molecular motion (Bro wnian motion) can cause transp ort of matter within a system. Within a v olume of w ater, molecules freely diuse in all directions. The w ater abundan t in biological systems is also sub ject to suc h sto c hastic motion. The prop erties of the surrounding tissue can aect the magnitude of diusion, and the directional prop erties as w ell. Tissue can form a barrier to diusion, restricting molecular motion. Within an orien ted structure, suc h as a bundle of axonal b ers, diusion can b e highly anisotropic. The white matter of the brain and spinal cord is c haracterized b y man y suc h bundles. The directional prop erties of diusion can b e c haracterized b y a tensor. The diusion tensor, D is a symmetric, p ositiv e-denite 3 3 matrix. W e will mak e use of the eigen v alues and eigen v ectors of this tensor, sorting the eigen v alues ( 1 ; 2 ; 3 ) from largest to smallest, and lab elling the corresp onding unit eigenv ectors ( e 1 ; e 2 ; e 3 ). The eigen v alues represen t the magnitude of diusion in the direction of their corresp onding eigen v ector. F or isotropic diusion 1 = 2 = 3 A

PAGE 19

9 p opular represen tation for describing anisotropic diusion is the diusion ellipsoid. This ellipsoid is the image of the unit sphere under the transformation dened b y the tensor, D The eigen v ectors of D form an orthogonal basis, represen ting the orien tation of the ellipsoid. The length of eac h axis of the ellipsoid is the corresp onding eigen v ector. F or isotropic diusion, the diusion ellipsoid is a sphere. 2.1.2 Ov erview of MR Imaging In this section w e will presen t a brief o v erview of the MRI acquisition pro cess. A detailed treatmen t of the sub ject w as done b y Haac k e et al. [23]. The protons in the n uclei of atoms align their axis of spin with the direction of an applied magnetic eld. The magnetic eld also induces a w obble, kno wn as precession, in the spin of the protons. This frequency the resonan t frequency is prop ortional to the strength of the applied eld. F or protons the resonance frequency lies in the RF range. In the MRI mac hine, a eld B 0 is applied through out the imaging pro cess. The direction of this eld denes the axial direction of the image. Protons will absorb energy from an RF pulse of the resonance frequency and tip a w a y from the direction induced b y B 0 The amoun t of tip is prop ortional to the pulse duration. The RF pulse also causes the protons to precess in phase with eac h other. This pulse is called the B 1 eld. When the B 1 transmitter is turned o, the absorb ed energy at the resonan t frequency is re-emitted b y the protons. This o ccurs as the spins, tipp ed b y B 1 return to their previous B 0 alignmen t. The time constan t asso ciated with this

PAGE 20

10 exp onen tial pro cess is kno wn as the T 1 relaxation time. The protons precessions also dephase exp onen tially with time constan t T 2 The nal image con trast is inruenced b y strength, width and rep etition time of the RF pulses in the B 1 signal. By spatially v arying the in tensities of B 0 and B 1 p osition information is enco ded. F or instance, sp ecially designed magnets add a gradien t eld to B 0 This causes the proton resonance frequency to b e a function of axial p osition. The frequency of B 1 can then b e c hosen to tip protons within a c hosen slice. T o enco de x; y p osition within a slice, t w o additional additional gradien t elds are emplo y ed. The rst gradien t, G y is pulsed, causing a phase v ariation, just as in T 2 relaxation. The phase v ariation is a function of p osition in the y direction. A p erp endicular gradien t, G x is then applied, c hanging resonance frequencies in the x direction. A 2D F ourier transform reconstructs the image of eac h slice from the data in the spatial frequency domain. 2.1.3 DT-MRI Acquisition By carefully designing gradien t pulse sequences, the measured signal from protons in w ater molecules undergoing diusion can b e atten uated. The rst gradien t pulse induces a kno wn phase shift in proton precession. After some dela y a second gradien t pulse is applied, inducing the opp osite phase shift. Protons whic h ha v e not mo v ed b et w een the t w o gradien t pulses are returned to their previous phase. Protons b elonging to molecules whic h ha v e c hanged lo cation ha v e some net c hange in phase, c hanging their T 2 relaxation time. Con v en tional MRI images of

PAGE 21

11 white matter in the brain suggest a material of homogeneous comp osition. The brous nature of white matter is apparen t in DT-MRI ho w ev er. 2.1.4 Stejsk al-T anner Equation The in tensit y of the receiv ed signal, S at eac h v o xel dep ends on the prop erties of the diusion-enco ding gradien t, and the apparen t diusion tensor, D at that lo cation. The Stejsk al-T anner equation relates all of these quan tities ln ( S S 0 ) = 3 X i =1 3 X j =1 b i;j D i;j (2.1) In 2.1, S 0 is the in tensit y with no diusion-enco ding gradien t presen t, and b is a matrix c haracterizing the gradien t pulse sequence. The ph ysics b ehind 2.1 is b ey ond the scop e of this thesis. The motiv ated reader ma y refer to the w ork of Haac k e et al. [23] for more details. This imaging pro cess m ust b e p erformed with 7 noncoplanar gradien t directions in order to fully generate a diusion tensor image. Multiple samples, usually 3 or 4 are tak en for eac h gradien t direction. 26666664 ln S 1 ... ln S 28 37777775 = 26666664 1 b 1xx b 1y y b 1z z 2 b 1xy 2 b 1y z 2 b 1xz ... ... ... ... ... ... ... 1 b 28xx b 28y y b 28z z 2 b 28xy 2 b 28y z 2 b 28xz 37777775 266666666666666666666664 ln S 0 D xx D y y D z z D xy D y z D xz 377777777777777777777775 (2.2)

PAGE 22

12 The o v erconstrained linear system, 2.1.4 is solv ed for S 0 and the elemen ts of the symmetric tensor D b y a least squares linear regression. 2.2 PDE Based Scalar-V alued Image Restoration Image data smo othing or denoising is a fundamen tal problem in image pro cessing. Image denoising (noise remo v al) is a tec hnique that enhances images b y attempting to rev erse the eects of degradations o ccurring during acquisition or transmission. Image noise mak es it dicult to p erform other pro cessing tasks suc h as edge-detection, segmen tation, or in our case, b er trac king. F or this reason, denoising is the rst step in most image analyses. The goal of image denoising is to reco v er an unkno wn, ideal image giv en some observ ed image. In our case, w e wish to reco v er the smo oth S v alues from whic h the DT image is calculated. The most common degradation source is the noise from the image acquisition system and is commonly mo deled b y additiv e Gaussian random noise. In the follo wing, w e will briery review represen tativ e sc hemes, sp ecically partial dieren tial equation (PDE) based metho ds that lend themselv es to fast n umerical implemen tations. Image denoising can b e form ulated using v ariational principles whic h in turn require solutions to PDEs. Recen tly there has b een a rurry of activit y on the PDE-based smo othing sc hemes. F or scalar v alued image smo othing using nonlinear diusion lters with scalar diusivit y co ecien t, w e refer the reader to the follo wing articles and references therein [1, 11 12, 15, 28, 31 32 34 43 45]. Anisotropic diusion lters that use a tensorial diusivit y parameter

PAGE 23

13 w ere in tro duced in W eic k ert [45]. These lters can b e tailored to enhance image structures (edges, parallel lines, curv es etc.) that o ccur in preferred directions. 2.2.1 Linear Filtering Linear lters are a simple and ecien t means of remo ving noise from images. One suc h lter ma y b e implemen ted b y the pro cess of isotropic diusion. This diusion pro cess is go v erned b y the heat equation @ I ( x; t ) @ t = r 2 I ( x; t ) I ( x; 0) = I 0 ( x ) (2.3) In the same w a y that a heated plate will seek an equilibrium state of a smo oth temp erature gradien t, so will an image ev olv ed according to 2.3 smo oth out discon tin uities in in tensit y It can b e sho wn that isotropic diusion is equiv alen t to con v olving the image with a gaussian k ernel. T aking the F ourier transform of 2.3 with resp ect to x and dening F ( I ( x; t )) = U ( w ; t ). @ U ( w ; t ) @ t = w 2 U ( w ; t ) U ( w ; 0) = U 0 ( w ) (2.4) The solution to 2.4 is U ( w ; t ) = U 0 ( w ) e w 2 t (2.5)

PAGE 24

14 The solution to 2.3 is then the in v erse F ourier transform of 2.5 I ( x; t ) = I 0 1 p 2 t e x 2 2 2 t 2 t = 2 t (2.6) Clearly the linear diusion pro cess con tin ues un til a I ( x; t ) b ecomes a constan tv alued image. By 2.6 w e can consider the eect of 2.3 as t 1 to equiv alen t to con v olving the original image with gaussian k ernels of ev er-increasing v ariance. Since the gaussian k ernel is separable, this t yp e of lter is simple to implemen t for images of arbitrary dimension. Although the images pro duced b y this simple ltering tec hnique sho w a reduction in the high frequency noise, there is also the un w elcome eect of blurred edges and lost details. The linear lter still has applications in image resampling, and generating scale-space image represen tations, suc h as image p yramids. 2.2.2 Nonlinear Filtering The nonlinear lters describ ed in this c hapter w ere designed to o v ercome the in ter-region blurring eect of linear lters. Nonlinear lters cannot b e mo delled with gaussian con v olution. The median lter is a nonlinear lter with a simple implemen tation. A median lter replaces eac h in tensit y v alue in an image with the median of the neigh b oring v alues. The size of this neigh b orho o d determines the amoun t of smo othing.

PAGE 25

15 There are n umerous mo dels of nonlinear diusion for image smo othing. T o implemen t the lters presen ted in the rest of this section w e m ust n umerically solv e a PDE. 2.2.2.1 P erona-Malik An anisotropic diusion pro cess whic h inhibits blurring at edges w as form ulated b y P erona and Malik [34] b y mo difying 2.3. @ I ( x; t ) @ t = div ( c ( x ) r I ( x; t )) I ( x; 0) = I 0 ( x ) (2.7) When comparing 2.7 to 2.3 recall that r 2 I r r I div ( r I ). By using a diusion co ecien t whic h is a decreasing function of jr I ( x; t ) j suc h as c ( x ) = 1 1 + jr I ( x; t ) j (2.8) w e can inhibit diusion at lo cations of high image gradien t, presumed to b e edges. By adding a reactiv e term to 2.7 w e can actually enhance edges. The presence of I ( x; t ) in 2.8 mak es 2.7 nonlinear. The diusion co ecien t serv es only to slo w diusion at edges. The steadystate ( t 1 ) solution of 2.7 is still a constan t v alued image. A constrain t can b e added to imp ose the condition that the smo othed image b e \close" to the original image. By in tro ducing a reaction term to the diusion equation w e can imp ose a

PAGE 26

16 data delit y requiremen t @ I ( x; t ) @ t = div ( c ( x ) r I ( x; t )) + ( I 0 ( x ) I ( x; t )) I ( x; 0) = I 0 ( x ) (2.9) w e p enalize solutions that dier from the initial image. The parameter con trols the degree of smo othing in the nal image. The steady-state solution of 2.9 is non trivial, so w e no longer need a stopping time. 2.2.2.2 T ensor anisotropic diusion An alternativ e to merely slo wing diusion at edges is to align the direction of diusion to b e parallel with edges in the image. W eic k ert [45 ] con trols diusion direction b y replacing the scalar c ( x ) in 2.7 with a diusion tensor-v alued function. @ I ( x; t ) @ t = div ( D ( x ) r I ( x; t )) I ( x; 0) = I 0 ( x ) (2.10) The function, D ( x ) pro duces tensors suc h that the unit eigen v ector, e 1 is parallel with r I and e 2 is p erp endicular to r I The eigen v ectors 1 = g ( jr I j ) and 2 = 1 are suggested to generate anisotropic tensors whose represen tativ e ellipsoid has the ma jor axis parallel with image edges. In this w a y diusion near edges still o ccurs, mostly along the edge. By limiting smo othing across the edge, edge lo cation and in tensit y can b e preserv ed.

PAGE 27

17 2.2.2.3 ALM Alv arez, Lions, and Morel [1] prop ose a nonlinear parab olic PDE of the form @ I ( x; t ) @ t = g ( r I ) jr I j div( r I ( x; t ) jr I ( x; t ) j ) I ( x; 0) = I 0 ( x ) (2.11) The eect is that I is smo othed in a direction orthogonal to the gradien t at eac h p oin t. This is b est considered in a lev el-set framew ork. Em b edding I in R 3 b y considering t to b e the third dimension, and ( t = 0) to b e the iso con tours of I w e ha v e @ ( x ( t ) ; t ) @ t + r ( x ( t ) ; t ) x 0 ( t ) = 0 (2.12) Recall the denitions of lev el-set normal N = r jr j and curv ature, = div ( r jr j ). The normal comp onen t of the v elo cit y will b e v n = r jr j x 0 ( t ) (2.13) So 2.12 b ecomes @ ( x ( t ) ; t ) @ t = v n jr j (2.14) Letting v n b e prop ortional to w e obtain 2.11. Curv e-shortening ro w of the iso con tours of smo othes the image I 2.2.3 V ariational F orm ulation The problem of image denoising is ill-p osed, as are man y in v erse problems. Since dieren t pixel in tensities could result in the same v alue when corrupted b y noise, the solution to the denoising problem is not unique. The tec hnique of

PAGE 28

18 Tikhono v regularization in v olv es incorp orating some additional kno wledge ab out the ideal image, similar to the prior distribution in v olv ed in a Ba y esian analysis. In Tikhono v regularization, the solution space is restricted b y p osing the problem as a minimization problem. The function (image) whic h minimizes some stabilizing functional can b e pro v en to b e unique for appropriate c hoice of functional. W e will consider problems of the form min I ( E s ( I ) + E d ( I )) (2.15) The functionals, E in v olv ed in image smo othing often corresp ond to ph ysical mo dels, and usually represen t energy F or instance, the energy asso ciated with data delit y is E d = 1 2 X i k i j I 0 I j 2 (2.16) A ph ysical analogy for the data constrain t is a n um b er of springs b et w een the initial image, I 0 and the ideal image, I The second energy term, E s represen ts the in ternal p oten tial energy of the surface describ ed b y the image. The regularization parameter, is in tro duced to con trol the amoun t of smo othing. 2.2.3.1 Mem brane spline A rst-order stabilizer, stretc hing energy (arc-length), is one p ossible stabilizer E ms = Z Z R 2 I 2 x + I 2 y dx dy (2.17)

PAGE 29

19 The Euler-Lagrange condition for the minimization of 2.17 is r 2 I = 0, the steadystate solution to 2.3. In fact, most of the lters considered so far can b e obtained from a v ariational principle. By setting E s = E ms in 2.15 w e obtain a mem brane spline solution. In this case, con trols the "tension" of the spline surface. 2.2.3.2 Thin-Plate spline Using a second-order stabilizer b ending energy (curv ature) E tps = Z Z R 2 I 2 xx + 2 I 2 xy + I 2 y y dx dy (2.18) as a stabilizer w e obtain a thin-plate spline solution. F or this spline, con trols the "stiness" of the spline surface. The Euler-Lagrange condition for the minimization of 2.18 is r 4 I = 0. 2.2.4 T otal V ariation Image Restoration Another side-constrain t that can b e used as an energy functional is called the total v ariation (TV) norm. This norm represen ts oscillation. In our case, image noise is considered to b e "wrinkling" of the surface describ ed b y the image. Minimizing the TV norm pro duces v ery smo oth images while p ermitting sharp discon tin uities b et w een regions [40]. The TV norm is form ulated b y 2.19. TV n; 1 ( I ( x )) = Z n jr I ( x ) j d x ; n R n (2.19) The TV norm is essen tially the L 1 norm of the image gradien t. A t discon tin uities, the w eak deriv ativ e D I ( x ) to calculate the TV norm. F or a piecewise con tin uous function, the TV norm is then the sum of the TV norm

PAGE 30

20 of eac h con tin uous piece plus the sum of the absolute v alues of the "jumps". F or example, functions f 1 and f 2 in 2.2 ha v e the same TV norm. The oscillatory function, f 1, has a higher TV than the function with a discon tin uit y f 3. Since Figure 2.2: TV( f 1) > TV( f 2) = TV( f 3) most images consist of piecewise smo oth regions separated b y discon tin uities (edges), this is a useful mo del for image denoising. The Euler-Lagrange condition for the minimization of 2.19 written in gradien t-descen t form is @ I ( x; t ) @ t = div ( r I ( x; t ) jr I ( x; t ) j ) I ( x; 0) = I 0 ( x ) (2.20) An example of image denoising b y TV norm minimization is sho wn in 2.3. The restored image on the righ t is c haracterized b y a high degree of in tra-region smo othing, and edges ha v e also b een preserv ed.

PAGE 31

21 Figure 2.3: Noisy image (left) and restored (righ t) b y TV norm minimization. 2.3 PDE Based V ector-V alued Image Restoration The main application of v ector-v alued image restoration has b een the restoration of color images. The simplest implemen tation of v ector-v alued smo othing is uncoupled smo othing. By treating eac h comp onen t of the v ector eld as an indep enden t scalar eld, w e can pro ceed b y smo othing c hannel-b y-c hannel using one of the metho ds presen ted in section 2.2 for scalar image smo othing. This can, ho w ev er, result in a loss of correlation b et w een c hannels as edges in eac h c hannel ma y mo v e indep enden tly due to diusion. T o prev en t this, there m ust b e coupling b et w een the c hannels. There are man y other PDE-based image smo othing tec hniques whic h w e will not co v er here, but will refer the reader to W eic k ert [45 ] and Caselles et al.[10].

PAGE 32

22 2.3.1 V ector-V alued Diusion Whitak er and Gerig in tro duced anisotropic v ector-v alued diusion whic h w as a direct extension of the w ork of P erona and Malik [34 ]. The equations @ I n ( x; t ) @ t = div ( c ( I ) r I n ( x; t )) I n ( x; 0) = I n; 0 ( x ) (2.21) are coupled through the function C and can b e written in v ector form as @ I ( x; t ) @ t = r ( c ( I ) r I ( x; t )) I ( x; 0) = I 0 ( x ) (2.22) 2.3.2 Riemannian Metric Based Anisotropic Diusion Sapiro et al. [41] in tro duced a selectiv e smo othing tec hnique where the selection term is not simply based on the gradien t of the v ector v alued image but based on the eigen v alues of the Riemannian metric of the underlying manifold. In 2 dimensions, the underlying manifold can b e though t of as the parametric surface describ ed b y the image I ( x ). By geometrically smo othing this surface, w e also smo oth the v ector-v alued image. The en tries of the metric tensor, G are dened b y g i;j = @ I @ x i @ I @ x j (2.23) The largest and smallest eigen v alues of G ( + ) and their corresp onding eigen v ectors ( + ) describ e surface prop erties of the Riemannian manifold. The

PAGE 33

23 degree and direction of maximal c hange are giv en b y + and + Smo othing is ac hiev ed b y ev olution in the direction of minimal c hange, that is, along edges in the v ector-v alued image @ I ( x; t ) @ t = g ( + ; ) @ 2 I ( x; t ) @ I ( x; 0) = I 0 ( x ) (2.24) 2.3.3 Beltrami Flo w More recen tly Kimmel et al. [25] presen ted a v ery general ro w called the Beltrami ro w as a general framew ork for image smo othing and sho w that most ro w-based smo othing sc hemes ma y b e view ed as sp ecial cases in their framew ork. @ I i ( x; t ) @ t = 1 p j G j 2 X =1 @ @ x ( p j G j 2 X =1 G @ I i @ x ) I ( x; 0) = I 0 ( x ) (2.25) where G is the metric tensor, and j G j is the determinan t of the metric tensor. The Beltrami ro w can b e though t of as nonlinear diusion on a manifold, where j G j acts as an edge-detecting diusion co ecien t. 2.3.4 Color T otal V ariation Blomgren and Chan [8] in tro duced the T V n;m norm for v ector v alued images. TV n;m ( I ( x )) = vuut m X i =1 [TV n; 1 ( I i )] 2 (2.26)

PAGE 34

24 F or m = 1, 2.26 reduces to the scalar TV norm 2.19. The Euler-Lagrange condition for the minimization of 2.26, written in gradien t descen t form is @ I i ( x; t ) @ t = TV n; 1 ( I i ) TV n;m ( I ) r ( r I i jr I i j ) I ( x; 0) = I 0 ( x ) (2.27) This w as sho wn to b e quite eectiv e for color images, preserving edges in the color space while atten uating noise. Ho w ev er, for m uc h larger dimensional data sets (m=7) as in the w ork prop osed here, the Color TV metho d b ecomes computationally v ery in tensiv e and th us ma y not b e the preferred metho d in suc h applications.

PAGE 35

CHAPTER 3 DT-MRI IMA GE RESTORA TION 3.1 Noise Mo del The degradation asso ciated with the measuremen t of the S v alues is mo delled b y an additiv e gaussian pro cess. Let ^ S ( X ) b e the v ector v alued image that w e w an t to smo oth where, X = ( x; y ; z ) and let S ( X ) b e the unkno wn smo oth appro ximation of the data that w e w an t to estimate. W e ha v e ^ S ( X ) = S ( X ) + Although w e can consider the noise to b e on the comp onen ts of the tensor, D, the form of this noise is no longer so simple. In fact, ev en computing D from S, using the Stejsk al-T anner relation ma y b e meaningless. Substituting the noise mo del in to 2.1 ln S + = ln S 0 3 X i =1 3 X j =1 b i;j D i;j (3.1) The ph ysical principles go v erning diusion do not allo w for negativ e v alues of S, as evidenced b y the dep endence of 3.1 on the natural logarithms of S and S 0 Lo w-in tensit y S measuremen ts ma y b e o v erwhelmed b y noise. Ho w ev er, since is a zero-mean distribution, these measuremen ts ma y b ecome negativ e. Clearly w e cannot calculate D in this case. Nor should w e propagate the noise mo del through the Sk esk al-T anner equation in order to write D as a function of W e prop ose to instead smo oth the v ector of S measuremen ts b efore an y further analysis. 25

PAGE 36

26 3.2 Restoration F orm ulation Smo othing the ra w v ector v alued image data is p osed as a v ariational principle in v olving a rst order smo othness constrain t on the solution to the smo othing problem. W e prop ose a w eigh ted TV-norm minimization for smo othing the v ector v alued image S The v ariational principle for estimating a smo oth S ( X ) is giv en b y min S E ( S ) = Z n [ g ( + ; ) 7 X i = 1 jr S i j + 2 7 X i = 1 j S i ^ S i j 2 ] dx (3.2) where, n is the image domain and is a regularization factor. The rst term here is the regularization constrain t on the solution to ha v e a certain degree of smo othness. The second term in the v ariational principle mak es the solution faithful to the data to a certain degree. The w eigh ting term in this case g ( s ) = 1 = (1 + s ) where s = F A is the fractional anisotrop y dened b y Basser et al. [2]. This selection criteria preserv es the dominan t anisotropic direction while smo othing the rest of the data. Note that since w e are only in terested in the b er tracts corresp onding to the streamlines of the dominan t anisotropic direction, it is apt to c ho ose suc h a selectiv e term. Here w e ha v e used a dieren t TV norm than the one used b y Blomgren and Chan [8]. The T V n;m norm is an L 2 norm of the v ector of T V n; 1 norms for eac h c hannel. W e use the L 1 norm.

PAGE 37

27 The gradien t descen t of the ab o v e minimization is giv en b y @ S i @ t = div g ( + ; ) r S i kr S i k (S i ^ S i ) i = 1 ; :::; 7 @ S i @ n j @ n R + = 0 and S ( x ; t = 0) = ^ S ( x ) (3.3) The deriv ation of equation 3.3 from 3.2 is presen ted in section 8.2.5. The use of a mo died TV norm in 3.2 results in a lo oser coupling b et w een c hannels than the use of the true T V n;m norm w ould ha v e. This reduces the n umerical complexit y of 3.3 and mak es solution for large data set feasible. Note that the TV n;m norm app ears in the gradien t descen t solution 2.27 of the minimization problem. Consider that our data sets consist of 7 images, corresp onding to gradien t directions. Eac h of these images consists of sev eral samples (usually 3 or 4) corresp onding to dieren t gradien t strengths. Calculating the TV n;m norm b y n umerically in tegrating o v er the 3-dimensional data set at eac h step of an iterativ e pro cess w ould ha v e b een prohibitiv ely exp ensiv e.

PAGE 38

CHAPTER 4 NEUR ONAL FIBER TRA CKING 4.1 Ov erview of Neuronal Fib er T rac king W ater in the brain preferen tially diuses along white matter b ers. By trac king the direction of fastest diusion, as measured b y MRI, non-in v asiv e b er trac king of the brain can b e accomplished. Fib ers trac ks ma yb e constructed b y rep eatedly stepping in the direction of fastest diusion. The direction along whic h the diusion is dominan t corresp onds to the direction of eigen v ector corresp onding to the largest eigen v alue of the tensor D In Con turo et al., [18], b er trac ks w ere constructed b y follo wing the dominan t eigen v ector in 0 : 5 mm steps un til a predened measure of anisotrop y fell b elo w some threshold. This usually o ccurred in grey matter. The tensor, D w as calculated at eac h step from in terp olated DT-MRI data. This trac king sc heme is primarily based on heuristics and is not grounded in w ell founded mathematical principles. Mori et al. [30] ac hiev ed b er trac king b y using sev eral heuristics. The trac king algorithm starts from a v o xel cen ter and pro ceeds in the direction of the ma jor axis of the diusion ellipsoid. When the edge of the v o xel is reac hed, the direction is c hanged to that of the neigh b oring v o xel. T rac king stops when a measure of adjacen t b er alignmen t crosses a giv en threshold. The sc heme ho w ev er is resolution dep enden t since the MRI data only rerects a v erage axonal orien tation 28

PAGE 39

29 within a v o xel. Small b ers adjacen t to eac h other ma y not b e distinguished. Another problem o ccurs with branc hing b ers. This metho d will only trac k one path. In this situation, m ultiple p oin ts within a bundle ma y b e indep enden tly trac k ed. W estin et al. [46 ] rep orted tec hniques for pro cessing DT-MRI data using tensor a v eraging. Diusion tensor a v eraging is an in teresting concept but do es not address the issue of estimating a smo oth tensor from the giv en noisy v ectorv alued data. More recen tly P oup on et al. [38 ] dev elop ed a Ba y esian form ulation of the b er tract mapping problem. Prior to mapping the b ers, they use robust regression to estimate the diusion tensor from the ra w v ector v alued image data. No image selectiv e smo othing is p erformed in their w ork prior to application of the robust regression for estimating the diusion tensors. Coulon et al. [20] determined b ers tracts after smo othing the eigen v alues and v ectors. Once again, this sc heme is faced with the same problem i.e., the eigen v ector and the eigen v alues are computed from a noisy tensor eld and hence ma y not b e meaningful at sev eral lo cations in the eld. P ark er et al. [33] also presen ted a follo w up article on b er tract mapping wherein, they use the idea of the fast marc hing metho d from Sethian et al. [42 ] for gro wing seeds initialized in the smo othed dominan t eigen v ector eld. Batc helor et al. [7] rep orted an in teresting b er tract mapping sc heme wherein, they pro duce a map indicating the probabilit y of a b er passing through eac h lo cation in the eld. Ho w ev er, they do not address the issue of computing the diusion tensor from a noisy set of v ector v alued images. This is a v ery imp ortan t issue and should not b e o v erlo ok ed as is demonstrated in the preliminary results of this prop osal.

PAGE 40

30 Giv en the dominan t eigen v ector eld of the diusion tensor in 3D, trac king the b ers (space curv es) tangen tial to this v ector eld is equiv alen t to nding the stream lines/in tegral curv es of this v ector eld. Finding in tegral curv es of v ector elds is a w ell researc hed problem in the eld of Fluid Mec hanics [17]. The simplest solution w ould b e to n umerically in tegrate the giv en v ector eld using a stable n umerical in tegration sc heme suc h as a fourth order Runge-Kutta in tegrator [39]. Ho w ev er, this ma y not yield a regularized in tegral curv e. In the con text of the b er tract mapping, t w o sub-tasks are in v olv ed namely (1) estimating a denoised diusion tensor from noisy v ector-v alued image measuremen ts and (2) estimating the streamlines of the dominan t eigen v ectors of the diusion tensor. The denoising in v olv es a w eigh ted TV-norm minimization of the v ector v alued data S follo w ed b y a linear least squares estimation of D from the smo othed S and then an estimation of regularized stream lines of the dominan t eigen v ectors of D 4.2 F orm ulation Diusion at eac h p oin t can b e c haracterized b y a diusion ellipsoid. The ellipsoid axes are the eigen v ectors of the diusion tensor. The radii are the corresp onding eigen v alues. The shap e of the ellipsoid rerects the isotrop y of diusion. Nearly spherical diusion ellipsoid represen t areas of free w ater, where diusion is unimp eded. Areas of white-matter b er bundles ha v e elongated ellipsoids, since w ater diusion is restricted in directions p erp endicular to b er

PAGE 41

31 direction. The phenomenon w as quan tied b y Basser and Pierpaoli [2]. F A = r 3 2 s ( 1 ) 2 + ( 2 ) 2 + ( 3 ) 2 21 + 22 + 23 (4.1) F or the stream line estimation problem, w e p ose the problem in a v ariational framew ork incorp orating smo othness constrain ts whic h regularize the stream lines/in tegral curv e. The v ariational principle form ulation leads to a PDE whic h can b e solv ed using ecien t n umerical tec hniques. The v ariational principle min c E ( p ) = min c Z 1 0 j c 0 ( p ) j + 2 j c 0 ( p ) v ( c ( p )) j 2 dp (4.2) where, c ( p ) = ( x ( p ) ; y ( p ) ; z ( p )) T is the in tegral curv e w e w an t to estimate and p 2 [0 ; 1] is the parameterization of the curv e, v ( c ( p )) is the v ector eld v restricted to the curv e c ( p ). The rst term of 4.2 is a smo othness constrain t. By minimizing arc-length w e p enalize spurious oscillation in the curv e. The second term pro vides data delit y : w e wish for the b er to b e tangen t to the v ector eld at ev ery p oin t along the curv e. The parameter, con trols the degree of smo othness of the solution. The gradien t descen t of (4.2) @ c @ t = j c 0 ( p ) j k n + [ c 00 ( p ) V ( c ( p )) c 0 ( p )] + V T ( c ( p ))( c 0 ( p ) v ( c ( p ))) (4.3)

PAGE 42

32 where k is the curv ature of the space curv e, is a regularization parameter and V T = The transp ose of V V = J v = 0BBBBBB@ v 1 x v 1 y v 1 z v 2 x v 2 y v 2 z v 3 x v 3 y v 3 z 1CCCCCCA (4.4) The initial condition, c ( p; t = 0) is pro vided b y an ordinary streamline in tegration routine.

PAGE 43

CHAPTER 5 NUMERICAL METHODS 5.1 Data Denoising The nonlinear PDE for denoising the ra w data is div ( g r s jr s j ) ( s s 0 ) = 0 (5.1) where g : R 3 R is the selectiv e smo othing and is the constan t regularization parameter. 5.1.1 Fixed-P oin t Lagged-Diusivit y Equation 5.1 is nonlinear due to the presence of jr s j in the denominator of the rst term. W e linearize 5.1 b y using the metho d of \lagged-diusivit y" presen ted b y Chan and Mulet [13 ]. By considering jr s j to b e a constan t for eac h iteration, and using the v alue from the previous iteration w e can instead solv e 1 jr s t j ( r g r s t + g r 2 s t +1 ) + ( s t +1 s 0 ) = 0 (5.2) Here the sup erscript denotes iteration n um b er. First, rewrite 5.2 with all of the s t terms on the righ t-hand side r 2 s t +1 + jr s t j g s t +1 = jr s t j s 0 + r g r s t g (5.3) 33

PAGE 44

34 5.1.2 Discretized Equations T o write 5.3 as a linear system ( A s t + 1 = y ), discretize the laplacian and gradien t terms. Using cen tral dierences for the laplacian w e ha v e r 2 s t +1 = s t +1 x 1 ;y ;z + s t +1 x;y 1 ;z + s t +1 x;y ;z 1 6 s t +1 x;y ;z + s t +1 x +1 ;y ;z + s t +1 x;y +1 ;z + s t +1 x;y ;z +1 (5.4) Dene the standard forw ard, bac kw ard and cen tral dierences to b e +x s = s x +1 ;y ;z s x;y ;z +y s = s x;y +1 ;z s x;y ;z +z s = s x;y ;z +1 s x;y ;z x s = s x;y ;z s x 1 ;y ;z y s = s x;y ;z s x;y 1 ;z z s = s x;y ;z s x;y ;z 1 x s = 1 2 ( s x +1 ;y ;z s x 1 ;y ;z ) y s = 1 2 ( s x;y +1 ;z s x 1 ;y ;z ) z s = 1 2 ( s x;y ;z +1 s x 1 ;y ;z ) (5.5) W e can rewrite 5.3 in discrete form using the denitions in 5.5 s x 1 ;y ;z s x;y 1 ;z s x;y ;z 1 + (6 + p ( x s t ) 2 +( y s t ) 2 +( z s t ) 2 g ) s x;y ;z s x +1 ;y ;z s x;y +1 ;z s x;y ;z +1 = 1 g ( s 0 p ( x s t ) 2 + ( y s t ) 2 + ( z s t ) 2 + x g x s t + y g y s t + z g z s t ) (5.6)

PAGE 45

35 This results in a banded-diagonal linear system with 7 nonzero co ecien ts p er ro w. 266666666664 6 + jr s t j 0 g 0 1 : : : 1 : : : 1 : : : : : : : : : 1 6 + jr s t j 1 g 1 1 : : : 1 : : : 1 : : : : : : 0 1 6 + jr s t j 3 g 3 1 : : : 1 : : : 1 : : : : : : . . . . . . : : : . : : : . 377777777775 266666666664 s t +1 0 s t +1 1 ... s t +1 n 3 377777777775 = 266666666664 f t 0 f t 1 ... f t n 3 377777777775 (5.7) where the righ t-hand side of 5.6 has b een replaced with f t n The matrix in equation 5.7 is symmetric and diagonally dominan t. W e ha v e successfully used conjugate gradien t descen t to solv e this system. The solution of 5.7 represen ts one xed-p oin t iteration. This iteration is con tin ued un til j s t s t +1 j < c where c is a small constan t. 5.2 Fib er Regularization The PDE for b er regularization is @ c @ t = j c 0 ( p ) j k n + [ c 00 ( p ) V ( c ( p )) c 0 ( p )] + V T ( c ( p ))( c 0 ( p ) v ( c ( p ))) (5.8) where k is the curv ature of the space curv e, is a regularization parameter and V T = The transp ose of V V = J v = 0BBBBBB@ v 1 x v 1 y v 1 z v 2 x v 2 y v 2 z v 3 x v 3 y v 3 z 1CCCCCCA (5.9)

PAGE 46

36 Using the denitions of tangen t, T curv ature, k and normal, n T = c 0 j c 0 j ; k = j T 0 j j c 0 j ; n = T 0 j T 0 j (5.10) w e ha v e @ c @ t = ( c 0 j c 0 j ) 0 + ( c 00 + V T ( c 0 v ) Vc 0 ) (5.11) 5.2.1 Lagged-Diusivit y W e linearize 5.11 using the concept of lagged-diusivit y as w e did in the data denoising case. ( c 0 j c 0 j ) 0 = ( c t + t 0 j c 0t j ) 0 = c t + t 00 j c 0t j (5.12) W e can simplify notation b y calling the co ecien ts of c 00t + t and c 0t + t resp ectiv ely t ( p ) and t ( p ) and calling the additiv e constan t K t @ c @ t = t ( p ) c 00t + t + t ( p ) c 0t + t + K t (5.13) 5.2.2 Crank-Nic holson Metho d W e solv e 5.13 using the Crank-Nic holson metho d. This metho d ac hiev es stabilit y b y using a v eraged dierences to estimate deriv ativ es. F or the second deriv ativ e term w e use @ 2 c @ p 2 = 1 2 p 2 ( c t ( p + p ) 2 c t ( p ) + c t ( p p ) + c t + t ( p + p ) 2 c t + t ( p ) + c t + t ( p p )) (5.14)

PAGE 47

37 The a v eraged dierence for the rst deriv ativ e is @ c @ p = 1 2 p ( c t ( p + p ) c t ( p p ) + c t + t ( p + p ) c t + t ( p p )) (5.15) In nite-dierence form c t + t ( p ) c t ( p ) t = t ( p ) 2 p 2 ( c t + t ( p + p ) 2 c t + t ( p ) + c t + t ( p p )) + t ( p ) 2 p ( c t ( p + p ) c t ( p p ) + t ( p ) 2 p 2 ( c t ( p + p ) 2 c t ( p ) + c t ( p p )) + t ( p ) 2 p ( c t + t ( p + p ) c t + t ( p p ) + K t ( p ) (5.16) Equation 5.16 can b e expressed as a linear system to b e solv ed for c t + t A t c t + t = M t c t + K t (5.17) This tridiagonal linear system can b e solv ed b y Crout's factorization, an optimized LU factorization for this t yp e of matrix.

PAGE 48

CHAPTER 6 NEUR ONAL FIBER VISUALIZA TION 6.1 Rendered Ellipsoids A v ery simple visualization strategy is to simply render the diusion ellipsoid at a subset of data p oin ts. Since a 3D eld of ellipsoids w ould o cclude eac h other, this visualization is usually done for 2D slices of data. Additionally only ellipsoids on a sparse grid can b e rendered in order for eac h ellipsoid to b e discerned. This t yp e of visualization can easily b ecome visually cluttered and con v ey so little information as to b e useless. Laidla w[26 ] ho w ev er, successfully incorp orated ellipsoids in a la y ered visualization approac h. 6.2 Stream tub es Stream tub es are a three-dimensional analogue to streamlines. In fact, the stream tub e is computed b y using a streamline as the cen terline of the tub e. W e can use the stream tub e diameter to enco de some additional information ab out the tensor eld b eing visualized, suc h as F A v alue. Previously Laidla w [27 ] has applied the stream tub e visualization approac h to DT-MRI. 6.3 Line In tegral Con v olution It is also p ossible to visualize the 3D v ector eld corresp onding to the dominan t eigen v alues of the diusion tensor using other visualization metho ds suc h as the line in tegral con v olution tec hnique in tro duced b y Cabral et al. [9 ] { a concept 38

PAGE 49

39 explored in this w ork as w ell. The adv an tage of this visualization tec hnique is that it is w ell suited for visualizing high densit y v ector elds and do es not dep end on the resolution of the v ector eld moreo v er, it also has the adv an tage of b eing able to deal with branc hing structures that cause singularities in the v ector eld. Since Figure 6.1: LIC visualization of syn thetic eld. the b er direction is parallel to the dominan t eigen v ector of the diusion tensor, w e can calculate b er paths as in tegral curv es of the dominan t eigen v ector eld. The stopping criterion is based on F A v alue. When F A falls b elo w 0.17 w e consider the diusion to b e nearly isotropic and stop trac king the b er at this p oin t. Once the diusion tensor has b een robustly estimated, the principal diusion direction can b e calculated b y nding the eigen v ector corresp onding to the dominan t eigen v alue of this tensor. The b er tracts ma y b e mapp ed b y visualizing the streamlines through the eld of eigen v ectors. LIC is a texture-based v ector eld visualization metho d. The tec hnique generates in tensit y v alues b y con v olving a noise texture with a curvilinear k ernel

PAGE 50

40 aligned with the streamline through eac h pixel, suc h as b y I ( x 0 ) = Z s 0 + L s 0 L T ( ( s )) k ( s 0 s ) ds (6.1) where I ( x 0 ) is the in tensit y of the LIC texture at pixel x 0 k is a lter k ernel of width 2 L T is the input noise texture, and is the streamline through p oin t x 0 The streamline, can b e found b y n umerical in tegration, giv en the discrete eld of eigen v ectors. The result is a texture with highly correlated v alues b et w een nearb y pixels on the same streamline, and con trasting v alues for pixels not sharing a streamline. In our case, an F A v alue b elo w a certain threshold can b e a stopping criterion for the in tegration since the diusion eld ceases to ha v e a principal direction for lo w F A v alues. Stalling and Hege [44] ac hiev e signican t computational sa vings b y lev eraging the correlation b et w een adjacen t p oin ts on the same streamline. F or a constan t v alued k ernel, k the in tensit y v alue at I ( ( s + ds )) can b e quic kly estimated b y I ( ( s )) + where is a small error term whic h can b e quic kly computed. Previously Chaing et al. [16] ha v e used LIC to visualize b ers from diusion tensor images of the m y o cardium. 6.4 P articles The LIC and stream tub e tec hniques presen ted in the previous sections are time-consuming op eration. All of the streamline are completely traced b efore an image can b e displa y ed. F or in teractiv e displa y of b ers, w e use a particle

PAGE 51

41 based visualization tec hnique. The particles are analogous to smok e in tro duced in to a wind-tunnel to visualize streamlines. Rather than sim ulating the actual diusion pro cess, the particles simply adv ect through a v elo cit y eld describ ed b y the dominan t eigen v ector of the diusion tensor at eac h p oin t. By seeding a few streamlines within a region of in terest, and p erforming a single step of n umerical in tegration at a time, in teractiv e frame-rates can b e ac hiev ed. Eac h particle is implemen ted as a small textured quadrilateral whic h is alw a ys orien ted to face the view er. W e v ary the size and color of this quad as a means of visualizing the F A v alue at the particle's lo cation in the tensor eld. W e adapt this tec hnique to tensor eld visualization b y incorp orating the F A v alue at eac h eld lo cation in to the LIC texture. By mo dulating the image in tensit y with an increasing function of F A, w e highligh t the areas of white matter, and are able to resolv e where the b ers ev en tually trac k in to grey matter. Unlik e LIC and stream tub e tracing, the particle visualization requires no prepro cessing time. The other tec hniques require completely in tegrating man y p erhaps thousands, of streamlines. P article-based tec hniques allo w immediate visualization.

PAGE 52

CHAPTER 7 EXPERIMENT AL RESUL TS 7.1 Data Denoising In all of the exp erimen ts, w e rst smo oth the sev en 3D directional images using the no v el selectiv e smo othing tec hnique outlined in section 3. F ollo wing this, the diusion tensor is estimated from the smo othed data using a standard least squares tec hnique. The results of F A calculation from the smo othed data, and from ra w data are presen ted in Figures 7.1 and 7.2 Figure 7.1: F A image of coronal slice of ra w rat brain data. Figure 7.2: F A results for smo othed data. 42

PAGE 53

43 7.2 Fib er T rac king 7.2.1 LIC Figures 7.3 and 7.4 depict the computed b er tracts for the reconstructed rat brain data. The in tensit y of the LIC texture has b e mo dulated with the F A image to emphasize the most anisotropic region of eac h image. Figure 7.3: LIC b er tracts in coronal slice of smo othed rat brain data. Figure 7.4: LIC b er tracts in axial slice of smo othed rat brain data.

PAGE 54

44 7.2.2 Stream tub es W e will no w compare our computed stream tub es with a ruoroscopic image. Fib ers are eviden t in in ruoroscopic images as high in tensit y treelik e structures. The ruoro image sho wn on the left of 7.5 sho ws a kno wn anatomical feature, the b er crossings in the corticospinal tract at the base of the brain. The stream tub e map from the same region is sho wn on the righ t of 7.5. The stream tub e map w as Figure 7.5: Comparison of ruoroscopic image (left) with extracted stream tub es (righ t). generated b y starting streamline in tegration at eac h p oin t of a sparse grid within the data if F A > 0 : 3. T rac king stopp ed at F A < 0 : 17. The stream tub e radius is a function of F A. 7.2.3 P articles Figure 7.8 sho ws a b er tracing sequence obtained b y seeding b er in the corpus callosum region of the rat brain data. The brigh ter particles signify high F A v alue. Dimmer particles can b e seen tracing in to grey matter.

PAGE 55

45 Figure 7.6: Axial view of stream tub es in corpus callosum.

PAGE 56

46 Figure 7.7: Details of corpus callosum.

PAGE 57

47 Figure 7.8: Sequence from Fib er visualization (2 seconds b et w een images).

PAGE 58

CHAPTER 8 CONCLUSIONS AND FUTURE W ORK 8.1 Conclusions In this pap er, w e presen ted a new w eigh ted TV-norm minimization form ulation for smo othing v ector-v alued data sp ecically tuned to computation of smo oth diusion tensor MR images. The smo othed v ector v alued data w as then used to compute a diusion tensor image using standard least squares tec hnique. Fib er tracts w ere estimated using the dominan t eigen v ector eld obtained from the diusion tensor image. Finally results of b er tract mapping of a rat spinal cord, and a rat brain w ere depicted using LIC, stream tub e and particles. The b er tracts are quite accurate when insp ected visually and corresp ond w ell with kno wn anatomical structures, suc h as the corpus callosum. 8.2 F uture W ork There are sev eral areas of impro v emen t and directions for future w ork in v olving this researc h. 8.2.1 Noise Mo del F or Measured Data The S v ector that w as smo othed in this w ork is not the actual measured data. These S v alues are actually the magnitude of F ourier transformed data. A noise mo del expressed in terms of the v oltage lev els observ ed during the imaging pro cess w ould b e fundamen tally correct and ma y lead to more accurate b er maps. 48

PAGE 59

49 8.2.2 Quan titativ e V alidation Ho w ev er, quan titativ e v alidation of the computed b er tracts is essen tial. This can b e ac hiev ed b y registration of our three-dimensional b ers with t w odimensional ruoroscop y images. This will b e the fo cus of our future eorts. 8.2.3 Robust Regression The time in tensiv e nature of DT-MRI data collection limits the n um b er of samples tak en in eac h direction. Outliers in the collected data ma y ha v e a signican t impact on the least-squares calculation of the diusion tensor. The use of a robust regression for calculating the diusion tensor should b e explored. 8.2.4 Non-T ensor Mo del of Diusion An alternativ e to the tensor mo del of diusion is to c haracterize diusion with a spherical distribution of diusion magnitudes. This requires measuremen ts to b e made in man y more directions, but has the adv an tage of handling b er crossings o ccurring within a single v o xel. 8.2.5 Automatic Region-Of-In terest Extraction By registering a new data set with an atlas, the user could select the regionof-in terest b y anatomical name. Curren tly the region-of-in terest for seeding b er tracts is done b y hand. This metho d is time consuming and p oten tially inaccurate.

PAGE 60

APPENDIX DERIV A TION OF EULER-LA GRANGE CONDITIONS Here w e will review the calculus of v ariations, and presen t the deriv ation the Euler-Lagrange equations for the image denoising and b er regularization PDEs. Ev ans [22] pro vides a thorough treatmen t of v ariational calculus. In v ariational calculus w e seek to obtain minima or maxima of expressions that dep end on functions instead of parameters. The expressions are referred to as functionals. Consider a general functional, I dep ending on function f ( x; y ; z ) I [ f ] = Z Z Z F ( f ; f x ; f y ; f z ; x; y ; z ) dx dy dz (1) W e seek the function f whic h minimizes the expression I Supp ose the function f minimizes the functional I Consider constructing a v ariation of f using parameter and test function This v ariation is giv en b y f + Consider a new function of one v ariable, i ( ) = I [ f + ]. Extrema of this function o ccur when di d = 0 (2) Since w e kno w that f minimizes I w e can sa y di d j =0 = 0 (3) is a condition for the minimization of I 50

PAGE 61

51 As an example, consider the v ariation of 1 i ( ) = I [ f + ] = Z Z Z F ( f + ; f x + x ; f y + y ; f z + z ; x; y ; z ) dx dy dz (4) Dieren tiate with resp ect to i 0 ( ) = Z Z Z @ F @ ( f + ) + @ F @ ( f x + x ) x + @ F @ ( f y + y ) y + @ F @ ( f z + z ) z dx dy dz (5) Ev aluate this deriv ativ e at = 0 i 0 (0) = Z Z Z @ F @ f + @ F @ f x x + @ F @ f y y + @ F @ f z z dx dy dz = 0 (6) In tegrating b y parts, and imp osing the condition that the test function, equals 0 on the b oundary w e get the result @ F @ f d dx ( @ F @ f x ) d dy ( @ F @ f y ) d dz ( @ F @ f z ) = 0 (7) This is the condition for minimization of the functional in 1 When second deriv ativ es are in v olv ed, suc h as the functional I [ f ] = Z Z Z F ( f ; f x ; f x x; x; y ; z ) dx dy dz (8) Then the in tegration b y parts step m ust b e p erformed t wice in order to write i 0 (0) in terms of This leads to the Euler-Lagrange condition @ F @ f d dx ( @ F @ f x ) + d 2 dx 2 ( @ F @ f xx ) = 0 (9)

PAGE 62

52 Data denoising W e can rewrite 3.2 as min s Z Z Z g q s 2x + s 2y + s 2z + 2 ( s s 0 ) 2 dx dy dz (10) Using the general result deriv ed in 7, the Euler-Lagrange condition for minimizing 10 is ( s s 0 ) d dx ( g s x jr s j ) d dy ( g s y jr s j ) d dz ( g s z jr s j ) = 0 (11) This can b e rewritten in v ector notation as ( s s 0 ) div ( g r s jr s j ) = 0 (12) Fib er regularization The rst v ariation of 4.2 with v ariation and test function is I [ c + ] = Z 1 0 j c 0 ( p ) + 0 j + 2 j c 0 ( p ) + 0 v ( c + ( p )) j 2 dp (13) Dieren tiating with resp ect to giv es @ I [ c + ] @ = Z 1 0 c 0 ( p ) + 0 j c 0 ( p ) + 0 j 0 + ( c 0 ( p ) + 0 v ( c + ))( 0 r v ( c + ) ) (14) The Euler-Lagrange condition is giv en b y @ I @ j =0 = Z 1 0 c 0 ( p ) j c 0 ( p ) j 0 + ( c 0 ( p ) v ( c ( p ))( 0 r v ( c ( p )) ) = 0 (15) In tegrating b y parts w e obtain c 0 ( p ) j c 0 ( p ) j 0 + [ c 00 ( p ) V ( c ( p )) c 0 ( p )] + V T ( c ( p ))( c 0 ( p ) v ( c ( p ))) (16)

PAGE 63

53 where V is the Jacobian of v dened as V = J v = 0BBBBBB@ v 1 x v 1 y v 1 z v 2 x v 2 y v 2 z v 3 x v 3 y v 3 z 1CCCCCCA (17)

PAGE 64

REFERENCES [1] L. Alv arez, P L. Lions, and J. M. Morel, \Image selectiv e smo othing and edge detection b y nonlinear diusion. ii," SIAM J. Numer. A nal. v ol. 29, no. 3, pp. 845{866, June 1992. [2] P J. Basser and C. Pierpaoli \Microstructural and ph ysiological features of tissue elucidated b y quan titativ e-diusion-tensor MRI," J. Magn. R eson. B v ol. 110, 209-219, 1996. [3] P J. Basser, J. Mattiello and D. LeBihan "MR diusion tensor sp ectroscop y and imaging" Biophys. J. v ol. 66, 259-267, 1994a. [4] P J. Basser, J. Mattiello and D. LeBihan "Estimate of the eectiv e selfdiusion tensor from the NMR spin ec ho," J. Magn. R eson. B v ol. 103, 247-254, 1994b. [5] P J. Basser "Inferring microstructural features and the ph ysiological state of tissues from diusion w eigh ted images" NMR in Biome d. v ol. 8, 333-344, 1995. [6] P J. Basser, S. P a jevic, C. Pierpaoli, J. Duda and A. Aldroubi, \In viv o b er tractograph y using DT-MRI data," Magnetic R esonanc e in Me dicine v ol. 44, 2000, pp. 625-632. [7] P G. Batc helor, D. L. G. Hill, F. Calaman te and D. A tkinson, \Study of connectivit y in the brain using the full diusion tensor from MRI," in Pr o c e e dings of the 17th Intl. Conf. on Information Pr o c essing in Me dic al Imaging 2001, Da vis, CA, pp. 121-133. [8] P Blomgren and T. F. Chan,"Color TV: total v ariation metho ds for restration of v ector-v alued images," IEEE T r ansaction on Image Pr o c essing v ol. 7, no. 3, pp. 304-309, Marc h 1998. [9] B. Cabral, \Imaging v ector elds using line in tegral con v olution," Computer Gr aphics Pr o c e e dings pp. 263-270, 1993. 54

PAGE 65

55 [10] V. Caselles, J. M. Morel, G. Sapiro and A. T annen baum, IEEE T r ans. on Image Pr o c essing sp ecial issue on PDEs and geometry-driv en diusion in image pro cessing and analysis, v ol. 7, No. 3, 1998. [11] F. Catte, P L. Lions, J. M. Morel, and T. Coll, \Image selectiv e smo othing and edge detection b y nonlinear diusion," SIAM Journal on Numeric al A nalysis v ol. 29, pp. 182{193, 1992. [12] P Charb onnier, L. Blanc-F eraud, G. Aub ert, and M. Barlaud, \Tw o deterministic half-quadratic regularization algorithms for computed imaging,," in in Pr o c. of the IEEE Intl. Conf. on Image Pr o c essing (ICIP), 1994, v ol. 2, pp. 168{172, IEEE Computer So ciet y Press. [13] T. Chan and P Mulet, \On the con v ergence of the lagged diusivit y xed p oin t metho d in total v ariation image restoration," SIAM Journal on Numeric al A nalysis v ol. 36, no. 2, pp. 354{367, 1999. [14] C. Chefd'hotel, D. Tsc h ump erl, R. Deric he, O. F augeras, "Constrained ro ws of matrix-v alued functions : Application to diusion tensor regularization," in Pro c. of the Eur op e an Confer enc e on Computer Vision 2002 (ECCV'02), Cop enhagen, Denmark, Ma y 2002, pp. 251-265. [15] Y. Chen, B. C. V em uri and L. W ang,\Image denoising and segmen tation via nonlinear diusion," Computers and Mathematics with Applic ations v ol. 39, No. 5/6, pp. 131-149, 2000. [16] P .-J. Chiang, B. Da vis and E. Hsu, "Line-in tegral con v olution reconstruction of tissue b er arc hitecture obtained b y MR diusion tensor imaging," BMES A nnual Me eting 2000. [17] A. Chorin, Computational Fluid Me chanics, Sele cte d p ap ers Academic Press, NY, 1989. [18] T. E. Con turo, R. C. McKinstry E. Akbudak, and B. H. Robinson "Enco ding of anisotropic diusion with tetrahedral gradien ts: A general mathematical diusion formalism and exp erimen tal results" Magn. R eson. Me d. v ol. 35, 399-412, 1996 [19] T. E. Con turo, N. F. Lori, T. S. Cull, E. Akbudak, A. Z. Sn yder, J. S. Shimon y R. C. McKinstry H. Burton and M. E. Raic hle "T rac king neuronal b er path w a ys in the living h uman brain" Pr o c. Natl. A c ad. Sci. USA v ol. 96, 10422-10427, 1999.

PAGE 66

56 [20] O. Coulon, D. C. Alexander and S. R. Arridge, \A regularization sc heme for diusion tensor magnetic resonance images," in Pr o c e e dings of the 17th Intl. Conf. on Information Pr o c essing in Me dic al Imaging Da vis, CA, pp. 92-105, 2001. [21] F. Cric k and E. Jones "Bac kw ardness of h uman neuroanatom y ," Natur e v ol. 361, 109-110, 1993. [22] L.C.Ev ans, Partial Dier ential Equations American Mathematical So ciet y Pro vidence RI, 1998. [23] E. Haac k e, R. Bro wn, M. Thompson, R. V enk atesan. Magnetic R esonanc e Imaging{Physic al Princip als and Se quenc e Design John Wiley and Sons, NY, 1999. [24] D. K. Jones, A. Simmons, S. C. R. Williams and M. A. Horseld, "Nonin v asiv e assessmen t of axonal b er connectivit y in the h uman brain via diusion tensor MRI" Magn. R eson. Me d. v ol. 42, 37-41, 1999. [25] R. Kimmel, N. So c hen, and R. Malladi, \Images as em b edding maps and minimal surfaces:mo vies, color and v olumetric medical images," in Pr o c. of the IEEE Conf. on Computer Vision and Pattern R e c o gnition pp. 350{355, June 1997. [26] D. H. Laidla w, E. T. Ahrens, D. Kremers, M. J. Av alos, \Mouse spinal cord diusion tensor visualization using concepts from pain ting," Siggraph T ec hnical Slide Set, 1998. [27] D. H. Laidla w, S. Zhang, C. Curry D. Morris, \Stream tub es and streamsurfaces for visualizing diusion tensor MRI v olume images," in Pr o c. of the IEEE Visualization Octob er 2000. [28] R. Malladi and J. A. Sethian, \A unied approac h to noise remo v al, image enhancemen t and shap e reco v ery ," IEEE T r ans. on Image Pr o c essing v ol. 5, no. 11, pp. 1554{1568, 1996. [29] N. Makris, A. J. W orth, A. G. Sorensen, G. M. P apadimitriou, O. W u, T. G. Reese, V. J. W edeen, T. L. Da vis, J. W. Stages, V. S. Ca viness, E. Kaplan, B. R. Rosen, D. N. P andy a, and D. N. Kennedy "Morphometry of in viv o h uman white matter asso ciation path w a ys with diusion-w eigh ted magnetic resonance imaging," A nn. Neur ol. v ol. 42, pp. 951-962, 1999.

PAGE 67

57 [30] S. Mori, B. J. Crain, V. P Chac k o and P C. M. v an Zijl, "Three-dimensional trac king of axonal pro jections in the brain b y magnetic resonance imaging" A nn. Neur ol. v ol. 45, pp. 265-269, 1999. [31] M. Nitzb erg and T. Shiota, \Nonlinear image ltering with edge and corner enhancemen t," IEEE T r ansactions on Pattern A nalysis and Machine Intel ligenc e v ol. 14, no. 8, pp. 826{832, 1992. [32] P J. Olv er, G. Sapiro, and A. T annen baum, \In v arian t geometric ev olutions of surfaces and v olumetric smo othing," SIAM J. Appl. Math. v ol. 57, pp. 176{194, 1997. [33] G. J. M. P ark er, C. A. M. Wheeler-Kingh tshott and G. J. Bak er, \Distributed anatomical brain connectivit y deriv ed from diusion tensor imaging," in Pr o c e e dings of the 17th Intl. Conf. on Information Pr o c essing in Me dic al Imaging Da vis, CA, pp. 106-120, 2001. [34] P P erona and J. Malik, \Scale-space and edge detection using anisotropic diusion," IEEE T r ans. P AMI v ol. 12, no. 7, pp. 629{639, 1990. [35] C. Pierpaoli, and P J. Basser "T o w ard a quan titativ e assessmen t of diusion anisotrop y ," Magn. R eson. Me d. v ol. 36, 893-906, 1996. [36] C. Pierpaoli, P Jezzard, P J. Basser, A. Barnett and G. Di Chiro "Diusion tensor MR imaging of the h uman brain," R adiolo gy v ol. 201, pp. 637-648, 1996. [37] L. Lapidus and G. F. Pinder, Numeric al solution of p artial dier ential e quations in scienc e and engine ering John Wiley and Sons, NY, 1982. [38] C. P oup on, C. A. Clark, V. F rouin, J. Regis, I. Blo c h, D. LeBihan, and J. Mangin, \Regularization of diusion-based direction maps for the trac king of brain white matter fascicles," Neur oImage v ol. 12, 184-195, 2000. [39] W. H. Press, B. P Flannery S. A. T euk olsky and W. T. V etterling, Numeric al r e cip es in C: The art of scientic c omputing Cam bridge Univ ersit y Press, Cam bridge, England, second edition, 1992. [40] L. Rudin, S. Osher, E. F atemi, \Nonlinear total v ariation based noise reduction algorithms," Physic a D v ol. 60, pp. 259-268, 1992.

PAGE 68

58 [41] G. Sapiro and D. L. Ringac h, \Anisotropic diusion of m ultiv alued images with applications to color ltering," IEEE T r ans. on Image Pr o c essing v ol. 5, pp. 1582{1586, 1996. [42] J. Sethian, L evel-set metho ds Cam bridge Univ ersit y Press, Cam bridge, 1998. [43] J. Shah, \A common framew ork for curv e ev olution, segmen tation and anisotropic diusion," in IEEE Conf. on Computer Vision and Pattern R e c o gnition pp. 136-142, 1996. [44] D. Stalling and H. Hege, \F ast and indep enden t line in tegral con v olution," in A CM SIGGRAPH pp. 249-256, 1995. [45] J. W eic k ert, \A review of nonlinear diusion ltering,," in Sc ale-sp ac e the ory in c omputer vision, (Eds.) B. ter Haar Romney L.Florac k, J. Ko enderink, and M. Viergev er, Eds., v ol. 1252, of L e ctur e Notes in Computer Scienc e, pp. 3{28, Springer-V erlag, 1997. [46] C. F. W estin, S. E. Maier, B. Khidir, P Ev erett, F. A. Jolesz and R. Kikinis, \Image pro cessing for diusion tensor magnetic resonance imaging," in Pr o c. of the Se c ond Intl. Conf. on Me dic al Image Computing and Computer Assiste d Interventions (MICCAI), pp. 441-452, 1999. [47] R. Xue, P C. M. v an Zijl, B. J. Crain, H. Solaiy appan, and S. Mori "In viv o three-dimensional reconstruction of rat brain axonal pro jections b y diusion tensor imaging," Magn. R eson. Me d. v ol. 42, pp. 1123-1127, 1999.

PAGE 69

BIOGRAPHICAL SKETCH Tim McGra w w as b orn in Key W est, Florida. He receiv ed his Bac helor of Science degree from the Mec hanical Engineering Departmen t at the Univ ersit y of Florida. He will receiv e his Master of Science degree in computer and information sciences from the Univ ersit y of Florida in August 2002. His researc h in terests include medical imaging, image pro cessing, computer vision and computer graphics. 59


Permanent Link: http://ufdc.ufl.edu/UFE0000573/00001

Material Information

Title: Neuronal fiber tracking in Dt-Mri
Physical Description: Mixed Material
Creator: McGraw, Tim E. ( Author, Primary )
Publication Date: 2002
Copyright Date: 2002

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0000573:00001

Permanent Link: http://ufdc.ufl.edu/UFE0000573/00001

Material Information

Title: Neuronal fiber tracking in Dt-Mri
Physical Description: Mixed Material
Creator: McGraw, Tim E. ( Author, Primary )
Publication Date: 2002
Copyright Date: 2002

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0000573:00001


This item has the following downloads:


Full Text











NEURONAL FIBER TRACKING IN DT-MRI


By

TIM E. MCGRAW















A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE

UNIVERSITY OF FLORIDA


2002


































Copyright 2002

by

Tim E. McGraw
















For my wife, Jo, and my mother, Patti.















ACKNOWLEDGMENTS


I would like to express thanks to everyone who made this research possible.

First, thanks go to Dr. Baba C. Vemuri for serving as the chairman of my thesis

committee, and being available for consultation and advice. Thanks go, also, to

Dr. Jorg Peters for serving on my committee and imparting the graphics knowledge

required for the visualization aspect of this project.

Also, many thanks go to the members of the Computer Vision, Graphics and

Medical Imaging (CVGMI) group, especially Zhizhou Wang, Jundong Liu and Fei

Wang, for their assistance in preparing this work.

I am grateful to the faculty members of the UF Mathematics Department, Dr.

Yunmei Chen, for serving on my committee, and Dr. Murali Rao, for advisement

and encouragement.

I would also like to acknowledge the McKnight Brain Institute for providing

data, analysis and validation of results. Thanks go to Dr. Tom Mareci, Dr. Steve

Blackband, Dr. Paul Reier, and Everen Ozarslan.

This research was funded in part by the NIH grant RO1-NS42075.















TABLE OF CONTENTS

page

ACKNOWLEDGMENTS ................... ...... iv

LIST OF FIGURES ..................... ......... vii

KEY TO ABBREVIATIONS ................... ....... viii

ABSTRACT ....................... ........... ix

CHAPTERS

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

1.1 Overview ............ ................. 1
1.2 Contributions .. .. .. ... .. .. .. .. ... .. .. .. 2
1.3 Outline of thesis ........... ............. 3

2 BACKGROUND OVERVIEW .......... ....... .... 5

2.1 DT-MRI Data Acquisition ........ ........... 5
2.1.1 Overview of Diffusion ........ ......... 8
2.1.2 Overview of MR Imaging ..... ....... 9
2.1.3 DT-MRI Acquisition ......... ......... 10
2.1.4 Stejskal-Tanner Equation ...... ........ 11
2.2 PDE Based Scalar-Valued Image Restoration ........ 12
2.2.1 Linear Filtering ........ ......... ... 13
2.2.2 Nonlinear Filtering ........ .......... 14
2.2.2.1 Perona-M alik ....... ......... 15
2.2.2.2 Tensor anisotropic diffusion .......... 16
2.2.2.3 ALM .............. ... 17
2.2.3 Variational Formulation ..... . . ..... 17
2.2.3.1 Membrane spline . . ..... 18
2.2.3.2 Thin-Plate spline .............. .. 19
2.2.4 Total Variation Image Restoration . .... 19
2.3 PDE Based Vector-Valued Image Restoration . ... 21
2.3.1 Vector-Valued Diffusion . . . 22
2.3.2 Riemannian Metric Based Anisotropic Diffusion . 22
2.3.3 Beltrami Flow ..... ........... ...... 23
2.3.4 Color Total Variation ............. .. 23









3 DT-MRI IMAGE RESTORATION .................. 25

3.1 Noise M odel ................... .... 25
3.2 Restoration Formulation ........ ............ 26

4 NEURONAL FIBER TRACKING ....... .......... 28

4.1 Overview of Neuronal Fiber Tracking ...... ....... 28
4.2 Formulation ................... .... 30

5 NUMERICAL METHODS .......... ............. 33

5.1 Data Denoising ......................... 33
5.1.1 Fixed-Point Lagged-Diffusivity ...... ..... 33
5.1.2 Discretized Equations ........ .......... 34
5.2 Fiber Regularization. ................ ..... 35
5.2.1 Lagged-Diffusivity .. ........ 36
5.2.2 Crank-Nicholson Method ............... .. 36

6 NEURONAL FIBER VISUALIZATION .............. .. 38

6.1 Rendered Ellipsoids .................. .. 38
6.2 Streamtubes .................. ........ .. 38
6.3 Line Integral Convolution ................... .. 38
6.4 Particles .................. .......... .. 40

7 EXPERIMENTAL RESULTS .................. ... 42

7.1 Data Denoising .................. ... .. 42
7.2 Fiber Tracking .................. ....... .. 43
7.2.1 LIC . . . . . ..... 43
7.2.2 Streamtubes .................. .. 44
7.2.3 Particles .................. ..... .. 44

8 CONCLUSIONS AND FUTURE WORK . . ..... 48

8.1 Conclusions .................. ...... .. .. 48
8.2 Future Work ....... ..... ........ 48
8.2.1 Noise Model For Measured Data . . ... 48
8.2.2 Quantitative Validation ................ .. 49
8.2.3 Robust Regression ............ .. .. .. 49
8.2.4 Non-Tensor Model of Diffusion . . ..... 49
8.2.5 Automatic Region-Of-Interest Extraction ...... .. 49

APPENDIX DERIVATION OF EULER-LAGRANGE CONDITIONS 50

REFERENCES ................... .......... .... .. 54

BIOGRAPHICAL SKETCH. .................. ........ .. 59















LIST OF FIGURES

Figure page

2.1 Diffusion Ellipsoid ............................. 8

2.2 TV(fl) > TV(f2) = TV(f3) .................. ..... 20

2.3 Noisy image (left) and restored (right) by TV norm minimization. .. 21

6.1 LIC visualization of synthetic field. .................... .. 39

7.1 FA image of coronal slice of raw rat brain data. . ..... 42

7.2 FA results for smoothed data. ................ ..... 42

7.3 LIC fiber tracts in coronal slice of smoothed rat brain data. . 43

7.4 LIC fiber tracts in axial slice of smoothed rat brain data. ...... ..43

7.5 Comparison of fluoroscopic image (left) with extracted streamtubes
(right). . . . . . . . . .. 44

7.6 Axial view of streamtubes in corpus callosum. ........... ..45

7.7 Details of corpus callosum. .................. ..... 46

7.8 Sequence from Fiber visualization (2 seconds between images). . 47















KEY TO ABBREVIATIONS


CNS: central nervous system

DT: diffusion tensor

DTI: diffusion tensor imaging

FA: fractional anisotropy

LIC: line integral convolution

MRI: magnetic resonance imaging

PDE: partial differential equation

RF: radio frequency

ROI: region of interest

TV: total variation















Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Science



NEURONAL FIBER TRACKING IN DT-MRI



By

Tim E. McGraw

December 2002

Chair: Baba C. Vemuri
Major Department: Computer and Information Sciences and Engineering

Diffusion tensor imaging (DTI) can provide the fundamental information

required for viewing structural connectivity. However, robust and accurate acqui-

sition and processing algorithms are needed to accurately map the nerve connec-

tivity. In this thesis, we present a novel algorithm for extracting and visualizing

the fiber tracts in the CNS, specifically the spinal cord. The automatic fiber tract

mapping problem will be solved in two phases, namely a data smoothing phase and

a fiber tract mapping phase. In the former, smoothing is achieved via a weighted

TV-norm minimization, which strives to smooth while retaining all relevant detail.

For the fiber tract mapping, a smooth 3D vector field indicating the dominant

anisotropic direction at each spatial location is computed from the smoothed data.

Neuronal fibers are traced by calculating the integral curves of this vector field.









Results are expressed using three modes of visualization. Line integral convo-

lution (LIC) produces an oriented texture which shows fiber pathways in a planar

slice of the data. A streamtube map is generated to present a three-dimensional

view of fiber tracts. Additional information, such as degree of anisotropy, can be

encoded in the tube radius, or by using color. A particle system form of visualiza-

tion is also presented. This mode of display allows for interactive exploration of

fiber connectivity with no additional preprocessing.















CHAPTER 1
INTRODUCTION


1.1 Overview


The understanding of neurological function, and medical diagnosis of disease

and injury require knowledge of the structural organization of the brain. Magnetic

Resonance Imaging provides a non-invasive means of studying anatomical con-

nectivity. This connectivity, in the form of axonal nerve fiber bundles, can only

be measured indirectly. The presence of these fibers must be inferred from the

behavior of water molecules near the tissue.

Diffusion tensor magnetic resonance imaging (DT-MRI) is a method of mea-

suring the rate of water diffusion in biological structures. A diffusion tensor at each

location on a regular lattice describes a volumetric average of the directional prop-

erties of diffusion within each voxel. The observation that diffusion is anisotropic

in areas of white-matter fiber bundles allows tracking of fibers through the lattice.

However, acquisition noise corrupts the data measurements. Voltage variations in

the receiving coil of the MRI machine due to thermal noise are a major source of

signal degradation. The noise is accepted to be zero-mean, and additive in nature.

Integrating a fiber path over a noisy field is an ill-posed problem. The goals of the

work presented in this thesis are to trace fibers through an ideal, unknown field

given only noisy observations, and to visualize the fibers in a useful way.









The fiber tracts are estimated in a two stage process, (a) smoothing of

the raw directional images acquired for varying magnetic field strengths and

estimating the diffusion tensor, D, field over the 3D image lattice, (b) computing

the dominant eigenvector field from this regularized D and estimating regularized

streamlines/integral curves as the desired fiber tracts.

Computed fibers must be visualized in a way that is useful for anatomical

study or medical diagnosis. Existing vector field visualization techniques, such

as streamlines/streamtubes can be used to visualize fibers. These techniques can

also be modified to convey more information about diffusion by incorporating

quantities derived from the tensor. Here, we adapt line integral convolution (LIC),

streamtubes, and particles to the task of tensor field visualization.

We will review the literature on the raw data smoothing, streamline genera-

tion, and visualization techniques where appropriate.

1.2 Contributions

Although is well known that the DT-MRI data requires denoising prior to

fiber tracking, previous work has concentrated on smoothing the vector field of

dominant diffusion directions, or even the diffusion tensor itself. Employing these

methods requires that the diffusion tensor, and possibly i, L-. -\ 4--t.' i of the tensor

be computed from noisy data. Results of these calculations may be meaningless

due to properties of the noise, and the diffusion tensor calculation. We propose to

smooth the data before these calculation to provide for more ]li\ -iP .'lly meaningful

results.









A new application of particle-based visualization is presented. This technique

has previously been used in fluid mechanics for visualization of velocity fields. Here

we will adapt the technique to tensor field visualization. This technique has no

associated preprocessing time, and allows real-time interactive visualization of the

reconstructed data. Many other techniques require extensive computing time for

streamline integration prior to visualization.

1.3 Outline of thesis

The rest of this thesis is organized as follows. In Chapter 2 we will discuss the

diffusion process and give a short introduction to the DT-MRI data acquisition

process. An overview of PDE based image denoising will follow. Some linear

and nonlinear techniques for scalar and vector-valued images will be reviewed.

Variational techniques for formulating image denoising, and the TV norm will be

presented.

Chapter 3 will describe our vector-valued image restoration technique, a

weighted TV norm minimization. The noise model used to formulate this technique

is discussed, as well as the motivation for modifying the TV norm usually used for

vector-valued images.

In Chapter 4 we present the process of tracking neuronal fiber bundles through

the restored DT-MRI data. Previous methods of fiber tracking will be discussed.

A quantification of diffusion tensor anisotropy will be given. In this chapter

streamline regularization will be presented as a way of enforcing a smoothness

constraint on the fibers.









Chapter 5 will show the linearization and discretization of the PDEs involved

in data denoising and fiber regularization. The numerical methods used to solve

linear equations will be described.

In Chapter 6 we will present visualization techniques for DT-MRI data. These

include rendered ellipsoids, streamtubes, LIC and particles.

In Chapter 7 experimental results will be reported. The data denoising results

for spinal cord and brain data sets will be presented. Fiber tracking results using

several visualization techniques will be presented.

Conclusions and possible directions for future work will be discussed in

Chapter 8.















CHAPTER 2
BACKGROUND OVERVIEW


2.1 DT-MRI Data Acquisition

Fundamental advances in understanding living biological systems require

detailed knowledge of structural and functional organization. This is particularly

important in the nervous system where anatomical connections determine the

information pathways and how this information is processed. Our current under-

standing of the nervous system is incomplete because of a lack of fundamental

structural information [21] necessary to understand function. In addition, under-

standing fundamental structural relationships is essential to the development and

application of therapies to treat pathological conditions (e.g. disease or injury).

However, most imaging methods give only an anatomically isolated represen-

tation of living tissue because the images do not contain connectivity information.

Such information would allow the identification and correlation of system elements

responding during function. For example in brain trauma, the relationship between

anatomy and behavior will only become apparent when we are able to discrimi-

nate the afferent nerve fiber pathways transmitting the sensation from a stimulus

to the brain or the efferent pathways transmitting impulses from the brain area

controlling behavior.









Research during the last few decades has shown that the central nervous

system is able to adapt to challenges and recover some function. However, the

structural basis for this adaptive ability is not well understood. For the entire

central nervous system, understanding and treating evolving pathology, such as

brain trauma, depends on a detailed understanding of the anatomical connectivity

changes and how they relate to function.

Recently MR measurements have been developed to measure the tensor of

diffusion. This provides a complete characterization of the restricted motion of

water through the tissue that can be used to infer tissue structure and hence

fiber tracts. In a series of papers, Basser and colleagues [2, 3, 4, 5, 6, 36] have

discussed in detail general methods of acquiring and processing the complete

apparent-diffusion-tensor of MR measured translational self-diffusion. They

showed that directly measured diffusion tensors could be recast in a rotationally

invariant form and reduced to parametric images that represent the average

rate of diffusion (tensor trace), diffusion anisotropy (relationship of eigenvalues),

and how the diffusion ellipsoid (eigenvalues and eigenvectors) can be related to

the laboratory reference frame. The parametric images [5, 35] of volume ratio,

fractional anisotropy, and lattice anisotropy index represent scalar measures of

diffusion that are independent of the lab reference frame and subject orientation.

Therefore, these measures can be used to characterize the tissue pathology, e.g.,

ischemia, independent of the specific frame of reference used to acquire the images.

The development of diffusion tensor acquisition, processing, and analysis methods

provides the framework for creating fiber tract maps based on this complete









diffusion tensor analysis [19, 24, 29, 30]. This has been used to produce fiber tract

maps in rat brains [30, 47] and to map fiber tracts in the human brain [24] then the

first steps were taken to relate this structural connectivity to function [19].

The directional properties of diffusion can be characterized by a diffusion

tensor, a 3 x 3 symmetric matrix of real values. In order to calculate the 6 indepen-

dent components of the tensor, the subject is imaged in seven different directions

with several magnetic field strengths. The relationship S = So exp(- EC byjDj)

allows the diffusion tensor, D, and the T-2 weighted image So to be calculated

given the samples S and field gradient strength, b. Previous work has concentrated

on smoothing the field of eigenvectors of D. More recent work [14] has formulated

regularization techniques for the tensor field itself, even constraining the resulting

tensors to be positive-definite. Here we will take the approach of smoothing the

observed vector-valued image S prior to calculating D. This decision to smooth in

this manner will be justified in Chapter 3.

In summary, the anisotropy of water translational diffusion can be used

to visualize structure in the brain and provides the basis for a new method of

visualizing nerve fiber tracts. Initial results have been very encouraging and

-.', -4 that this approach to fiber mapping may be applied to a wide range of

studies in living subjects. However, it is essential to optimize the acquisition and

processing algorithms for fiber tract mapping and validate the results relative to

known measures of fiber tracts.






















Figure 2.1: Diffusion Ellipsoid


2.1.1 Overview of Diffusion

Random molecular motion (Brownian motion) can cause transport of matter

within a system. Within a volume of water, molecules freely diffuse in all direc-

tions. The water abundant in biological systems is also subject to such stochastic

motion. The properties of the surrounding tissue can affect the magnitude of diffu-

sion, and the directional properties as well. Tissue can form a barrier to diffusion,

restricting molecular motion.

Within an oriented structure, such as a bundle of axonal fibers, diffusion

can be highly anisotropic. The white matter of the brain and spinal cord is

characterized by many such bundles.

The directional properties of diffusion can be characterized by a tensor. The

diffusion tensor, D, is a symmetric, positive-definite 3 x 3 matrix. We will make

use of the eigenvalues and eigenvectors of this tensor, sorting the eigenvalues

(A1, A2, A3) from largest to smallest, and labelling the corresponding unit eigen-

vectors (el, e2, e3). The eigenvalues represent the magnitude of diffusion in the

direction of their corresponding eigenvector. For isotropic diffusion A = A2 =A3 A









popular representation for describing anisotropic diffusion is the diffusion ellipsoid.

This ellipsoid is the image of the unit sphere under the transformation defined

by the tensor, D. The eigenvectors of D form an orthogonal basis, representing

the orientation of the ellipsoid. The length of each axis of the ellipsoid is the

corresponding eigenvector. For isotropic diffusion, the diffusion ellipsoid is a sphere.

2.1.2 Overview of MR Imaging

In this section we will present a brief overview of the MRI acquisition process.

A detailed treatment of the subject was done by Haacke et al. [23].

The protons in the nuclei of atoms align their axis of spin with the direction

of an applied magnetic field. The magnetic field also induces a wobble, known

as precession, in the spin of the protons. This frequency, the resonant frequency,

is proportional to the strength of the applied field. For protons the resonance

frequency lies in the RF range.

In the MRI machine, a field Bo is applied through out the imaging process.

The direction of this field defines the axial direction of the image. Protons will

absorb energy from an RF pulse of the resonance frequency, and tip away from the

direction induced by Bo. The amount of tip is proportional to the pulse duration.

The RF pulse also causes the protons to process in phase with each other. This

pulse is called the B1 field.

When the B1 transmitter is turned off, the absorbed energy at the resonant

frequency is re-emitted by the protons. This occurs as the spins, tipped by B1

return to their previous Bo alignment. The time constant associated with this









exponential process is known as the Ti relaxation time. The protons precessions

also dephase exponentially with time constant T2. The final image contrast is

influenced by strength, width and repetition time of the RF pulses in the B1 signal.

By spatially varying the intensities of Bo and B1, position information is

encoded. For instance, specially designed magnets add a gradient field to Bo. This

causes the proton resonance frequency to be a function of axial position. The

frequency of B1 can then be chosen to tip protons within a chosen slice.

To encode x, y position within a slice, two additional additional gradient fields

are employed. The first gradient, G, is pulsed, causing a phase variation, just as in

T2 relaxation. The phase variation is a function of position in the y direction. A

perpendicular gradient, Gx is then applied, changing resonance frequencies in the

x direction. A 2D Fourier transform reconstructs the image of each slice from the

data in the spatial frequency domain.

2.1.3 DT-MRI Acquisition

By carefully designing gradient pulse sequences, the measured signal from

protons in water molecules undergoing diffusion can be attenuated. The first

gradient pulse induces a known phase shift in proton precession. After some delay,

a second gradient pulse is applied, inducing the opposite phase shift. Protons which

have not moved between the two gradient pulses are returned to their previous

phase. Protons belonging to molecules which have changed location have some net

change in phase, changing their T2 relaxation time. Conventional MRI images of









white matter in the brain -,.-' -t a material of homogeneous composition. The

fibrous nature of white matter is apparent in DT-MRI however.

2.1.4 Stejskal-Tanner Equation


The intensity of the received signal, S, at each voxel depends on the properties

of the diffusion-encoding gradient, and the apparent diffusion tensor, D at that

location. The Stejskal-Tanner equation relates all of these quantities

3 3
In () D (2.1)
i=1 j=1

In 2.1, So is the intensity with no diffusion-encoding gradient present, and b is

a matrix characterizing the gradient pulse sequence. The ]iih\'-i -~ behind 2.1 is

beyond the scope of this thesis. The motivated reader may refer to the work of

Haacke et al. [23] for more details.

This imaging process must be performed with 7 noncoplanar gradient direc-

tions in order to fully generate a diffusion tensor image. Multiple samples, usually

3 or 4 are taken for each gradient direction.


In So

Dxx

In S1 1 -bl -bl -b' -2b' -2b' -2bJl D


= D (2.2)

ln S28 1 -b -b 2 -b 2 -2b28 -2b28 -2b28 D

Dy

Dxz









The overconstrained linear system, 2.1.4 is solved for So and the elements of the

symmetric tensor D by a least squares linear regression.

2.2 PDE Based Scalar-Valued Image Restoration


Image data smoothing or denoising is a fundamental problem in image

processing. Image denoising (noise removal) is a technique that enhances images

by attempting to reverse the effects of degradations occurring during acquisition

or transmission. Image noise makes it difficult to perform other processing tasks

such as edge-detection, segmentation, or in our case, fiber tracking. For this reason,

denoising is the first step in most image analyses.

The goal of image denoising is to recover an unknown, ideal image given some

observed image. In our case, we wish to recover the smooth S values from which

the DT image is calculated. The most common degradation source is the noise

from the image acquisition system and is commonly modeled by additive Gaussian

random noise. In the following, we will briefly review representative schemes,

specifically, partial differential equation (PDE) based methods that lend themselves

to fast numerical implementations.

Image denoising can be formulated using variational principles which in

turn require solutions to PDEs. Recently, there has been a flurry of activity on

the PDE-based smoothing schemes. For scalar valued image smoothing using

nonlinear diffusion filters with scalar diffusivity coefficient, we refer the reader

to the following articles and references therein [1, 11, 12, 15, 28, 31, 32, 34,

43, 45]. Anisotropic diffusion filters that use a tensorial diffusivity parameter









were introduced in Weickert [45]. These filters can be tailored to enhance image

structures (edges, parallel lines, curves etc.) that occur in preferred directions.

2.2.1 Linear Filtering

Linear filters are a simple and efficient means of removing noise from images.

One such filter may be implemented by the process of isotropic diffusion. This

diffusion process is governed by the heat equation


dI(:::, t)
(x,t) = V21(X, t)
at

I(x, 0) = Io(x) (2.3)


In the same way that a heated plate will seek an equilibrium state of a smooth

temperature gradient, so will an image evolved according to 2.3 smooth out

discontinuities in intensity. It can be shown that isotropic diffusion is equivalent to

convolving the image with a gaussian kernel. Taking the Fourier transform of 2.3

with respect to x, and defining S(I(x, t)) = U(w, t).

U(w,t) 2U(w,t)
&t

U(w, ) = Uo(w) (2.4)


The solution to 2.4 is


U(w, t) = Uo(w)e-w,2 t


(2.5)









The solution to 2.3 is then the inverse Fourier transform of 2.5

1 -4
I(x, t) = Jo e


2 = 2t (2.6)


Clearly, the linear diffusion process continues until a I(x, t) becomes a constant-

valued image. By 2.6 we can consider the effect of 2.3 as t -+ oo to equivalent to

convolving the original image with gaussian kernels of ever-increasing variance.

Since the gaussian kernel is separable, this type of filter is simple to implement

for images of arbitrary dimension. Although the images produced by this simple

filtering technique show a reduction in the high frequency noise, there is also

the unwelcome effect of blurred edges and lost details. The linear filter still has

applications in image i -.,uIiliIL. and generating scale-space image representations,

such as image pyramids.

2.2.2 Nonlinear Filtering

The nonlinear filters described in this chapter were designed to overcome the

inter-region blurring effect of linear filters. Nonlinear filters cannot be modelled

with gaussian convolution. The median filter is a nonlinear filter with a simple

implementation. A median filter replaces each intensity value in an image with the

median of the neighboring values. The size of this neighborhood determines the

amount of smoothing.









There are numerous models of nonlinear diffusion for image smoothing. To

implement the filters presented in the rest of this section we must numerically solve

a PDE.

2.2.2.1 Perona-Malik

An anisotropic diffusion process which inhibits blurring at edges was formu-

lated by Perona and Malik [34] by modifying 2.3.

0 =(, div(c(x)VI(x, t))


I(x, 0) = Io(x) (2.7)


When comparing 2.7 to 2.3 recall that V21 = V VI div(VI). By using a

diffusion coefficient which is a decreasing function of IVI(x, t)l, such as

1
c(x) = (2.8)
1 + VI(x, t)|

we can inhibit diffusion at locations of high image gradient, presumed to be edges.

By adding a reactive term to 2.7 we can actually enhance edges. The presence of

I(x, t) in 2.8 makes 2.7 nonlinear.

The diffusion coefficient serves only to slow diffusion at edges. The steady-

state (t -- oo) solution of 2.7 is still a constant valued image. A constraint can be

added to impose the condition that the smoothed image be "close" to the original

image. By introducing a reaction term to the diffusion equation we can impose a









data fidelity requirement

OI(x t)
=, t div(c(x)VI(x, t)) + y(lo(x) I(x, t))

I(x, 0) = 1o(X) (2.9)


we penalize solutions that differ from the initial image. The parameter /[ controls

the degree of smoothing in the final image. The steady-state solution of 2.9 is

nontrivial, so we no longer need a stopping time.

2.2.2.2 Tensor anisotropic diffusion

An alternative to merely slowing diffusion at edges is to align the direction

of diffusion to be parallel with edges in the image. Weickert [45] controls diffusion

direction by replacing the scalar c(x) in 2.7 with a diffusion tensor-valued function.

( =, t) div(D(x)VI(x, t))
Ot

I(x, 0) = Io(x) (2.10)


The function, D(x) produces tensors such that the unit eigenvector, el is parallel

with VI and e2 is perpendicular to VI. The eigenvectors A = g(|VII) and A2 = 1

are -,-. -fI .1 to generate anisotropic tensors whose representative ellipsoid has the

major axis parallel with image edges. In this way, diffusion near edges still occurs,

mostly along the edge. By limiting smoothing across the edge, edge location and

intensity can be preserved.









2.2.2.3 ALM

Alvarez, Lions, and Morel [1] propose a nonlinear parabolic PDE of the form

dI(x, t) VI(, t)
=g(VI)|VI div( )
Ot IVI(z, t)|

I(x, 0)= o() (2.11)


The effect is that I is smoothed in a direction orthogonal to the gradient at each

point. This is best considered in a level-set framework. Embedding I in IR3 by

considering t to be the third dimension, and 4(t = 0) to be the isocontours of I we

have

St) + V (x(t), t) '(t) = 0 (2.12)
&t

Recall the definitions of level-set normal ,N = --and curvature, = div( ).

The normal component of the velocity will be


V,, V .'(t) (2.13)


So 2.12 becomes

o((,t) lVl (2.14)


Letting v, be proportional to K we obtain 2.11. Curve-shortening flow of the

isocontours of 4 smoothes the image I.

2.2.3 Variational Formulation

The problem of image denoising is ill-posed, as are many inverse problems.

Since different pixel intensities could result in the same value when corrupted

by noise, the solution to the denoising problem is not unique. The technique of









Tikhonov regularization involves incorporating some additional knowledge about

the ideal image, similar to the prior distribution involved in a B.,\' -i.,l analysis. In

Tikhonov regularization, the solution space is restricted by posing the problem as

a minimization problem. The function (image) which minimizes some stabilizing

functional can be proven to be unique for appropriate choice of functional.

We will consider problems of the form


min(pE,(I) + Ed(I)) (2.15)
I

The functionals, E involved in image smoothing often correspond to ]11i\--i. ,i

models, and usually represent energy. For instance, the energy associated with data

fidelity is

Ed = klIo 2 (2.16)


A ]11\--1i. analogy for the data constraint is a number of springs between the

initial image, Io and the ideal image, I.

The second energy term, E,, represents the internal potential energy of the

surface described by the image. The regularization parameter, [t, is introduced to

control the amount of smoothing.

2.2.3.1 Membrane spline

A first-order stabilizer, stretching energy (arc-length), is one possible stabilizer


E,= Ij + + I ldxdy (2.17)
J JR2









The Euler-Lagrange condition for the minimization of 2.17 is V21 = 0, the steady-

state solution to 2.3. In fact, most of the filters considered so far can be obtained

from a variational principle. By setting E, = Ems in 2.15 we obtain a membrane

spline solution. In this case, p/ controls the "t. i,-i. ii of the spline surface.

2.2.3.2 Thin-Plate spline

Using a second-order stabilizer bending energy (curvature)


Et,, =J I + 2I + I2 dx dy (2.18)
R2

as a stabilizer we obtain a thin-plate spline solution. For this spline, p/ controls the

"stiffness" of the spline surface. The Euler-Lagrange condition for the minimization

of 2.18 is V41 = 0.

2.2.4 Total Variation Image Restoration

Another side-constraint that can be used as an energy functional is called

the total variation (TV) norm. This norm represents oscillation. In our case,

image noise is considered to be "wrinkling" of the surface described by the image.

Minimizing the TV norm produces very smooth images while permitting sharp

discontinuities between regions [40]. The TV norm is formulated by 2.19.


TV,,i(I(x))= VI(x) dx, Q C T (2.19)


The TV norm is essentially the L1 norm of the image gradient.

At discontinuities, the weak derivative DI(x) to calculate the TV norm. For

a piecewise continuous function, the TV norm is then the sum of the TV norm










of each continuous piece plus the sum of the absolute values of the "jumps". For

example, functions fl and f2 in 2.2 have the same TV norm. The oscillatory

function, fl, has a higher TV than the function with a discontinuity, f3. Since

f2" If3



S/ J2






a /




Figure 2.2: TV(f 1) > TV(f 2) = TV(f 3)


most images consist of piecewise smooth regions separated by discontinuities

(edges), this is a useful model for image denoising.

The Euler-Lagrange condition for the minimization of 2.19 written in




dI(x, t) VI(x, t)
//"
























Oit VI(x, t)
/ f3


I--
a h










Figure 2.2: TV() > TV(.202) TV(3)











An example of image denoising by TV norm minimization is shown in 2.3.
The rest images consist on the right ise smooth regions separated by a high degree of intra-region

(edges), this is a useful model for image denoising.
The Euler-Lagrange condition for the minimization of 2.19 written in

gradient-descent form is


I(x,t) d VI(,t)
dt Iv(I,t))

I(x,0) =Io(x) (2.20)


An example of image denoising by TV norm minimization is shown in 2.3.

The restored image on the right is characterized by a high degree of intra-region

smoothing, and edges have also been preserved.
























Figure 2.3: Noisy image (left) and restored (right) by TV norm minimization.


2.3 PDE Based Vector-Valued Image Restoration

The main application of vector-valued image restoration has been the restora-

tion of color images. The simplest implementation of vector-valued smoothing

is uncoupled smoothing. By treating each component of the vector field as an

independent scalar field, we can proceed by smoothing channel-by-channel using

one of the methods presented in section 2.2 for scalar image smoothing. This can,

however, result in a loss of correlation between channels as edges in each channel

may move independently due to diffusion. To prevent this, there must be coupling

between the channels.

There are many other PDE-based image smoothing techniques which we will

not cover here, but will refer the reader to Weickert [45] and Caselles et al.[10].









2.3.1 Vector-Valued Diffusion

Whitaker and Gerig introduced anisotropic vector-valued diffusion which was a

direct extension of the work of Perona and Malik [34]. The equations


=, t) div(c(I)Vl,(x, t))
Ot

In(x,0) = I,o(x) (2.21)


are coupled through the function C, and can be written in vector form as

dI(z,t)
( =, V (c(I)VI(x, t))
Ot

I(x, 0) =Io(z) (2.22)


2.3.2 Riemannian Metric Based Anisotropic Diffusion

Sapiro et al. [41] introduced a selective smoothing technique where the

selection term is not simply based on the gradient of the vector valued image but

based on the eigen values of the Riemannian metric of the underlying manifold. In

2 dimensions, the underlying manifold can be thought of as the parametric surface

described by the image I(x). By geometrically smoothing this surface, we also

smooth the vector-valued image.

The entries of the metric tensor, G, are defined by

0I 0I
gij = (2.23)


The largest and smallest eigenvalues of G, (A+, A_) and their corresponding

eigenvectors ( +,(_) describe surface properties of the Riemannian manifold. The









degree and direction of maximal change are given by A+ and (+. Smoothing is

achieved by evolution in the direction of minimal change, that is, along edges in the

vector-valued image

I(x, t) g 2I(z, t)
Oi-t g(A+, A-)

I(x, 0) =Io(z) (2.24)


2.3.3 Beltrami Flow

More recently, Kimmel et al. [25] presented a very general flow called the

Beltrami flow as a general framework for image smoothing and show that most

flow-based smoothing schemes may be viewed as special cases in their framework.




OIj (, t) 1 2 0 ( G 1,


I(x, 0) =Io(z) (2.25)


where G is the metric tensor, and IGI is the determinant of the metric tensor.

The Beltrami flow can be thought of as nonlinear diffusion on a manifold,

where IGI acts as an edge-detecting diffusion coefficient.

2.3.4 Color Total Variation

Blomgren and Chan [8] introduced the TVn,m norm for vector valued images.

mTV Ix)) ] (.6
TVnm(I(x)) = [TV1,i(I()]2 (2.26)
i -1







24

For m = 1, 2.26 reduces to the scalar TV norm 2.19. The Euler-Lagrange condition

for the minimization of 2.26, written in gradient descent form is

01i (x,,t) _TV ,I (1,) V -],i
Ot TV,,m(1) VI

I(x, 0) = Io(x) (2.27)


This was shown to be quite effective for color images, preserving edges in the

color space while attenuating noise. However, for much larger dimensional data

sets (m=7) as in the work proposed here, the Color TV method becomes com-

putationally very intensive and thus may not be the preferred method in such

applications.















CHAPTER 3
DT-MRI IMAGE RESTORATION


3.1 Noise Model


The degradation associated with the measurement of the S values is modelled

by an additive gaussian process. Let S(X) be the vector valued image that we

want to smooth where, X = (x, y, z) and let S(X) be the unknown smooth

approximation of the data that we want to estimate. We have S(X) = S(X) + p.

Although we can consider the noise to be on the components of the tensor, D, the

form of this noise is no longer so simple. In fact, even computing D from S, using

the Stejskal-Tanner relation may be meaningless. Substituting the noise model into

2.1
3 3
lnS+ = lnSo -b,jD,j (3.1)
i=1 j=1
The ]1il-. ,1l principles governing diffusion do not allow for negative values of S, as

evidenced by the dependence of 3.1 on the natural logarithms of S and So.

Low-intensity S measurements may be overwhelmed by noise. However,

since p is a zero-mean distribution, these measurements may become negative.

Clearly, we cannot calculate D in this case. Nor should we propagate the noise

model through the Skeskal-Tanner equation in order to write D as a function of

pq. We propose to instead smooth the vector of S measurements before any further

analysis.









3.2 Restoration Formulation


Smoothing the raw vector valued image data is posed as a variational principle

involving a first order smoothness constraint on the solution to the smoothing

problem. We propose a weighted TV-norm minimization for smoothing the vector

valued image S.

The variational principle for estimating a smooth S(X) is given by

7 7
minS(S) = [g(A+, A_) VSi + Si il2]dx (3.2)
s 2n 2
i=1 ii

where, 2 is the image domain and p is a regularization factor. The first term

here is the regularization constraint on the solution to have a certain degree of

smoothness. The second term in the variational principle makes the solution

faithful to the data to a certain degree. The weighting term in this case g(s) =

1/(1 + s) where s = FA is the fractional anisotropy defined by Basser et al. [2].

This selection criteria preserves the dominant anisotropic direction while smoothing

the rest of the data. Note that since we are only interested in the fiber tracts

corresponding to the streamlines of the dominant anisotropic direction, it is apt to

choose such a selective term.

Here we have used a different TV norm than the one used by Blomgren and

Chan [8]. The TV,,, norm is an La norm of the vector of TV,i1 norms for each

channel. We use the L1 norm.









The gradient descent of the above minimization is given by

OS, g(A+,A_)VSi.
= div VS -/[ (Si S) i = 1, ..., 7
Ot ||VS ||ll

OS
-On\9x+ = 0 and S(x, t = 0) = S(x) (3.3)


The derivation of equation 3.3 from 3.2 is presented in section 8.2.5.

The use of a modified TV norm in 3.2 results in a looser coupling between

channels than the use of the true TVn,m norm would have. This reduces the

numerical complexity of 3.3 and makes solution for large data set feasible.

Note that the TVn,m norm appears in the gradient descent solution 2.27

of the minimization problem. Consider that our data sets consist of 7 images,

corresponding to gradient directions. Each of these images consists of several

samples (usually 3 or 4) corresponding to different gradient strengths. Calculating

the TV,,m norm by numerically integrating over the 3-dimensional data set at each

step of an iterative process would have been prohibitively expensive.















CHAPTER 4
NEURONAL FIBER TRACKING


4.1 Overview of Neuronal Fiber Tracking

Water in the brain preferentially diffuses along white matter fibers. By

tracking the direction of fastest diffusion, as measured by MRI, non-invasive fiber

tracking of the brain can be accomplished. Fibers tracks maybe constructed by

repeatedly stepping in the direction of fastest diffusion. The direction along which

the diffusion is dominant corresponds to the direction of eigen vector corresponding

to the largest eigenvalue of the tensor D. In Conturo et al., [18], fiber tracks

were constructed by following the dominant eigenvector in 0.5 mm steps until a

predefined measure of anisotropy fell below some threshold. This usually occurred

in grey matter. The tensor, D, was calculated at each step from interpolated

DT-MRI data. This tracking scheme is primarily based on heuristics and is not

grounded in well founded mathematical principles.

Mori et al. [30] achieved fiber tracking by using several heuristics. The

tracking algorithm starts from a voxel center and proceeds in the direction of

the major axis of the diffusion ellipsoid. When the edge of the voxel is reached,

the direction is changed to that of the neighboring voxel. Tracking stops when a

measure of adjacent fiber alignment crosses a given threshold. The scheme however

is resolution dependent since the MRI data only reflects average axonal orientation









within a voxel. Small fibers adjacent to each other may not be distinguished.

Another problem occurs with branching fibers. This method will only track one

path. In this situation, multiple points within a bundle may be independently

tracked. Westin et al. [46] reported techniques for processing DT-MRI data using

tensor averaging. Diffusion tensor averaging is an interesting concept but does

not address the issue of estimating a smooth tensor from the given noisy vector-

valued data. More recently, Poupon et al. [38] developed a B.\. -i,11 formulation

of the fiber tract mapping problem. Prior to mapping the fibers, they use robust

regression to estimate the diffusion tensor from the raw vector valued image data.

No image selective smoothing is performed in their work prior to application

of the robust regression for estimating the diffusion tensors. Coulon et al. [20]

determined fibers tracts after smoothing the eigen values and vectors. Once again,

this scheme is faced with the same problem i.e., the eigen vector and the eigen

values are computed from a noisy tensor field and hence may not be meaningful at

several locations in the field. Parker et al. [33] also presented a follow up article on

fiber tract mapping wherein, they use the idea of the fast marching method from

Sethian et al. [42] for growing seeds initialized in the smoothed dominant eigen

vector field. Batchelor et al. [7] reported an interesting fiber tract mapping scheme

wherein, they produce a map indicating the probability of a fiber passing through

each location in the field. However, they do not address the issue of computing the

diffusion tensor from a noisy set of vector valued images. This is a very important

issue and should not be overlooked as is demonstrated in the preliminary results of

this proposal.









Given the dominant eigen vector field of the diffusion tensor in 3D, tracking

the fibers (space curves) tangential to this vector field is equivalent to finding the

stream lines/integral curves of this vector field. Finding integral curves of vector

fields is a well researched problem in the field of Fluid Mechanics [17]. The simplest

solution would be to numerically integrate the given vector field using a stable

numerical integration scheme such as a fourth order Runge-Kutta integrator [39].

However, this may not yield a regularized integral curve.

In the context of the fiber tract mapping, two sub-tasks are involved namely,

(1) estimating a denoised diffusion tensor from noisy vector-valued image mea-

surements and (2) estimating the streamlines of the dominant eigen vectors of the

diffusion tensor. The denoising involves a weighted TV-norm minimization of the

vector valued data S followed by a linear least squares estimation of D from the

smoothed S and then an estimation of regularized stream lines of the dominant

eigen vectors of D.

4.2 Formulation

Diffusion at each point can be characterized by a diffusion ellipsoid. The

ellipsoid axes are the eigenvectors of the diffusion tensor. The radii are the

corresponding eigenvalues. The shape of the ellipsoid reflects the isotropy of

diffusion. Nearly spherical diffusion ellipsoid represent areas of free water, where

diffusion is unimpeded. Areas of white-matter fiber bundles have elongated

ellipsoids, since water diffusion is restricted in directions perpendicular to fiber









direction. The phenomenon was quantified by Basser and Pierpaoli [2].


FA 3 /(A )2 (A2 -X2 (A3 -) (4.1)
2 Al + A 2+A A3+

For the stream line estimation problem, we pose the problem in a variational

framework incorporating smoothness constraints which regularize the stream

lines/integral curve. The variational principle formulation leads to a PDE which

can be solved using efficient numerical techniques.

The variational principle


minE(p) = min c'(p) + -c'(p) v(c(p))2dp (4.2)
c c J

where, c(p) = (x(p), y(p), z(p))T is the integral curve we want to estimate and

p E [0, 1] is the parameterization of the curve, v(c(p)) is the vector field v restricted

to the curve c(p). The first term of 4.2 is a smoothness constraint. By minimizing

arc-length we penalize spurious oscillation in the curve. The second term provides

data fidelity : we wish for the fiber to be tangent to the vector field at every

point along the curve. The parameter, 3 controls the degree of smoothness of the

solution.

The gradient descent of (4.2)

Dc
S= c'(p)kn + 3[c"(p) V(c(p))c'(p)] +

P3V(c(p))(c'(p) v(c(p))) (4.3)







32

where k is the curvature of the space curve, /3 is a regularization parameter and

VT = The transpose of V


Vlx Vly Vlz

V =Jv= (4.4)


V 3x C'

The initial condition, c(p, t = 0) is provided by an ordinary streamline integration

routine.














CHAPTER 5
NUMERICAL METHODS


5.1 Data Denoising

The nonlinear PDE for denoising the raw data is

Vs
div(g ) (s so) = 0 (5.1)
|Vs

where g : R3 -- IR is the selective smoothing and p is the constant regularization

parameter.

5.1.1 Fixed-Point Lagged-Diffusivity

Equation 5.1 is nonlinear due to the presence of |Vs in the denominator of the

first term. We linearize 5.1 by using the method of "lagged-diffusivity" presented

by Chan and Mulet [13]. By considering |Vs| to be a constant for each iteration,

and using the value from the previous iteration we can instead solve

1
(Vg Vst + gV2st+1) + (st+l sO) = 0 (5.2)


Here the superscript denotes iteration number. First, rewrite 5.2 with all of the st

terms on the right-hand side


V2t+l +P 1 = V + Vg V (5.3)
g g









5.1.2 Discretized Equations


To write 5.3 as a linear system (As + 1 = y), discretize the laplacian and

gradient terms. Using central differences for the laplacian we have

V2t+1 t+1 + t+1 .+1 f+l 1 t+1 t+1 t+1l
V2 +1 += S l,yz 1 xl,z 51,z_1 6 z x+l,y,z x,y+l,z + 5x,y,z+l (5.4)


Define the standard forward, backward and central differences to be


aS = Sx+1,y,z Sx,y,z

Ss = sx,y+l,z Sx,y,z

A+ = Sx,y,z+l Sxy,tz

A S = Sx,y,z SX-1,yz

A 's = Sx,y,z Sx,y-l,z

A S = Sxy,z -- Sx,y,z-1

aS = S(+1,,yz Sw-1,y,z)

AyS = I(Sx,y+l,z Sx-1,y,z)

1
Azs = (SX,y,z+l Sx-1,y,z) (5.5)


We can rewrite 5.3 in discrete form using the definitions in 5.5


-S_-1,y,z Sx,y-l,z SX,y,z-1 + (6 + p )SX, Y,

-Sx+1,y,z S,y+ l,z Sx,y,z+l

= l(s (A( )2 + (A0^)2 + (AS)2 + AgA2S + A yAS s + A ) (5.6)
9 is 56










This results in a banded-diagonal linear system with 7 nonzero coefficients per row.


6+ 1 -1 ... -1 ... ... ... s f

1 1 ... 1 ... t+1 f
91 (5.7)

0 -1 6+w3 1 ... ... 1 ...
93
+ t




where the right-hand side of 5.6 has been replaced with f$. The matrix in equation

5.7 is symmetric and diagonally dominant. We have successfully used conjugate

gradient descent to solve this system.

The solution of 5.7 represents one fixed-point iteration. This iteration is

continued until Ist sJ+1 < c, where c is a small constant.

5.2 Fiber Regularization


The PDE for fiber regularization is

dc
-= c'(p)lkn + p[c"(p) V(c(p))c'(p)] +
Ot

pVT(c(p))(c'(p) v(c(p))) (5.8)


where k is the curvature of the space curve, p is a regularization parameter and

VT = The transpose of V


VUlx Vy Ulz

V =Jv= ,. ,. (5.9)


UV3x V3z









Using the definitions of tangent, T, curvature, k, and normal, n

c' IT'I T'
T = k =c ,n =T' (5.10)
|C'| k IC'I T '|

we have

Oc c'
( C' + p(c" + VT(c v) VI') (5.11)
Ot |c'|


5.2.1 Lagged-Diffusivity

We linearize 5.11 using the concept of lagged-diffusivity, as we did in the data

denoising case.

c (t+At Ct+At(
( ) = ( ) (5.12)
IC/t

We can simplify notation by calling the coefficients of c" and c+', respectively

at(p) and Lh(p) and calling the additive constant Kt.

Oc
= at(p)c"+At + (p)ct+At + Kt (5.13)


5.2.2 Crank-Nicholson Method

We solve 5.13 using the Crank-Nicholson method. This method achieves

stability by using averaged differences to estimate derivatives. For the second

derivative term we use


2 (ct(p + Ap) 2ct(p) + ct(p Ap) +

ct+At(p + Ap) 2ct+t(p) + ct+At(p A))


(5.14)









The averaged difference for the first derivative is

ac 1
= 2 (ct (p + Ap) c(p Ap) + ct+At (p + Ap) ct+At (p Ap)) (5.15)

In finite-difference form

cAT c) p J(ct+At(p + Ap) 2ct+t(p) + ct+t(p Ap)) +


2(ct(p + Ap) ct(p Ap) +

2p(c (p + Ap) 2ct(p) + ct(p Ap)) +

S(ct+At (p + Ap) ct+At (p Ap) + Kt(p) (5.16)


Equation 5.16 can be expressed as a linear system to be solved for ct+t.


Atct+At = Mtc + K, (5.17)


This tridiagonal linear system can be solved by Grout's factorization, an optimized

LU factorization for this type of matrix.















CHAPTER 6
NEURONAL FIBER VISUALIZATION


6.1 Rendered Ellipsoids

A very simple visualization strategy is to simply render the diffusion ellipsoid

at a subset of data points. Since a 3D field of ellipsoids would occlude each other,

this visualization is usually done for 2D slices of data. Additionally, only ellipsoids

on a sparse grid can be rendered in order for each ellipsoid to be discerned. This

type of visualization can easily become visually cluttered and convey so little

information as to be useless. Laidlaw[26] however, successfully incorporated

ellipsoids in a layered visualization approach.

6.2 Streamtubes

Streamtubes are a three-dimensional analogue to streamlines. In fact, the

streamtube is computed by using a streamline as the centerline of the tube. We

can use the streamtube diameter to encode some additional information about the

tensor field being visualized, such as FA value. Previously, Laidlaw [27] has applied

the streamtube visualization approach to DT-MRI.

6.3 Line Integral Convolution

It is also possible to visualize the 3D vector field corresponding to the domi-

nant eigenvalues of the diffusion tensor using other visualization methods such as

the line integral convolution technique introduced by Cabral et al. [9] -a concept







39

explored in this work as well. The advantage of this visualization technique is that

it is well suited for visualizing high density vector fields and does not depend on

the resolution of the vector field moreover, it also has the advantage of being able

to deal with branching structures that cause singularities in the vector field. Since













Figure 6.1: LIC visualization of synthetic field.


the fiber direction is parallel to the dominant eigenvector of the diffusion tensor, we

can calculate fiber paths as integral curves of the dominant eigenvector field. The

stopping criterion is based on FA value. When FA falls below 0.17 we consider the

diffusion to be nearly isotropic and stop tracking the fiber at this point.

Once the diffusion tensor has been robustly estimated, the principal diffusion

direction can be calculated by finding the eigen vector corresponding to the

dominant eigen value of this tensor. The fiber tracts may be mapped by visualizing

the streamlines through the field of eigen vectors.

LIC is a texture-based vector field visualization method. The technique

generates intensity values by convolving a noise texture with a curvilinear kernel









aligned with the streamline through each pixel, such as by


I(xo) = T(a(s))k(so s)ds (6.1)
Jso-L

where I(xo) is the intensity of the LIC texture at pixel xo, k is a filter kernel of

width 2L, T is the input noise texture, and a is the streamline through point Xo.

The streamline, a can be found by numerical integration, given the discrete field of

eigen vectors.

The result is a texture with highly correlated values between nearby pixels on

the same streamline, and contrasting values for pixels not sharing a streamline. In

our case, an FA value below a certain threshold can be a stopping criterion for the

integration since the diffusion field ceases to have a principal direction for low FA

values.

Stalling and Hege [44] achieve significant computational savings by leveraging

the correlation between adjacent points on the same streamline. For a constant

valued kernel, k, the intensity value at I(r(s + ds)) can be quickly estimated by

I(o(s)) + e, where e is a small error term which can be quickly computed.

Previously, Chaing et al. [16] have used LIC to visualize fibers from diffusion

tensor images of the myocardium.

6.4 Particles


The LIC and streamtube techniques presented in the previous sections are

time-consuming operation. All of the streamline are completely traced before

an image can be displayed. For interactive display of fibers, we use a particle









based visualization technique. The particles are analogous to smoke introduced

into a wind-tunnel to visualize streamlines. Rather than simulating the actual

diffusion process, the particles simply advect through a velocity field described by

the dominant eigenvector of the diffusion tensor at each point. By seeding a few

streamlines within a region of interest, and performing a single step of numerical

integration at a time, interactive frame-rates can be achieved.

Each particle is implemented as a small textured quadrilateral which is always

oriented to face the viewer. We vary the size and color of this quad as a means of

visualizing the FA value at the particle's location in the tensor field.

We adapt this technique to tensor field visualization by incorporating the

FA value at each field location in to the LIC texture. By modulating the image

intensity with an increasing function of FA, we highlight the areas of white matter,

and are able to resolve where the fibers eventually track into grey matter.

Unlike LIC and streamtube tracing, the particle visualization requires no

preprocessing time. The other techniques require completely integrating many,

perhaps thousands, of streamlines. Particle-based techniques allow immediate

visualization.















CHAPTER 7
EXPERIMENTAL RESULTS


7.1 Data Denoising

In all of the experiments, we first smooth the seven 3D directional images

using the novel selective smoothing technique outlined in section 3. Following this,

the diffusion tensor is estimated from the smoothed data using a standard least

squares technique. The results of FA calculation from the smoothed data, and from

raw data are presented in Figures 7.1 and 7.2












Figure 7.1: FA image of coronal slice of raw rat brain data.













Figure 7.2: FA results for smoothed data.









7.2 Fiber Tracking

7.2.1 LIC

Figures 7.3 and 7.4 depict the computed fiber tracts for the reconstructed rat

brain data. The intensity of the LIC texture has be modulated with the FA image

to emphasize the most anisotropic region of each image.


Figure 7.3: LIC fiber tracts in coronal slice of smoothed rat brain data.


Figure 7.4: LIC fiber tracts in axial slice of smoothed rat brain data.









7.2.2 Streamtubes

We will now compare our computed streamtubes with a fluoroscopic image.

Fibers are evident in in fluoroscopic images as high intensity treelike structures.

The fluoro image shown on the left of 7.5 shows a known anatomical feature, the

fiber crossings in the corticospinal tract at the base of the brain. The streamtube

map from the same region is shown on the right of 7.5. The streamtube map was
















Figure 7.5: Comparison of fluoroscopic image (left) with extracted streamtubes
(right).


generated by starting streamline integration at each point of a sparse grid within

the data if FA > 0.3. Tracking stopped at FA < 0.17. The streamtube radius is a

function of FA.

7.2.3 Particles

Figure 7.8 shows a fiber tracing sequence obtained by seeding fiber in the

corpus callosum region of the rat brain data. The brighter particles signify high FA

value. Dimmer particles can be seen tracing into grey matter.










































Figure 7.6: Axial view of streamtubes in corpus callosum.

























































Figure 7.7: Details of corpus callosum.










































Figure 7.8: Sequence from Fiber visualization (2 seconds between images).















CHAPTER 8
CONCLUSIONS AND FUTURE WORK


8.1 Conclusions

In this paper, we presented a new weighted TV-norm minimization formulation

for smoothing vector-valued data specifically tuned to computation of smooth

diffusion tensor MR images. The smoothed vector valued data was then used to

compute a diffusion tensor image using standard least squares technique. Fiber

tracts were estimated using the dominant eigen vector field obtained from the

diffusion tensor image. Finally, results of fiber tract mapping of a rat spinal cord,

and a rat brain were depicted using LIC, streamtube and particles. The fiber

tracts are quite accurate when inspected visually, and correspond well with known

anatomical structures, such as the corpus callosum.

8.2 Future Work


There are several areas of improvement and directions for future work involv-

ing this research.

8.2.1 Noise Model For Measured Data

The S vector that was smoothed in this work is not the actual measured data.

These S values are actually the magnitude of Fourier transformed data. A noise

model expressed in terms of the voltage levels observed during the imaging process

would be fundamentally correct and may lead to more accurate fiber maps.









8.2.2 Quantitative Validation

However, quantitative validation of the computed fiber tracts is essential.

This can be achieved by registration of our three-dimensional fibers with two-

dimensional fluoroscopy images. This will be the focus of our future efforts.

8.2.3 Robust Regression

The time intensive nature of DT-MRI data collection limits the number

of samples taken in each direction. Outliers in the collected data may have a

significant impact on the least-squares calculation of the diffusion tensor. The use

of a robust regression for calculating the diffusion tensor should be explored.

8.2.4 Non-Tensor Model of Diffusion

An alternative to the tensor model of diffusion is to characterize diffusion with

a spherical distribution of diffusion magnitudes. This requires measurements to be

made in many more directions, but has the advantage of handling fiber crossings

occurring within a single voxel.

8.2.5 Automatic Region-Of-Interest Extraction

By registering a new data set with an atlas, the user could select the region-

of-interest by anatomical name. Currently the region-of-interest for seeding fiber

tracts is done by hand. This method is time consuming and potentially inaccurate.














APPENDIX
DERIVATION OF EULER-LAGRANGE CONDITIONS


Here we will review the calculus of variations, and present the derivation the

Euler-Lagrange equations for the image denoising and fiber regularization PDEs.

Evans [22] provides a thorough treatment of variational calculus.

In variational calculus we seek to obtain minima or maxima of expressions

that depend on functions instead of parameters. The expressions are referred to as

functionals. Consider a general functional, I depending on function f(x, y, z)


I[f] = J J F(f, ff, f, x,y,z)dxdydz (1)


We seek the function f which minimizes the expression I.

Suppose the function f minimizes the functional I. Consider constructing a

variation of f, using parameter e and test function pr. This variation is given by

f + er/. Consider a new function of one variable, i(e) = I[f + erq]. Extrema of this

function occur when
di
= o (2)

Since we know that f minimizes I, we can say


o = 0 (3)


is a condition for the minimization of I.






51

As an example, consider the variation of 1


i(e) = I[f + ] F(f + cq, f + crq, ff + er, f + erz, x, y, z)dxdydz (4)

Differentiate with respect to e


i'(0) = Jf


OF
O(f + Er)


OF OF OF
(f + 'x) x + qf + ,y) + + /zf~ r dx dy dz
9(f, + e) 9(fy + eJq) + (f, + eo)


Evaluate this derivative at e = 0


i'(0) =


OF
Of F


OF
+ Of'
6 fm


OF OF
+ jrr + rz dx dy dz
aj' (Ji


Integrating by parts, and imposing the condition that the test function, tr, equals 0

on the boundary we get the result


OF d OF
Of dxz Of


d (OF
dy Ofy


OF
Of)


This is the condition for minimization of the functional in 1 When second deriva-

tives are involved, such as the functional


I[f = F(f, f,, fx, x,y, z)dxdydz


Then the integration by parts step must be performed twice in order to write i'(0)

in terms of tr. This leads to the Euler-Lagrange condition


OF d OF
Of dx Of,


d2 OF
dX2 Of) 0








Data denoising. We can rewrite 3.2 as


min fg s + S S+ s + (s


so)2dx dy dz


Using the general result derived in 7, the Euler-Lagrange condition for minimizing

10 is


d sx
(s so) (g
dx |Vs|


d sy
dy SVs
dy Ws


d g
dz Vs|


(11)


This can be rewritten in vector notation as


Vs
y(s so) div(g)Vs
IVS


Fiber regularization. The first variation of 4.2 with variation e and test

function r is


I[c + eC] c = Ic'(p) + + 1c'(p) + Fe'
S

Differentiating with respect to e gives


v(c + eq(p))12dp


(12)


(13)


OI[c + erq]
Oe


c'(p) + eW'1


The Euler-Lagrange condition is given by


Ic'(p) I + 0(c'(p)
Jo c'(p)i


Vv(c(p)). )


0 (15)


Integrating by parts we obtain


c'(p) '
C + P[c"(p)
Ic'(p)l


V(c(p))c'(p)] + PVT(c(p))(c'(p)


(10)


v(c + Ce))(7'


OI 0
e 1 0


Vv(c + cq). ) (14)


v(c(p)))


(16)







53

where V is the 1. .,1L,; of v, (d flu, 1 as


Vix Vl-y Vz

(17)
V = J= v2; U2z


3I 3















REFERENCES


[1] L. Alvarez, P. L. Lions, and J. M. Morel, "Image selective smoothing and edge
detection by nonlinear diffusion. ii," SIAM J. Numer. Anal., vol. 29, no. 3, pp.
845-866, June 1992.

[2] P. J. Basser and C. Pierpaoli ili. rostructural and ]pl\ --i.l. -i ., features of
tissue elucidated by quantitative-diffusion-tensor MRI," J. Magn. Reson. B
vol. 110, 209-219, 1996.

[3] P. J. Basser, J. Mattiello and D. LeBihan "MR diffusion tensor spectroscopy
and imaging" B;i..I,,, J. vol. 66, 259-267, 1994a.

[4] P. J. Basser, J. Mattiello and D. LeBihan "Estimate of the effective self-
diffusion tensor from the NMR spin echo," J. Magn. Reson. B vol. 103,
247-254, 1994b.

[5] P. J. Basser "Inferring microstructural features and the ]1i\i, ,-1 i.i, .di state of
tissues from diffusion weighted images" NMR in Biomed. vol. 8, 333-344, 1995.

[6] P. J. Basser, S. Pajevic, C. Pierpaoli, J. Duda and A. Aldroubi, "In vivo fiber
tractography using DT-MRI data," Magnetic Resonance in Medicine, vol. 44,
2000, pp. 625-632.

[7] P. G. Batchelor, D. L. G. Hill, F. Calamante and D. Atkinson, "Study
of connectivity in the brain using the full diffusion tensor from MRI," in
Proceedings of the 17th Intl. Conf. on Information Processing in Medical
Imaging, 2001, Davis, CA, pp. 121-133.

[8] P. Blomgren and T. F. Chan,"Color TV: total variation methods for restoration
of vector-valued images," IEEE Transaction on Image Processing, vol. 7, no. 3,
pp. 304-309, March 1998.

[9] B. Cabral, "Imaging vector fields using line integral convolution," Computer
Graphics Proceedings, pp. 263-270, 1993.









[10] V. Caselles, J. M. Morel, G. Sapiro and A. Tannenbaum,IEEE Trans. on
Image Processing, special issue on PDEs and geometry-driven diffusion in
image processing and analysis, vol. 7, No. 3, 1998.

[11] F. Catte, P. L. Lions, J. M. Morel, and T. Coll, "Image selective smoothing
and edge detection by nonlinear diffusion," SIAM Journal on Numerical
A.l.i,li,, vol. 29, pp. 182-193, 1992.

[12] P. Charbonnier, L. Blanc-Feraud, G. Aubert, and M. Barlaud, "Two deter-
ministic half-quadratic regularization algorithms for computed imaging,," in
in Proc. of the IEEE Intl. Conf. on Image Processing (ICIP), 1994, vol. 2, pp.
168-172, IEEE Computer Society Press.

[13] T. Chan and P. Mulet, "On the convergence of the lagged diffusivity fixed
point method in total variation image restoration," SIAM Journal on Numeri-
cal Aira,i1,. vol. 36, no. 2, pp. 354-367, 1999.

[14] C. Chefd'hotel, D. Tschumperl, R. Deriche, O. Faugeras, "Constrained flows
of matrix-valued functions : Application to diffusion tensor regularization,"
in Proc. of the European Conference on Computer Vision 2002 (ECCV'02),
Copenhagen, Denmark, May 2002, pp. 251-265.

[15] Y. Chen, B. C. Vemuri and L. Wang, "Image denoising and segmentation via
nonlinear diffusion," Computers and Mathematics with Applications, vol. 39,
No. 5/6, pp. 131-149, 2000.

[16] P.-J. Chiang, B. Davis and E. Hsu, "Line-integral convolution reconstruction
of tissue fiber architecture obtained by MR diffusion tensor imaging," BMES
Annual Meeting, 2000.

[17] A. Chorin, Computational Fluid Mechanics, Selected papers, Academic Press,
NY, 1989.

[18] T. E. Conturo, R. C. McKinstry, E. Akbudak, and B. H. Robinson "Encoding
of anisotropic diffusion with tetrahedral gradients: A general mathematical
diffusion formalism and experimental i. -tlt- Magn. Reson. Med. vol. 35,
399-412, 1996

[19] 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 '.,In Proc. Natl. Acad. Sci. USA vol. 96,
10422-10427, 1999.









[20] O. Coulon, D. C. Alexander and S. R. Arridge, "A regularization scheme for
diffusion tensor magnetic resonance images," in Pii'... i;,.1'i of the 17th Intl.
Conf. on Information Processing in Medical Imaging, Davis, CA, pp. 92-105,
2001.

[21] F. Crick and E. Jones "Backwardness of human neuroanatomy," Nature vol.
361, 109-110, 1993.

[22] L.C.Evans, Partial Differential Equations, American Mathematical Society,
Providence RI, 1998.

[23] E. Haacke, R. Brown, M. Thompson, R. Venkatesan. Magnetic Resonance
Imaging-Phi,,~ i,1 Principals and Sequence Design, John Wiley and Sons, NY,
1999.

[24] D. K. Jones, A. Simmons, S. C. R. Williams and M. A. Horsfield, "Non-
invasive assessment of axonal fiber connectivity in the human brain via
diffusion tensor MRI" Magn. Reson. Med., vol. 42, 37-41, 1999.

[25] R. Kimmel, N. Sochen, and R. Malladi, "Images as embedding maps and
minimal surfaces:movies, color and volumetric medical images," in Proc. of the
IEEE Conf. on Computer Vision and Pattern Recognition, pp. 350-355, June
1997.

[26] D. H. Laidlaw, E. T. Ahrens, D. Kremers, M. J. Avalos, "Mouse spinal
cord diffusion tensor visualization using concepts from painting," Siggraph
Technical Slide Set, 1998.

[27] D. H. Laidlaw, S. Zhang, C. Curry, D. Morris, "Streamtubes and streamsur-
faces for visualizing diffusion tensor MRI volume images," in Proc. of the IEEE
Visualization, October 2000.

[28] R. Malladi and J. A. Sethian, "A unified approach to noise removal, image
enhancement and shape recovery," IEEE Trans. on Image Processing, vol. 5,
no. 11, pp. 1554-1568, 1996.

[29] N. Makris, A. J. Worth, A. G. Sorensen, G. M. Papadimitriou, O. Wu, T. G.
Reese, V. J. Wedeen, T. L. Davis, J. W. Stages, V. S. Caviness, E. Kaplan, B.
R. Rosen, D. N. Pandya, and D. N. Kennedy M Ir'phometry of in vivo human
white matter association pnthiv-s-s with diffusion-weighted magnetic resonance
imaging," Ann. Neurol., vol. 42, pp. 951-962, 1999.









[30] S. Mori, B. J. Crain, V. P. 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. 265-269, 1999.

[31] M. Nitzberg and T. Shiota, "Nonlinear image filtering with edge and cor-
ner enhancement," IEEE Transactions on Pattern Aiil,i,J,; and Machine
Intelligence, vol. 14, no. 8, pp. 826-832, 1992.

[32] P. J. Olver, G. Sapiro, and A. Tannenbaum, "Invariant geometric evolutions
of surfaces and volumetric smoothing," SIAM J. Appl. Math., vol. 57, pp.
176-194, 1997.

[33] G. J. M. Parker, C. A. M. Wheeler-Kinghtshott and G. J. Baker, "Distributed
anatomical brain connectivity derived from diffusion tensor imaging," in
Proceedings of the 17th Intl. Conf. on Information Processing in Medical
Imaging, Davis, CA, pp. 106-120, 2001.

[34] P. Perona and J. Malik, "Scale-space and edge detection using anisotropic
diffusion," IEEE Trans. PAMI, vol. 12, no. 7, pp. 629-639, 1990.

[35] C. Pierpaoli, and P. J. Basser "Toward a quantitative assessment of diffusion
anisotropy," Magn. Reson. Med. vol. 36, 893-906, 1996.

[36] C. Pierpaoli, P. Jezzard, P. J. Basser, A. Barnett and G. Di Chiro "Diffusion
tensor MR imaging of the human brain," RiIl..I i vol. 201, pp. 637-648, 1996.

[37] L. Lapidus and G. F. Pinder, Numerical solution of partial differential equa-
tions in science and engineering, John Wiley and Sons, NY, 1982.

[38] C. Poupon, C. A. Clark, V. Frouin, J. Regis, I. Bloch, D. LeBihan, and J.
Mangin, "Regularization of diffusion-based direction maps for the tracking of
brain white matter fascicles," Neurolmage, vol. 12, 184-195, 2000.

[39] W. H. Press, B. P. Flannery, S. A. Teukolsky and W. T. Vetterling, Numerical
recipes in C: The art of scientific computing. Cambridge University Press,
Cambridge, England, second edition, 1992.

[40] L. Rudin, S. Osher, E. Fatemi, "Nonlinear total variation based noise reduc-
tion algorithms," Phi,,;i ,. D, vol. 60, pp. 259-268, 1992.









[41] G. Sapiro and D. L. Ringach, "Anisotropic diffusion of multivalued images
with applications to color filtering," IEEE Trans. on Image Processing, vol. 5,
pp. 1582-1586, 1996.

[42] J. Sethian, Level-set methods, Cambridge University Press, Cambridge, 1998.

[43] J. Shah, "A common framework for curve evolution, segmentation and
anisotropic diffusion," in IEEE Conf. on Computer Vision and Pattern
Recognition, pp. 136-142, 1996.

[44] D. Stalling and H. Hege, "Fast and independent line integral convolution," in
ACM SIGGRAPH, pp. 249-256, 1995.

[45] J. Weickert, "A review of nonlinear diffusion filtering,," in Scale-space theory
in computer vision,, (Eds.) B. ter Haar Romney, L.Florack, J. Koenderink, and
M. Viergever, Eds., vol. 1252, of Lecture Notes in Computer Science,, pp. 3-28,
Springer-Verlag, 1997.

[46] C. F. Westin, S. E. Maier, B. Khidir, P. Everett, F. A. Jolesz and R. Kikinis,
"Image processing for diffusion tensor magnetic resonance imaging," in Proc.
of the Second Intl. Conf. on Medical Image Computing and Computer Assisted
Interventions(MICCAI), pp. 441-452, 1999.

[47] R. Xue, P. C. M. van Zijl, B. J. Crain, H. Solaiyappan, and S. Mori "In vivo
three-dimensional reconstruction of rat brain axonal projections by diffusion
tensor imaging," Magn. Reson. Med. vol. 42, pp. 1123-1127, 1999.















BIOGRAPHICAL SKETCH


Tim McGraw was born in Key West, Florida. He received his Bachelor of

Science degree from the Mechanical Engineering Department at the University of

Florida. He will receive his Master of Science degree in computer and information

sciences from the University of Florida in August 2002. His research interests

include medical imaging, image processing, computer vision and computer graphics.