University of Southern California Dissertations and Theses

Diffusion MRI white matter tractography: estimation of multiple fibers per voxel using independent component analysis

(USC Thesis Other)

Diffusion MRI white matter tractography: estimation of multiple fibers per voxel using independent component analysis

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content DIFFUSION MRI WHITE MATTER TRACTOGRAPHY: ESTIMATION OF
MULTIPLE FIBERS PER VOXEL USING INDEPENDENT COMPONENT
ANALYSIS
by
Chi Wah (Alec) Wong
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(BIOMEDICAL ENGINEERING)
August 2009
Copyright 2009 Chi Wah (Alec) Wong
Dedication
My thesis work is dedicated to
the Trinity God,
my loving family Frederick, Beatrice, Agatha Wong, and
my significant other Angela Lee
ii
Acknowledgements
Professor Manbir Singh has been a supportive advisor to me throughout my graduate studies. He
gave me the freedom to pursue independent work while providing invaluable supervision to make
sure I was on the right track. He is a student-oriented professor in that he is willing to spend
time interacting and discussing topics with his students. I truly appreciate his effort in helping
me develop critical thinking skills.
I am thankful to my dissertation committee (Professor K. Kirk Shung, Professor Richard Leahy
and Professor Krishna Nayak) for their comments on my thesis. A special thanks to Professor
Nayak, my first teacher at USC for introducing me into the exciting field of MRI. I also appreciate
his time discussing my work in detail.
I am so blessed to have friends and colleagues (Witaya Sungkarat, Jeongwon Jeong, Darryl
Hwang, Sinchai Tsao, Bryce Wilkins, Namgyun Lee, Vaibhav Bora, Aarti and Armita) to study
alongside with me throughout these years and for providing a supporting and serious yet fun
research environment. My memories with them are unforgettable. Special thanks to Darryl and
Bryce for the proofreading of my dissertation.
I would also like to express my sincere appreciation to my buddies at church. Their uncondi-
tional love and prayers have enabled me to grow throughout these years. They have helped me
through thick and thin and encouraged me to keep moving forward. Finally, I am thankful to my
Lord Jesus Christ. Throughout these years, He has faithfully showed me the way so I could grow
and bear fruit. All glory to God!
iii
Table of Contents
Dedication ii
Acknowledgements iii
List Of Tables vi
List Of Figures vii
Abstract xiii
Chapter 1: Introduction 1
1.1 Briefhistory. .. .. ... . . . . ... .. .. ... .. .. ... .. .. ... .. .. 1
1.2 Motivation . .. .. ... . . . . . .. .. .. ... .. .. ... .. .. ... .. .. 1
1.3 OriginalContributions .. .. .. ... . . . . . .. .. .. ... .. .. ... .. .. 3
1.4 Outline ... .. .. ... .. .. ... . . . . . .. .. .. ... . . . . ... .. .. 3
Chapter 2: Overview of Diffusion MRI and tractography 5
2.1 OverviewofMRI .. ... . . . . . .. .. .. ... .. .. ... .. .. ... . . . . 5
2.1.1 PrinciplesofMRI .. .. ... . . . . . .. .. .. ... . . . . . .. .. .. 5
2.1.2 Offresonanceandspinecho . . . . . . ... .. .. ... .. .. ... . . . . 8
2.1.3 EchoPlanarImaging. . . . .. .. .. ... . . . . . .. .. .. ... .. .. 8
2.2 Thehumanbrain .. ... .. .. ... . . . . . .. .. .. ... . . . . ... .. .. 10
2.3 OverviewofDiffusionMRI . . . . .. .. .. ... .. .. ... . . . . ... .. .. 12
2.3.1 DiffusionTensorImaging(DTI). . . . . .. .. .. ... .. .. ... .. .. 12
2.3.2 Eddycurrentcorrection . ... . . . . ... .. .. ... .. .. ... .. .. 14
2.3.3 Reconstructionofmultiplefiberorientations.. .. ... .. .. ... .. .. 15
2.3.3.1 Multi-compartmentmodel . ... .. .. ... . . . . ... .. .. 15
2.3.3.2 Diffusion spectrum imaging/Q-space imaging . . . . . . . . . . . . 16
2.3.3.3 Spherical de-convolution methods and Q-ball imaging . . . . . . . 17
2.4 DiffusionMRItractography .. .. . . . . . . . . .. .. .. . . . . . . . . .. .. . 19
2.4.1 Streamlinetractography . . .. .. .. ... .. .. ... .. .. ... . . . . 19
2.4.2 Probabilistic tractography . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Chapter 3: Overview of Independent Component Analysis (ICA) 28
3.1 Relationship with Principal Component Analysis (PCA) . . . . . . . . . . . . . . . 28
3.2 Maximizingnon-Gaussianity. .. ... . . . . ... .. .. ... .. .. ... .. .. 30
iv
Chapter 4: Estimating multiple fibers per voxel using ICA and statistical analysis 34
4.1 FormulationofICAindiffusionMRI .. .. .. . . . . . . . . .. .. .. . . . . . . 34
4.2 ICAprerequisitesindiffusionMRI .. .. .. ... . . . . ... .. .. ... .. .. 35
4.3 Methodology .. .. ... . . . . . .. .. .. ... .. .. ... .. .. ... .. .. 36
4.3.1 Filtering of Cerebral Spinal Fluid (CSF) with FA and MD using k-means
segmentation ... . . . . ... .. .. ... . . . . ... .. .. ... . . . . 37
4.3.2 ICA-basedmultiplefibersestimation . ... . . . . . .. .. .. ... .. .. 38
4.3.3 An example of source fiber direction recovery by ICA . . . . . . . . . . . . 41
4.3.4 Multiplefiberstreamlinetractography ... .. .. ... .. .. ... . . . . 42
Chapter 5: Monte Carlo simulation studies 44
5.1 Setup and generation of synthetic diffusion signal . . . . . . . . . . . . . . . . . . . 44
5.2 Estimatingorientationsofindividualfiberpervoxel . .. ... .. .. ... .. .. 47
5.2.1 case 1: number of original sources equals number of estimated sources with-
outCSF . . . . .. .. .. ... . . . . ... .. .. ... .. .. ... .. .. 48
5.2.2 case 2: number of original sources equals number of estimated sources with
CSF . .. .. ... . . . . . .. .. .. ... .. .. ... . . . . ... .. .. 49
5.2.3 case 3: number of original sources equals number of estimated sources with
partialzeroweighting .. ... . . . . . .. .. .. ... .. .. ... .. .. 49
5.2.4 case 4: number of original sources does not equal number of estimated sources 58
5.3 Estimatingnumberoffiberdirectionspervoxel .. .. .. ... .. .. ... .. .. 58
Chapter 6: Simulation study with synthetic phantom 65
6.1 Orthogonal crossing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
6.1.1 Fiberreconstruction . . . . .. .. .. ... . . . . . .. .. .. ... .. .. 65
6.1.2 Tractographyresults . .. ... . . . . . .. .. .. ... . . . . ... .. .. 69
6.2 60degreescrossing . . .. .. .. ... . . . . . .. .. .. ... .. .. ... . . . . 72
6.2.1 Fiberreconstruction . . . . .. .. .. ... . . . . . .. .. .. ... .. .. 77
6.2.2 Tractographyresults . .. ... . . . . . .. .. .. ... . . . . ... .. .. 81
Chapter 7: Experimental human studies of ICA-based tractography 89
7.1 Example1. . . . . . . .. .. .. ... . . . . . .. .. .. ... . . . . ... .. .. 90
7.1.1 Mapfornumberoffiberdirections .. ... . . . . ... .. .. ... .. .. 90
7.1.2 Frontal-occipital . .. .. ... . . . . ... .. .. ... .. .. ... . . . . 90
7.1.3 Cingulum . . . .. .. .. ... . . . . . .. .. .. ... .. .. ... . . . . 95
7.2 Example2. . . . . . . .. .. .. ... . . . . . .. .. .. ... . . . . ... .. .. 95
7.2.1 Mapfornumberoffiberdirections .. ... . . . . ... .. .. ... .. .. 95
7.2.2 Frontal-occipital . .. .. ... . . . . ... .. .. ... .. .. ... . . . . 95
7.2.3 Cingulum . . . .. .. .. ... . . . . . .. .. .. ... .. .. ... . . . . 102
Chapter 8: Conclusions and future research work 106
8.1 Conclusion . .. .. ... . . . . ... .. .. ... . . . . ... .. .. ... .. .. 106
8.2 Futureresearchwork ... .. .. ... . . . . ... .. .. ... . . . . . .. .. .. 107
References 108
v
List Of Tables
4.1 Mean and standard deviation of negentropy for diffusion MRI data from voxels
with known single fiber direction. Because there are only 25 data points per voxel,
bootstrapping is used to increase the data size in calculating negentropy in each
iteration. .. .. .. ... . . . . . .. .. .. ... .. .. ... .. .. ... .. .. 35
4.2 The cluster centers from k-means segmentation with lowest FA and higher MD
obtainedfromfivehealthysubjects. .. .. .. . . . . . . . . .. .. .. . . . . . . 38
5.1 The mean and standard deviation (s.d.) of the three eigen values obtained from cor-
pus callosum of 12 healthy adult subjects. The statistics are then used to randomly
generate eigen values in creating synthetic tensors. . . . . . . . . . . . . . . . . . . 45
5.2 The mean and standard deviation (s.d.) of the three eigen values obtained from
CSFof12healthyadultsubjects.... . . . . ... .. .. ... .. .. ... . . . . 45
vi
List Of Figures
2.1 The principle of spin echo: the drawing is presented in transverse (xy)plane. (a)
shows the spins right after 90
o
excitation pulse. (b) shows that the spins dephase
at time t due to difference in precessional frequency. A 180
o
pulse is applied to flip
all spins with respect to x-axis resulting in the spins in (c). After another time t,
thespinsrephaseasshownin(d). . .. .. .. . . . . . . . . .. .. .. . . . . . . 9
2.2 Spin echo pulse sequence (left) and its corresponding k−space trajectory (right). . 10
2.3 EPI sequence (left) and its corresponding k−space trajectory (right). . . . . . . . . 10
2.4 The construction of the brain from a high resolution T
1
MRIimage. . . . .. .. . 23
2.5 Using SPM to segment a T
1
images into grey matter, white matter and cerebral-
spinal fluid. The probability maps shown are adjusted for dynamic contrast for
easyvisualization. . ... .. .. ... .. .. ... .. .. ... .. .. ... .. .. 24
2.6 Dual-spin echo Pulse sequence for DTI using spin echo. Eddy current effect is
reduced using two 1809
o
pulseswithreversedgradients. . ... .. .. ... . . . . 24
2.7 Example for a color-coded FA map. Red, green and blue are used to represent the
three orthogonal directions. Higher FA implies the diffusion is more anisotropic. . 25
2.8 Map for Mean Diffusivity (MD). Diffusion is more restricted with lower MD value. 26
2.9 Tractography result for Frontal-Occipital tracts for a healthy adult subject with FA
threshold = 0.15 for tensor and vector interpolations. Vector interpolation trac-
tography results with more tracts because in the tracking process, a tract proceeds
as long as at least one of the eight neighbors carry FA larger than threshold. For
tensor interpolation tractography, the tract proceeds only if the interpolated tensor
carryFAlargerthanthreshold. . ... .. .. ... .. .. ... .. .. ... . . . . 27
4.1 Non-Gaussianity property of voxels with single fiber showing the histogram of ne-
gentropy of the corpus callosum. The red curve shows the corresponding Gaussian
distribution with the same mean and standard deviation. . . . . . . . . . . . . . . 36
4.2 Whole brain streamline tractography results with FA threshold of 0.5 and passing
through the genu/splenium of corpus callosum. . . . . . . . . . . . . . . . . . . . . 37
vii
4.3 Flow diagram of the proposed methodology for ICA-based multiple fibers estima-
tionandtractography. . . . . . . ... .. .. ... .. .. ... .. .. ... .. .. 38
4.4 10 neighborhood configuration of ICA source estimation. The 8-connected voxels
inthesamesliceandtheimmediateneighboraboveandbelowthevoxelofinterest
are used to obtain multiple measurements which are assumed to be the different
mixtureofsources. . ... . . . . ... .. .. ... .. .. ... .. .. ... .. .. 41
4.5 An example shows how ICA recover the source diffusion signal of one fiber direction.
Row 1: original diffusion signal from voxels with known single fiber direction. Row
2: estimated sources obtained from output of FastICA. Row 3: ICA estimated
sources scaled between 0.1 and 0.9 such that PCA can be used to estimate one
fiber direction for each source. The correlation coefficients between original sources
and scaled ICA estimated sources 1 and 2 are 0.9698 and 0.9798 respectively. . . . 43
5.1 The cone used to generated 3 fibers with identical inter-fiber angles. The radius of
the circle and the height of the isosceles triangle are adjusted to vary the inter-fiber
angle. . ... .. .. ... . . . . . .. .. .. ... .. .. ... . . . . ... .. .. 46
5.2 Simulation result with 2 source fibers and 2 estimated fibers and SNR=15. It is
observed that the error angle for ICA is around 5 degrees less than using Behrens
model. In ICA, Over 80% of trials have angle error less than 15 degrees with inter-
fiber angle larger than 50 degrees. In Behrens model, the percentage is less than
80%. . . ... .. .. ... . . . . . .. .. .. ... .. .. ... .. .. ... . . . . 50
5.3 Simulation result with 3 source fibers and 3 estimated fibers and SNR=15. It is
observed that the error angle for ICA is around 15-20 degrees with inter-fiber angle
between 30 and 70 degrees, which is less than Behrens model. The error angle is
higherthanthatofthe2-fibercase... .. .. ... .. .. ... .. .. ... .. .. 51
5.4 Simulation result with 2 source fibers and 2 estimated fibers and SNR=30. It is
observed that the error angle for ICA is less than 5 degrees while Behrens model
has error angle around 10 degree range with inter-fiber angle larger than 50 degrees.
In ICA, Over 80% of trials have angle error less than 15 degrees with inter-fiber
anglelargerthan40degrees. . . ... .. .. ... . . . . ... .. .. ... .. .. 52
5.5 Simulation result with 3 source fibers and 3 estimated fibers and SNR=30. It is
observed that the error angle for ICA is around 10-15 degrees with inter-fiber angle
between 30 and 70 degrees, which is less than Behrens model. The error angle is
higher than that of the 2-fiber case. The green curve shows the result in that in
each trial, the fiber with the largest error is eliminated. It is observed that the
mean error angles significantly drop. This shows that ICA can reliably recover 2
ofthe3sources. . . . .. .. .. ... .. .. ... .. .. ... .. .. ... .. .. 53
5.6 Simulation result with 2 source fibers and 2 estimated fibers and SNR=45. It is
observed that the result is similar to that of SNR=30, except with a few degrees
ofimprovement.. . . ... .. .. ... . . . . ... .. .. ... . . . . . .. .. .. 54
5.7 Simulation result with 3 source fibers and 3 estimated fibers and SNR=45. It is
observed that the result is similar to that of SNR=30, except with a few degrees
ofimprovement.. . . ... .. .. ... . . . . ... .. .. ... . . . . . .. .. .. 55
viii
5.8 Simulation result with a CSF component added into synthetic diffusion data and
ICA is used to recover the components. SNR=30 is used in this study. The number
of source fibers is identical to the number of estimated fibers. The CSF component
is not counted as one fiber source. It is observed that the performance is similar to
thepreviouscasewherethereisnoCSFcomponent. . . . ... .. .. ... .. .. 56
5.9 Simulation result with zero weighting assigned to three voxels for each fiber source.
SNR=30 is used in this study. This is to simulate the scenario where a fiber is not
spanning the entire neighborhood. It is observed that the result is better than the
case without zero weighting. It is because zero weighting enhances the weighting
contrastamongdifferentsources. . . . . . . . . . .. .. .. . . . . . . . . .. .. . 57
5.10 Simulation result where there is one fiber source, and with two and three estimated
sources. SNR=30 is used. It is observed that the original fiber is estimated with
lessthan5degreeerrorinbothcases. . . . . ... .. .. ... . . . . ... .. .. 59
5.11 Simulation result where there are two fiber sources, and with two and three esti-
mated sources. SNR=30 is used. It is observed that in both cases the error angle
is very similar except that in 3-fiber case,theextrafiberiserroneous. ... .. .. 60
5.12 Simulation result where there are three fiber sources, and with two and three esti-
mated sources. SNR=30 is used. It is observed that in both cases the error angle
islargerthan10degrees.. .. .. ... .. .. ... .. .. ... . . . . ... .. .. 60
5.13 Samples for 1 source fiber. It is observed that the original fiber source is reliably
estimated. However, extra estimatedfibersarenoise-like. ... .. .. ... .. .. 61
5.14 Samples for 2 source fibers. Similarly, it is observed that the original fiber sources
arereliablyestimated. . .. .. ... . . . . ... .. .. ... . . . . ... .. .. 61
5.15 Samples for 3 source fibers. It is observed that the original fiber sources are reliably
estimated in the 3 estimated fibers case. In the case where only two fibers are
estimated, they lie in-between the threeoriginalfibersources... .. .. . . . . . . 62
5.16 Number of fiber directions obtained using F-test with p-value 0.1 and SNR=30. It
is observed that the true positive percentage for single fiber is above 80%, while
for 2-fiber case, the percentage is above 80% with inter fiber angle larger than 30
degrees. The true positive percentage for three fibers is below 50%. However, it
is observed that the true positive percentage for multiple fibers is between 70-80%
with inter-fiber angle larger than 50%. This implies that a significant amount of
threefibervoxelsareclassifiedastwofibers.. ... .. .. ... .. .. ... .. .. 63
6.1 T
2
-weighted slice obtained from the b
0
diffusion MRI set with SNR=30. . . . . . . 66
6.2 FA map obtained from the diffusion MRI set for one slice. A mask is drawn to
includeonlythewhitemattervoxels.. . . . . ... .. .. ... .. .. ... .. .. 66
6.3 Color index for representing 1, 2 and 3 fiber directions per voxel. . . . . . . . . . . 67
ix
6.4 Map showing number of fiber directions per voxel with SNR=30 and p−value=
0.005. It is observed that F-test tends to overestimate the number of fibers in this
case. . . ... .. .. ... . . . . . .. .. .. ... .. .. ... .. .. ... . . . . 67
6.5 Map showing number of fiber directions per voxel with SNR=30 and p−value=
0.001. Compared to previous case, the single fiber braches are more clearly identified. 68
6.6 Map showing number of fiber directions per voxel with SNR=30 and p−value=
0.0005. It is observed that the lower the p-value, F-test tends to select model with
fewer fiber directions. p−value= 0.001 seems to offer the best compromise and it
isadoptedinsubsequentstudies. . .. .. .. ... . . . . ... .. .. ... .. .. 68
6.7 Map showing number of fiber directions per voxel with SNR=15 and p−value=
0.001. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
6.8 Reconstructed fiber directions with 90 degrees crossing and SNR=30. It is observed
that PCA failed to resolve the fiber crossing directions. In Behrens model, although
the fiber crossing is resolved, it is observed that the model did not converge for many
voxels. The background fiber direction is not properly recovered. In ICA, the fiber
directionsarereliablyrecovered. ... . . . . ... .. .. ... .. .. ... .. .. 70
6.9 Reconstructed fiber directions with 90 degrees crossing and SNR=15. Observations
are similar to previous case, except that the error is slightly worse than SNR=30. . 71
6.10 Tractography results with 90 degrees crossing and SNR=30. The tracts that pass
through the crossing region are shown. It is observed that in PCA, the tract fails to
propagate through the crossing region. In Behrens model, there is some tilting in
the tracts even though they could propagate through the crossing region. It seems
thatICAoffersthebesttractographyresult.. ... .. .. ... . . . . ... .. .. 73
6.11 Tractography results with 90 degrees crossing and SNR=15. The tracts that pass
through the crossing region are shown. Similar results are observed as in SNR=30,
exceptthatthereexistsslightlymoreerrorintracts. . . . ... .. .. ... .. .. 74
6.12 Tractography results with 90 degrees crossing and SNR=30. The tracts that origi-
nate from one branch are shown. It is observed that the tracts in PCA are diverted
to the wrong branch. In ICA and Behrens model, the tracts could pass through
thecrossingregion. . ... .. .. ... . . . . . .. .. .. ... .. .. ... .. .. 75
6.13 Tractography results with 90 degrees crossing and SNR=15. The tracts that origi-
nate from one branch are shown. Similar observations are obtained as in SNR=30
case. . . ... .. .. ... . . . . . .. .. .. ... .. .. ... .. .. ... . . . . 76
6.14 T
2
-weighted slice obtained from the b
0
diffusion MRI set with SNR=30. . . . . . . 77
6.15 FA map obtained from the diffusion MRI set for one slice. A mask is drawn to
includeonlythewhitemattervoxels.. . . . . ... .. .. ... .. .. ... .. .. 78
6.16 Map showing number of fiber directions per voxel with SNR=30 and p−value=
0.005. It is observed that F-test tends to overestimate the number of fibers in this
case. . . ... .. .. ... . . . . . .. .. .. ... .. .. ... .. .. ... . . . . 79
x
6.17 Map showing number of fiber directions per voxel with SNR=30 and p−value=
0.001. It is observed that the branches are more clearly idenitified than previous
case. . . ... .. .. ... . . . . . .. .. .. ... .. .. ... .. .. ... . . . . 79
6.18 Map showing number of fiber directions per voxel with SNR=30 and p−value=
0.0005. It is observed that the lower the p-value, F-test tends to select model with
lower fiber directions. p−value= 0.001 seems to offer the best compromise and it
isadoptedinsubsequentstudies. . .. .. .. ... . . . . ... .. .. ... .. .. 80
6.19 Map showing number of fiber directions per voxel with SNR=15 and p−value= 0.001. 80
6.20 Reconstructed fiber directions with 60 degrees crossing and SNR=30. It is observed
that in ICA, the individual 60 degree crossing fiber directions are recovered. In PCA
and Behrens model, the crossing directions are incorrect. . . . . . . . . . . . . . . . 82
6.21 Reconstructed fiber directions with 60 degrees crossing and SNR=15. . . . . . . . . 83
6.22 Tractography results with 60 degrees crossing and SNR=30. The tracts that pass
through the crossing region are shown. In this study, only ICA is capable of recov-
eringtheweakertractinblue... ... .. .. ... .. .. ... .. .. ... .. .. 84
6.23 Tractography results with 60 degrees crossing and SNR=15. The tracts that pass
throughthecrossingregionareshown.. .. .. . . . . . . . . .. .. .. . . . . . . 85
6.24 Tractography results with 60 degrees crossing and SNR=30. The tracts that origi-
nate from one branch are shown. It is observed that the stronger tract dominates
in the tractography result. Behrens and PCA failed to allow the weaker tract to
propagate through the fiber crossing region. . . . . . . . . . . . . . . . . . . . . . . 87
6.25 Tractography results with 60 degrees crossing and SNR=15. The tracts that origi-
natefromonebranchareshown. ... . . . . ... .. .. ... .. .. ... . . . . 88
7.1 Indexinformationforourexperimentalstudy. ... .. .. ... .. .. ... .. .. 90
7.2 Map of number of fiber directions with p−value= 0.001. It is observed that there
arevariousregionswithmultiplefibers. . .. ... . . . . ... .. .. ... .. .. 91
7.3 Corresponding FA map. It is observed that the high FA value voxels correspond
well with single fiber voxels including the corpus callosum. . . . . . . . . . . . . . . 92
7.4 Reconstructed fiber directions. Fiber directions are drawn within the ROI shown
in the FA map. It is known that corpus callosum, frontal occipital and a few other
minors tracts cross in this region. It is observed that ICA is capable of reducing
the partial volume problems between fibers. In ICA result, it is observed that the
Frontal-Occipital vectors are recovered as pointed by the arrow. In PCA case, the
tractisdominatedbythecorpuscallosum... ... . . . . . .. .. .. ... . . . . 93
7.5 Tractography results showing tracts passing through both ROIs. It is observed that
in PCA the frontal occipital tracts failed to propagate into the occipital area. In
ICA,thereissignificantimprovement. . . . . ... .. .. ... . . . . ... .. .. 94
xi
7.6 Tractography results showing tracts passing through both ROIs. . . . . . . . . . . 94
7.7 Reconstructed fiber directions for cingulum in sagittal view. Fiber directions are
drawn within the ROI shown in FA map. It is known that corpus callosum and the
cingulate gyrus tracts cross in this region. In ICA result, it is observed that the
individual fiber directions are recovered as pointed by arrows. In PCA, the fiber is
lyingbetweenthetwodirections. . .. .. .. ... . . . . ... .. .. ... .. .. 96
7.8 Reconstructed fiber directions for cingulum in coronal view. Similar observation is
obtained in that the individual directions between corpus callosum and cingulate
gyrusarerecovered. ... . . . . . .. .. .. ... . . . . ... .. .. ... .. .. 97
7.9 Tractography results showing tracts passing through both ROIs in axial view. It
is observed that in PCA, there is a significant amount of tract diverted into the
corpus callosum where the problem disappears in ICA. . . . . . . . . . . . . . . . . 98
7.10 Tractography results showing tracts passing through both ROIs in sagittal view.
There is a significant reduction in partial volume problems between corpus callosum
and cingulate gyrus. It is observed that in ICA, the cingulate could propagate into
thehippocampus.. . . .. .. .. ... .. .. ... .. .. ... .. .. ... .. .. 98
7.11 Map of number of fiber directions with p−value= 0.001. . . . . . . . . . . . . . . . 99
7.12 Corresponding FA map showing that high FA value voxels correspond well with
single fiber voxels including the corpus callosum. . . . . . . . . . . . . . . . . . . . 100
7.13 Reconstructed fiber directions. From the pointed arrows, it is observed that ICA
is capable of recovering multiple fiber directions in fiber crossing regions. . . . . . . 101
7.14 Tractography results showing tracts passing through both ROIs. It is observed that
the partial volume problems between corpus callosum and frontal occipital tracts
issignificantlyreducedusingICA. . . . . . . ... .. .. ... . . . . ... .. .. 102
7.15 Tractography results showing tracts passing through both ROIs, plus a filter re-
movingtractspassingthroughthemiddle-line. .. .. .. ... .. .. ... .. .. 103
7.16 Reconstructed fiber directions for cingulum in sagittal view. Fiber directions are
drawn within the ROI showin in the FA map. It is known that corpus callosum
and the cingulate gyrus tracts cross in this region. Similar to the first example, it
is observed that the individual fiber directions are recovered as pointed by arrows.
InPCA,thefiberislyingbetweenthetwodirections. . . ... .. .. ... .. .. 104
7.17 Tractography results showing tracts passing through both ROIs in sagittal view.
It is observed that the cingulate gyrus tract propagates towards the hippocampus
whilePCAfailedtodoso. .. .. ... .. .. ... .. .. ... . . . . ... .. .. 105
xii
Abstract
Brownian motion is a phenomenon describing the random movement of particles. Water molecules,
the most abundant particle inside the body, diffuse in all directions with equal probability in free
space. In the presence of organized structure, diffusion is restricted. Diffusion Tensor Imaging
(DTI) was developed to model water diffusion in tissues. In the human brain, myelinated neuronal
fiber bundles hinder water diffusion. By knowing the amount of diffusion along different directions,
the underlying fiber bundle directions can be inferred.
DTI is implemented by adding pulsed magnetic field gradients to a T
2
-weighted spin echo pulse
sequence. A diffusion signal is obtained for each voxel by combining data using different gradient
directions. Apparent Diffusion Coefficient (ADC) profile is calculated to show the amount of
diffusion in different directions. A rank-2 tensor (ellipsoid) is then used to model the diffusion
profile. From the tensor, the direction of the fiber bundle can be estimated. This technique is
generally regarded as Principal Component Analysis (PCA).
However, a rank-2 tensor model fails to account for the possibility of multiple fibers in a
single voxel. When multiple fiber bundle directions exist within a voxel, the PCA estimated
orientation lies somewhere in-between the multiple fiber directions, and is an incomplete and
incorrect representation of the underlying fibers in the voxel.
A major problem of existing multiple fiber reconstruction methods is that they require High
Angular Resolution Diffusion Imaging (HARDI) which usually takes 60 to 100+ gradient directions
with multiple averages, requiring at least 15 minutes per scan.
xiii
In this thesis, a novel method is proposed using Independent Component Analysis (ICA) to
estimate the orientations of multiple fiber bundles in a single voxel. Diffusion data used with
this technique is obtained using clinical pulse sequences and are acquired in around 3-7 minutes.
The novelty of the approach is to use the non-Gaussian property of diffusion signals arising from
individual fiber directions. White matter neuronal fiber bundles are in orders of millimeters,
which are larger than a voxel, is utilized. It is very likely that a fiber bundle will span the width
of several voxels and this neighborhood information is used to estimate multiple fibers.
First, the pre-requisites of ICA were examined. By calculating the negentropy of the diffusion
data from voxels known to contain only a single fiber direction (corpus callosum) and bootstrap-
ping, we found that the diffusion signal caused by single fiber direction is non-Gaussian.
Monte Carlo simulation studies were carried out to estimate the accuracy of recovering the
number of fiber directions per voxel, as well as the accuracy of individual recovered fiber directions.
Results showed that the average absolute error angle was less than 5 degrees with inter-fiber angle
of above 50 degrees in the case of two fiber directions per voxel at SNR=30 and SNR=45.
By incorporating multiple fiber streamline tractography, the method was validated on synthetic
phantom data, which revealed significant improvement over both DTI and an alternative multi-
compartment approach. The approach was then applied to experimental human data.
xiv
Chapter 1
Introduction
1.1 Brief history
The Bloch equation was first developed [Blo53] to describe the behavior of nuclear magnetization in
the presence of a magnetic field. A measurable voltage is generated by electro-magnetic induction
due to the precession of nuclear spins around the main field in a direction perpendicular to the
main field. This voltage (signal) can be measured by a coil (the receiver coil). Stejskal and
Tanner modified the Bloch equations to account for anisotropic restricted diffusion by including
a pulsed magnetic field gradient into the spin echoes pulse sequence. Tanner developed a method
[Tan78] to calculate a scalar representing the effective diffusion that is averaged over the echo
time TE in a voxel. Basser et al [BML94] defined an effective diffusion tensor showing that the
equation developed by Stejskal can be simplified into a linear relationship known as Diffusion
Tensor Imaging (DTI).
1.2 Motivation
Diffusion Magnetic Resonance Imaging (D-MRI) offers a unique way to infer tissue microstructure
[BML94, BP96, AHK
+
00]. In particular, D-MRI is capable of recovering information on neuronal
fiber directions in white matter of the human brain through water diffusion. A rank-2 tensor is
1
commonly used to model the diffusion profile of a fiber in DTI. However, this model becomes
inaccurate when two or more fiber directions cross within a voxel [AHL
+
01].
Several strategies have been reported in the literature to resolve multiple fiber directions
inside a voxel. One alternative is to adopt High Angular Resolution Diffusion Imaging (HARDI).
Several methods were proposed in the literature to process the HARDI data [TRW
+
02, TRWW03,
TCGC04, Ale05c, And05]. Although HARDI based approaches allow reconstruction of multiple-
fiber directions per voxel, the major drawback is a significant increase in scanning time. In
particular, [TRW
+
02, TRWW03, TCGC04, And05] proposed several methods requiring 40 mins
to over one hour to finish the acquisition. Such long scan times cause motion artifacts and affect
the comfort-level of patients. Hence, it is desirable to have alternative methods to reconstruct
multiple fiber directions per voxel reliably using clinical MRI sequences that could acquire data
typically under 10 minutes.
In a previous study [KKJS05], an approach was proposed called the “modified mixture density
ICA” which tuned its nonlinearity functions to minimize errors between the prefixed nonlinearity
functions and the probability density functions of intermediate sources. It was assumed that there
were one or more neighboring voxels in the same imaging plane containing the same fiber tracts
but with different volume fractions. In estimating sources, 2 connected voxels including the voxel
of interest were used. The voxels were chosen to be in L shape and four configurations are chosen.
A two compartment model was then used to evaluate which configuration provided the best fit.
Even though the feasibility of this ICA model was demonstrated with preliminary human data
obtained using six gradient directions and four nonzero b−values, its application has been limited
due to its computational complexity and slow convergence. The method requires an accurate
estimation of the probability density function. In addition, no whole brain tractography results
were conducted to evaluate the accuracy of the method.
2
1.3 Original Contributions
In this thesis, the original contributions are as follows:
• Investigated statistical properties of diffusion profile with single fiber voxels.
• Used ICA to estimate 1-3 fiber directions per voxel.
• Used statistical hypothesis testing (F-test) for model selection among 1-3 fiber directions
using ICA estimated fiber directions.
• Extended single-fiber tractography to multiple fiber tractography.
• Validated ICA approach with Monte Carlo simulation and phantom simulation studies at
various Signal-to-Noise Ratios (SNR).
• Compared tractography results between ICA approach and conventional PCA-based DTI
using experimental human data.
• Implemented tractography algorithms in C and ICA estimation algorithms in Matlab
1.4 Outline
This dissertation is divided into two parts. Part I is an overview of the various topics pertinent to
the research work described in this dissertation. Chapter 2 provides a summary of Diffusion MRI
and tractography including techniques about acquiring diffusion weighted images, calculating
tensors and whole brain tractography. Common problems including partial volume effect and
eddy current correction are also discussed. Chapter 3 provides an introduction to Independent
Component Analysis (ICA), a statistical approach for separating the independent latent sources
from different measurements.
3
DTI fails to work when there is more than one fiber direction per voxel. This problem leads
to Part II of the dissertation, which contains the major original contributions. Chapter 4 formu-
lates the multiple-fiber reconstruction problem using negentropy-based ICA statistical analysis.
Negentropy is a measure describing the non-Gaussianity of a signal. ICA prerequisites are also
investigated in Chapter 4. Monte Carlo simluation is used in Chapter 5 to evaluate the accuracy
of our methods. In Chapter 6, a common DTI dataset developed by the ISMRM Perfusion and
Diffusion Study Group and multiple fiber streamline tractography is used to validate our method.
DTI data from a phantom with crossing fibers are used to compare our approach with conventional
DTI. An experimental human study is then described in Chapter 7.
4
Chapter 2
Overview of Diffusion MRI and tractography
2.1 Overview of MRI
2.1.1 Principles of MRI
Nuclei within atoms carry magnetic resonance property when they possess odd number of protons
or odd number of neutrons [Nis]. More specifically, each nucleons (protons or neutrons) carries
a nuclear spin giving rise to a magnetic moment. The charged nature of protons make them
contribute more to the magnetic moment than the neutral neutrons. When there is no net
external magnetic field in a volume, the magnetic moments of individual protons (also referred to
as spins in MRI) are oriented in random directions and hence, the net magnetic moment of the
volume is zero. On the other hand, in the presence of an external magnetic field (usually defined
as z or longitudinal direction), the spins will either align along or in the exact opposite direction
of the external magnet field. The difference between the parallel and antiparallel populations
creates a net magnetic moment in the volume. The nuclear spins will resonate at the Larmor
frequency ω, with its relationship with the external magnetic field B
0
as follows:
ω = γB
0
(2.1)
5
where γ is the gyromagnetic ratio (a constant that is unique for each element). In living
tissue, hydrogen H
1
is the most abundant element (part of water molecules). Hence, the protons
within water molecules is the most studied in MRI. The gyromagnetic ratio for a proton within
a hydrogen atom is
γ
2π
=42.58 MHz/Tesla.
Next, a radiofrequency RF electro-magnetic signal B
1
with frequency equal to the Larmor
frequency, is applied in the xy or transverse plane such that the spins are excited from their
equilibrium position. The spins tip away from the z direction and precess with respect to it.
More specifically, the B
1
field provides a torque to the spins such that the magnetic moments
flip from the z direction to the xy plane. The angle of flipping depends on the strength of B
1
and its duration. When the B
1
field is off, the magnetic moments will gradually flip back to the
equilibrium position (z−direction). The time constant characterizing the speed of retoring the
magnetic moment in z direction is called T
1
, while the time constant characterizing the speed of
the decay of the magnetic moment in the xy plane is called T
2
. Different tissue carries different
T
1
and T
2
values, which is used to create contrast in MRI.
A receiving coil is set up along the xy plane such that any changes of magnetization in the xy
plane are detected. Changes in magnetization in the xy plane induce an electromotive force which
becomes alternating current in the receiving coil. This generated signal is called Free Induction
Decay (FID). This signal is processed and an MR image is reconstructed. This general behavior
of magnetization M is characterized by the Bloch equation as follows:
dM
dt
= M× γB−
M
x
i + M
y
j
T
2
−
(M
z
− M
0
)k
T
1
(2.2)
A limitation of the above setup is that excitation and relaxation are applied to the entire vol-
ume. To achieve localization, a linear gradient field G is added to B
0
such that the magnetization
strength is spatially varying. For instance, if we want to enable a varying magnetization strength
in z direction, the gradient field is B
0
+ G
z
z. With this setup, the resonance frequency of spins
6
become a function of space in the z direction. By applying a B
1
field with a particular range of
RF frequencies, only those spins that are of frequency corresponding to B
1
are affected. Thus
spins located only at a certain spatial position along the z direction are excited. This is known as
selective excitation. In other words, this corresponds to exciting a particular slice in z direction
where all spins in the xy plane of that slice are excited.
In addition to applying linear gradients during selective excitation, gradients are also applied
in the phase encoding (y) and readout (x) directions. For instance, if we denote the linear gradient
in x and y directions as G
x
and G
y
respectively, and the magnetization of a slice as m(x, y), the
signal from free induction decay generated after application of linear gradients are as follows:
s(t
x
,t
y
)=

