Registration, Atlas Generation, and Statistical
Analysis of High Angular Resolution Diffusion
Imaging based on Riemannian Structure of
Orientation Distribution Functions
Jia Du
Department of Bioengineering
National University of Singapore
A DISSERTATION SUBMITTED IN PARTIAL FULFILLMENT OF
THE REQUIREMENTS OF
PhilosophiæDoctor (PhD)
June 18, 2013
1. Reviewer:
2. Reviewer:
3. Reviewer:
ii
Declaration
I hereby declare that the thesis is my original work and it has been written
by me in its entirety. I have duly acknowledged all the sources of
information which have been used in the thesis.
This thesis has also not been submitted for any degree in any university
previously.
Jia Du,
June 18, 2013
Acknowledgements
I would like to express my gratitude to my advisor, Dr. Anqi Qiu, for
giving me professional guidance, unending encouragement and full research
support. I appreciate all her contributions of time and ideas over the course
of my Ph.D. studies. I am also deeply grateful to our collaborator, Dr.
Alvina Goh, in the Department of Mathematics, National University of
Singapore. I thank her for the technical support and constructive advice
throughout this work.

and Unlabeled Points, IEEE International Conference on Complex Medical Engi-
neering, 2011. (oral presentation).
• Jia Du
, A. Pasha Hosseinbor, Moo K. Chung, Andrew L. Alexander, Anqi Qiu,
Diffeomorphic Metric Mapping of Hybrid Diffusion Imaging based on BFOR
iii
Signal Basis, Information Processing in Medical Imaging (IPMI), 2013.
• Jia Du
, Alvina Goh, Anqi Qiu, Bayesian Atlas Estimation from High Angular
Resolution Diffusion Imaging (HARDI), Geometric Science of Information (GSI),
2013. (oral presentation).
Conference Abstracts
• Jia Du
, Anqi Qiu, Whole Brain Diffeomorphic Mapping via the Integration of
Sulcal Curves, Cortical Surfaces, and Images, Organization for Human Brain
Mapping (OHBM), 2011, Quebec City (Trainee Abstract Travel Award, Interactive
Poster, highlighting as top ranked abstracts).
• Jia Du
, Alvina Goh, Anqi Qiu, Bayesian Atlas Estimation for High Angular
Resolution Diffusion Imaging, Organization for Human Brain Mapping (OHBM),
2012, Beijing.
• Jia Du
, Alvina Goh, Anqi Qiu, Large Deformation Diffeomorphic Metric Map-
ping of High Angular Resolution Diffusion Imaging, Organization for Human
Brain Mapping (OHBM), 2012, Beijing.
iv
Contents
List of Tables xi
List of Figures xiii
List of Symbols xix

3.5.1
Comparison of LDDMM-FA, LDDMM-DTI and LDDMM-ODF
47
3.5.2 Comparison with existing ODF registration algorithm 52
3.5.3
Computational complexity of LDDMM-FA, LDDMM-ODF
and LDDMM-Raffelt 54
3.6 Summary 55
4 Bayesian Estimation of White Matter Atlas from High Angular Resolution
Diffusion Imaging 59
4.1 General Framework of Bayesian HARDI Atlas Estimation 61
4.2
The Shape Prior of the Atlas and the Distribution of Random Diffeo-
morphisms 63
4.3 The Conditional Likelihood of the ODF Data 65
4.4 Expectation-Maximization Algorithm 67
4.4.1 Derivation of Update Equations of σ
2
and m
0
inEM 70
4.5 Results 71
4.5.1 HARDI Atlas Generation 73
4.5.2
Convergence and Effects of Hyperatlas Choice of the HARDI
Atlas Estimation 74
4.5.3 Aging HARDI atlases 81
4.5.4 Comparison with existing method 82
4.6 Summary 83
5 Geodesic Regression of Orientation Distribution Functions with its Appli-

diffusion with an orientation distribution function (ODF) that can capture
multiple axonal orientations at a voxel. For both scientific and clinical
applications, it is necessary to develop methods to represent, compare and
make correct inferences from the rich information provided by HARDI
data. However, the main challenge arises from the complexity of HARDI
data as the existing analysis frameworks based on scalar images or DTI are
unable to handle such data.
The main contribution of this thesis is providing an HARDI-based analysis
framework for the studies of white matter similarities and differences
across large populations. Under a unified Riemannian manifold of ODF,
the framework includes three components: registration, atlas generation
and statistical analysis. Firstly, we propose a novel ODF-based registration
algorithm, which seeks an optimal diffeomorphism between ODFs of two
subjects in a spatial volume domain and at the same time, locally reorients
an ODF in a manner such that it remains consistent with the surrounding
anatomical structure. Next, we develop a Bayesian probabilistic model
to estimate the ODF atlas for a specified population, which serves as a
common space to eliminate statistical bias introduced during the atlas
selection. In addition, to perform statistical inference on ODFs in the
common space, we develop an algorithm that allows for geodesic regression
on directly the manifold of ODFs and thus, avoid the loss of potential
information during dimension reduction or feature extraction from ODFs.
Finally, we apply this framework to examine the effects of normal aging in
a large group of healthy subjects spanning the adult age range. The results
show that the ODF-based framework is able to detect age-related changes
in the white matter regions where fibers cross, and thus, offer new insights
into the understanding of white matter microstructure deterioration during
normal aging.
CONTENTS
x

