|UFDC Home||myUFDC Home | Help|
This item has the following downloads:
NEURONAL FIBER TRACKING IN DT-MRI
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
Tim E. McGraw
For my wife, Jo, and my mother, Patti.
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
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
ACKNOWLEDGMENTS ................... ...... iv
LIST OF FIGURES ..................... ......... vii
KEY TO ABBREVIATIONS ................... ....... viii
ABSTRACT ....................... ........... ix
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
220.127.116.11 Perona-M alik ....... ......... 15
18.104.22.168 Tensor anisotropic diffusion .......... 16
22.214.171.124 ALM .............. ... 17
2.2.3 Variational Formulation ..... . . ..... 17
126.96.36.199 Membrane spline . . ..... 18
188.8.131.52 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
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
Tim E. McGraw
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.
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.
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
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
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
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
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  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
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  then the
first steps were taken to relate this structural connectivity to function .
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  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. .
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
In () D (2.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.  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 S1 1 -bl -bl -b' -2b' -2b' -2bJl D
= D (2.2)
ln S28 1 -b -b 2 -b 2 -2b28 -2b28 -2b28 D
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 . 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
(x,t) = V21(X, t)
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, ) = Uo(w) (2.4)
The solution to 2.4 is
U(w, t) = Uo(w)e-w,2 t
The solution to 2.3 is then the inverse Fourier transform of 2.5
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
An anisotropic diffusion process which inhibits blurring at edges was formu-
lated by Perona and Malik  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
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
=, 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.
184.108.40.206 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  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))
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.
Alvarez, Lions, and Morel  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
St) + V (x(t), t) '(t) = 0 (2.12)
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)
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
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.
220.127.116.11 Membrane spline
A first-order stabilizer, stretching energy (arc-length), is one possible stabilizer
E,= Ij + + I ldxdy (2.17)
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.
18.104.22.168 Thin-Plate spline
Using a second-order stabilizer bending energy (curvature)
Et,, =J I + 2I + I2 dx dy (2.18)
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 . 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
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)
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)
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  and Caselles et al..
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 . The equations
=, t) div(c(I)Vl,(x, t))
In(x,0) = I,o(x) (2.21)
are coupled through the function C, and can be written in vector form as
( =, V (c(I)VI(x, t))
I(x, 0) =Io(z) (2.22)
2.3.2 Riemannian Metric Based Anisotropic Diffusion
Sapiro et al.  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
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
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.  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  introduced the TVn,m norm for vector valued images.
mTV Ix)) ] (.6
TVnm(I(x)) = [TV1,i(I()]2 (2.26)
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
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
lnS+ = lnSo -b,jD,j (3.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
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
minS(S) = [g(A+, A_) VSi + Si il2]dx (3.2)
s 2n 2
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. .
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 . 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
= div VS -/[ (Si S) i = 1, ..., 7
Ot ||VS ||ll
-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.
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., , 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.  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.  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.  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. 
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.  also presented a follow up article on
fiber tract mapping wherein, they use the idea of the fast marching method from
Sethian et al.  for growing seeds initialized in the smoothed dominant eigen
vector field. Batchelor et al.  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
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 . 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 .
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.
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 .
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
The gradient descent of (4.2)
S= c'(p)kn + 3[c"(p) V(c(p))c'(p)] +
P3V(c(p))(c'(p) v(c(p))) (4.3)
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
5.1 Data Denoising
The nonlinear PDE for denoising the raw data is
div(g ) (s so) = 0 (5.1)
where g : R3 -- IR is the selective smoothing and p is the constant regularization
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 . By considering |Vs| to be a constant for each iteration,
and using the value from the previous iteration we can instead solve
(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)
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)
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
0 -1 6+w3 1 ... ... 1 ...
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
-= c'(p)lkn + p[c"(p) V(c(p))c'(p)] +
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)
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 '|
( C' + p(c" + VT(c v) VI') (5.11)
We linearize 5.11 using the concept of lagged-diffusivity, as we did in the data
c (t+At Ct+At(
( ) = ( ) (5.12)
We can simplify notation by calling the coefficients of c" and c+', respectively
at(p) and Lh(p) and calling the additive constant Kt.
= 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))
The averaged difference for the first derivative is
= 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.
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 however, successfully incorporated
ellipsoids in a layered visualization approach.
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  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.  -a concept
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)
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
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
Stalling and Hege  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.  have used LIC to visualize fibers from diffusion
tensor images of the myocardium.
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
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
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.
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
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.
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).
CONCLUSIONS AND FUTURE WORK
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.
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  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
= o (2)
Since we know that f minimizes I, we can say
o = 0 (3)
is a condition for the minimization of I.
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
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
+ jrr + rz dx dy dz
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
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,
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
(s so) (g
This can be rewritten in vector notation as
y(s so) div(g)Vs
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'
Differentiating with respect to e gives
v(c + eq(p))12dp
OI[c + erq]
c'(p) + eW'1
The Euler-Lagrange condition is given by
Ic'(p) I + 0(c'(p)
Integrating by parts we obtain
C + P[c"(p)
V(c(p))c'(p)] + PVT(c(p))(c'(p)
v(c + Ce))(7'
e 1 0
Vv(c + cq). ) (14)
where V is the 1. .,1L,; of v, (d flu, 1 as
Vix Vl-y Vz
V = J= v2; U2z
 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.
 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.
 P. J. Basser, J. Mattiello and D. LeBihan "MR diffusion tensor spectroscopy
and imaging" B;i..I,,, J. vol. 66, 259-267, 1994a.
 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,
 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.
 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.
 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.
 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.
 B. Cabral, "Imaging vector fields using line integral convolution," Computer
Graphics Proceedings, pp. 263-270, 1993.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 A. Chorin, Computational Fluid Mechanics, Selected papers, Academic Press,
 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,
 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,
 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,
 F. Crick and E. Jones "Backwardness of human neuroanatomy," Nature vol.
361, 109-110, 1993.
 L.C.Evans, Partial Differential Equations, American Mathematical Society,
Providence RI, 1998.
 E. Haacke, R. Brown, M. Thompson, R. Venkatesan. Magnetic Resonance
Imaging-Phi,,~ i,1 Principals and Sequence Design, John Wiley and Sons, NY,
 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.
 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
 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.
 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.
 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.
 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.
 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.
 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.
 P. J. Olver, G. Sapiro, and A. Tannenbaum, "Invariant geometric evolutions
of surfaces and volumetric smoothing," SIAM J. Appl. Math., vol. 57, pp.
 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.
 P. Perona and J. Malik, "Scale-space and edge detection using anisotropic
diffusion," IEEE Trans. PAMI, vol. 12, no. 7, pp. 629-639, 1990.
 C. Pierpaoli, and P. J. Basser "Toward a quantitative assessment of diffusion
anisotropy," Magn. Reson. Med. vol. 36, 893-906, 1996.
 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.
 L. Lapidus and G. F. Pinder, Numerical solution of partial differential equa-
tions in science and engineering, John Wiley and Sons, NY, 1982.
 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.
 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.
 L. Rudin, S. Osher, E. Fatemi, "Nonlinear total variation based noise reduc-
tion algorithms," Phi,,;i ,. D, vol. 60, pp. 259-268, 1992.
 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.
 J. Sethian, Level-set methods, Cambridge University Press, Cambridge, 1998.
 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.
 D. Stalling and H. Hege, "Fast and independent line integral convolution," in
ACM SIGGRAPH, pp. 249-256, 1995.
 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,
 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.
 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.
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.