x

y
m(x, y)e
−iγGxxt
e
−iγGyyt
dxdy (2.3)
The equation shows the 2D Fourier transformation relationship between magnetization and
the signal. By varying G
x
,G
y
and t, we could collect magnetic signals corresponding to different
spatial frequencies and estimate m(x, y) using inverse Fourier Transform.
In general, the gradient could be time-varying and the signal equation becomes as follows:
s(t
x
,t
y
)=

x

y
m(x, y)e
−iγ

t
0
Gx(τ)dτx
e
−iγ

t
0
Gy(τ)dτy
dxdy (2.4)
By letting
k
x
=
γ
2π

t
0
G
x
(τ)dτ and k
y
=
γ
2π

t
0
G
y
(τ)dτ (2.5)
we could easily visualize the relationship between the magnetization m and signal s in k-space.
Pulse sequence design is to develop protocols specifying how to perform selective excitation, as
well as how to navigate the k-space by varying the gradients to acquire spatial frequency data for
different applications.
7
2.1.2 Off resonance and spin echo
So far in the discussion we assume the external magnetic field B
0
is perfectly uniform. In reality,
the MRI magnet is not homogenous. In addition, because of electronic shielding from electrons
orbiting outside the nucleus, the effective magnetic field experienced by the nucleus is reduced
by a small amount depending on the molecular configuration, a phenomenon known as chemical
shift. Moreover, there also exist differences in magnetic susceptibility in various regions in the
volume. These are most severe near the boundaries between materials of different densities.
Magnetic field inhomogeneities create spatial inhomogeneity in resonance frequencies, which
result in phase shift and de-phasing among spins. Phase shift in k-space results in position shift
in the reconstructed image, which distorts the original image especially when the volume carries
different molecules with different chemical shifts. De-phasing among spins will cause the T
2
decay
to become faster. When magnetic field inhomogeneities are included, the transverse decay is called
a T
∗
2
decay.
The spin echo pulse sequence was developed to compensate for the phase shift and de-phasing
artifacts caused by B
0
inhomogeneity and chemical shift by applying a 180
o
excitation. Figure
2.1 shows the principle in generating spin echo. All de-phased spins are re-phased at the echo
time and the readout signal goes back to the T
2
decay envelop. Figure 2.2 shows an example of
the pulse sequence for spin echo and its corresponding k−space trajectory.
2.1.3 Echo Planar Imaging
Single-shot Echo Planar Imaging (EPI) is an efficient imaging method because it scans through an
entire slice in one Free Induction Decay (FID). In multi-shot imaging, each FID scans only part
of the slice and each FID carries different phase errors. In single-shot EPI, since the entire slice
is scanned in one shot, the phase error problem is prevented. A rapid scanning protocol can also
prevent bulk motion artifacts. Figure 2.3 shows an example pulse sequence and its corresponding
8
t
180
o
t
x
y
1
2
3
3
2
1
(a) (b)
(c) (d)
Figure 2.1: The principle of spin echo: the drawing is presented in transverse (xy)plane. (a)
shows the spins right after 90
o
excitation pulse. (b) shows that the spins dephase at time t due to
difference in precessional frequency. A 180
o
pulse is applied to flip all spins with respect to x-axis
resulting in the spins in (c). After another time t, the spins rephase as shown in (d).
k-space trajectory for the spin-echo version of EPI. The main difference with spin-echo sequence
in Figure 2.2 is that instead of acquiring one line in k−space per echo, EPI acquires the entire
k−space in one echo.
A major drawback for single-shot EPI is its hardware requirement. Single-shot EPI requires
high powered gradients for very fast rising time. There is also more eddy current effect associated
with rapid change in gradients. However, with advances in MRI hardware, image quality for
single-shot EPI has improved significantly in the past few years and now it is one of the most
popular pulse sequences. Single-shot EPI is widely used in functional MRI and Diffusion Tensor
Imaging (DTI).
9
RF
G
z
G
y
G
x
90
o
180
o
tt
Slice selection
Phase encoding
Readout
Echo
k
x
k
y
12 3
1
2
3
Figure 2.2: Spin echo pulse sequence (left) and its corresponding k−space trajectory (right).
RF
G
z
G
y
G
x
90
o
180
o
tt
Slice selection
Phase encoding
Readout
Echo
k
x
k
y
Figure 2.3: EPI sequence (left) and its corresponding k−space trajectory (right).
2.2 The human brain
Thebraincanbepartitionedinto cerebrum, cerebellum, and brain stem. The cerebrum is divided
into two hemispheres connected by thick bundles of neuronal fibers called the corpus callosum.
The two hemispheres differ in how they process information. The left hemisphere focuses on
details, while the right hemisphere focuses on broad background [oA08]. The hemispheres are
further divided into frontal lobe, parietal lobe, occipital lobe and temporal lobe. Each of them
has different roles in processing information from outside world, controlling voluntary movement
and regulating cognitive functions.
10
The brain stem connects the spinal cord to the brain. It relays information between the brain
and the rest of the body. It also controls involuntary functions including heart rate, blood pressure
and breathing. The cerebellum resides on top of the brain stem right under occipital lobe. It is
responsible for maintaining balance and coordination of our body through receiving signal from
eyes, ears, muscles and joints. It is also responsible for motor learning process.
The limbic system lies in the center of the brain linking the brain stem and cerebral cortex
together. It mainly consists of the amygdala, hippocampus, fornix, cingulate gyrus, hypothalamus
and the thalamus. Figure 2.4 shows the basic construction of the brain.
In cellular level, the brain is made up of billions of neurons. Each neuron is composed of a cell
body including a nucleus, an axon and dendrites. Axon and dendrites extend from the cell body.
Axons transmit signals among neurons. The transmission distance is often very long and hence,
axons are myelinated which is made up of glial cells. It acts as an insulating layer to speed up
signal transmission. Dendrites receive signal from other axons.
Hundreds of myelinated axons group together in parallel forming neuronal fiber bundles, trans-
mitting information from one part of the brain to another. For instance, the corpus callosum
bundle connect the two hemisphere as described above.
Myelinated axons appear in white while the neuronal cell bodies, dendrites and blood vessels
appear in grey-brown. The human brain can be simplified as being composed of grey matter, white
matter and Cerebral-spinal fluid (CSF). Figure 2.5 shows an example of a T
1
-weighted image slice
of the human brain. Statistical Parametric Mapping software (SPM) is used to segment the T
1
images into grey matter, white matter and CSF. The objective of this thesis is to reconstruct
axonal fiber bundle connectivity in white matter.
11
2.3 Overview of Diffusion MRI
2.3.1 Diffusion Tensor Imaging (DTI)
Stejskal and Tanner modified the Bloch equations to account for anisotropic restricted diffusion
by inserting a pulsed magnetic field gradient into the spin echoes pulse sequence. The signal
attenuation formula can be written as follows:
ln[S/S
o
]=−γ
2
G
2
δ
2
(Δ− δ/3)D =−bD (2.6)
where S is the signal readout with gradient applied, S
0
is the signal without gradient, G is the
applied gradient field, D is the effective diffusion constant, and b is the diffusion weighted matrix
which is a function of the gradient pulse field duration and strength.
Figure 2.6 shows the pulse sequence with applied gradient fields. It is essentially an EPI
pulse sequence, except that a dephasing gradient and a rephasing gradient are added to the
pulse sequence. In a diffusion-free environment, when a dephasing gradient is applied, the spin
precession frequency varies spatially along the gradient direction. After that, when a rephasing
gradient (the exact same opposite effect as the dephasing gradient) is applied, the spins rephase.
However, when there is water diffusion, the spins’ location are altered when the rephasing gradient
is applied. Hence, the spins cannot be perfectly rephased. When the direction of molecular
motion is aligned with the gradient direction, there is the most phase incoherence while in the
case when the direction of molecular motion is perpendicular to the gradient direction, all spins
are perfectly re-phased. As a result, by applying gradients in different directions, we could obtain
diffusion weighted signals, which carry information about the extent of water diffusion in different
directions.
Tanner developed a method [Tan78] to calculate a scalar representing the effective diffusion
that is averaged over the echo time TE in a voxel. Basser et al [BML94] defined an effective
12
diffusion tensor showing that the equation developed by Stejskal [Ste65] can be simplified into a
linear relationship, which is known as Diffusion Tensor Imaging (DTI). A rank-2 tensor is used to
model the effective diffusion in one voxel. Thus, for one voxel, equation 2.6 becomes as follows:
ln[S/S
o
]=−Db (2.7)
where D and b are two vectors defined as follows:
D =[D
xx
,D
yy
,D
zz
,D
xy
,D
xz
,D
yz
]
b=[b
x
,b
y
,b
z
, 2
√
b
x

b
y
, 2
√
b
x
√
b
z
, 2

b
y
√
b
z
]
(2.8)
It is a linear equation in which S, S
0
are measured signals, while b is given. b is designed
based on the number of gradient directions, the strength and duration of the gradient fields. D
could then be estimated by linear regression. For instance, in our typical DTI setting, there are
25 gradient directions and hence, in the matrix formulation, ln[S/S
o
]isofsize1× 25, D is of size
1× 6, and b is of size 6× 25.
A rank-2 tensor in 3 dimensions can be characterized by its three eigen values and the three
corresponding eigen vectors. Various parameters can be generated based on these parameters.
For instance, Fractional Anisotropy (FA) is used to describe the anisotropy of the tensor, which
is defined as follows:
FA =