white matter similarities and differences across large populations. . . . 7
2.1
Illustration of the manifold of square-root ODFs. (picture token and
modified from [1])) 16
2.2 Flow equation. 19
2.3
Mapping one shape to another via the group action of diffeomorphic
transformation. On the right, the diffeomorphisim
φ
is shown on the
square grid. 19
2.4 Geodesic specified by initial momentum. 21
3.1 The role of Chapter 3 in the ODF-based analysis framework. 24
3.2
Illustration of affine transformation on square-root ODFs. (Similar to
the shape of ODF, the colors of ODF also indices the relative values of
ODF in each direction, where blue stands for low ODF value and red
for high value.) 25
xiii
LIST OF FIGURES
3.3
Examples of local affine transformations on an isotropic ODF in the
first row, an ODF with a single orientation fiber in the middle row,
and an ODF with crossing fibers in the bottom row. From left to
right, three types of affine transformations,
A
, on the ODFs are demon-
strated: in panel (a), a rotation with angle
θ
z

shape, while the image on panel (c) shows different ODFs but the
same circular image shape. Panels (d) and (e) show the deformations
(grid) and the corresponding momenta (arrows), calculated using
∇
φ
1
E
in Eq.
(3.19)
, for mapping the image on panel (a) to panels (b) and
(c), respectively. Panels (f) and (g) show the deformations and the
corresponding momenta, calculated using the gradient in our previous
work [
2
], for mapping the image on panel (a) to panels (b) and (c),
respectively. 45
3.6
Comparison between the LDDMM-ODF and LDDMM-DTI algorithms.
Panels (a, b) respectively show the template and target HARDI and their
enlarged images, where the ODF or diffusion tensor at each location
contains two crossing fibers with equal orientation distribution. Panel (c)
illustrates the template HARDI image transformed via the deformation
given in panel (d), the result of the LDDMM-ODF algorithm. Panel (e)
illustrates no deformation found via the LDDMM-DTI algorithm and
thus the template HARDI image remains. 46
xiv
LIST OF FIGURES
3.7
Panels (a-d) respectively show the maps of mean symmetrized Kullback–
Leibler (sKL) divergence of the ODFs between the template and the sub-

(blue), LDDMM-FA (cyan), LDDMM-DTI (yellow), or LDDMM-ODF
(red) are applied. 57
3.12
Dice overlap ratios averaged over all 25 subjects deformed by affine
(blue), LDDMM-FA (cyan), LDDMM-DTI (yellow), or LDDMM-ODF
(red). 57
xv
LIST OF FIGURES
3.13
Spatial distribution of mean squared error of spherical harmonics coef-
ficients across all subjects. For each voxel, the red (MSE=0.1) indicts
the difference between the template and the deformed subjects is large,
while the blue (MSE=0) indicts the two corresponding ODFs are equal. 58
3.14 Illustration of three fiber tract masks, CST-left, CST-right and CC. . . 58
3.15
Mean squared error of spherical harmonics coefficients (MSE of SH) in
the region where the CC and CST intersects each other. 58
4.1 The role of Chapter 4 in the ODF-based analysis framework. 61
4.2
Illustration of the general framework of the Bayesian HARDI atlas
estimation. 63
4.3
The evolution of
ψ
0
(s,x)
over the optimization of the atlas estimation.
Panels from left to right show
ψ
0

ODF
metric square between the
two hyperatlases (Figure 4.6 (a, c)) and between the two estimated
atlases (Figure 4.6 (b, d)) are respectively shown in the dashed and
solid lines. For each voxel, the red (metric=0.1) indicts the difference
between the template and the deformed subjects is large, while the blue
(metric=0) indicts the two corresponding ODFs are equal. 80
4.8
Comparison of HARDI atlases respectively generated from young and
old adults. In each row, the last three columns show three zoom-
in regions for branching and crossing bundles corresponding to the
anatomy given on the first panel. 82
4.9 Comparison between Bayesian and averaged atlases 84
5.1 The role of Chapter 5 in the ODF-based analysis framework. 86
5.2 Geodesic regression on manifold Ψ. 88
5.3 Illustration of synthetic ODFs for single (a) and crossing fibers (b). . . 96
5.4
Illustration for synthetic
√
ODF
data, regression result and ground
truth under four levels of noise (
M =0.1, 0.5, 1.0, 2.0
): In each panel,
each column shows the ODFs at
x
i
=0, 0.2, 0.4, 0.6, 0.8, 1
.The first
five rows illustrate the synthetic ODFs, while the next row shows the

