UFDC Home  myUFDC Home  Help 



Full Text  
NEURONAL FIBER TRACKING IN DTMRI 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 RO1NS42075. 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 DTMRI Data Acquisition ........ ........... 5 2.1.1 Overview of Diffusion ........ ......... 8 2.1.2 Overview of MR Imaging ..... ....... 9 2.1.3 DTMRI Acquisition ......... ......... 10 2.1.4 StejskalTanner Equation ...... ........ 11 2.2 PDE Based ScalarValued Image Restoration ........ 12 2.2.1 Linear Filtering ........ ......... ... 13 2.2.2 Nonlinear Filtering ........ .......... 14 2.2.2.1 PeronaM 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 ThinPlate spline .............. .. 19 2.2.4 Total Variation Image Restoration . .... 19 2.3 PDE Based VectorValued Image Restoration . ... 21 2.3.1 VectorValued 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 DTMRI 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 FixedPoint LaggedDiffusivity ...... ..... 33 5.1.2 Discretized Equations ........ .......... 34 5.2 Fiber Regularization. ................ ..... 35 5.2.1 LaggedDiffusivity .. ........ 36 5.2.2 CrankNicholson 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 NonTensor Model of Diffusion . . ..... 49 8.2.5 Automatic RegionOfInterest Extraction ...... .. 49 APPENDIX DERIVATION OF EULERLAGRANGE 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 DTMRI 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 TVnorm 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 threedimensional 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 noninvasive 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 (DTMRI) 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 whitematter 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 zeromean, and additive in nature. Integrating a fiber path over a noisy field is an illposed 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 DTMRI 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. \ 4t.' 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 particlebased 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 realtime 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 DTMRI data acquisition process. An overview of PDE based image denoising will follow. Some linear and nonlinear techniques for scalar and vectorvalued images will be reviewed. Variational techniques for formulating image denoising, and the TV norm will be presented. Chapter 3 will describe our vectorvalued 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 vectorvalued images. In Chapter 4 we present the process of tracking neuronal fiber bundles through the restored DTMRI 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 DTMRI 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 DTMRI 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 apparentdiffusiontensor of MR measured translational selfdiffusion. 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 T2 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 positivedefinite. Here we will take the approach of smoothing the observed vectorvalued 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, positivedefinite 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 reemitted 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 DTMRI 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 DTMRI however. 2.1.4 StejskalTanner Equation The intensity of the received signal, S, at each voxel depends on the properties of the diffusionencoding gradient, and the apparent diffusion tensor, D at that location. The StejskalTanner 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 diffusionencoding 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 ScalarValued 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 edgedetection, 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 PDEbased 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)ew,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 everincreasing 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 scalespace image representations, such as image pyramids. 2.2.2 Nonlinear Filtering The nonlinear filters described in this chapter were designed to overcome the interregion 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 PeronaMalik 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 steadystate 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 tensorvalued 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 levelset 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 levelset 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. Curveshortening flow of the isocontours of 4 smoothes the image I. 2.2.3 Variational Formulation The problem of image denoising is illposed, 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 firstorder stabilizer, stretching energy (arclength), is one possible stabilizer E,= Ij + + I ldxdy (2.17) J JR2 The EulerLagrange 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 ThinPlate spline Using a secondorder stabilizer bending energy (curvature) Et,, =J I + 2I + I2 dx dy (2.18) R2 as a stabilizer we obtain a thinplate spline solution. For this spline, p/ controls the "stiffness" of the spline surface. The EulerLagrange condition for the minimization of 2.18 is V41 = 0. 2.2.4 Total Variation Image Restoration Another sideconstraint 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 EulerLagrange 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 intraregion (edges), this is a useful model for image denoising. The EulerLagrange condition for the minimization of 2.19 written in gradientdescent 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 intraregion smoothing, and edges have also been preserved. Figure 2.3: Noisy image (left) and restored (right) by TV norm minimization. 2.3 PDE Based VectorValued Image Restoration The main application of vectorvalued image restoration has been the restora tion of color images. The simplest implementation of vectorvalued smoothing is uncoupled smoothing. By treating each component of the vector field as an independent scalar field, we can proceed by smoothing channelbychannel 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 PDEbased 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 VectorValued Diffusion Whitaker and Gerig introduced anisotropic vectorvalued 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 vectorvalued 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 vectorvalued image I(x, t) g 2I(z, t) Oit 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 flowbased 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 edgedetecting 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 EulerLagrange 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 DTMRI 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 StejskalTanner 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. Lowintensity S measurements may be overwhelmed by noise. However, since p is a zeromean distribution, these measurements may become negative. Clearly, we cannot calculate D in this case. Nor should we propagate the noise model through the SkeskalTanner 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 TVnorm 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 3dimensional 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, noninvasive 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 DTMRI 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 DTMRI 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 RungeKutta integrator [39]. However, this may not yield a regularized integral curve. In the context of the fiber tract mapping, two subtasks are involved namely, (1) estimating a denoised diffusion tensor from noisy vectorvalued image mea surements and (2) estimating the streamlines of the dominant eigen vectors of the diffusion tensor. The denoising involves a weighted TVnorm 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 whitematter 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 arclength 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 FixedPoint LaggedDiffusivity 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 "laggeddiffusivity" 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 righthand 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 SX1,yz A 's = Sx,y,z Sx,yl,z A S = Sxy,z  Sx,y,z1 aS = S(+1,,yz Sw1,y,z) AyS = I(Sx,y+l,z Sx1,y,z) 1 Azs = (SX,y,z+l Sx1,y,z) (5.5) We can rewrite 5.3 in discrete form using the definitions in 5.5 S_1,y,z Sx,yl,z SX,y,z1 + (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 bandeddiagonal 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 righthand 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 fixedpoint 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 LaggedDiffusivity We linearize 5.11 using the concept of laggeddiffusivity, 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 CrankNicholson Method We solve 5.13 using the CrankNicholson 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 finitedifference 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 threedimensional 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 DTMRI. 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 texturebased 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) JsoL 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 timeconsuming 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 windtunnel 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 framerates 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. Particlebased 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 TVnorm minimization formulation for smoothing vectorvalued 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 threedimensional 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 DTMRI data collection limits the number of samples taken in each direction. Outliers in the collected data may have a significant impact on the leastsquares calculation of the diffusion tensor. The use of a robust regression for calculating the diffusion tensor should be explored. 8.2.4 NonTensor 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 RegionOfInterest Extraction By registering a new data set with an atlas, the user could select the region ofinterest by anatomical name. Currently the regionofinterest for seeding fiber tracts is done by hand. This method is time consuming and potentially inaccurate. APPENDIX DERIVATION OF EULERLAGRANGE CONDITIONS Here we will review the calculus of variations, and present the derivation the EulerLagrange 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 EulerLagrange 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 EulerLagrange 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 EulerLagrange 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 Vly 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. 845866, June 1992. [2] P. J. Basser and C. Pierpaoli ili. rostructural and ]pl\ i.l. i ., features of tissue elucidated by quantitativediffusiontensor MRI," J. Magn. Reson. B vol. 110, 209219, 1996. [3] P. J. Basser, J. Mattiello and D. LeBihan "MR diffusion tensor spectroscopy and imaging" B;i..I,,, J. vol. 66, 259267, 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, 247254, 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, 333344, 1995. [6] P. J. Basser, S. Pajevic, C. Pierpaoli, J. Duda and A. Aldroubi, "In vivo fiber tractography using DTMRI data," Magnetic Resonance in Medicine, vol. 44, 2000, pp. 625632. [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. 121133. [8] P. Blomgren and T. F. Chan,"Color TV: total variation methods for restoration of vectorvalued images," IEEE Transaction on Image Processing, vol. 7, no. 3, pp. 304309, March 1998. [9] B. Cabral, "Imaging vector fields using line integral convolution," Computer Graphics Proceedings, pp. 263270, 1993. [10] V. Caselles, J. M. Morel, G. Sapiro and A. Tannenbaum,IEEE Trans. on Image Processing, special issue on PDEs and geometrydriven 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. 182193, 1992. [12] P. Charbonnier, L. BlancFeraud, G. Aubert, and M. Barlaud, "Two deter ministic halfquadratic regularization algorithms for computed imaging,," in in Proc. of the IEEE Intl. Conf. on Image Processing (ICIP), 1994, vol. 2, pp. 168172, 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. 354367, 1999. [14] C. Chefd'hotel, D. Tschumperl, R. Deriche, O. Faugeras, "Constrained flows of matrixvalued functions : Application to diffusion tensor regularization," in Proc. of the European Conference on Computer Vision 2002 (ECCV'02), Copenhagen, Denmark, May 2002, pp. 251265. [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. 131149, 2000. [16] P.J. Chiang, B. Davis and E. Hsu, "Lineintegral 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, 399412, 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, 1042210427, 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. 92105, 2001. [21] F. Crick and E. Jones "Backwardness of human neuroanatomy," Nature vol. 361, 109110, 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 ImagingPhi,,~ 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, 3741, 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. 350355, 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. 15541568, 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 pnthivss with diffusionweighted magnetic resonance imaging," Ann. Neurol., vol. 42, pp. 951962, 1999. [30] S. Mori, B. J. Crain, V. P. Chacko and P. C. M. van Zijl, "Threedimensional tracking of axonal projections in the brain by magnetic resonance imaging" Ann. Neurol., vol. 45, pp. 265269, 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. 826832, 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. 176194, 1997. [33] G. J. M. Parker, C. A. M. WheelerKinghtshott 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. 106120, 2001. [34] P. Perona and J. Malik, "Scalespace and edge detection using anisotropic diffusion," IEEE Trans. PAMI, vol. 12, no. 7, pp. 629639, 1990. [35] C. Pierpaoli, and P. J. Basser "Toward a quantitative assessment of diffusion anisotropy," Magn. Reson. Med. vol. 36, 893906, 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. 637648, 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 diffusionbased direction maps for the tracking of brain white matter fascicles," Neurolmage, vol. 12, 184195, 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. 259268, 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. 15821586, 1996. [42] J. Sethian, Levelset 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. 136142, 1996. [44] D. Stalling and H. Hege, "Fast and independent line integral convolution," in ACM SIGGRAPH, pp. 249256, 1995. [45] J. Weickert, "A review of nonlinear diffusion filtering,," in Scalespace 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. 328, SpringerVerlag, 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. 441452, 1999. [47] R. Xue, P. C. M. van Zijl, B. J. Crain, H. Solaiyappan, and S. Mori "In vivo threedimensional reconstruction of rat brain axonal projections by diffusion tensor imaging," Magn. Reson. Med. vol. 42, pp. 11231127, 1999. 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. 