(λ
1
− λ
2
)
2
+(λ
2
− λ
3
)
2
+(λ
3
− λ
1
)
2
2(λ
2
1
+ λ
2
2
+ λ
2
3
)
(2.9)
The range of FA is 0−1, the higher the FA value, the more anisotropic the tensor. Furthermore,
color information could be encoded into the FA map showing the direction of anisotropic diffusion.
Figure 2.7(b) shows the color sphere definition, as well as an example of color-coded FA map.
Besides FA, Mean Diffusivity (MD) is another diffusion measure commonly used in DTI:
13
MD =
λ
1
+ λ
2
+ λ
3
3
(2.10)
MD tells us the diffusivity of a voxel. In particular, Cerebral Spinal Fluid (CSF) carries the
highest MD because of its free diffusion characteristic. White and grey matters carry lower MD
values. Figure 2.8 shows an example of a MD map.
Linear anisotropy (CL) and planar anisotropy (CP) are among other anisotropic measures
proposed in [AHK
+
00] defined as follows:
CL =
λ1−λ2
λ1+λ2+λ3
CP =
λ2−λ3
λ1+λ2+λ3
(2.11)
2.3.2 Eddy current correction
Due to the rapid gradient switching in EPI, eddy current geometric distortion remains a major
problem in DTI. These distortions depend on gradient directions and strength. The following
summarize the major cause of eddy current artifacts.
• For B
0
induced eddy current, the entire volume experiences a constant phase shift, which
in turn becomes a bulk shifting in the MR image in the same direction as the applied B
0
.
• For phase encode direction induced eddy current, each phase encode experiences different
phase shift in the phase encode direction, which in turn becomes a scaling effect in the MR
image in the same direction.
• For readout direction induced eddy current, each phase encode experiences different phase
shift in the readout direction, which in turn becomes a shearing effect in the MR image in
the same direction.
There have been various methods proposed to reduce eddy current effects. In general, they
can be classified into three categories.
14
• Pulse sequences are developed such that eddy current effect is minimized or reversed [ALA97,
PMP
+
00, RHWW03, TCS08].
• K-space calibration to obtain an eddy current distortion map, which is then used to correct
image distortion [JBP98, Hor99, CPGC99, BA00, PSBM05, CGS06].
• Post-processing techniques to co-register Diffusion-Weighted images to a reference image
with minimal distortion [HM96, Bas99, Bas01, BKKT04, RBB
+
04, SLC
+
04, AS05, ZHK
+
06].
2.3.3 Reconstruction of multiple fiber orientations
2.3.3.1 Multi-compartment model
Several models have been used previously to study the multi-fiber problem. Tuch et. al. modeled
the intra-voxel fiber crossing with a discrete mixture of Gaussian diffusion processes in slow
exchange [TRW
+
02], which is written as follows:
S/S
o
=

j
f
j
e
−bDj
(2.12)
The model is solved using a gradient descent method. Correlation between the reconstructed
signal and source signal is used for model selection to identify the number of fiber directions.
Kreher et. al. proposed a similar model [KSM
+
05] using two anisotropy tensors and one
isotropic tensor. In their case, F-test is used for model selection. In addition, Hosey et. al. applied
such multi-tensor model [HWA05] to estimate the number of fiber directions using probabilistic
model selection besides reconstructing fibers.
Behrens et. al. [BWJ
+
03] proposed a “stick and ball” model in which a stick was used to
represent one fiber while a ball was used to represent isotropic diffusion. The model is written in
Equation 4.4.
15
S
S
0
=(1−
m

i=1
f
i
)e
−bd
+
m

i=1
f
i
e
−bdr
T
RiAR
T
i
r
(2.13)
where b and r is the b-value and gradient direction associated with data acquisition, d is the
diffusivity, f is the compartment weighting, and RAR
T
is the anisotropy diffusion tensor such
that
A =
⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 10 0
00 0
00 0
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (2.14)
The advantage of this model is that the number of parameters was significantly reduced com-
pared to multi-tensor models without the second and third eigen vectors and eigen values. How-
ever, the major drawback of multi-compartment methods is that they are only capable of recov-
ering two fiber directions per voxel. When there are three or more fibers per voxel, the model
becomes unstable and results are not recorded in the literature.
2.3.3.2 Diffusion spectrum imaging/Q-space imaging
Diffusion Spectrum Imaging (DSI) or Q-space Imaging (QSI) applies the Fourier relationship
between the MRI diffusion signal and the diffusion probability density function p(x) without
assumptions. The relationship is written as follows:
S(q)= S
0

p(x)e
iqx
dx (2.15)
where q is the strength, orientation and duration of pulsed magnetic field gradients, S(q) is
the diffusion signal acquired at a specific q.
16
In DSI, HARDI is used where diffusion signals are acquired along a large number of wave
numbers q (> 100) to reconstruct the pdf [WHT
+
05]. The orientation distribution function
(ODF) can then be reconstructed as follows:
φ(ˆ x)=

∞
0
p(αˆ x)dα (2.16)
where ˆ x is a unit vector and α is a scalar. The ODF is the projection of p(x) on a unit sphere.
The strength is this approach is that, no model or assumption is necessary and the exact
probability density function of displacement of particles in different directions is recovered.
Liu et. al. developed a generalized diffusion tensor framework [LBAM04, LBM05]. The
advantages of higher order tensors is that high-order moments of the PDF are obtained. This
is significant when mutliple fibers exist per voxel and the underlying diffusion is non-Gaussian.
They showed that the rank-2 diffusion tensor approach was a second order approximation of this
framework. Its direct relationship with q-space imaging was also derived. The drawback of q-
space or related methods is that a large number of gradient direction acquisitions are required
which is not practical for clinical purposes.
2.3.3.3 Spherical de-convolution methods and Q-ball imaging
Tournier et. al. assumed that the diffusion characteristics of all fiber populations in the brain are
identical [TCGC04]. As a result, the diffusion signal (S) becomes the convolution between the
fiber orientation density function (ODF) and the diffusion profile (DF) corresponding to single
fiber direction written as follows:
S(θ,φ)= ODF(θ,φ)⊗ DF(θ) (2.17)
If n
th
order spherical harmonics decomposition was applied to each of the components, the
convolution relationship became a multiplication written as follows:
17
S
n
= DF
n
× ODF
n
(2.18)
ODF could then be recovered by inverting the diffusion profile matrix of a single fiber direction
and multiplying that with the diffusion signal matrix. Furthermore, an additional constraint that
the fiber orientation distribution be non-negative is proposed in [TYC
+
08]. This technique is
known as constrained spherical deconvolution.
On the other hand, Anderson developed a de-convolution kernel based on axially symmetric
tensors known as FORE-CAST [And05]. Similarly, the developed convolution relationship is
expressed as spherical harmonics.
Tuch showed that the ODF is the spherical radon transform of the diffusion signal evaluated at
a particular radius [TRWW03, Tuc04], which equals the projection of the PDF along a zero-order
Bessel function centered at the origin. This method is known as the Q-ball imaging. In contrast
to q-space imaging where the diffusion signals are acquired on 3-D q-space, q-ball imaging only
requires diffusion signal to be acquired on a sphere at a particular radius, which reduces the
scanning time.
Alexander et. al. proposed a maximum-entropy parameterization of the PDF projected on a
sphere known as the Persistent Angular Structure (PAS) [JA03, Ale05c]. The Fourier Transform of
this new PDF is fit against the diffusion signal acquired on the sphere using Levenberg-Marquardt
optimization.
Various kernel functions are developed in the family of de-convolution methods [JVO
+
07,
HMH
+
06]. De-convolution methods allow the continuous ODF to be successfully recovered with
less gradient directions than DSI. Usually, the actual fiber directions are identified by identifying
local maxima in the ODF and choosing the ones with highest values.
18
2.4 Diffusion MRI tractography
The principal eigen vector of the estimated tensor lies parallel to the local fiber tract direction
in white matter. Tractography is conducted to link the fiber directions together in an attempt
to reconstruct fiber tract information. There are various approaches to tractography. Streamline
tractography was proposed [BPP
+
00] using the Frenet equation. Park et. al. proposed fast
marching tractography with tractography results in [PSB
+
02]. Probabilistic Index of Connectivity
(PICo) is proposed in [PHWK03] to calculate the probability of connection between points or
regions. Tensor Deflection method is proposed in [LWT
+
03] utilizing the entire tensor to do
the tractography. Bootstrap white matter tractography is proposed in [LA05]. Among all, it
seems that streamline tractography and probability-based tractography are the most popular
methods. In our work, streamline tractography is adopted due to its simplicity. The method is
also convenient for visualization because different fiber tract bundles can be clearly shown in the
same brain.
2.4.1 Streamline tractography
Basser et. al. proposed that a white matter fiber tract could be represented by a vector which
constructs a curve in 3D space [BPP
+
00]. The unit tangent vector can be set equal to the principal
eigen vector of the estimated tensor at any position. As a result, the following relationship is
obtained:
dr(s)
ds
= ξ(r(s)) (2.19)
where r(s) is the vector at the current position, ξ is the principal eigen vector.
The differential equation can be solved by the Euler’s method or the Runge-Kutta method.
Runge-Kutta method has the advantage of its fourth-order accuracy while Euler method has the
19
advantage of its simplicity. For a given position s
0
, the next location s
1
for tracking is obtained
as follows:
r(s
1
)= r(s
0
)+ αξ(r(s
0
)) (2.20)
where α is the step size of the tractography.
Each tract begins at a point known as the seed point. In single-fiber streamline tractography,
every voxel in the brain is allocated a seed point. The Euler tracking process begins at each
seed point and two sub-tracks are created. In particular, one sub-track goes along a particular
direction and the other sub-track goes along the exact opposite direction. The two sub-tracks are
joined to create a final track corresponding to the seed point.
Fractional Anisotropy (FA) is usually used as a parameter in streamline tractography. In
general, the higher the FA, the more likely there will be a white matter fiber tract in that voxel.
If the FA is too low, the signal to noise ratio is too low for reliable estimation of the tensor. The
FA threshold for tractography is usually set to 0.1to0.2 in which tractography will stop if the
FA threshold of the current position is less than threshold value.
During tractography, an angular threshold is maintained at each step and within 1cm. In our
case, a 45 degrees angular threshold is chosen because tracts generally do not bend sharply. These
two adjustments are to make sure a tract does not loop in itself perpetually. The step size can
be chosen. Usually, step size of 0.2mm is used. If the step size is too low, the computational
complexity as well as storage requirement is increased. If the step size is too large, the tracking
error accumulates quickly.
In DTI, tensors and eigen vectors are only calculated at discrete positions in voxels. During
real-time tractography, at any location in the brain, the principal eigen vector is obtained by
interpolating among its eight neighbors. In general, it can be done by using tensor interpolation
or vector interpolation.
20
In tensor interpolation tractography, at any location, the principal eigen vector is calculated
from the interpolated tensor among its eight neighbors. The tracking process keeps going until
one of the above conditions are violated (FA and angular thresholds). In vector interpolation
tractography, at any given location, the principal eigen vector is calculated from the eigen vector
among its eight neighbors. The tracking process keeps going until all the FA of its eight neighbors
are below the threshold.
Figure 2.9 shows the results of tractography using tensor and vector interpolations for frontal
occipital tract of a healthy adult subject. A ”blue” brain is overlayed on the tracts to visualize
the location of the tracts within the brain. The frontal occipital tract is obtained by filtering out
the tracts passing through both the ROIs shown in green. It is observed that vector interpolation
recovers more tracts than tensor interpolation. It is because in our current implementation, the
FA threshold criterion for vector interpolation is more relaxed than that of tensor interpolation.
Vector interpolation only requires at least one neighboring voxel to exceed the FA threshold; in
tensor interpolation, it requires the interpolated tensor to exceed the FA threshold. It is observed
that the tracts with vector interpolation are marginally noisier than that of tensor interpolation.
2.4.2 Probabilistic tractography
Streamline tractography is a deterministic approach, in which fiber directions are generated by
interpolation of neighboring voxels. Such methods do not consider the uncertainty of the prin-
cipal fiber directions due to noise. In probabilistic streamline tractography or simply regarded
as probabilistic tractography, uncertainty is added into the calculation of fiber direction in each
step. Such uncertainty is characterized by a local probability density function (pdf), in which
the direction with highest probability align with the calculated principal fiber direction as in
streamline tractography. In other words, if the local pdf is a delta function, probabilistic trac-
tography is reduced into streamline tractography. Various pdfs are proposed in the literature
[PHWK03, BWJ
+
03, BBJ
+
07].
21
In general, the procedure of probabilistic tractography is summarized as follows:
1. Define a seed point
2. Construct a local pdf to represent the uncertainty of fiber directions
3. Randomly select a sample of fiber direction from the local pdf
4. Move in the selected direction
5. Repeat until stopping criteria is reached, which is identical to those in streamline tractog-
raphy
For each seed point, the above process is repeated many times (∼ 10000). A connectivity
probability map is generated by counting how many tracts passed through each voxel. The map
is then normalized to form the final map.
22
Parietal
lobe
Occipital
lobe
Frontal
lobe
Temporal
lobe
Cerebral
cortex
Cerebellum
Brain
stem
(a) Four lobes of the brain
Occipital
lobe
Parietal
lobe
Frontal
lobe
Thalamus
Brain stem
Cerebellum
Hypothalamus
(b) Structure of the brain
Figure 2.4: The construction of the brain from a high resolution T
1
MRI image.
23
(a) T
1
(b) grey matter (c) white matter (d) CSF
Figure 2.5: Using SPM to segment a T
1
images into grey matter, white matter and cerebral-spinal
fluid. The probability maps shown are adjusted for dynamic contrast for easy visualization.
Slice selection
Phase encoding
Readout
RF
G
z
G
y
G
x
G
90
o
180
o
180
o
Figure 2.6: Dual-spin echo Pulse sequence for DTI using spin echo. Eddy current effect is reduced
using two 1809
o
pulses with reversed gradients.
24
(a) Sphere shows color definition of di-
rections
(b) Map for color-coded FA
Figure 2.7: Example for a color-coded FA map. Red, green and blue are used to represent the
three orthogonal directions. Higher FA implies the diffusion is more anisotropic.
25
Figure 2.8: Map for Mean Diffusivity (MD). Diffusion is more restricted with lower MD value.
26
(a) Tensor interpolation
(b) Vector interpolation
Figure 2.9: Tractography result for Frontal-Occipital tracts for a healthy adult subject with FA
threshold = 0.15 for tensor and vector interpolations. Vector interpolation tractography results
with more tracts because in the tracking process, a tract proceeds as long as at least one of the
eight neighbors carry FA larger than threshold. For tensor interpolation tractography, the tract
proceeds only if the interpolated tensor carry FA larger than threshold.
27
Chapter 3
Overview of Independent Component Analysis (ICA)
In the statistical model of ICA [AH01], there are n observations modeled as linear combinations
of k sources:
x
i
=
k