ODF
data at four noise
levels (
M =0.1, 0.5, 1.0, 2.0
) against the number of observations
n
.
Three types of metric between ODFs: geodesic distance; L2 norm
of spherical harmonics coefficients and symmetric Kullback-Leibler
divergence are calculated, one for each row. Under one type of ODF
metrics of that row, the mean variance of synthetic data,
MSE(ψ
i
,
¯
ψ)
;
the mean squared residuals of the geodesic regression,
MSE(ψ
i
,
ˆ
ψ
i
)
;
and
R
2
of the geodesic regression are shown in each column respectively.

direction of the ODFs. 109
5.11
Interpretation for the regions between corpus callosum and corticospinal
tracts: Panel (a) shows the T1 image of the selected slice and ROIs.
Panel (b) illustrates the ODFs in the selected ROI and the peaks of
each ODF, where blue and red lines indicate the diffusion directions of
CC and CST based on their ODF values. Panels (c)-(h) show the ODF
values against the age for corresponding voxels labeled in Panel (b),
where blue denotes the ODF values in the diffusion direction of CC and
red denotes the ODF values in the diffusion direction of CST. For both
diffusion directions, the crosses (or dots) and lines represent the ODF
values for observation data in the chosen direction and the projections
of estimated regression geodesics on the chosen direction respectively. 110
xviii
List of Symbols
DWI Diffusion-Weighted Imaging
DTI Diffusion Tensor Imaging
HARDI High Angular Resolution Diffusion Imaging
PDF Probability Density Function
ODF Orientation Distribution Function
√
ODF square-root Orientation Distribution Function
EAP Ensemble Average Propagator
LDDMM Large Deform Diffeomorphic Metric Mapping
GRF Gaussian Random Field
FA Fractional Anisotropy
GFA Generalized Fractional Anisotropy
MSE Mean Squared Error
SH Spherical Harmonics
sKL symmetrized Kullback-Leibler divergence

magnetic resonance imaging (DW-MRI) technique has exploited the property that water
molecules move faster along neural axons than against them. By measuring water diffu-
sion in the brain, the location and trajectories of axons can be visualized and the axonal
pathways can be reconstructed using DW-MRI. Several techniques may be used to re-
construct the local orientation of brain tissue from DW-MRI data. A classical method is
known as Diffusion Tensor Imaging (DTI) [
3
], which characterizes the diffusivity profile
of water molecules in brain tissue by a single oriented 3D Gaussian probability distribu-
tion function (PDF). In DTI, the diffusivity profile is often represented mathematically
by a symmetric positive definite (SPD) tensor field
D : R
3
→ SPD(3) ⊂ R
3×3
that
1
1. INTRODUCTION
measures the extent of diffusion in any direction
v ∈ R
3
as
v

Dv
. The geometry of
SPD(3)
is well-studied and several metrics for comparing tensors have been proposed
[
4

D
1
(x)
and
D
2
(x)
at voxel
x
, it is no longer true that
D
1
(x) ≈ D
2
(φ(x))
and each tensor must be reoriented in such a way that it remains consistent with the
surrounding anatomical structure. There exist several approaches for reorientation
that are used in DTI [
12
]. For instance, the Finite Strain (FS) scheme decomposes
an affine transformation matrix
A
into
A = RS
, where
R
is the rigid rotation and
S
is the deformation, and reorients the tensor
D

the population that well characterizes the in-vivo white matter anatomy of the human
2
1.1 Motivation
brain. To this end, there are several works on DTI-based atlas generation. Park et al.
constructed an atlas with full tensor information in [
19
]. Since then there have been
several different approaches to DTI atlas construction, either using scalar registration
[
20
,
21
], multi-channel methods [
19
] or by directly optimizing tensor similarity [
22
].
Since then, there have been works that leveraged the DTI template information and some
of these include several exploratory works that provided anatomical validation [
23
] and
anatomical labeling of fiber tracts [
24
,
25
], and evaluation. The comprehensive work
by Mori et al. [
24
] provides a three-dimensional and two-dimensional in-vivo atlas of
various white matter tracts in the human brain based on DTI and has become an essential