j=1
a
ij
s
j
for all i=1, ..., n (3.1)
where all a
i
are real and all s
j
are statistically independent. s
j
are latent variables which
means they cannot be directly measured and are called independent components. All we have is
the observations x
i
. Using vector matrix notation, it is written as follows:
X = AS (3.2)
where X is the observation random vector of size n, A is the mixing matrix and S is the source
latent matrix. Our objective is to estimate the original sources S from observations X.Inthis
thesis, the upper case of indexed variables correspond to the matrix representation.
3.1 Relationship with Principal Component Analysis (PCA)
If Principal Component Analysis (PCA) is used to solve the above problem, a covariance matrix is
first obtained from X. Next, singular value decomposition (SVD) or eigen value decomposition is
28
used to recover the individual random sources s
n
that best represent the data. Since the recovered
sources are uncorrelated, the following is obtained:
E{s1,s2} = E{s1}E{s2} (3.3)
In ICA, on the other hand, we require the sources to be statistically independent. In order
for two sources to be statistically independent, for any two functions h
1
and h
2
, the following is
required:
E{h
1
(s
1
),h
1
(s
2
)} = E{h
1
(s
1
)}E{h
2
(s
2
)} (3.4)
As a result, independence implies uncorrelatedness but not the reverse. Hence, It is reasonable
to claim that ICA is more than PCA. PCA is half-way ICA [AH01]. Suppose Z is the new
observation matrix that is obtained from whitening the original observation matrix X. Whitening
is a process where the observation matrix is de-meaned and PCA is applied to extract components
that are orthogonal and best represent the data. The following is obtained:
Z = VX = VAS (3.5)
Now, let’s consider an orthogonal transformation of U of Z:
K = UZ (3.6)
The following is then obtained:
E{KK
T
} = E{UKK
T
U
T
} = UIU
T
= I (3.7)
29
which implies that K is uncorrelated as well. Hence, it is not sufficient to tell if the independent
components are given by Z or K if we just consider the uncorrelatedness or orthogonality alone.
However, whitening or conducting PCA before ICA has a significant advantage in implementing
ICA as the estimation of independent components now becomes more efficient. This is because an
n×n orthogonal mixing matrix contains n(n−1)/2 degrees of freedom, which reduces significantly
the number of parameters to estimate compared to estimating n
2
parameters for an exhaustive
search in the mixing matrix.
3.2 Maximizing non-Gaussianity
From Equation 3.2, independent components can be estimated as follows:
Y = BX (3.8)
where Y is the estimated S and B is the estimation mixing matrix, ideally we want B = A
−1
.
Non-Gaussianity, the most important element in ICA estimation, is then taken into account.
First of all, it is required that the independent components S be non-Gaussian. Let’s re-write the
above equation as follows:
y = b
T
X = bAS = QS (3.9)
where y is one estimated component and b is the corresponding mixing matrix.
From the central limit theorem, any mixture of random independent components is always
more Gaussian than each of the individual components. Hence, from Equation 3.9, y would be
more Gaussian than any of the components in S. It becomes least Gaussian when there exists a
mixing matrix such that only one of the components among S is non-zero. Thus, by maximizing
30
the non-Gaussianity of b
T
X,we obtain a Y with only one non-zero component, i.e., we recover
one component of the mixture.
There are various metrics to measure non-Gaussianity. A simple one is kurtosis, which is
defined as follows:
kurt(y)= E{y
4
}− 3(E{y
2
})
2
(3.10)
The major drawback of kurtosis is that it is very sensitive to outliers. As a result, another
consistent metric, negentropy, is used to determine non-Gaussianity of the signal, which is defined
as follows:
J(y)= H(y
gauss
)− H(y) (3.11)
where H is the differential entropy of a random vector y and y
gauss
is a Gaussian random
variable with the same covariance matrix as y. Negentropy has the property that it is always non-
negative and zero only when y has a Gaussian distribution. [AH01] shows that the negentropy
can be approximated as follows:
J(y)∝ [E{G(y)}− E{G(v)}]
2
(3.12)
G can be any non-quadratic function and v is a Gaussian variable with the same mean and
variance as y.Someexamplesof G that are used:
G(y)=tanha
1
y
G(y)= yexp(−
a2y
2
2
)
(3.13)
where a
1
is some constant between 1 to 2.
31
A simple example is given here to illustrate how the above descriptions can be used to recover
the independent components by maximizing non-Gaussianity.
Let there be two latent sources s
1
and s
2
. Let the negentropy for s
1
and s
2
be neg
s1
and neg
s2
respectively. Negentropy measures are non-negative and become larger with more non-Gaussian
source. Let there be two observations x
1
and x
2
, which are mixtures of the sources as written in
Equation 3.2. Let the negentropy of x
1
and x
2
be neg
x1
and neg
x2
respectively.
By central limit theorem,
neg
x1
≤ neg
s1
neg
x2
≤ neg
s2
(3.14)
Now, the observations, which is our only available information, contain a mixture of the
components as written in Equation 3.8. Let y
1
and y
2
be the mixture of x
1
and x
2
with negentropy
values neg
y1
and neg
y2
respectively. B is the corresponding mixing matrix. From the relationships
in Equation 3.9, we have:
neg
y1
≤ neg
s1
neg
y2
≤ neg
s2
(3.15)
By central limit theorem, the only case in which neg
y1
= neg
s1
and neg
y2
= neg
s2
is that
y
1
= ks
1
and y
2
= ks
2
,where k is a scalar. In this case, the non-Gaussianity of Y is maximized.
As a result, by varying the mixing matrix B such that the non-Gaussianity of Y is maximized,
the sources are recovered up to a scaling factor.
A fast fixed-point ICA algorithm is developed based on negentropy [AH01]. The advantage
of this ICA implementation is that the convergence is cubic, which is much faster than other
optimization approaches like gradient descent. In addition, the algorithm is automatic requiring
no settings for the learning rate. For the case of recovering one component, the algorithm is
summarized as follows [AH01]:
32
1. Randomly select a mixing matrix b
0
2. Calculate b
k+1
= E{xg((b
kT
x}− E{g

(b
kT
x)}b
k
3. Orthogonalize b
k+1
4. Repeat until convergence
where b
T
means the matrix transpose of b, g can be any non-quadratic function and g

is the
corresponding Gaussian variable with the same mean and variance as b
kT
x. For more than one
unit, in every iteration, the weighting matrix B is decorrelated using Gram-Schmidt decorrelation
or similar schemes. In our case, symmetric orthogonalization was used under the assumption that
all signals carry equal ranking and are estimated at the same time [AH01].
Besides the requirement that sources must possess non-Gaussian statistics, the sources should
also be statistically independent. However, in practical situations, ICA works well when the
sources are not perfectly independent [AH01].
In ICA, the variance of the independent components is not determined. From the mixing
equation, we can always add any scalar multiplier to one of the sources and divide the same scalar
in the corresponding mixing constant:
x
i
=
k

j=1
(
1
α
a
ij
)(s
j
x
j
), for all i=1, ..., n (3.16)
The formulation of ICA in diffusion MRI to estimate multiple fibers per voxel is addressed in
Chapter 4.
33
Chapter 4
Estimating multiple fibers per voxel using ICA and
statistical analysis
4.1 Formulation of ICA in diffusion MRI
As described inEquation2.6, S is the diffusion signal measured in a voxel at given b,where b is
the diffusion weighted matrix of one gradient direction. Here we normalize the diffusion signal
with the MRI signal without diffusion weighting and let that be a measurement x:
S(b)
S(b=0)
= x
i
(4.1)
If g gradient directions are used per voxel, the measurement vector x
i
is with dimension 1×g.
In our work, using the statistical model of ICA [AH01], we obtain n measurements x
1
...x
n
from
n different voxels that are modeled as linear combinations of k latent sources s
1
...s
k
.Eachofthe
k sources corresponds to a diffusion source signal with one fiber direction. In other words,
x
i
=
k

j=1
a
ij
s
j
for all i=1, ..., n (4.2)
where a
ij
denotes the mixing fraction. In vector form,
34
mean of negentropy variance of negentropy
genu of corpus callosum 2.6926 0.0000272
splenium of corpus callosum 2.595 0.000121
Table 4.1: Mean and standard deviation of negentropy for diffusion MRI data from voxels with
known single fiber direction. Because there are only 25 data points per voxel, bootstrapping is
used to increase the data size in calculating negentropy in each iteration.
X = AS (4.3)
The n measurements are taken from one current voxel and its n− 1 neighbors. Here, we
assume that a fiber bundle spans several voxels. This implies that neighboring voxels carry
similar diffusion signal with the same mixed fibers as the current voxel, but with a different
mixture weighting.
4.2 ICA prerequisites in diffusion MRI
In Chapter 3, ICA prerequisites are discussed. In this section, we examine them in the context of
diffusion MRI. For the requirement of source non-Gaussianity, we investigated the statistics of 3T
diffusion data from the genu and splenium of corpus callosum since these structures are composed
mostly of single fiber. These structures were identified by isolating voxels where the FA> 0.5in
12 healthy human subjects. Next, bootstrapping was used to obtain the statistics of the signal.
25 gradient directions were applied yielding 25 data points per voxel. In one step, 40 voxels were
chosen at random giving 25× 40 = 1000 data points. The negentropy was then computed for
these points, and the process was repeated 100 times. Table 4.1 shows the mean and variance of
the negentropy of the genu and splenium of corpus callosum.
Since the negentropy of the combined statistics obtained from bootstrapping method has very
low variance, we conclude that the diffusion signal corresponding to single fibers in the corpus
callosum is non-Gaussian. We extend this assumption to all diffusion signals carrying single
35
fiber direction. Figure 4.1 shows the probability density function (pdf) generated from one step,
showing the non-Gaussian property of the genu of the corpus callosum and the splenium of the
corpus callosum. The red curve shows the corresponding Gaussian distribution with the same
mean and variance. Figure 4.2 shows an example of tracts with FA threshold of 0.5.
0 0.2 0.4 0.6 0.8 1
0
0.005
0.01
0.015
0.02
0.025
0.03
S(b
1000
)/S(b
0
)
p(x)
(a) genu corpus callosum
0 0.2 0.4 0.6 0.8 1
0
0.005
0.01
0.015
0.02
0.025
S(b
1000
)/S(b
0
)
p(x)
(b) splenium corpus callosum
Figure 4.1: Non-Gaussianity property of voxels with single fiber showing the histogram of negen-
tropy of the corpus callosum. The red curve shows the corresponding Gaussian distribution with
the same mean and standard deviation.
4.3 Methodology
Figure 4.3 shows our proposed ICA multiple fiber reconstruction algorithm with streamline trac-
tography. First, diffusion weighted data with 25 gradient directions are acquired at GE 3T scanner
using double spin echo EPI sequence with TR= 8000ms and TE= 86.1ms. 25 gradient directions
with b = 1000s/mm
2
and two averages are used, along with an acquisition without the diffusion
gradients. Voxel size is 2× 2× 4mm
3
. Total scanning time is around 7 minutes per subject.
Next, DTI processing techniques are applied to calculate a rank-2 tensor per voxel and its associ-
ated values for Fractional Anisotropy (FA) and Mean Diffusivity (MD). In particular, voxels with
FA < 0.05 are very isotropic and likely to be pure CSF or noise. Those voxels are eliminated.
After that, we propose a series of processing methods described below to achieve multiple fiber
tractography.
36
Figure 4.2: Whole brain streamline tractography results with FA threshold of 0.5 and passing
through the genu/splenium of corpus callosum.
4.3.1 Filtering of Cerebral Spinal Fluid (CSF) with FA and MD using
k-means segmentation
K-means segmentation is conducted to segment a brain volume into three clusters with FA and
MD as features. We extract the centroid with the lowest FA and highest MD values, which is
likely to correspond to CSF. This step is not redundant with the previous one. It is because we
want to keep grey matter or multiple fibers crossing voxels with FA > 0.05 and low diffusivity
but remove CSF voxels which may also have FA > 0.05. The process is repeated on five healthy
subjects and results are listed in Figure 4.2.
We observe that the centroids are consistent among five subjects. The mean FA is 0.1 and
mean MD is 0.0014. The thresholds are then used on all subjects for CSF filtering.
37
Acquire diffusion-
weighted data with 25
gradient directions
Calculate fractional
anisotropy (FA) and mean
diffiusivity (MD) per voxel
CSF filtering using
k-means
segmentation with
FA and MD as
features
Remove noisy
voxels with a global
filter of FA<0.05
Estimate
orientation of
fibers assuming
1 (PCA), 2 and
3 (ICA)
Multiple-fiber
streamline
tractography
Apply model to
select among 1, 2
or 3 fibers using
statistical
hypothesis testing
Figure 4.3: Flow diagram of the proposed methodology for ICA-based multiple fibers estimation
and tractography.
FA MD
Subject 1 0.0992 0.0014
Subject 2 0.104 0.0013
Subject 3 0.0966 0.0013
Subject 4 0.0952 0.0015
Subject 5 0.0992 0.0014
Mean 0.0988 0.0014
Table 4.2: The cluster centers from k-means segmentation with lowest FA and higher MD obtained
from five healthy subjects.
4.3.2 ICA-based multiple fibers estimation
The model for the forward problem, i.e., to estimate measured data from the mixture of sources,
was based on the multi-compartment model proposed by Behrens et. al. [BWJ
+
03], which is
written in Equation 4.4.
S
S
0
=(1−
m

i=1
f
i
)e
−bd
+
m

i=1
f
i
e
−bdr
T
RiAR
T
i
r
(4.4)
where b and r is the b-value and gradient direction associated with data acquisition, d is the
diffusivity, f is the compartment weighting, and RAR
T
is the anisotropy diffusion tensor such
that
38
A =
⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 10 0
00 0
00 0
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (4.5)
The term r
T
RAR
T
r is the correlation scalar between the gradient direction and the actual
fiber direction. In particular, if both directions align, the correlation becomes 1 and hence, the
diffusion signal S is minimum. If both directions are orthogonal to each other, the correlation
becomes 0 and hence, the diffusion signal S is maximum.
The first term corresponds to compartment of isotropy diffusion (CSF), while the second
term corresponds to compartments of anisotropy diffusion (fibers). In our work, CSF filtering
is conducted before multiple fiber estimation. The equation can be approximated without the
isotropic diffusion compartment:
S
S
0
=
m

i=1
f
i
e
−bdr
T
RiAR
T
i
r
,where
m

i=1
f
i
=1 (4.6)
Diffusivity d is constant throughout all compartments in a voxel. It is observed that there
exists ambiguity between d and f
i
. Specifically, let’s examine one compartment:
f
i
e
−bdr
T
RiAR
T
i
r
(4.7)
Assume there is an error in estimating d which is written as follows:
f
i
e
−b(d+Δd)r
T
RiAR
T
i
r
(4.8)
This error can be absorbed into the weighting factor as follows:
f
i
e
−b(d+Δd)r
T
RiAR
T
i
r
= f
i
e
−bΔdr
T
RiAR
T
i
r
e
−bdr
T
RiAR
T
i
r
= f

i
e
−bdr
T
RiAR
T
i
r
(4.9)
39
where f

i
= f
i
e
−bΔdr
T
RiAR
T
i
r
which is a function of fiber/gradient directions and the estimation
error of d. The Behrens model is incomplete in modeling the underlying diffusion process in a
voxel. Even though CSF filtering is conducted before multiple fibers estimation, there is some
residue in each voxel. As a result, the weighting factor constraint

m
i=1
f
i
= 1 is relaxed. Hence,
the simplified Behrens model is re-written as:
S
S
0
=
m

i=1
f
i
e
−bdr
T
RiAR
T
i
r
,where0≤
m

i=1
f
i
≤ 1 (4.10)
which is adopted in our work for statistical analysis.
Kreher et. al. in [KSM
+
05] proposed a F-test procedure to choose the better model between 1
and 2 fibers. However, they did not incorporate 3-fiber case because of difficulty in convergence.
In our work, we extend their framework into 3-fiber case.
Assuming m=1, 2 and 3, F-test is based on the residue errors conducted to classify the voxels
into single and multiple fiber orientations, The F-test formulation is written as:
F(A, B)=
SSEA−SSEB
nB−nA
SSEB
S−n
b
(4.11)
where n is the number of variables in the model and S is the number of data points. In our
case it is 25, corresponding to 25 gradient directions. Here, model B represents more number of
fiber directions than model A. The null hypothesis is that both models are equally suitable given
a confidence level p, and the F-test is carried out as follows:
1. Specify a confidence interval (e.g. p=0.001)
2. Set the null hypothesis that models with m=1 and m = 2 are equal in variance
3. If the value obtained from F-test is smaller than the threshold of the p-value T
p
, classify the
voxel as single fiber orientation, else
4. Set the null hypothesis that models with m=2 and m = 3 are equal in variance
40
5. If the value obtained from F-test is smaller than the threshold of the p-value T
p
, classify the
voxel as two fiber orientations, else
6. Classify the voxel as three fiber orientations
The p-value can be varied to adjust sensitivity of voxel classification. Specifically, the larger
the p-value, the more likely for the F-test to select a higher-order model.
For each model, ICA is used to estimate the fiber directions R
i
. In our ICA configuration, we
used 10 neighbors of the current voxel to obtain 11 measurements as shown in Figure 4.4. Each
measurement is assumed to contain the same set of diffusion signals corresponding to different
fibers in a voxel, but with different mixing ratios.
Figure 4.4: 10 neighborhood configuration of ICA source estimation. The 8-connected voxels in
the same slice and the immediate neighbor above and below the voxel of interest are used to
obtain multiple measurements which are assumed to be the different mixture of sources.
As a result, f
i
become the remaining variables of the model. Hence, instead of non-linear
optimization, least square estimation can be used, which guarantees global minimum solution.
4.3.3 An example of source fiber direction recovery by ICA
Figure 4.5 shows an example of how ICA recovered the original sources. Normalized diffusion
data were taken from two voxels, each known to contain a single fiber. The first voxel was taken
from the corpus callosum and the second voxel was taken from the spinal cortical projection. The
41
original sources were mixed with random mixing ratios to create 11 measurements. ICA was then
used to estimate the source. The first row shows the original diffusion signal of the two sources.
The second row shows the ICA estimated sources and the third row shows the ICA estimated
sources scaled into 0.1 to 0.9. The range is chosen to ensure the sources are within the range of
normalized diffusion signal without getting too close to 0 or 1.
The correlation coefficients between original sources and scaled ICA estimated sources 1 and
2 are 0.9698 and 0.9798 respectively. In addition, the errors of estimated fiber direction for source
1 and 2 are 5.5 degrees and 2.14 degrees respectively. The shape of the signals between original
sources and ICA estimated sources are recovered well, and ICA is capable to separate the two fiber
orientations. In subsequent chapters where simulation studies are conducted, extensive analysis
in reconstructing fiber directions are shown.
4.3.4 Multiple fiber streamline tractography
The streamline tractography described in Chapter 2.4.1 is valid only for a single fiber per voxel, be-
cause the diffusion profile is modeled by a rank-2 tensor. To incorporate multiple fiber directions,
the conventional streamline tractography methodis modified as follows:
1. If a voxel carries single fiber, PCA is used to estimate the fiber direction.
2. if a voxel carries two or three fibers, the fiber directions calculated from ICA are used.
3. If there are n fiber orientations in a voxel, n seed points are allocated in that voxel. Each
seed point corresponds to one particular fiber direction.
4. When a tract navigates into a region in which at least one of the neighbors carry multiple
fibers, the fiber that is closest to the current fiber direction is chosen for each of such
neighbors.
5. Vector interpolation is always used in calculating the next tract direction.
42
5 10 15 20 25
0
0.5
1
Source 1
5 10 15 20 25
0
0.5
1
Source 2
5 10 15 20 25
4
6
8
5 10 15 20 25
4
6
8
5 10 15 20 25
0
0.5
1
5 10 15 20 25
0
0.5
1

Figure 4.5: An example shows how ICA recover the source diffusion signal of one fiber direction.
Row 1: original diffusion signal from voxels with known single fiber direction. Row 2: estimated
sources obtained from output of FastICA. Row 3: ICA estimated sources scaled between 0.1 and
0.9 such that PCA can be used to estimate one fiber direction for each source. The correlation
coefficients between original sources and scaled ICA estimated sources 1 and 2 are 0.9698 and
0.9798 respectively.
6. Instead of using all 8 neighbors, m neighbors with the least angles to the current direction
are used (m=1 to 8). The reason is that some of the neighbors might have a large angle
with the current direction and it might be erroneous to include them to calculate the new
direction. If too few neighbors are used, the tractography result will be noisier. In our case,
6 neighbors are chosen. The resulting direction is normalized into unit vector before moving
forward with specified step size.
In subsequent chapters, this new scheme is used to compare single and multiple fiber tractog-
raphy.
43
Chapter 5
Monte Carlo simulation studies
5.1 Setup and generation of synthetic diffusion signal
In this chapter, results of Monte Carlo simulation studies conducted to evaluate the accuracy
of the ICA-based multiple fibers estimation approach are presented. Synthetic diffusion signals
generated by a Gaussian mixture model with two and three fiber directions were considered as
follows:
S/S
o
=

j
f
j
e
−bDj
(5.1)
The eigen values of the rank-2 tensor D were randomly generated according to mean and
standard deviation of the eigen values of the corpus callosums of 12 healthy adult subjects. It is
known that the genu and splenium portions of the corpus callosum contain mostly a single fiber
direction per voxel. ROIs were first drawn for each subject and those voxels within those ROIs
and with Fractional Anisotropy values larger than 0.6 were used for the calculations. The mean
and standard deviation of the three eigen values for corpus callosum and CSF are shown in Table
5.1 and 5.2 respectively.
44
×0.001 in cm
2
s
−1
mean s.d.
1
st
eigen value 1.68 0.18
2
nd
eigen value 0.37 0.07
3
rd
eigen value 0.275 0.075
Table 5.1: The mean and standard deviation (s.d.) of the three eigen values obtained from corpus
callosum of 12 healthy adult subjects. The statistics are then used to randomly generate eigen
values in creating synthetic tensors.
×0.001 in cm
2
s
−1
mean s.d.
1
st
eigen value 3.64 0.24
2
nd
eigen value 3.14 0.14
3
rd
eigen value 3.11 0.21
Table 5.2: The mean and standard deviation (s.d.) of the three eigen values obtained from CSF
of 12 healthy adult subjects.
For each synthetic signal, random complex noise (Rician noise) was added with independent
real and imaginary parts, each with distribution N(0,σ
2
)and σ =1/SNR. 25 gradient directions
were used with b = 1000
s
mm
2
leading to matrix b with dimension 25× 6.
Each study was divided into eight sections with inter-fiber angles ranging from 10 degrees to
90 degrees. Each section corresponds to a specific range of inter-fiber angle in degrees:
[10− 20], [20− 30], [30− 40], [40− 50], [50− 60], [60− 70], [70− 80], [80− 90]
In each section, an inter-fiber angle is randomly selected from the specified range. 1000 trials
were conducted for each section. For each trial, after an inter-fiber angle is randomly selected,
the synthetic diffusion signal is created based on the Gaussian mixture model. In particular, for
the case of 3 synthetic fibers, the fiber directions are generated using a cone with a circular base.
The main axis of the cone passes through the center of the base at right angle to its plane. By
adjusting the length of the axis and the radius of the circle, the fiber orientations are created upon
the surface of the cone. In this manner, the inter-fiber angle between all 3 fibers are identical.
Figure 5.1 shows the configuration of the cone.
The neighborhood configuration described in Figure 4.4 is used for ICA estimation. 11 syn-
thetic diffusion signals ”measurements” are generated. The program was run on Intel Xeon
Quad-core 2.66GHz workstation.
45
−1 −0.5 0 0.5 1
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
y
z
(a) Side view of the cone. It is a isosceles triangle.
−1 −0.5 0 0.5 1
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
x
y
(b) Base of the cone. It is a circle.
(c) 3D view of the cone
Figure 5.1: The cone used to generated 3 fibers with identical inter-fiber angles. The radius of
the circle and the height of the isosceles triangle are adjusted to vary the inter-fiber angle.
46
5.2 Estimating orientations of individual fiber per voxel
In this section, the accuracy of estimating fiber orientations are examined. Four test cases are
investigated:
Case 1: The number of source fibers and estimated fibers is identical.
Case 2: The number of source fibers and estimated fibers is identical. Cerebral Spinal Fluid (CSF)
component with isotropic diffusion is added into the mixture. The CSF component is not
counted as a fiber.
Case 3: The number of source fibers and estimated fibers is identical. For each source fiber, 3 of the
10 neighbors are set to zero weighting. This case is to simulate when a fiber does not span
the entire 10-voxel neighborhood.
Case 4: The number of source fibers and estimated fibers is different.
For each case, the mean and standard deviation of the error angle is calculated as follows:
1. For each trial, there are n
src
source fibers and n
est
estimated fibers.
2. For each of the source fibers, a closest match with the smallest angle is chosen among the
estimated fiber.
3. n
src
angles are obtained and they are counted as n
src
samples for that trial.
4. Combine all the results in 1000 trials and obtain 1000*n
src
data points.
5. Calculate mean and standard deviation.
The percentage of all recorded data points in which the absolute error angle is less than 15
degrees is also calculated. This is similar to the study in [Ale05c], where they use a correlation
coefficient of 0.95, which corresponds to an error angle of 18.2 degrees.
47
In each study, the ICA approach was compared to orientations estimated from Behrens model
using multi-parameter estimation but without isotropic diffusion as described in Equation 4.4.
The diffusivity d is fixed to
1
b
,where b is the b− value of diffusion data to reduce the number of
variables in the model.
5.2.1 case 1: number of original sources equals number of estimated
sources without CSF
Figures 5.2 to 5.7 show the simulation results for 2 and 3-fiber cases with SNR=15, 30 and 45
with equal number of source fibers and estimated fibers. There is no CSF (isotropic diffusion)
component in this study. For 2-fiber case, there is siginificant improvement from SNR=15 to
SNR=30, but not from SNR=30 to SNR=45. For SNR=30 and SNR=45, the error angle is
around 8 degrees with inter-fiber angle 40-50 degrees and is less than 5 degrees with inter fiber
angle 50 degrees and above. Over 90 percentage of trials with error angle less than 15 degrees are
observed with inter fiber angle 40 degrees and above for SNR=30 and SNR=45.
For 3-fiber case, the mean error and the standard deviation becomes larger. For instance, for
SNR=30, the mean error for inter-fiber angle between 80 and 90 is 20 degrees and the standard
deviation is around 15 degrees. The best performance region is with inter-fiber angles between
30-70 degrees, where the mean error angle is between 10-15 degrees. The worse performance at
inter-fiber angle larger than 70 degrees is probably because of its isotropic nature and it becomes
more susceptible to noise.
In Figure 5.5, there is an extra plot where among the three estimated fiber directions, the two
best estimates are retained and the worst one is eliminated. It shows that the accuracy improves
significantly. This shows that in 3-fiber case, with inter-fiber angle above 40 degrees, ICA is able
to recover two of the fiber directions accurately (with mean error less than 10 degrees), while the
third direction has a higher error.
48
When compared to the Behrens model, ICA performs better for 2-fiber case in all inter-fiber
angles. However, in the case of three fibers, ICA performs worse than Behrens model by around
5 degrees with inter-fiber angle larger than 70 degrees. The computational speed for ICA is 0.005
second per voxel while the computational speed for Behrens model is 0.5 second per voxel.
5.2.2 case 2: number of original sources equals number of estimated
sources with CSF
In this study, a CSF component is added into the synthetic diffusion signal. The CSF component is
not counted towards the number of source fibers. ICA is then used to estimate the fiber directions.
It is observed that there is no significant difference with the case without CSF.
5.2.3 case 3: number of original sources equals number of estimated
sources with partial zero weighting
In case 3, we simulate the scenario when each fiber cross the voxel of interest (the center voxel),
but not all of its neighbors. This is the case when the voxel is at the boundary of a fiber bundle.
In this study, three of the 10 neighbors are set to zero weighting for each fiber. CSF component
is not added.
Figure 5.9 shows some interesting results. The performance with zero weighting is better than
case 1 where all weightings are non-zero. The percentage of trials with error angle less than 15
degrees with zero weighting is above 95% with inter-fiber angle between 40 and 60 degrees, while
it is less than 80% without zero weighting.
One possible explanation is that zero weighting enhances the mixing contrasts of the sources
enabling ICA to obtain better results in maximizing non-Gaussianity of measurements.
49
10 20 30 40 50 60 70 80 90
0
5
10
15
20
25
30
35
40
45
Error angle (degree)
Inter Fiber Angle (degree)

Behrens model
ICA
(a)
10 20 30 40 50 60 70 80 90
0
20
40
60
80
100
% of trials with error angle < 15
0
Inter Fiber Angle (degree)

Behrens model
ICA
(b)
Figure 5.2: Simulation result with 2 source fibers and 2 estimated fibers and SNR=15. It is
observed that the error angle for ICA is around 5 degrees less than using Behrens model. In
ICA, Over 80% of trials have angle error less than 15 degrees with inter-fiber angle larger than
50 degrees. In Behrens model, the percentage is less than 80%.
50
10 20 30 40 50 60 70 80 90
0
5
10
15
20
25
30
35
40
45
Error angle (degree)
Inter Fiber Angle (degree)

Behrens model
ICA
(a)
10 20 30 40 50 60 70 80 90
0
20
40
60
80
100
% of trials with error angle < 15
0
Inter Fiber Angle (degree)

Behrens model
ICA
(b)
Figure 5.3: Simulation result with 3 source fibers and 3 estimated fibers and SNR=15. It is
observed that the error angle for ICA is around 15-20 degrees with inter-fiber angle between 30
and 70 degrees, which is less than Behrens model. The error angle is higher than that of the
2-fiber case.
51
10 20 30 40 50 60 70 80 90
0
5
10
15
20
25
30
35
40
45
Error angle (degree)
Inter Fiber Angle (degree)

ICA
Behrens model
(a)
10 20 30 40 50 60 70 80 90
0
20
40
60
80
100
% of trials with error angle < 15
0
Inter Fiber Angle (degree)

Behrens model
ICA
(b)
Figure 5.4: Simulation result with 2 source fibers and 2 estimated fibers and SNR=30. It is
observed that the error angle for ICA is less than 5 degrees while Behrensmodel haserrorangle
around 10 degree range with inter-fiber angle larger than 50 degrees. In ICA, Over 80% of trials
have angle error less than 15 degrees with inter-fiber angle larger than 40 degrees.
52
10 20 30 40 50 60 70 80 90
0
5
10
15
20
25
30
35
40
45
Error angle (degree)
Inter Fiber Angle (degree)

Behrens model
ICA
ICA, but ignore the fiber with largest error
(a)
10 20 30 40 50 60 70 80 90
0
20
40
60
80
100
% of trials with error angle < 15
0
Inter Fiber Angle (degree)

Behrens model
ICA
(b)
Figure 5.5: Simulation result with 3 source fibers and 3 estimated fibers and SNR=30. It is
observed that the error angle for ICA is around 10-15 degrees with inter-fiber angle between 30
and 70 degrees, which is less than Behrens model. The error angle is higher than that of the
2-fiber case. The green curve shows the result in that in each trial, the fiber with the largest error
is eliminated. It is observed that the mean error angles significantly drop. This shows that ICA
can reliably recover 2 of the 3 sources.
53
10 20 30 40 50 60 70 80 90
0
5
10
15
20
25
30
35
40
45
Error angle (degree)
Inter Fiber Angle (degree)

Behrens model
ICA
(a)
10 20 30 40 50 60 70 80 90
0
20
40
60
80
100
% of trials with error angle < 15
0
Inter Fiber Angle (degree)

Behrens model
ICA
(b)
Figure 5.6: Simulation result with 2 source fibers and 2 estimated fibers and SNR=45. It is
observed that the result is similar to that of SNR=30, except with a few degrees of improvement.
54
10 20 30 40 50 60 70 80 90
0
5
10
15
20
25
30
35
40
45
Error angle (degree)
Inter Fiber Angle (degree)

Behrens model
ICA
(a)
10 20 30 40 50 60 70 80 90
0
20
40
60
80
100
% of trials with error angle < 15
0
Inter Fiber Angle (degree)

Behrens model
ICA
(b)
Figure 5.7: Simulation result with 3 source fibers and 3 estimated fibers and SNR=45. It is
observed that the result is similar to that of SNR=30, except with a few degrees of improvement.
55
10 20 30 40 50 60 70 80 90
0
5
10
15
20
25
30
35
40
45
Error angle (degree)
Inter Fiber Angle (degree)

2 fibers
3 fibers
(a)
10 20 30 40 50 60 70 80 90
0
10
20
30
40
50
60
70
80
90
100
% of trials with error angle < 15
0
Inter Fiber Angle (degree)

2 fibers
3 fibers
(b)
Figure 5.8: Simulation result with a CSF component added into synthetic diffusion data and ICA
is used to recover the components. SNR=30 is used in this study. The number of source fibers
is identical to the number of estimated fibers. The CSF component is not counted as one fiber
source. It is observed that the performance is similar to the previous case where there is no CSF
component.
56
10 20 30 40 50 60 70 80 90
0
5
10
15
20
25
30
35
40
45
Error angle (degree)
Inter Fiber Angle (degree)

2 sources
3 sources
(a)
10 20 30 40 50 60 70 80 90
0
10
20
30
40
50
60
70
80
90
100
% of trials with error angle < 15
0
Inter Fiber Angle (degree)

2 sources
3 sources
(b)
Figure 5.9: Simulation result with zero weighting assigned to three voxels for each fiber source.
SNR=30 is used in this study. This is to simulate the scenario where a fiber is not spanning the
entire neighborhood. It is observed that the result is better than the case without zero weighting.
It is because zero weighting enhances the weighting contrast among different sources.
57
5.2.4 case 4: number of original sources does not equal number of
estimated sources
Case 4, the scenario where the estimation of number of fiber directions is incorrect. In addition
to numerical figures, visual results are given showing how mis-estimation of the number of fibers
affects the estimated fiber directions.
The results are summarized as follows:
• Figure 5.10 shows the case with one source fiber. If there are 2 or 3 estimated fibers, the
original fiber could be recovered with the extra components being noise as shown in Figure
5.13.
• Figure 5.11 shows the case with two source fibers. If there are 3 estimated fibers, the original
2 fibers could be recovered with the third component being noise as shown in Figure 5.14.
• Figure 5.12 shows the case with three source fibers. If there are 2 estimated fibers, the fibers
somewhat lie between the three fibers as shown in Figure 5.15.
Please note that if there is one estimated source, PCA is used to calculate the single fiber
direction. In conclusion, if the number of estimated sources is larger than original sources, the
original sources are still recovered with the extra components being noise. If the number of
estimated sources is smaller than original sources, the estimated directions lie in-between the
original fiber directions.
5.3 Estimating number of fiber directions per voxel
The accuracy of estimating the number of fiber directions is evaluated in this section based on
the statistical hypothesis testing procedure described in Chapter 4. Figure 5.16 shows the results
with p-value=0.1. SNR=30 is used.
58
10 20 30 40 50 60 70 80 90
0
5
10
15
20
25
30
35
40
45
Error angle (degree)
Inter Fiber Angle (degree)

2 estimated sources
3 estimated sources
Figure 5.10: Simulation result where there is one fiber source, and with two and three estimated
sources. SNR=30 is used. It is observed that the original fiber is estimated with less than 5 degree
error in both cases.
59
10 20 30 40 50 60 70 80 90
0
5
10
15
20
25
30
35
40
45
Error angle (degree)
Inter Fiber Angle (degree)

2 estimated sources
3 estimated sources
Figure 5.11: Simulation result where there are two fiber sources, and with two and three estimated
sources. SNR=30 is used. It is observed that in both cases the error angle is very similar except
that in 3-fiber case, the extra fiber is erroneous.
10 20 30 40 50 60 70 80 90
0
5
10
15
20
25
30
35
40
45
Error angle (degree)
Inter Fiber Angle (degree)

2 estimated sources
3 estimated sources
Figure 5.12: Simulation result where there are three fiber sources, and with two and three esti-
mated sources. SNR=30 is used. It is observed that in both cases the error angle is larger than
10 degrees.
60
(a) 2 estimated fibers (b) 3 estimated fibers
Figure 5.13: Samples for 1 source fiber. It is observed that the original fiber source is reliably
estimated. However, extra estimated fibers are noise-like.
(a) 2 estimated fibers (b) 3 estimated fibers
Figure 5.14: Samples for 2 source fibers. Similarly, it is observed that the original fiber sources
are reliably estimated.
61
(a) 2 estimated fibers (b) 3 estimated fibers
Figure 5.15: Samples for 3 source fibers. It is observed that the original fiber sources are reliably
estimated in the 3 estimated fibers case. In the case where only two fibers are estimated, they lie
in-between the three original fiber sources.
62
10 20 30 40 50 60 70 80 90
0
20
40
60
80
100
True positive in %
Inter Fiber Angle (degree)

two fibers
three fibers
multiple fibers
Figure 5.16: Number of fiber directions obtained using F-test with p-value 0.1 and SNR=30. It
is observed that the true positive percentage for single fiber is above 80%, while for 2-fiber case,
the percentage is above 80% with inter fiber angle larger than 30 degrees. The true positive
percentage for three fibers is below 50%. However, it is observed that the true positive percentage
for multiple fibers is between 70-80% with inter-fiber angle larger than 50%. This implies that a
significant amount of three fiber voxels are classified as two fibers.
63
It is observed that the accuracy (true positive percentage) of estimating 1-2 fibers is around
or above 80% with inter-fiber angle above 30 degrees. The accuracy for three fibers is below
50%. However, it is observed that the true positive percentage for estimating multiple fibers
per voxel is around or above 70%. This shows that a significant amount of 3 fiber voxels are
being estimated as 2. However, in the phantom simulation study in the next Chapter, our F-
test method works reasonably well and our approach reveals significant improvement over DTI
and multi-compartment approaches. Experimental results with human data also shows that our
approach is capable of recovering crossing fibers that cannot be separated by PCA-based DTI.
64
Chapter 6
Simulation study with synthetic phantom
In this Chapter, a Common DTI dataset created by [DG] is used to validate our approach for
multiple fiber tractography. In particular, the orthogonal crossing and 60 degrees crossing con-
figurations are used. In this dataset, the diffusion MRI data of the phantoms are acquired with
30 gradient directions using an algorithm analogous to electrostatic repulsion [JHS99]. b−value
of 1000 is used with TE=90ms. A T
2
-weighted image is also used. The T
2
of fiber tracts and
background are 65ms and 95ms respectively. The background is a simulation for grey matter
mixed with the tracts. In our studies, ROIs are drawn such that only those voxels with white
matter information are captured. Data with SNR=15 and SNR=30 are used.
6.1 Orthogonal crossing
Figure 6.1 shows one b
0
slices with 90 degree crossing.
6.1.1 Fiber reconstruction
Figure 6.2 shows a slice of FA map. There are two white matter tracts with the horizontal tract
having a lower FA value than the vertical one. The surrounding low FA regions are grey matter.
There exists mixing between white and grey matters in the boundaries.
65
Figure 6.1: T
2
-weighted slice obtained from the b
0
diffusion MRI set with SNR=30.
Figure 6.2: FA map obtained from the diffusion MRI set for one slice. A mask is drawn to include
only the white matter voxels.
66
One fiber per voxel
Two fibers per voxel
Three fibers per voxel
Figure 6.3: Color index for representing 1, 2 and 3 fiber directions per voxel.
Figure 6.4: Map showing number of fiber directions per voxel with SNR=30 and p−value= 0.005.
It is observed that F-test tends to overestimate the number of fibers in this case.
Figure 6.3 shows the index used to represent different numbers of fiber directions per voxel,
which are used in all studies.
Figures 6.4 to 6.6 show the map of number of fiber directions with two p-values at SNR=30.
It is observed that a lower p-value causes the F-test to choose a model with less fiber directions.
In both results, the fiber crossing region is cleared defined, showing two fiber directions per voxel.
Figure 6.7 shows the map of number of fiber directions at SNR=15. In the following studies,
p=0.001 is chosen for fiber reconstruction and tractography. The effects of overestimation and
underestimation of number of fiber directions, which have been described in the previous Chapter,
are examined here on tractography.
67
Figure 6.5: Map showing number of fiber directions per voxel with SNR=30 and p−value= 0.001.
Compared to previous case, the single fiber braches are more clearly identified.
Figure 6.6: Map showing number of fiber directions per voxel with SNR=30 and p−value= 0.0005.
It is observed that the lower the p-value, F-test tends to select model with fewer fiber directions.
p−value= 0.001 seems to offer the best compromise and it is adopted in subsequent studies.
Figure 6.7: Map showing number of fiber directions per voxel with SNR=15 and p−value= 0.001.
68
Figure 6.8 shows the estimated fiber directions vectors for different voxels with SNR=30. From
the Figures, it is trivial to realize there are two fiber bundles with high FA (white matter) and one
background fiber direction with low FA (grey matter). Observations are summarized as follows:
• In PCA approach, all fibers are modeled with a rank-2 tensor. When a voxel contains more
than one fiber, the estimated direction lies between the multiple directions, or the fiber
bundle with higher FA dominates.
• In Behrens model, all voxels with 3 fiber directions fail to converge. Such voxels are assigned
2 fiber directions and the corresponding results are also shown. It is observed that the
method could estimate the correct number of fiber directions for the 90 degrees crossing
between the two white matter bundle, but not the case for white and grey matter crossing.
• In the ICA approach, the directions are correctly estimated if the number of fibers per voxel
is correct. In voxels with overestimated number of fibers, the extra components are random.
In voxels with underestimation, the components lie between the source fiber directions.
Figure 6.9 shows the estimated fiber directions vectors for different voxels with SNR=15. The
interpretation is similar to the case with SNR=30, except that the estimated fiber directions are
visibly noisier.
6.1.2 Tractography results
Figures 6.10 and 6.11 show the tractography results with SNR=30 and SNR=15 respectively.
Tracts that successfully originates from one branch and end at the other end are displayed. The
tracts are obtained by filtering with two ROIs. Tracts in different colors are filtered with ROIs in
the corresponding color. The observations are summarized as follows:
• In PCA approach, it is observed that the horizontal white matter tract (with lower FA than
the other white matter tract) failed to pass through the crossing region. The vertical tract
is thinner compared to the ICA approach.
69
(a) PCA (b) ICA
(c) Behrens model (d) Behrens model with unconverged voxels
assigned 2 fibers
Figure 6.8: Reconstructed fiber directions with 90 degrees crossing and SNR=30. It is observed
that PCA failed to resolve the fiber crossing directions. In Behrens model, although the fiber cross-
ing is resolved, it is observed that the model did not converge for many voxels. The background
fiber direction is not properly recovered. In ICA, the fiber directions are reliably recovered.
70
(a) PCA (b) ICA
(c) Behrens model (d) Behrens model with unconverged voxels
assigned 2 fibers
Figure 6.9: Reconstructed fiber directions with 90 degrees crossing and SNR=15. Observations
are similar to previous case, except that the error is slightly worse than SNR=30.
71
• In Behrens model, the two white matter tracts are recovered. However, the horizontal tract
is not perfectly horizontal. There is an observable tilt from the horizontal direction.
• In the ICA approach, both tracts are recovered. In particular, the vertical tract is the widest
among all methods and it spans into regions that are mixed with grey matter. The tilting is
significantly less than using Behrens model. When comparing ICA approach using SNR=15
and SNR=30, the tracts for SNR=30 is slightly thicker and the tilting in SNR=15 is worse.
Figures 6.12 and 6.13 show the tractography results with SNR=30 and SNR=15 respectively.
In this case, we investigated how tracts propagate into the crossing area from one branch. Obser-
vations are summarized as follows:
• In PCA approach, it is observed that the horizontal white matter tract is diverted into the
vertical branch. This is due to the partial volume problems of the two 90-degree tracts.
Similarly, a significant portion of the vertical tract is diverted into horizontal branch.
• In Behrens model, majority of the two white matter tracts could pass through the fiber
crossing region. However, there is a significant amount of distortion in tract direction inside
the crossing region.
• In the ICA approach, majority of the tracts could pass through the fiber crossing region. Due
to the multiple fiber tractography algorithm, some of the white matter tracts are diverted
into grey matter ones.
6.2 60 degrees crossing
Figure 6.14 shows one b
0
slices with 60 degree crossing.
72
(a) PCA
(b) ICA (c) Behrens model
Figure 6.10: Tractography results with 90 degrees crossing and SNR=30. The tracts that pass
through the crossing region are shown. It is observed that in PCA, the tract fails to propagate
through the crossing region. In Behrens model, there is some tilting in the tracts even though
they could propagate through the crossing region. It seems that ICA offers the best tractography
result.
73
(a) PCA
(b) ICA (c) Behrens model
Figure 6.11: Tractography results with 90 degrees crossing and SNR=15. The tracts that pass
through the crossing region are shown. Similar results are observed as in SNR=30, except that
there exists slightly more error in tracts.
74
(a) PCA
(b) ICA (c) Behrens model
Figure 6.12: Tractography results with 90 degrees crossing and SNR=30. The tracts that originate
from one branch are shown. It is observed that the tracts in PCA are diverted to the wrong branch.
In ICA and Behrens model, the tracts could pass through the crossing region.
75
(a) PCA
(b) ICA (c) Behrens model
Figure 6.13: Tractography results with 90 degrees crossing and SNR=15. The tracts that originate
from one branch are shown. Similar observations are obtained as in SNR=30 case.
76
Figure 6.14: T
2
-weighted slice obtained from the b
0
diffusion MRI set with SNR=30.
6.2.1 Fiber reconstruction
Figure 6.15 shows a slice of FA map. There are two white matter tracts crossing at 60 degrees.
One of the branches is thicker with a higher FA than the other branch. Similar to the 90 degree
crossing case, there exists mixing between white and grey matters in the boundaries.
Figures 6.16 to 6.18 show the map of number of fiber directions with two p-values at SNR=30.
It is observed that a lower p-value causes the F-test to choose a model with less fiber directions.
Figure 6.7 shows the map of number of fiber directions at SNR=15. The effect of noise in this
scenario causes the crossing region to disappear in that, the voxels in the crossing region are
classified as single fiber. If the fiber map with SNR=15 is used, it is obvious that no tract would
be able to pass through the fiber crossing region successfully. As a result, in the case for SNR=15,
the number of fiber map of SNR=30 is used, while the fiber directions are reconstructed using data
of SNR=15. The rationale is to investigate the effect of noise in the accuracy of estimating fiber
direction. In the following studies, p=0.001 is chosen for fiber reconstruction and tractography.
77
Figure 6.15: FA map obtained from the diffusion MRI set for one slice. A mask is drawn to
include only the white matter voxels.
The effects of overestimation and underestimation of number of fiber directions are described in
previous Chapter and in this chapter, such effects on tractography are examined.
Figure 6.8 shows the estimated fiber directions vectors for different voxels with SNR=30. From
the Figures, it is trivial to realize there are two fiber bundles with high FA (white matter) and one
background fiber direction with low FA (grey matter). Observations are summarized as follows:
• In PCA approach, all fibers are modeled with a rank-2 tensor. When a voxel contains more
than one fiber, the estimated direction lies between the multiple directions. Specifically, all
fibers in the crossing regions are exactly the mean of the two branching directions.
• In Behrens model, all voxels with 3 fiber directions fail to converge. Such voxels are assigned
2 fiber directions and the corresponding results are also shown. It is observed that the
method failed to reconstruct the crossing fibers.
78
Figure 6.16: Map showing number of fiber directions per voxel with SNR=30 and p−value= 0.005.
It is observed that F-test tends to overestimate the number of fibers in this case.
Figure 6.17: Map showing number of fiber directions per voxel with SNR=30 and p−value= 0.001.
It is observed that the branches are more clearly idenitified than previous case.
79
Figure 6.18: Map showing number of fiber directions per voxel with SNR=30 and p−value=
0.0005. It is observed that the lower the p-value, F-test tends to select model with lower fiber
directions. p−value= 0.001 seems to offer the best compromise and it is adopted in subsequent
studies.
Figure 6.19: Map showing number of fiber directions per voxel with SNR=15 and p−value= 0.001.
80
• In the ICA approach, the 60 degree crossing directions are observed even though the crossing
region is narrow. There are a few voxels with overestimated number of fibers (red). In the
case of SNR=15, there is a noticeable amount of error in the central crossing region.
6.2.2 Tractography results
Figures 6.22 and 6.23 show the tractography results with SNR=30 and SNR=15 respectively.
Tracts that successfully originate from one branch and ends at the other end are displayed. The
tracts are obtained by filtering with two ROIs. Tract in different colors are filtered with ROIs in
the corresponding color. The observations are summarized as follows:
• In PCA approach, it is observed that none of the tracts are correctly recovered. For the red
(thicker/higher FA) tract, there is significant amount of distortion due to the partial volume
problems.
• In Behrens model, only the red tract is recovered. The blue (thinner/lower FA) tract fails
to pass through the crossing region, which corresponds to the result of the reconstructed
fibers as shown in Figure 6.20.
• In the ICA approach, both tracts are recovered. There is some distortion in the fiber crossing
region. Although the thinner tract with lower FA is recovered, it is thinner than expected.
When comparing ICA approach using SNR=15 and SNR=30, the tracts for SNR=30 is
slightly thicker.
Figures 6.24 and 6.25 show the tractography results with SNR=30 and SNR=15 respectively.
In this case, we investigated how tracts propagate into the crossing area from the branch with
tract of lower FA. Observations are summarized as follows:
• In PCA approach, it is observed that the entire tract is diverted into the wrong branch.
• In Behrens model, the same happens. There is no tract that could pass the crossing region.
81
(a) PCA (b) ICA
(c) Behrens model (d) Behrens model with unconverged voxels
assigned 2 fibers
Figure 6.20: Reconstructed fiber directions with 60 degrees crossing and SNR=30. It is observed
that in ICA, the individual 60 degree crossing fiber directions are recovered. In PCA and Behrens
model, the crossing directions are incorrect.
82
(a) PCA (b) ICA
(c) Behrens model (d) Behrens model with unconverged voxels
assigned 2 fibers
Figure 6.21: Reconstructed fiber directions with 60 degrees crossing and SNR=15.
83
(a) PCA
(b) ICA (c) Behrens model
Figure 6.22: Tractography results with 60 degrees crossing and SNR=30. The tracts that pass
through the crossing region are shown. In this study, only ICA is capable of recovering the weaker
tract in blue.
84
(a) PCA
(b) ICA (c) Behrens model
Figure 6.23: Tractography results with 60 degrees crossing and SNR=15. The tracts that pass
through the crossing region are shown.
85
• In the ICA approach, portion of the tracts could pass through the fiber crossing region.
A major portion of the tract is still diverted into the wrong branch because of the severe
partial volume problems.
86
(a) PCA
(b) ICA (c) Behrens model
Figure 6.24: Tractography results with 60 degrees crossing and SNR=30. The tracts that originate
from one branch are shown. It is observed that the stronger tract dominates in the tractography
result. Behrens and PCA failed to allow the weaker tract to propagate through the fiber crossing
region.
87
(a) PCA
(b) ICA (c) Behrens model
Figure 6.25: Tractography results with 60 degrees crossing and SNR=15. The tracts that originate
from one branch are shown.
88
Chapter 7
Experimental human studies of ICA-based tractography
A direct consequence of multi-fiber per voxel tractography is the reduction of partial volume
errors. Experimental Diffusion MRI data acquired from healthy human subjects were processed
using our ICA approach and compared to the PCA-based DTI approach to show a reduction in the
partial volume problems with ICA based tractography. Results from two subjects are presented
here. Comparison of the ICA approach to Behrens model is omitted in this study because of
its computational complexity (0.5 second per voxel versus 0.005 second per voxel in ICA). The
algorithm described in 4.3 is used. Frontal-occipital tracts and the cingulum are investigated as
examples where partial volume errors coupd play a major role. Figure 7.1 shows the index used to
representing different numbers of fiber directions per voxel, as well as the color-coded directional
information used in drawing fiber vectors.
In this study, diffusion MRI data are acquired using a 3T GE MRI scanner with 25 gradient
directions with b=1000s/mm
2
, 2 averages, 28 slices with voxel size 2× 2× 4mm
3
, TR=8000ms,
and TE=86.1ms. One set of data is also acquired with b=0.
89
One fiber per voxel
Two fibers per voxel
Three fibers per voxel
(a) color index for 1-
3fiberdirections
(b) color sphere for showing directional
information for fiber vectors
Figure 7.1: Index information for our experimental study.
7.1 Example 1
7.1.1 Map for number of fiber directions
Figure 7.2 and 7.3 show the map of number of fiber directions and the corresponding FA map. It
is observed that the single fiber per voxel (red) matches well with high FA voxels.
7.1.2 Frontal-occipital
Figure 7.4 shows the reconstructed fiber directions from the ROI shown in the FA map. It is
observed that the partial volume effects due to crossing fibers are mostly resolved using our ICA
approach. Highlighted regions are pointed with arrows.
Figures 7.5 and 7.6 show the tractography results for the frontal-occipital tracts. The cor-
responding filtering ROIs are shown in green. Only tracts that pass through all specified ROIs
are retained. It is observed that there are severe partial volume problems between the ventricle,
corpus callosum and frontal-occipital in PCA approach. Most of the tracts are diverted into the
corpus callosum on the left side of the brain. The problem is greatly reduced in our ICA approach.
It is observed that most tracts propagate to the occipital region.
90
Figure 7.2: Map of number of fiber directions with p−value= 0.001. It is observed that there are
various regions with multiple fibers.
91
Figure 7.3: Corresponding FA map. It is observed that the high FA value voxels correspond well
with single fiber voxels including the corpus callosum.
92
(a) axial slice 13 of the FA map
(b) PCA (c) ICA
Figure 7.4: Reconstructed fiber directions. Fiber directions are drawn within the ROI shown in
the FA map. It is known that corpus callosum, frontal occipital and a few other minors tracts
cross in this region. It is observed that ICA is capable of reducing the partial volume problems
between fibers. In ICA result, it is observed that the Frontal-Occipital vectors are recovered as
pointed by the arrow. In PCA case, the tract is dominated by the corpus callosum.
93
(a) PCA (b) ICA
Figure 7.5: Tractography results showing tracts passing through both ROIs. It is observed that
in PCA the frontal occipital tracts failed to propagate into the occipital area. In ICA, there is
significant improvement.
(a) PCA (b) ICA
Figure 7.6: Tractography results showing tracts passing through both ROIs.
94
7.1.3 Cingulum
Figures 7.7 and 7.8 show the reconstructed fiber vectors from the ROI shown in the FA map. It
is observed that the partial volume effects due to crossing fibers between the cingulum and other
tracts including the corpus callosum are mostly resolved using our ICA approach. Highlighted
regions are pointed with arrows.
Figures 7.10 and 7.10 show the tractography results of the tracts passing through both specified
ROIs. It is observed that there are significant partial volume problems with the corpus callosum
in PCA approach, seen in the axial view. In our ICA approach, the partial volume problems
disappears. The cingulum reconstruction is significantly better using our ICA approach in that
the cingulum extends into the hippocampus. The cingulum reconstructed using PCA approach
is cut pre-maturely.
7.2 Example 2
7.2.1 Map for number of fiber directions
Figure 7.2 and 7.3 show the map of number of fiber directions and the corresponding FA map. It
is observed that the single fiber per voxel (red) matches well with high FA voxels.
7.2.2 Frontal-occipital
Figure 7.13 shows the reconstructed fiber directions from the ROI shown in the FA map. It is
observed that the partial volume effects due to crossing fibers are mostly resolved using our ICA
approach. Highlighted regions are pointed with arrows.
Figures 7.14 and 7.15 show the results for the frontal-occipital tracts for the second subject.
The partial volume effect with the ventricle and corpus callosum in PCA approach are significantly
95
(a) sagittal slice 62 of FA map
(b) PCA
(c) ICA
Figure 7.7: Reconstructed fiber directions for cingulum in sagittal view. Fiber directions are
drawn within the ROI shown in FA map. It is known that corpus callosum and the cingulate
gyrus tracts cross in this region. In ICA result, it is observed that the individual fiber directions
are recovered as pointed by arrows. In PCA, the fiber is lying between the two directions.
96
(a) coronal slice 75 of FA map
(b) PCA
(c) ICA
Figure 7.8: Reconstructed fiber directions for cingulum in coronal view. Similar observation
is obtained in that the individual directions between corpus callosum and cingulate gyrus are
recovered.
97
(a) PCA (b) ICA
Figure 7.9: Tractography results showing tracts passing through both ROIs in axial view. It is
observed that in PCA, there is a significant amount of tract diverted into the corpus callosum
where the problem disappears in ICA.
(a) PCA (b) ICA
Figure 7.10: Tractography results showing tracts passing through both ROIs in sagittal view.
There is a significant reduction in partial volume problems between corpus callosum and cingulate
gyrus. It is observed that in ICA, the cingulate could propagate into the hippocampus.
98
Figure 7.11: Map of number of fiber directions with p−value= 0.001.
99
Figure 7.12: Corresponding FA map showing that high FA value voxels correspond well with
single fiber voxels including the corpus callosum.
100
(a) axial slice 13 of FA map
(b) PCA (c) ICA
Figure 7.13: Reconstructed fiber directions. From the pointed arrows, it is observed that ICA is
capable of recovering multiple fiber directions in fiber crossing regions.
101
(a) PCA (b) ICA
Figure 7.14: Tractography results showing tracts passing through both ROIs. It is observed that
the partial volume problems between corpus callosum and frontal occipital tracts is significantly
reduced using ICA.
worse than our ICA approach. In our ICA approach, there are very few tracts that are diverted
into the corpus callosum. The branching at the right frontal region is better using ICA.
7.2.3 Cingulum
Figure 7.7 shows the reconstructed fiber vectors from the ROI shown in the FA map. It is
observed that the partial volume effects due to crossing fibers between the cingulum and other
tracts including the corpus callosum are mostly resolved using our ICA approach. Highlighted
regions are pointed with arrows.
Figure 7.17 shows the tractography results of the cingulum. Similarly, the tracts reconstructed
using ICA are significantly better than PCA in that the tracts extend into the hippocampus.
102
(a) PCA (b) ICA
Figure 7.15: Tractography results showing tracts passing through both ROIs, plus a filter removing
tracts passing through the middle-line.
103
(a) sagittal slice 70 of FA map
(b) PCA
(c) ICA
Figure 7.16: Reconstructed fiber directions for cingulum in sagittal view. Fiber directions are
drawn within the ROI showin in the FA map. It is known that corpus callosum and the cingulate
gyrus tracts cross in this region. Similar to the first example, it is observed that the individual
fiber directions are recovered as pointed by arrows. In PCA, the fiber is lying between the two
directions.
104
(a) PCA (b) ICA
Figure 7.17: Tractography results showing tracts passing through both ROIs in sagittal view. It
is observed that the cingulate gyrus tract propagates towards the hippocampus while PCA failed
to do so.
105
Chapter 8
Conclusions and future research work
8.1 Conclusion
In this thesis, a solution for the Diffusion MRI fiber crossing problem is proposed using Indepen-
dent Component Analysis (ICA) and statistical hypothesis testing. One of the major strengths of
our approach is that routine clinically acquired data, where diffusion data are typically acquired
in less than 10 minutes, can be used.
A multi-compartment model proposed by Behrens et. al. [BWJ
+
03] is used to estimate the
number of fibers per voxel. Instead of estimating all model parameters by non-linear optimization,
ICA is used to estimate the directions of multiple fibers. Application of ICA requires us to
investigate neighbors of the voxel of interest. We assume that a fiber bundle spans several voxels.
This implies that neighboring voxels carry similar diffusion signal with the same mixed fibers as
the current voxel, but with a different mixture weighting. As a prerequisite of ICA is that the
”sources” must be non-Gaussian, we investigate the statistical properties of the diffusion profile
for single fiber voxels to establish non-Gaussianity.
By obtaining fiber directions through ICA, the Behrens model is then simplified into a linear
one in which least square estimation can be used to calculate residue in fitting the diffusion signal.
F-test is then used for model selection between 1, 2 or 3 fiber directions.
106
After fiber directions calculation, whole brain multiple fiber streamline tractography is con-
ducted. Various tracts of interest are obtained by filtering with different ROIs.
Monte Carlo simulation studies are conducted at SNR=15, 30 and 45. Results show that in
2-fibers case, the error angle is less than 5 degrees with SNR=30 or above with a inter-fiber angle
fiber larger than 40 degrees. For 3-fibers case, the error angel is less than 15 degrees with a
inter-fiber angle between 30 and 70 degrees. Various scenarios are taken into consideration:
• with the addition of CSF component,
• with zero weightings for different fibers, and
• when the number of source fibers is either more or less than the number of estimated fibers.
Our ICA approach is then validated using a common DTI dataset developed by ISMRM
diffusion and perfusion study group with diffusion phantom data at SNR=15 and 30. Orthogonal
and 60 degree fiber crossings are investigated. Results show that our approach works significantly
better than PCA-based DTI as well as multi-parameter estimation using Behrens model.
Finally, our approach is applied to experimental human brain diffusion MRI data and trac-
tography results are compared using PCA approach. Results show that our method is capable of
significantly reducing the partial volume problems caused by mixing among fibers and with CSF,
which is a direct consequence of being able to identify multiple fibers within a voxel.
8.2 Future research work
Possible extension of current work is highlighted in this section:
• Performance of ICA at higher number of gradient directions. How does an increase of
number of gradient directions affect the performance of ICA?
• Comparison with other multiple fiber approaches including Q-ball at HARDI settings.
107
• The p-value in our model selection process is a variable. By increasing the p-value, higher
order model is more likely to be selected. Effects of varying p-value on tractography results
are a natural extension to current work.
• Quantitation of multiple fibers per voxel. There are various approaches to quantify tensor
anisotropy and single fiber tractography. However, multiple fiber quantification is still an
open research topic. For instance, Diffusion Kurtosis Imaging (DKI) is a recent approach
to characterize the non-Gaussian property of diffusion per voxel.
• ICA-based multiple fiber tractography could be applied to quantitative analysis and group
studies to improve the statistical significance of differences between groups.
108
References
[AAC04] Liat Avram, Yaniv Assaf, and Yoram Cohen. The effect of rotational angle
and experimental parameters on the diffraction patterns and micro-structural
information obtained from q-space diffusion nmr: implication for diffusion in
white matter fibers. J Magn Reson, 169(1):30–38, Jul 2004.
[ABA02] D. C. Alexander, G. J. Barker, and S. R. Arridge. Detection and modeling of
non-gaussian apparent diffusion coefficient profiles in human brain data. Magn
Reson Med, 48(2):331–340, Aug 2002.
[ABBC
+
02] Y. Assaf,D. Ben-Bashat,J.Chapman,S. Peled,I.E.Biton, M.Kafri,Y. Segev,
T. Hendler, A. D. Korczyn, M. Graif, and Y. Cohen. High b-value q-space ana-
lyzed diffusion-weighted mri: application to multiple sclerosis. Magn Reson Med,
47(1):115–126, Jan 2002.
[ACH
+
02] Konstantinos Arfanakis, Dietmar Cordes, Victor M Haughton, John D Carew,
and M. Elizabeth Meyerand. Independent component analysis applied to diffusion
tensor mri. Magn Reson Med, 47(2):354–363, Feb 2002.
[AFRB04] Yaniv Assaf, Raisa Z Freidlin, Gustavo K Rohde, and Peter J Basser. New
modeling and experimental framework to characterize hindered and restricted
water diffusion in brain white matter. Magn Reson Med, 52(5):965–978, Nov
2004.
[AH01] Erkki Oja. Aapo Hyvarinen, Juha Karhunen. Independent Component Analysis.
Wiley Interscience. John Wiley and Sons Inc., 2001.
[AHA
+
01] J. L. Andersson, C. Hutton, J. Ashburner, R. Turner, and K. Friston. Modeling
geometric deformations in epi time series. Neuroimage, 13(5):903–919, May 2001.
[AHK
+
00] A. L. Alexander, K. Hasan, G. Kindlmann, D. L. Parker, and J. S. Tsuruda. A
geometric analysis of diffusion tensor measurements of the human brain. Magn
Reson Med, 44(2):283–291, Aug 2000.
[AHL
+
01] A. L. Alexander, K. M. Hasan, M. Lazar, J. S. Tsuruda, and D. L. Parker.
Analysis of partial volume effects in diffusion-tensor mri. Magn Reson Med,
45(5):770–780, May 2001.
[ALA97] Dennis L. Parker Andrew L. Alexander, Jay S. Tsuruda. Elimination of eddy
current artifacts in diffusion-weighted echo-planar images: The use of bipolar
gradients. Magnetic Resonance in Medicine, 38(6):1016–1021, 1997.
[Ale05a] D. C. Alexander. An introduction to computational diffusion MRI: the diffusion
tensor and beyond. Springers, 2005.
109
[Ale05b] Daniel C Alexander. Maximum entropy spherical deconvolution for diffusion mri.
Inf Process Med Imaging, 19:76–87, 2005.
[Ale05c] Daniel C Alexander. Multiple-fiber reconstruction algorithms for diffusion mri.
Ann N Y Acad Sci, 1064:113–133, Dec 2005.
[And01] A. W. Anderson. Theoretical analysis of the effects of noise on diffusion tensor
imaging. Magn Reson Med, 46(6):1174–1188, Dec 2001.
[And05] Adam W Anderson. Measurement of fiber orientation distributions using high
angular resolution diffusion imaging. Magn Reson Med, 54(5):1194–1206, Nov
2005.
[APBG01] D. C. Alexander, C. Pierpaoli, P. J. Basser, and J. C. Gee. Spatial transforma-
tions of diffusion tensor magnetic resonance images. IEEE Trans Med Imaging,
20(11):1131–1139, Nov 2001.
[AS02] Jesper L R Andersson and Stefan Skare. A model-based method for retrospec-
tive correction of geometric distortions in diffusion-weighted epi. Neuroimage,
16(1):177–199, May 2002.
[AS05] Siamak Ardekani and Usha Sinha. Geometric distortion correction of high-
resolution 3 t diffusion tensor brain images. Magn Reson Med, 54(5):1163–1171,
Nov 2005.
[BA00] M. E. Bastin and P. A. Armitage. On the use of water phantom images to cali-
brate and correct eddy current induced artefacts in mr diffusion tensor imaging.
Magn Reson Imaging, 18(6):681–687, Jul 2000.
[Bas99] M. E. Bastin. Correction of eddy current-induced artefacts in diffusion tensor
imaging using iterative cross-correlation. Magn Reson Imaging, 17(7):1011–1024,
Sep 1999.
[Bas01] M. E. Bastin. On the use of the flair technique to improve the correction of eddy
current induced artefacts in mr diffusion tensor imaging. Magn Reson Imaging,
19(7):937–950, Sep 2001.
[Bas02] Peter J Basser. Relationships between diffusion tensor and q-space mri. Magn
Reson Med, 47(2):392–397, Feb 2002.
[BBJ
+
07] T. E J Behrens, H. Johansen Berg, S. Jbabdi, M. F S Rushworth, and M. W.
Woolrich. Probabilistic diffusion tractography with multiple fibre orientations:
What can we gain? Neuroimage, 34(1):144–155, Jan 2007.
[BCM
+
08] Jeffrey I Berman, SungWon Chung, Pratik Mukherjee, Christopher P Hess, Eric T
Han, and Roland G Henry. Probabilistic streamline q-ball tractography using the
residual bootstrap. Neuroimage, 39(1):215–222, Jan 2008.
[Ber09] Jeffrey Berman. Diffusion mr tractography as a tool for surgical planning. Magn
Reson Imaging Clin N Am, 17(2):205–214, May 2009.
[BFW06] Saurav Basu, Thomas Fletcher, and Ross Whitaker. Rician noise removal in
diffusion tensor mri. Med Image Comput Comput Assist Interv Int Conf Med
Image Comput Comput Assist Interv, 9(Pt 1):117–125, 2006.
110
[BGKM
+
04] Naama Barnea-Goraly, Hower Kwon, Vinod Menon, Stephan Eliez, Linda Lotspe-
ich, and Allan L Reiss. White matter structure in autism: preliminary evidence
from diffusion tensor imaging. Biol Psychiatry, 55(3):323–326, Feb 2004.
[BHH
+
09] Angelos Barmpoutis, Min Sig Hwang, Dena Howland, John R Forder, and Baba C
Vemuri. Regularized positive-definite fourth order tensor field estimation from
dw-mri. Neuroimage, 45(1 Suppl):S153–S162, Mar 2009.
[BJ02] Peter J Basser and Derek K Jones. Diffusion-tensor mri: theory, experimental
design and data analysis - a technical review. NMR Biomed, 15(7-8):456–467,
2002.
[BJBW
+
03] T. E J Behrens, H. Johansen-Berg, M. W. Woolrich, S. M. Smith, C. A M
Wheeler-Kingshott, P. A. Boulby, G. J. Barker, E. L. Sillery, K. Sheehan, O. Cic-
carelli, A. J. Thompson, J. M. Brady, and P. M. Matthews. Non-invasive mapping
of connections between human thalamus and cortex using diffusion imaging. Nat
Neurosci, 6(7):750–757, Jul 2003.
[BKKT04] Nils Bodammer, Jrn Kaufmann, Martin Kanowski, and Claus Tempelmann.
Eddy current correction in diffusion-weighted imaging using pairs of images ac-
quired with opposite diffusion gradient polarity. Magn Reson Med, 51(1):188–193,
Jan 2004.
[BKLW06] Orjan Bergmann, Gordon Kindlmann, Arvid Lundervold, and Carl-Fredrik
Westin. Diffusion k-tensor estimation from q-ball imaging using discretized prin-
cipal axes. Med Image Comput Comput Assist Interv Int Conf Med Image Com-
put Comput Assist Interv, 9(Pt 2):268–275, 2006.
[Blo53] Felix Bloch. The principle of nuclear induction. Science, 118(3068):425–430, Oct
1953.
[BML94] P J Basser, J Mattiello, and D LeBihan. Mr diffusion tensor spectroscopy and
imaging. Biophys. J., 66(1):259–267, 1994.
[BP96] Peter J. Basser and Carlo Pierpaolib. Microstructural and physiological features
of tissues elucidated by quantitative-diffusion-tensor mri. Journal of Magnetic
Resonance, 111(3):209–219, 1996.
[BP00] P. J. Basser and S. Pajevic. Statistical artifacts in diffusion tensor mri (dt-mri)
caused by background noise. Magn Reson Med, 44(1):41–50, Jul 2000.
[BPP
+
00] P. J. Basser, S. Pajevic, C. Pierpaoli, J. Duda, and A. Aldroubi. In vivo fiber
tractography using dt-mri data. Magn Reson Med, 44(4):625–632, Oct 2000.
[BSAO
+
08] Amnon Bar-Shir, Liat Avram, Evren Ozarslan, Peter J Basser, and Yoram Cohen.
The effect of the diffusion time and pulse gradient duration ratio on the diffrac-
tion pattern and the structural information estimated from q-space diffusion mr:
experiments and simulations. J Magn Reson, 194(2):230–236, Oct 2008.
[BSC08a] Amnon Bar-Shir and Yoram Cohen. Crossing fibers, diffractions and nonhomo-
geneous magnetic field: correction of artifacts by bipolar gradient pulses. Magn
Reson Imaging, 26(6):801–808, Jul 2008.
[BSC08b] Amnon Bar-Shir and Yoram Cohen. The effect of the rotational angle on mr dif-
fusion indices in nerves: is the rms displacement of the slow-diffusing component
a good measure of fiber orientation? J Magn Reson, 190(1):33–42, Jan 2008.
111
[BWJ
+
03] T.EJBehrens,M.W.Woolrich,M.Jenkinson,H.Johansen-Berg,R.G.Nunes,
S. Clare, P. M. Matthews, J. M. Brady, andS.M.Smith. Characterizationand
propagation of uncertainty in diffusion-weighted mr imaging. Magn Reson Med,
50(5):1077–1088, Nov 2003.
[CGS06] Bin Chen, Hua Guo, and Allen W Song. Correction for direction-dependent
distortions in diffusion tensor imaging using matched magnetic field maps. Neu-
roimage, 30(1):121–129, Mar 2006.
[CH05] Bin Chen and Edward W Hsu. Noise removal in magnetic resonance diffusion
tensor imaging. Magn Reson Med, 54(2):393–401, Aug 2005.
[CKPB07] Lin-Ching Chang, Cheng Guan Koay, Carlo Pierpaoli, and Peter J Basser. Vari-
ance of estimated dti-derived parameters via first-order perturbation methods.
Magn Reson Med, 57(1):141–149, Jan 2007.
[CMW
+
06] Peng Cheng, Vincent A Magnotta, Dee Wu, Peg Nopoulos, David J Moser, Jane
Paulsen, Ricardo Jorge, and Nancy C Andreasen. Evaluation of the gtract diffu-
sion tensor tractography algorithm: a validation and reliability study. Neuroim-
age, 31(3):1075–1085, Jul 2006.
[CPGC99] F. Calamante, D. A. Porter, D. G. Gadian, and A. Connelly. Correction for eddy
current induced bo shifts in diffusion-weighted echo-planar imaging. Magn Reson
Med, 41(1):95–102, Jan 1999.
[Cra75] J. Crank. The mathematics of diffusion. Oxford, UK: Oxford University Press,
1975.
[CRMGIM
+
08] Erick Jorge Canales-Rodrguez, Lester Melie-Garca, Yasser Iturria-Medina, Ed-
uardo Martnez-Montes, Yasser Alemn-Gmez, and Ching-Po Lin. Inferring mul-
tiple maxima in intravoxel white matter fiber distribution. Magn Reson Med,
60(3):616–630, Sep 2008.
[CRMGIM09] Erick Jorge Canales-Rodrguez, Lester Melie-Garca, and Yasser Iturria-Medina.
Mathematical description of q-space in spherical coordinates: Exact q-ball imag-
ing. Magn Reson Med, Mar 2009.
[CSR
+
05] Jennifer S W Campbell, Kaleem Siddiqi, Vladimir V Rymar, Abbas F Sadikot,
and G. Bruce Pike. Flow-based fiber tracking with diffusion tensor and q-ball
data: validation and comparison to principal diffusion direction techniques. Neu-
roimage, 27(4):725–736, Oct 2005.
[CTP
+
03] O. Ciccarelli, A. T. Toosy, G. J M Parker, C. A M Wheeler-Kingshott, G. J.
Barker, D. H. Miller, and A. J. Thompson. Diffusion tractography based group
mapping of major white-matter pathways in the human brain. Neuroimage,
19(4):1545–1555, Aug 2003.
[CWC
+
06] Ming-Chung Chou, Ming-Long Wu, Cheng-Yu Chen, Chao-Ying Wang, Teng-Yi
Huang, Yi-Jui Liu, Chun-Jung Juan, and Hsiao-Wen Chung. Tensor deflection
(tend) tractography with adaptive subvoxel stepping. J Magn Reson Imaging,
24(2):451–458, Aug 2006.
[CWW07] M. Singh C. W. Wong. Estimating number of fiber directions per voxel for
multiple fiber dti tractography. Proc. Joint Annual Scientific Meeting, Int. Soc.
Magnetic Resonance In Medicine and European Society of Magnetic Resonance
in Medicine and Biology, 2007.
112
[CYC
+
09] Kuan-Hung Cho, Chun-Hung Yeh, Yi-Ping Chao, Jiun-Jie Wang, Jyh-Horng
Chen, and Ching-Po Lin. Potential in reducing scan times of hardi by accu-
rate correction of the cross-term in a hemispherical encoding scheme. JMagn
Reson Imaging, 29(6):1386–1394, Jun 2009.
[CYT
+
08] Kuan-Hung Cho, Chun-Hung Yeh, Jacques-Donald Tournier, Yi-Ping Chao, Jyh-
Horng Chen, and Ching-Po Lin. Evaluation of the accuracy and angular resolu-
tion of q-ball imaging. Neuroimage, 42(1):262–271, Aug 2008.
[DAFD07] Maxime Descoteaux, Elaine Angelino, Shaun Fitzgibbons, and Rachid Deriche.
Regularized, fast, and robust analytical q-ball imaging. Magn Reson Med,
58(3):497–510, Sep 2007.
[DDKA09] Maxime Descoteaux, Rachid Deriche, Thomas R Knsche, and Alfred Anwander.
Deterministic and probabilistic tractography based on complex fibre orientation
distributions. IEEE Trans Med Imaging, 28(2):269–286, Feb 2009.
[DG] ISMRM Diffusion and Perfusion Study Group. Common dti dataset.
[DSH
+
06] J. T. Devlin, E. L. Sillery, D. A. Hall, P. Hobden, T. E J Behrens, R. G. Nunes,
S. Clare, P. M. Matthews, D. R. Moore, and H. Johansen-Berg. Reliable iden-
tification of the auditory thalamus using multi-modal structural analyses. Neu-
roimage, 30(4):1112–1120, May 2006.
[FLT
+
03] Christopher G Filippi, Doris D M Lin, Apostolos J Tsiouris, Richard Watts,
A. Maurine Packard, Linda A Heier, and Aziz M Ulu? Diffusion-tensor mr
imaging in children with developmental delay: preliminary findings. Radiology,
229(1):44–50, Oct 2003.
[Fra01] L. R. Frank. Anisotropy in high angular resolution diffusion-weighted mri. Magn
Reson Med, 45(6):935–939, Jun 2001.
[Fra02] Lawrence R Frank. Characterization of anisotropy in high angular resolution
diffusion-weighted mri. Magn Reson Med, 47(6):1083–1099, Jun 2002.
[FVN07] Hubert M J Fonteijn, Frans A J Verstraten, and David G Norris. Probabilistic
inference on q-ball imaging data. IEEE Trans Med Imaging, 26(11):1515–1524,
Nov 2007.
[HM96] J. C. Haselgrove and J. R. Moore. Correction for distortion of echo-planar images
used to calculate the apparent diffusion coefficient. Magn Reson Med, 36(6):960–
964, Dec 1996.
[HME
+
09] Hamied A Haroon, David M Morris, Karl V Embleton, Daniel C Alexander,
and Geoffrey J M Parker. Using the model-based residual bootstrap to quantify
uncertainty in fiber orientations from q-ball analysis. IEEE Trans Med Imaging,
28(4):535–550, Apr 2009.
[HMH
+
06] Christopher P Hess, Pratik Mukherjee, Eric T Han, Duan Xu, and Daniel B Vi-
gneron. Q-ball reconstruction of multimodal fiber orientations using the spherical
harmonic basis. Magn Reson Med, 56(1):104–117, Jul 2006.
[Hor99] M. A. Horsfield. Mapping eddy current induced fields for the correction of
diffusion-weighted echo planar images. Magn Reson Imaging, 17(9):1335–1345,
Nov 1999.
113
[HTNS09] Nathan S Hageman, Arthur W Toga, Katherine L Narr, and David W Shattuck.
A diffusion tensor imaging tractography algorithm based on navier-stokes fluid
mechanics. IEEE Trans Med Imaging, 28(3):348–360, Mar 2009.
[HWA05] Tim Hosey, Guy Williams, and Richard Ansorge. Inference of multiple fiber
orientations in high angular resolution diffusion imaging. Magn Reson Med,
54(6):1480–1489, Dec 2005.
[JA03] Kalvis M Jansons and Daniel C Alexander. Persistent angular structure: new in-
sights from diffusion mri data. dummy version. Inf Process Med Imaging, 18:672–
683, Jul 2003.
[JA08] Ha-Kyu Jeong and Adam W Anderson. Characterizing fiber directional uncer-
tainty in diffusion tensor mri. Magn Reson Med, 60(6):1408–1421, Dec 2008.
[JB04] Derek K Jones and Peter J Basser. ”squashing peanuts and smashing pumpkins”:
how noise distorts diffusion-weighted mr data. Magn Reson Med, 52(5):979–993,
Nov 2004.
[JBP98] P. Jezzard, A. S. Barnett, and C. Pierpaoli. Characterization of and correction
for eddy current artifacts in echo planar diffusion imaging. Magn Reson Med,
39(5):801–812, May 1998.
[JHR
+
05] Jens H Jensen, Joseph A Helpern, Anita Ramani, Hanzhang Lu, and Kyle
Kaczynski. Diffusional kurtosis imaging: the quantification of non-gaussian water
diffusion by means of magnetic resonance imaging. Magn Reson Med, 53(6):1432–
1440, Jun 2005.
[JHS99] D. K. Jones, M. A. Horsfield, and A. Simmons. Optimal strategies for measuring
diffusion in anisotropic systems by magnetic resonance imaging. Magn Reson
Med, 42(3):515–525, Sep 1999.
[Jon03] Derek K Jones. Determining and visualizing uncertainty in estimates of fiber
orientation from diffusion tensor mri. Magn Reson Med, 49(1):7–12, Jan 2003.
[JP05] Derek K Jones and Carlo Pierpaoli. Confidence mapping in diffusion tensor
magnetic resonance imaging tractography using a bootstrap approach. Magn
Reson Med, 53(5):1143–1149, May 2005.
[JRHH08] M. R. Jayachandra, N. Rehbein, C. Herweh, and S. Heiland. Fiber tracking
of human brain using fourth-order tensor and high angular resolution diffusion
imaging. Magn Reson Med, 60(5):1207–1217, Nov 2008.
[JVO
+
07] Bing Jian, Baba C Vemuri, Evren Ozarslan, Paul R Carney, and Thomas H
Mareci. A novel tensor distribution model for the diffusion-weighted mr signal.
Neuroimage, 37(1):164–176, Aug 2007.
[KCWT08] Li-Wei Kuo, Jyh-Horng Chen, Van Jay Wedeen, and Wen-Yih Isaac Tseng. Opti-
mization of diffusion spectrum imaging and q-ball imaging on clinical mri system.
Neuroimage, 41(1):7–18, May 2008.
[KHR
+
94] M. D. King, J. Houseman, S. A. Roussel, N. van Bruggen, S. R. Williams, and
D. G. Gadian. q-space imaging of the brain. Magn Reson Med, 32(6):707–713,
Dec 1994.
114
[KKJ07] Timothy A Keller, Rajesh K Kana, and Marcel Adam Just. A developmental
study of the structural integrity of white matter in autism. Neuroreport, 18(1):23–
27, Jan 2007.
[KKJS05] Sungheon Kim, Sungheon Kim, Jeong-Won Jeong, and M. Singh. Estimation of
multiple fiber orientations from diffusion tensor mri using independent component
analysis. Nuclear Science, IEEE Transactions on, 52(1):266–273, 2005.
[KMK08] B. W. Kreher, I. Mader, and V. G. Kiselev. Gibbs tracking: a novel approach for
the reconstruction of neuronal pathways. Magn Reson Med, 60(4):953–963, Oct
2008.
[KSM
+
05] B. W. Kreher, J. F. Schneider, I. Mader, E. Martin, J. Hennig, and K. A. Il’yasov.
Multitensor approach for analysis and tracking of complex fiber configurations.
Magn Reson Med, 54(5):1216–1225, Nov 2005.
[KTNU00] J. Kybic, P. Thvenaz, A. Nirkko, and M. Unser. Unwarping of unidirectionally
distorted epi images. IEEE Trans Med Imaging, 19(2):80–93, Feb 2000.
[KWT07] Mark H Khachaturian, Jonathan J Wisco, and David S Tuch. Boosting the
sampling efficiency of q-ball imaging using multiple wavevector fusion. Magn
Reson Med, 57(2):289–296, Feb 2007.
[LA03] Mariana Lazar and Andrew L Alexander. An error analysis of white matter
tractography methods: synthetic diffusion tensor field simulations. Neuroimage,
20(2):1140–1153, Oct 2003.
[LA05] Mariana Lazar and Andrew L Alexander. Bootstrap white matter tractography
(boot-trac). Neuroimage, 24(2):524–532, Jan 2005.
[LBA
+
07] Jee Eun Lee, Erin D Bigler, Andrew L Alexander, Mariana Lazar, Molly B
DuBray, Moo K Chung, Michael Johnson, Jubel Morgan, Judith N Miller,
William M McMahon, Jeffrey Lu, Eun-Kee Jeong, and Janet E Lainhart. Diffu-
sion tensor imaging of white matter in the superior temporal gyrus and temporal
stem in autism. Neurosci Lett, 424(2):127–132, Sep 2007.
[LBAM04] Chunlei Liu, Roland Bammer, Burak Acar, and Michael E Moseley. Character-
izing non-gaussian diffusion by using generalized diffusion tensors. Magn Reson
Med, 51(5):924–937, May 2004.
[LBM05] Chunlei Liu, Roland Bammer, and Michael E Moseley. Limitations of apparent
diffusion coefficient-based models in characterizing non-gaussian diffusion. Magn
Reson Med, 54(2):419–428, Aug 2005.
[LCD
+
09] C. Lenglet, J. S W Campbell, M. Descoteaux, G. Haro, P. Savadjiev, D. Wasser-
mann, A. Anwander, R. Deriche, G. B. Pike, G. Sapiro, K. Siddiqi, and P. M.
Thompson. Mathematical methods for diffusion mri processing. Neuroimage,
45(1 Suppl):S111–S122, Mar 2009.
[LJRH06] Hanzhang Lu, Jens H Jensen, Anita Ramani, and Joseph A Helpern. Three-
dimensional characterization of non-gaussian water diffusion in humans using
diffusion kurtosis imaging. NMR Biomed, 19(2):236–247, Apr 2006.
[LJXH08] Mariana Lazar, Jens H Jensen, Liang Xuan, and Joseph A Helpern. Estimation
of the orientation distribution function from diffusional kurtosis imaging. Magn
Reson Med, 60(4):774–781, Oct 2008.
115
[LLA05] Mariana Lazar, Jong Hoon Lee, and Andrew L Alexander. Axial asymmetry
of water diffusion in brain white matter. Magn Reson Med, 54(4):860–867, Oct
2005.
[LWT
+
03] Mariana Lazar, David M Weinstein, Jay S Tsuruda, Khader M Hasan, Konstanti-
nos Arfanakis, M. Elizabeth Meyerand, Benham Badie, Howard A Rowley, Victor
Haughton, Aaron Field, and Andrew L Alexander. White matter tractography
using diffusion tensor deflection. Hum Brain Mapp, 18(4):306–321, Apr 2003.
[LZZ
+
09] A. D. Leow, S. Zhu, L. Zhan, K. McMahon, G. I. de Zubicaray, M. Meredith,
M. J. Wright, A. W. Toga, and P. M. Thompson. The tensor distribution function.
Magn Reson Med, 61(1):205–214, Jan 2009.
[MHX
+
08] Pratik Mukherjee, Christopher P Hess, Duan Xu, Eric T Han, Douglas A Kelley,
and Daniel B Vigneron. Development and initial evaluation of 7-t q-ball imaging
of the human brain. Magn Reson Imaging, 26(2):171–180, Feb 2008.
[MS07a] A. Shetty A. Rajagopalan D. Hwang C. W. Wong H. Chui N. Schuff M. Weiner
M. Singh, W. Sungkarat. Toward quantitation of whole-brain tractography group
comparisons with application to alzheimer disease. Proc. Joint Annual Scientific
Meeting, Int. Soc. Magnetic Resonance In Medicine and European Society of
Magnetic Resonance in Medicine and Biology, 2007.
[MS07b] C. W. Wong M. Singh. Whole-brain tractography incorporating ica based
crossing-fiber orientations. Proc. Joint Annual Scientific Meeting, Int. Soc. Mag-
netic Resonance In Medicine and European Society of Magnetic Resonance in
Medicine and Biology, Berlin, 2007.
[MS09] C. W. Wong M. Singh. Ica based multi-fiber tractography. Proc. 17-th Scientific
Meeting, Int. Soc. Magnetic Resonance In Medicine, 2009.
[Nis] Dwight G. Nishimura. Principles of magnetic resonance imaging.
[oA08] National Institute on Aging. Alzheimer’s disease: Unraveling the mystery,
September 2008.
[OM03] Evren Ozarslan and Thomas H Mareci. Generalized diffusion tensor imaging
and analytical relationships between diffusion tensor imaging and high angular
resolution diffusion imaging. Magn Reson Med, 50(5):955–965, Nov 2003.
[PA05] Geoffrey J M Parker and Daniel C Alexander. Probabilistic anatomical connec-
tivity derived from the microscopic persistent angular structure of cerebral tissue.
Philos Trans R Soc Lond B Biol Sci, 360(1457):893–902, May 2005.
[PB96] C. Pierpaoli and P. J. Basser. Toward a quantitative assessment of diffusion
anisotropy. Magn Reson Med, 36(6):893–906, Dec 1996.
[PB03] Sinisa Pajevic and Peter J Basser. Parametric and non-parametric statistical
analysis of dt-mri data. J Magn Reson, 161(1):1–14, Mar 2003.
[PCF
+
00] C. Poupon, C. A. Clark, V. Frouin, J. Rgis, I. Bloch, D. Le Bihan, and J. Mangin.
Regularization of diffusion-based direction maps for the tracking of brain white
matter fascicles. Neuroimage, 12(2):184–195, Aug 2000.
116
[PHWK03] Geoffrey J M Parker, Hamied A Haroon, and Claudia A M Wheeler-Kingshott.
A framework for a streamline-based probabilistic index of connectivity (pico)
using a structural interpretation of mri diffusion measurements. J Magn Reson
Imaging, 18(2):242–254, Aug 2003.
[PJB
+
96] C Pierpaoli, P Jezzard, PJ Basser, A Barnett, and G Di Chiro. Diffusion tensor
mr imaging of the human brain. Radiology, 201(3):637–648, 1996.
[PMP
+
00] N. G. Papadakis,K.M. Martin,J.D.Pickard,L.D.Hall,T. A. Carpenter,and
C. L. Huang. Gradient preemphasis calibration in diffusion-weighted echo-planar
imaging. Magn Reson Med, 44(4):616–624, Oct 2000.
[PPC
+
05] M. Perrin, C. Poupon, Y. Cointepas, B. Rieul, N. Golestani, C. Pallier, D. Rivire,
A. Constantinesco, D. Le Bihan, and J. F. Mangin. Fiber tracking in q-ball fields
using regularized particle trajectories. Inf Process Med Imaging, 19:52–63, 2005.
[PPR
+
05] Muriel Perrin, Cyril Poupon, Bernard Rieul, Patrick Leroux, Andr Constanti-
nesco, Jean-Franois Mangin, and Denis Lebihan. Validation of q-ball imaging
with a diffusion fibre-crossing phantom on a clinical scanner. Philos Trans R Soc
LondBBiol Sci, 360(1457):881–891, May 2005.
[PPR
+
07] C. Poupon, F. Poupon, A. Roche, Y. Cointepas, J. Dubois, and J. F. Mangin.
Real-time mr diffusion tensor and q-ball imaging using kalman filtering. Med
Image Comput Comput Assist Interv Int Conf Med Image Comput Comput Assist
Interv, 10(Pt 1):27–35, 2007.
[PRD
+
08] Cyril Poupon, Alexis Roche, Jessica Dubois, Jean-Franois Mangin, and Fabrice
Poupon. Real-time mr diffusion tensor and q-ball imaging using kalman filtering.
Med Image Anal, 12(5):527–534, Oct 2008.
[PRK
+
08] Cyril Poupon, Bernard Rieul, Irina Kezele, Muriel Perrin, Fabrice Poupon, and
Jean-Franois Mangin. New diffusion phantoms dedicated to the study and vali-
dation of high-angular-resolution diffusion imaging (hardi) models. Magn Reson
Med, 60(6):1276–1283, Dec 2008.
[PRP
+
08] V. Prckovska, A. F. Roebroeck, W. L P M Pullens, A. Vilanova, and B. M. ter
Haar Romeny. Optimal acquisition schemes in high angular resolution diffusion
weighted imaging. Med Image Comput Comput Assist Interv Int Conf Med Image
Comput Comput Assist Interv, 11(Pt 2):9–17, 2008.
[PSB
+
02] Geoffrey J M Parker, Klaas E Stephan, Gareth J Barker, James B Rowe, David G
MacManus, Claudia A M Wheeler-Kingshott, Olga Ciccarelli, Richard E Passing-
ham, Rachel L Spinks, Roger N Lemon, and Robert Turner. Initial demonstration
of in vivo tracing of axonal projections in the macaque brain and comparison with
the human brain using diffusion tensor imaging and fast marching tractography.
Neuroimage, 15(4):797–809, Apr 2002.
[PSBM05] Nikos G Papadakis, Theo Smponias, Jason Berwick, and John E W Mayhew.
k-space correction of eddy-current-induced distortions in diffusion-weighted echo-
planar imaging. Magn Reson Med, 53(5):1103–1111, May 2005.
[PWKB02] Geoffrey J M Parker, Claudia A M Wheeler-Kingshott, and Gareth J Barker.
Estimating distributed anatomical connectivity using fast marching methods and
diffusion tensor imaging. IEEE Trans Med Imaging, 21(5):505–512, May 2002.
117
[QRO
+
08] Arish A Qazi, Alireza Radmanesh, Lauren O’Donnell, Gordon Kindlmann,
Sharon Peled, Stephen Whalen, Carl-Fredrik Westin, and Alexandra J Golby.
Resolving crossings in the corticospinal tract by two-tensor streamline tractogra-
phy: Method and clinical assessment using fmri. Neuroimage, Jul 2008.
[RBB
+
04] G. K. Rohde, A. S. Barnett, P. J. Basser, S. Marenco, and C. Pierpaoli. Com-
prehensive approach for correction of motion and distortion in diffusion-weighted
mri. Magn Reson Med, 51(1):103–114, Jan 2004.
[RHWW03] T. G. Reese, O. Heid, R. M. Weisskoff, and V. J. Wedeen. Reduction of eddy-
current-induced distortion in diffusion mri using a twice-refocused spin echo.
Magn Reson Med, 49(1):177–182, Jan 2003.
[RVH05] Elliot A. Tanis Robert V Hogg. Probability and Statistical Inference.Prentice
Hall, 2005.
[Sin06] C.-W. Jeong J.-W. Singh, M. Wong. Simulation and experimental study of
multiple-fibers per voxel detection by ica in dti tractography. Nuclear Science
Symposium Conference Record, IEEE, 5:3176–3179, 2006.
[SJC
+
06] P. Staempfli, T. Jaermann, G. R. Crelier, S. Kollias, A. Valavanis, and P. Boe-
siger. Resolving fiber crossing using advanced fast marching tractography based
on diffusion tensor imaging. Neuroimage, 30(1):110–120, Mar 2006.
[SKWP06] Manbir Singh, Aditya Kwatra, Chi-Wah Wong, and Viktor Prasanna. Accelera-
tion of fiber tracking in dti tractography by reconfigurable computer hardware.
Conf Proc IEEE Eng Med Biol Soc, 1:4819–4822, 2006.
[SLC
+
04] Yuji Shen, David J Larkman, Serena Counsell, Ida M Pu, David Edwards, and
Joseph V Hajnal. Correction of high-order eddy current induced geometric dis-
tortion in diffusion-weighted echo-planar images. Magn Reson Med, 52(5):1184–
1189, Nov 2004.
[Ste65] E. O. Stejskal. Use of spin echoes in a pulsed magnetic-field gradient to study
anisotropic, restricted diffusion and flow. The Journal of Chemical Physics,
43(10):3597–3603, 1965.
[SW06] Manbir Singh and Chi-Wah Wong. Recovery of multiple fibers per voxel by ica
in dti tractography. Conf Proc IEEE Eng Med Biol Soc, 1:735–738, 2006.
[Tan78] J. E. Tanner. Transient diffusion in a system partitioned by permeable barriers.
application to nmr measurements with a pulsed field gradient. J. Chem. Phys.,
69(4):1748–1754, 1978.
[TCC07] J-Donald Tournier, Fernando Calamante, and Alan Connelly. Robust determi-
nation of the fibre orientation distribution in diffusion mri: non-negativity con-
strained super-resolved spherical deconvolution. Neuroimage, 35(4):1459–1472,
May 2007.
[TCGC04] J-Donald Tournier, Fernando Calamante, David G Gadian, and Alan Connelly.
Direct estimation of the fiber orientation density function from diffusion-weighted
mri data using spherical deconvolution. Neuroimage, 23(3):1176–1185, Nov 2004.
[TCP
+
04] Ahmed T Toosy, Olga Ciccarelli, Geoff J M Parker, Claudia A M Wheeler-
Kingshott, David H Miller, and Alan J Thompson. Characterizing function-
structure relationships in the human visual system with functional mri and dif-
fusion tensor imaging. Neuroimage, 21(4):1452–1463, Apr 2004.
118
[TCS08] Trong-Kha Truong, Bin Chen, and Allen W Song. Integrated sense dti with cor-
rection of susceptibility- and eddy current-induced geometric distortions. Neu-
roimage, 40(1):53–58, Mar 2008.
[TRW
+
02] David S Tuch, Timothy G Reese, Mette R Wiegell, Nikos Makris, John W Bel-
liveau, and Van J Wedeen. High angular resolution diffusion imaging reveals
intravoxel white matter fiber heterogeneity. Magn Reson Med, 48(4):577–582,
Oct 2002.
[TRWW03] David S Tuch, Timothy G Reese, Mette R Wiegell, and Van J Wedeen. Diffusion
mri of complex neural architecture. Neuron, 40(5):885–895, Dec 2003.
[Tuc04] David S Tuch. Q-ball imaging. Magn Reson Med, 52(6):1358–1372, Dec 2004.
[TWK
+
05] David S Tuch, Jonathan J Wisco, Mark H Khachaturian, Leeland B Ekstrom,
Rolf Ktter, and Wim Vanduffel. Q-ball imaging of macaque white matter archi-
tecture. Philos Trans R Soc Lond B Biol Sci, 360(1457):869–879, May 2005.
[TYC
+
08] J-Donald Tournier, Chun-Hung Yeh, Fernando Calamante, Kuan-Hung Cho,
Alan Connelly, and Ching-Po Lin. Resolving crossing fibres using constrained
spherical deconvolution: validation using diffusion-weighted imaging phantom
data. Neuroimage, 42(2):617–625, Aug 2008.
[WFA08] Yu-Chien Wu, Aaron S Field, and Andrew L Alexander. Computation of diffusion
function measures in q-space using magnetic resonance hybrid diffusion imaging.
IEEE Trans Med Imaging, 27(6):858–865, Jun 2008.
[WHT
+
05] Van J Wedeen, Patric Hagmann, Wen-Yih Isaac Tseng, Timothy G Reese, and
Robert M Weisskoff. Mapping complex tissue architecture with diffusion spec-
trum magnetic resonance imaging. Magn Reson Med, 54(6):1377–1386, Dec 2005.
[WS07] C. W. Wong and M. Singh. Estimating number of fiber directions per voxel for
ica dti tractography. Progress in Biomedical Optics and Imaging: Physiology,
Function and Structure from Medical Images, 2007.
[WTLW03] Mette R Wiegell, David S Tuch, Henrik B W Larsson, and Van J Wedeen. Auto-
matic segmentation of thalamic nuclei from diffusion tensor magnetic resonance
imaging. Neuroimage, 19(2 Pt 1):391–401, Jun 2003.
[WWS
+
08] V. J. Wedeen, R. P. Wang, J. D. Schmahmann, T. Benner, W. Y I Tseng, G. Dai,
D. N. Pandya, P. Hagmann, H. D’Arceuil, and A. J. de Crespigny. Diffusion
spectrum magnetic resonance imaging (dsi) tractography of crossing fibers. Neu-
roimage, 41(4):1267–1277, Jul 2008.
[XHB
+
08] Dongrong Xu, Xuejun Hao, Ravi Bansal, Kerstin J Plessen, and Bradley S Pe-
terson. Seamless warping of diffusion tensor fields. IEEE Trans Med Imaging,
27(3):285–299, Mar 2008.
[XMS
+
02] Dongrong Xu, Susumu Mori, Meiyappan Solaiyappan, Peter C M van Zijl, and
Christos Davatzikos. A framework for callosal fiber distribution analysis. Neu-
roimage, 17(3):1131–1143, Nov 2002.
[XMS
+
03] Dongrong Xu, Susumu Mori, Dinggang Shen, Peter C M van Zijl, and Christos
Davatzikos. Spatial normalization of diffusion tensor fields. Magn Reson Med,
50(1):175–182, Jul 2003.
119
[YCL
+
08] Chun-Hung Yeh, Kuan-Hung Cho, Hsuan-Cheng Lin, Jiun-Jie Wang, and Ching-
Po Lin. Reduced encoding diffusion spectrum imaging implemented with a bi-
gaussian model. IEEE Trans Med Imaging, 27(10):1415–1424, Oct 2008.
[ZHK
+
06] Jiancheng Zhuang, Jan Hrabe, Alayar Kangarlu, Dongrong Xu, Ravi Bansal,
Craig A Branch, and Bradley S Peterson. Correction of eddy-current distortions
in diffusion tensor images using the known directions and strengths of diffusion
gradients. J Magn Reson Imaging, 24(5):1188–1193, Nov 2006.
120

Abstract (if available)

Abstract Brownian motion is a phenomenon describing the random movement of particles. Water molecules, the most abundant particle inside the body, diffuse in all directions with equal probability in free space. In the presence of organized structure, diffusion is restricted. Diffusion Tensor Imaging (DTI) was developed to model water diffusion in tissues. In the human brain, myelinated neuronal fiber bundles hinder water diffusion. By knowing the amount of diffusion along different directions, the underlying fiber bundle directions can be inferred.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Diffusion tensor tractography: visualization and quantitation with applications to Alzheimer disease and traumatic brain injury

PDF

Theory and simulation of diffusion magnetic resonance imaging on brain's white matter

PDF

Unveiling the white matter microstructure in 22q11.2 deletion syndrome with diffusion magnetic resonance imaging

PDF

Diffusion MRI of the human brain: signal modeling and quantitative analysis

Asset Metadata

Creator Wong, Chi Wah (Alec) (author)

Core Title Diffusion MRI white matter tractography: estimation of multiple fibers per voxel using independent component analysis

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Biomedical Engineering

Publication Date 07/21/2009

Defense Date 06/10/2009

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag diffusion MRI,diffusion tensor imaging,fiber crossing,independent componenet analysis,magnetic resonance imaging,multiple fibers,OAI-PMH Harvest,tractography

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Singh, Manbir (committee chair), Leahy, Richard M. (committee member), Nayak, Krishna S. (committee member), Shung, Kirk K. (committee member)

Creator Email cwwong.alec@gmail.com,cwwong_alec@hotmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m2384

Unique identifier UC1230024

Identifier etd-Wong-3166 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-406386 (legacy record id),usctheses-m2384 (legacy record id)

Legacy Identifier etd-Wong-3166.pdf

Dmrecord 406386

Document Type Dissertation

Rights Wong, Chi Wah (Alec)

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu