University of Southern California Dissertations and Theses

Novel computational techniques for connectome analysis based on diffusion MRI

(USC Thesis Other)

Novel computational techniques for connectome analysis based on diffusion MRI

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Novel Computational Techniques for Connectome Analysis Based on Diffusion MRI
by
Xinyu Nie
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(ELECTRICAL ENGINEERING)
May 2024
Copyright 2024 Xinyu Nie

Acknowledgements
I would first express my sincere gratitude to my advisor, Professor Yonggang Shi. His deep knowledge of
diffusion MRI inspired me to build my research focus. His algorithms and programming experience helped
me realize my research projects smoothly. His ability to see the big picture hugely influenced me to explore
and think deeply about the research areas. His critical thinking helps me figure out my work’s weaknesses
and find solutions to improve them. He was always available and patient when I was in need of his help
and advice. His insights into the academic world inspired what I need to be a successful researcher.
I would like to thank Professors Richard M. Leahy, Cauligi S. Raghavendra, and Danny J.J. Wang for
joining my Ph.D. dissertation committee and for their constructive feedback and comments on my research
works. I thank Professors Paul M. Thompson and Anand A. Joshi for joining my Ph.D. qualifying exam
committee and for their suggestions on improving my research work and presentation skills. I am grateful
for Dr. John M. Ringman, who provided us with important datasets, without which we could not have
developed research for Alzheimer’s disease.
I was fortunate to join the Laboratory of NeuroImaging (LONI), which provides me with an excellent
environment for studying and working. Since joining LONI, Yihao Xia has been a close friend and colleague. We often communicated our research ideas and plans, and he gave me valuable suggestions for my
research and career plans. I want to thank Jialiang Ruan, who helped me build up the software development, which will be valuable in our future research. I very much enjoyed working with Yuan Li, and we
collaborated to develop several works more efficiently. I enjoyed conversations and communication with
ii

Jianwei Zhang and Yuan Li, and I was less stressed and had more time to share ideas with them. They also
helped me to organize my qualifying exam and defense. I would like to thank my friends, Jiakang Wang
and Bo Wu, who gave me suggestions and intuitions to build up my mathematical frameworks. I want to
thank all my colleagues and friends for their support and encouragement.
I want to sincerely thank my parents, who have supported me despite the long distance between China
and the USA. I would not have begun the Ph.D. program without their support. They always encouraged
me and provided me with advice during my Ph.D. study. I also want to thank my other family members
for their help and support of my parents.
iii

Table of Contents
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii
Chapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Objective of This Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 Organization of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.5 Relevant Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
Chapter 2: Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1 Principle of Diffusion MRI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Diffusion Tensor Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Multiple Fiber Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.4 Spherical Harmonics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.5 FOD with Compartments Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.6 Tractography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.7 Opitmal Transport and Wasserstein Space . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Chapter 3: Surface-based Probabilistic Fiber Tracking in Superficial White Matter . . . . . . . . . . 21
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2.1 Mesh Representation of SWM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2.2 FOD Projection onto SWM Surface . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2.3 Parallel Transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.2.4 Tracking Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3 Experiments and Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.3.1 U-fiber Tracking on HCP Subjects . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.3.2 Validation with 7T Post-mortem MRI . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.3.3 U-Fiber Connectivity Changes in ADAD . . . . . . . . . . . . . . . . . . . . . . . . 44
3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.4.1 Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
iv

3.4.2 More General SWM Connectivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.4.3 Code and Computation Speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.4.4 Clinical Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
Chapter 4: Flow-based Geometric Interpolation of Fiber Orientation Distribution Functions . . . . 51
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.2.1 FOD Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.2.2 Modeling Single Peak FOD Components as Flow of Vector Fields . . . . . . . . . . 55
4.2.3 Rotation Calculation for SPHARM-based FODs . . . . . . . . . . . . . . . . . . . . 55
4.2.4 Evaluation Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.3 Experiment Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
Chapter 5: A Geometric Framework for Processing Fiber Orientation Distribution Functions . . . . 62
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.2.1 Challenges of FOD Functions Processing . . . . . . . . . . . . . . . . . . . . . . . . 65
5.2.2 FOD Function Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.2.3 New Metric and Interpolation for FOD functions . . . . . . . . . . . . . . . . . . . 73
5.2.4 Numerical Computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.2.5 Fast Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.3 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.3.1 Synthetic Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.3.2 HCP Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.3.3 Application in Tractography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
Chapter 6: Topographic Filtering of Tractograms as Vector Field Flows . . . . . . . . . . . . . . . . 87
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
6.2.1 Vector Field Flows for Topographically Regular Tracts . . . . . . . . . . . . . . . . 89
6.2.2 Vector Field Computation from FODs . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.2.3 Measure for Topographic Regularity . . . . . . . . . . . . . . . . . . . . . . . . . . 93
6.3 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
Chapter 7: Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
v

List of Tables
3.1 Algorithm 1. Surface-based Probabilistic Tractography . . . . . . . . . . . . . . . . . 33
5.1 Distances between Probability distributions . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.2 Algorithm: Geometric Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.3 Mean Distances between the interpolated FOD functions and the ground truth . . . . . . . 81
vi

List of Figures
3.1 Accurate placement of the WM cortical surface in the dMRI space. (a) T1-weighted
MRI; (b) the extra-axonal diffusivity calculated from dMRI by our compartment models.
After the preprocessing by HCP-Pipeline, yellow arrows in (c) highlight the residual
misalignment of the WM cortical surface (cyan) and the intra-axonal tissue map calculated
by our compartment models. (d) shows the much-improved alignment of the WM cortex
(green) with the same tissue map after nonlinear registration between T1-weighted MRI
in (a) and the extra-axonal diffusivity map in (b). . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2 The projection of 3D FODs onto 2D triangles local coordinate systems. (a) The F OD3D
function is represented in a local spherical coordinate system defined on the center of
a triangle of the SWM mesh. The reference plane of the system is the triangle plane;
meanwhile, the zenith direction z is aligned with the normal direction. The x − axis
is aligned with an edge of the triangle. (b) The projected F OD2D using the Eq. 3.1.
The F OD2D function is parametrized by ϕ in the local coordinate system and can be
normalized to a probability distribution function on the unit circle S
1
. . . . . . . . . . . . . 25
3.3 The local coordinate system Oxyz differs from the physical coordinate system Ox0y
0
z
0
by a rotation R. (θ0, ϕ0) and (θ, ϕ) represent the same point p in the two systems. The
rotation R is decomposed as three successive rotations around three axes: (a) Ox0y
0
z
0
is rotated by an angle α around Oz0 → Ox1y
1
z
0
; (b) Ox1y
1
z
0
is rotated by a second
angle β around Oy1 → Ox2y
1
z
2
; (c) Ox2y
1
z
2
is rotated by a third angle γ around
Oz2 → Ox3y
3
z
2
, which is Oxyz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.4 For a vector up in the tangent space of point p on a Riemannian manifold, the Levi-Civita
connection parallel transports up along a geodesic γ(t) to a vector uq in the tangent space
of point q, which allows the comparison of vectors in the tangent spaces of nearby points. 31
3.5 The parallel transport of a vector v in the current triangle T1 to a neighboring triangle T2.
The parallel transport consists of unfolding the triangles to a common plane (a) → (b),
translating the red tangent vector v from T1 to T2 (b), and folding the triangles back to
their original positions (b) → (c). The transported vector vp can then be compared with
any sampling direction w in the triangle T2. . . . . . . . . . . . . . . . . . . . . . . . . . . 32
vii

3.6 Representative HCP examples of U-fibers between the pre-and post-central gyrus
reconstructed by the two surface-based methods (SP: surface-based probabilistic tracking;
SD: surface-based deterministic tracking) and the MRtrix software. For better visualization,
we downsampled all results to 1000 streamlines. The subject ID is plotted at the top of the
results from each HCP subject. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.7 Quantitative comparisons of the U-fibers reconstructed by the proposed surface-based
probabilistic (SP) fiber tracking, surface-based deterministic (SD) fiber tracking, and
volume-based fiber tracking in the MRtrix software on 484 HCP subjects. (a) Box plots of
the ‘well-U-connected’ measure of the three methods. Box plots of the ‘well-distributed’
measure on the precentral and postcentral gyrus are shown in (b) and (c), respectively. (d)
Box plots of the U-ratio measure across subjects. (e) Box plots of the Procrustes distance
are used to assess the topographic regularity of the reconstructed U-fibers using the three
methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.8 Across the HCP subjects, box plots of the number of touched sections in (a) pre- and (b)
post-central gyrus by MRtrix results with varying numbers of seeds from 3000 to 50,000
are shown. Box plots of the number of touched sections by results from the proposed
surface-based probabilistic (SP) method (30,000 seeds) are also shown for ease of comparison. 39
3.9 T2-weighted 7T MRIs of two post-mortem brains are shown in (a) Subject 1 (91-year-old
male) and (b) Subject 2 (97-year-old female). SWM with higher iron concentrations is
highlighted with blue arrows. Manually delineated SWM masks in the lateral frontal lobe
for subjects 1 and 2 are shown in (c) and (d), respectively. . . . . . . . . . . . . . . . . . . . 41
3.10 Overlay of reconstructed U-fibers with the SWM on T2-weighted post-mortem MRIs. Top
row: subject 1; Bottom row: subject 2. (a) and (d) show a zoomed view of the overlay
of the U-fibers with the MRI in the red box highlighted in (b) and (e), respectively. (c)
and (f) show a zoomed view of the overlay of the U-fibers with the MRI in the blue box
highlighted in (b) and (e), respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.11 (a) The mean distance from the U-fibers reconstructed by the proposed algorithm to
the SWM mask as a function of the deformation distance parameter δ. Boxplots of the
distance of the reconstructed U-fibers to the SWM masks in (b) subject 1 and (c) subject 2.
MRtrix: iFOD1 algorithm in the MRtrix software. SP: proposed surface-based probabilistic
tracking method with δ = 0.5mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.12 The reconstructed U-fibers in the parietal lobe of two ADAD patients are superimposed
over their SWM meshes. The left and right hemispheres of subject 1 are shown in (a) and
(b), respectively. The left and right hemispheres of subject 2 are shown in (c) and (d),
respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.13 The linear regression results between the CDR-SOB score, tau SUVR, and the mean
F OD2D peak value on the U-fibers of the parietal cortices of the ADAD patients. (a), (c)
are results from the left hemisphere and (b), (d) are results from the right hemisphere. . . . 47
viii

4.1 An illustrative example of FOD interpolations. (A) FODs of four neighboring voxels from
a bending fiber bundle were highlighted by red circles in (B) and (C). The blue arrow
shows the direction of the fiber bundle. (B) Interpolated FODs are generated by linear
interpolation, where artificial peaks are generated. (C) The interpolated result by our
proposed method correctly accounts for rotation and follows the bending geometry of the
fiber bundle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.2 A FOD function with two peak lobes is decomposed into two SPHARM-based FOD
functions, each with only a salient peak lobe. . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.3 A vector field within a tube centering at a voxel P0. (A) Tube T
k
p0
(yellow box) centering
at p0 along the peak direction vp0 of the k
th FOD component F
k
p0
(green). (B) The FODs
within the tube and the {F
k
pt
} are in green. (C) The vectors {vpt} are picked at voxels
within the tube, and these vectors are used to solve Eq. 4.3 and compute the supporting
vector field V
k
p0
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.4 In the flowchart, the original FODs are highlighted by red circles at the four corners. The
FOD functions are decomposed into single-peak components in the first row and locally
fitted by vector field flows. In the second row, we rotate the single-peak components at
four corners to align with the vectors at each target point, and then we compute their
weighted mean at each target point. Finally, we combine the single-peak components into
complete FODs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.5 One interpolated slice of the FODs. For FODs within the red box, a zoomed view of the
proposed interpolation, the linear interpolation, and the ground truth are plotted. . . . . . 58
4.6 Boxplots from quantitative comparisons using data from 40 HCP subjects. (A) Mean FAHM
of the ground-truth FODs, upsampled FODs from the proposed and linear interpolation.
(B) The relative error between up-sampled FODs and ground truth. (C) Procrustes
distance of CST bundles from fiber tracking based on linear interpolation (original) and
our up-sampled FODs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.7 Tractography of CST on original HCP data and upsampled data for three subjects, and for
each bundle, we downsample the number of streamlines to 1000 for visualization. . . . . . 61
5.1 A FOD function with three peak lobes plotted in spherical radiation pattern (a) and on the
unit sphere (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5.2 Probability distribution functions for the illustration of different distances. (a) shows
three disjoint distributions with single sharp peak lobes. (b) shows two single peak lobe
distributions and a distribution with two peak lobes. . . . . . . . . . . . . . . . . . . . . . . 66
5.3 Interpolations between Probability distribution functions based on different distances.
Two single peak lobe distributions are shown in (a). (b) to (e) show the interpolations
between them based on different metrics. One distribution with two peak lobes and
another with three are shown in (f). (g) to (j) show their interpolations. . . . . . . . . . . . 69
ix

5.4 (a) The FOD function f with peak lobes mixed. (b) Deconluluted function fs from the
FOD function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.5 The FOD function with multiple peak lobes is decomposed into five single-peak component
FOD functions, and the decomposed components are rescaled for better visualization. . . . 73
5.6 The flowchart of the geometric interpolation framework. The blue, purple, green, and
red arrows mean hierarchical clustering, peak lobe field generation, interpolation, and
summation. The FOD functions from neighboring voxels are decomposed and clustered
at the first row. Each peak lobe is used as the generator for a peak lobe field within the
same cluster. The same color border highlights the peak lobe field and its generator peak
lobe. Interpolation is computed within each peak lobe field from the second to the third
row. Finally, we combine interpolations from all peak lobe fields to get the complete FOD
functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.7 The interpolated FOD functions between the first peak lobe (first column in each figure)
and the second peak lobe (last column in each figure) using different methods. The
second peak lobes are rotated to perform different angle differences, shown on the left,
to the first peak lobes. ’Linear’ represents the linear interpolation, ’RF’ is the Riemannian
Framework, ’WB’ is the Wasserstein Barycenter, and ’FA’ is the Fast Approximation. . . . 79
5.8 The Full Area at Half Maximum (FAHM) of interpolated FOD peak lobes in the middle
columns of each figure in Fig. 5.7. ’Linear’ represents the linear interpolation, ’RF’ is
the Riemannian Framework, ’WB’ is the Wasserstein Barycenter, and ’FA’ is the Fast
Approximation, ’Baseline’ represents the baseline peak lobe. (a) and (b) are the 8
th and
16th FOD functions, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.9 The interpolated FOD functions of the synthetic half-cycle data are based on different
methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.10 The boxplots of the distances between the interpolated FOD functions and the ground
truth FOD functions in Fig. 5.9. ’Linear’ represents the linear interpolation, ’RF’ is the
Riemannian Framework interpolation, ’WB’ is the Wasserstein Barycenter, and ’FA’ is the
Fast Approximation. (a) and (b) are the 8
th and 16th FOD functions, respectively. . . . . . 81
5.11 One interpolated slice of the FODs. For FODs within the red box, a zoomed view of the
Ground truth, linear interpolation, Wasserstein Barycenter, and Fast Approximation. . . . 82
5.12 Boxplots of the FAHM of the interpolated FOD functions in (a) and the distance between
the interpolated FOD functions to the ground truth FOD functions in (b). ’Linear’
represents the linear interpolation, ’GR’ is the ground truth, ’WB’ is the Wasserstein
Barycenter, and ’FA’ is the Fast Approximation. . . . . . . . . . . . . . . . . . . . . . . . . 83
5.13 The tractography results on the synthetic half-circle data based on different interpolations. 84
5.14 Tractography of CST on original HCP data and upsampled data for three subjects, and for
each bundle, we downsample the number of streamlines to 200 for visualization. . . . . . . 85
x

5.15 Boxplots of the Procrustes distances for CSTs of 40 subjects, HCP represents the
reconstructed CST on the original HCP images, and Upsampled presents the reconstructed
CST on the Upsampled super-resolution images. . . . . . . . . . . . . . . . . . . . . . . . . 86
6.1 A vector field flow representation for the visual pathway bundle. (A) Fiber bundle from
FOD-based tractography. The red box is the ROI used for visualization in (B)-(D). (B) The
FODs in the ROI. (C) The multiple peaks of the FODs. (D) The vector field extracted from
the FOD peaks. Note that we use the QIT software ( http://cabeen.io/qitwiki) to visualize
all FOD peaks and vector fields as sticks in this work. . . . . . . . . . . . . . . . . . . . . . 89
6.2 Main steps in our topographic tract filtering algorithm. (A) Input tracts of the callosal
motor (CM) pathway from FOD-based tractography. (B) The field of MVD from FOD
peaks for the ROI highlighted by the yellow box in (A). (C) The PVF for the CM pathways
is computed by solving the MRF model. (D) The preserved fiber tracts after removing
30% of tracts with high VFD value. (E) For the ROI highlighted by the white box in (A),
the 10% of tracts with the highest VFD values are plotted as white tubes together with
the original input tracts. This helps visualize the irregularity of the filtered tracts by our
algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.3 Representative tract filtering results from 3 HCP subjects for the visual (left) and CM (left)
pathway. For each case, regions with outlier tracts were highlighted on the input tracts. . . 94
6.4 Box plots of the average G2
SD measures of the removed and preserved bundle of 30 HCP
subjects. (A) Visual pathway; (B) CM pathway. . . . . . . . . . . . . . . . . . . . . . . . . . 95
xi

Abstract
Diffusion MRI (dMRI) is an essential non-invasive imaging tool that allows us to study the microstructure
and connectivity of the human brain in vivo. Tractography, a critical technique for reconstructing fiber
connections based on dMRI, has become a critical tool in brain connectome research, but it is error-prone
and suffers from various limitations. In this thesis, we will develop novel computational techniques to
overcome several limitations in current tractography techniques and advance the state-of-the-art in human
connectome research.
A critical challenge for tractography in connectome research is reconstructing fiber connections in
the superficial white matter (SWM) distributed beneath the cortical surface. Indeed, complex structures
and the high curvature of the fibers in SWM lead to a false-negative reconstruction in the conventional
tractography algorithms. Another weakness of tractography is that it generates enormous anatomically invalid connections that reduce the accuracy of the connectivity analysis, even with advanced mathematical
models such as fiber orientation distribution (FOD). In this thesis, we propose a novel surface-based tractography algorithm that successfully reconstructs the fiber connectivity in the SWM. On the other hand,
we develop geometric frameworks to process and analyze the FOD functions and provide novel interpolations of the FOD that help reduce the implausible propagations in tractography. Finally, we proposed a
post-processing vector-flow-based filter to reduce the unrealistic streamlines in tractography.
To overcome the difficulties of traditional volume-based tractography in SWM, we propose a novel
surface-based framework for the probabilistic tracking of fibers on the triangular mesh representation of
xii

the SWM. By deriving a closed-form solution to transform the spherical harmonics (SPHARM) coefficients
of 3D fiber orientation distributions (FODs) to local coordinate systems on each triangle, we develop a
novel approach to project the FODs onto the tangent space of the SWM. After that, we utilize parallel
transport to realize the intrinsic propagation of streamlines on SWM following probabilistically sampled
fiber directions. Our intrinsic and surface-based method eliminates the need to perform the necessary but
challenging sharp turns in 3D compared with conventional volume-based tractography methods. We performed quantitative comparisons to demonstrate the proposed algorithm can more effectively reconstruct
the fiber connections in SWM than previous methods.
Interpolation plays a vital role in the tractography algorithms since the step size of the tracking is much
smaller than the resolution of the dMRI data. Previous interpolation methods of FODs often generate false
artificial information, leading to anatomically incorrect propagations in tractography. We propose a flowrotation-based interpolation framework based on the local geometry of the fiber bundles. Our approach
decomposes each FOD function with multiple peak lobes into single-peak components, and neighboring
single-peak components consistent in direction are modeled as vector field flows. The vector field flow
representation enables geometrically consistent interpolation of FOD functions in a rotation-based framework. Experiments demonstrate that our method produces anatomically meaningful FOD interpolation
and improves state-of-the-art tractography.
We developed an advanced geometric framework to improve the flow-rotation-based framework and
put the essential concepts in more efficient and mathematically rigorous frameworks. The decomposed
single-peak FOD functions are embedded and processed in Wasserstein space, a well-defined metric space
with excellent properties. We define a new metric to compare and analyze the FOD functions, and the
interpolation is computed as the Barycenter of the new metric. For efficient computation, we propose a
Fast Approximation, which is a generalization of the flow-rotation-based interpolation, for the Barycenter.
In experiments, we use synthetic and real datasets to demonstrate that the new geometric framework can
xiii

handle FOD functions representing complicated fiber geometry and provide accurate interpolations that
improve tractography performance.
To reduce the invalid streamlines in the post-processing of tractography, we proposed a filter based
on the vector flow representation, which is used to model the local fiber geometry in flow-rotation-based
interpolation. The vector flow representation is consistent with the topographic principle in brain connections, and it can be naturally extended to the global representation of the topographically regular fiber
bundles. Thus, we developed a framework that models the streamlines of fiber bundles as vector field flows.
We compute a fiber bundle’s principal vector field (PVF) by solving a Markov Random Field problem based
on the FOD peaks. We propose a metric named Vector Flow Deviation (VFD) for tractography by quantifying the consistency between the PVF and the tractography streamlines and applying the VFD to filter out
topographically irregular streamlines. In the experiments, we visually and quantitatively show that our
method successfully removes the fiber bundles’ outliers while significantly improving the topographical
regularity.
In this thesis, we developed four novel computational techniques for connectome analysis: a surfacebased tractography for reconstructing SWM connections, a flow-rotation-based FOD interpolation that
overcomes the difficulties in commonly used interpolations, a geometric framework to process and analyze
the FOD functions efficiently, and a filter that reduces the invalid connections in tractography. All four
techniques overcome critical difficulties in tractography and connectome analysis and outperform stateof-the-art techniques in visual and quantitative comparisons. We develop software tools for the proposed
algorithms and will distribute them publicly to the research community.
xiv

Chapter 1
Introduction
1.1 Motivation
Diffusion magnetic resonance imaging (dMRI) [1] is a vital imaging technique to non-invasively study
biological tissue microarchitecture based on the diffusivity of water molecules within the human body. The
dMRI is widely used to reveal the structural connectivity of the brain, the development and degeneration of
the nervous system, and pathological changes in brain tissue. In the last decades, advanced improvements
in image acquisition techniques significantly increased the spatial and angular resolution of dMRI [2, 3, 4]
and provided preconditions for advanced mathematical models such as fiber orientation distribution (FOD)
[5, 6, 7]. These models facilitated the development of the dMRI-based technique, tractography, which is the
only non-invasive tool for studying the brain’s white matter pathways in vivo. It has become the preferred
tool for reconstructing connections between brain regions and studying the brain’s underlying structure.
Various tractography algorithms have been developed for connectome research. However, the reliability
of tractography is a critical problem; studies [8] show that a large number of fiber connections are not
reconstructed by tractography (false negatives); meanwhile, it is prone to develop numerous anatomically
invalid streamlines (false positives) [9]. Critical limitations reduce the reliability of connectome research
and need improvement.
1

• Incomplete Fiber Connections in SWM. In recent years, more and more attention has been paid
to the fibers in the superficial white matter (SWM) [10]. SWM has approximately twice as many
fiber connections as the DWM [11] and is indispensable to constructing the entire cerebrum network [12]. However, as mentioned, tractography suffers from false negatives; in particular, it is hard
to provide satisfactory reconstructions of fibers within the SWM [10]. Indeed, tractography needs
strong regularizations in tracking procedures to generate anatomically realistic fibers since noise
and complex fiber geometry caused ambiguities and inaccuracies in dMRI data [9]. In traditional
tractography, the algorithms control the turn angle at each tracking step to avoid weird and sharp
turns of the streamlines. The most used regularization is the continuity of the tangent direction
of the tracking; meanwhile, some algorithms use polynomials to model the smooth propagation of
the tracking. However, the strong regularizations prevent numerous potential highly curved connections in the SWM, which lies directly under the cortical layer. Thus, current tractography can
only provide a small partial representation of connectivity in SWM. There is still a lack of advanced
tractography techniques to reconstruct the fibers in SWM effectively.
• False Positive Problem. Another critical limitation of tractography is that they are prone to generate numerous false positive streamlines [9]. Each streamline of the tractography starts at a seed
point and follows the local diffusion information on a step-by-step propagation. However, the local
diffusion data may suffer from noise and low angular resolution. They cannot represent the fiber
direction accurately, and these factors lead to deviations in the propagation of streamlines, which are
invalid connections. Moreover, the local diffusion information can represent complex fiber geometry such as fiber crossing, kissing, and fanning situations [9]; the tractography cannot distinguish
these fiber geometries without manually selecting priors. These complicated fiber geometries often
lead to anatomically unrealistic connections between brain regions and reduce the reliability of connectome research. To improve the reliability of tractography, efficient methods still need to alleviate
2

the errors in the propagation of the streamlines in the tracking procedure and reduce unrealistic
streamlines in post-processing.
• Lack of Framework for Processing FOD functions. The fiber orientation distribution (FOD),
which is an advanced model representing the complicated cross or bifurcate fiber geometry and
quantifying the fraction of fiber portions [5, 6], has become more popular in the last decade. However, fewer works have been proposed for analyzing the FOD data, and no proper framework exists
for processing FOD functions. Thus, people cannot compare and analyze the FOD data efficiently,
and general operators on data, such as interpolation, smoothing, and convolution, can not be defined. In particular, interpolation, the essential operation in image processing [13], has received rare
attention. The interpolation, which is necessary for FOD-based tractography, highly affects the performance of tractography. A unified framework is still necessary for the process and interpolation
of the FOD functions.
1.2 Objective of This Thesis
This thesis aims to develop a new tractography algorithm that efficiently reconstructs the fiber connections in SWM. Meanwhile, we will develop new frameworks to process and analyze the fiber orientation
distribution (FOD) functions and propose novel interpolations to reduce propagation errors in FOD-based
tractography algorithms. Finally, we expect a new filer to reduce the invalid streamlines in the postprocessing of tractography.
SWM Tractography. Previous works on fiber tracking in SWM primarily rely on existing tractography methods, and these methods can be characterized into two categories. In the first category, cortical and
volume ROIs were used to generate and select fiber tracts for specific bundles of U-fibers [14, 15, 16, 17].
In the second category, whole-brain tractography was applied, then they clustered U-fibers throughout
3

the cortex [18, 19, 20, 21]. One known limitation of existing volume-based tractography methods is that
they tend to have many false negatives in the reconstruction of highly curved fibers, particularly the fibers
in SWM. Previous tractography methods can, thus, typically only provide a partial representation of the
fiber connections in the SWM. Recently, increasing interest has been in advancing the tracking of SWM
fibers by leveraging the surface representation of the SWM. A surface-based seeding approach [22] was
proposed recently to enhance the reconstruction of fibers in SWM. However, it still relies on conventional
volume-based tractography for fiber tracking and cannot prevent the reconstructed tracts from deviating
away from the sheet-like structure of the SWM. To overcome the difficulties of traditional tractography
in reconstructing the fibers within the SWM layer, we will develop a novel surface-based probabilistic
tracking framework.
Frameworks for Processing FODs. The fiber orientation distribution (FOD) represents complicated
fiber configurations as a function on a sphere and has been used widely for fiber tractography [6]. However, no efficient framework exists for processing and analyzing the FOD functions, and critical operators
on data, such as interpolation, smoothing, and convolution, can not be developed for FOD functions. In
particular, interpolation, the essential operation in data processing, has received rare attention. Interpolation of FOD functions is essential in tractography since the step size of the tracking is much smaller than
the dMRI resolution, and improper interpolations lead to error propagation in the tracking procedure [23].
In popular FOD-based tractography, linear interpolation is commonly adopted for numerical efficiency.
Still, it often generates artificial directions and ignores rotations. Several rotation-based frameworks [24,
25] proposed alleviating linear interpolation problems by simultaneously averaging the shape and rotation of FODs. However, since only one rotation is used for the whole FOD function, these methods cannot
handle more general situations where individual FOD peaks experience different rotations; in particular,
different FOD peaks may belong to different fiber bundles. We will develop novel frameworks to process
4

and analyze the FOD functions and provide new interpolation techniques to overcome the difficulties in
previous works and reduce the errors in the propagation of tractography.
Tractography Filter. The tractography algorithms are prone to generate enormous invalid connections, reducing the reliability of the connectivity analysis [9]. Thus, numerous tractography filters are
proposed and used in post-processing tractography data. These filtering methods can be characterized
into several categories. Geometric-based methods develop distance-based metrics, which measure the
geometrical differences between streamlines, to evaluate and filter the streamlines [26, 27]. Registrationbased filters register the subject streamlines to a common space with atlas tractography, and they remove
the streamlines mismatch with the atlas [28, 29]. Some filtering algorithms, such as SIFT [30], are devoted to improving the coincidence between the tractography and the dMRI data to reduce the artifacts
induced by tractography algorithms. Recently, a new filter based on autoencoder learning was proposed
[31]; they train the autoencoder and represent the streamlines in the latent space, and then the undesired
streamlines are filtered using the nearest neighbor method. However, outliers recognized by the distancebased metrics, registration, and latent space are not the only invalid streamlines, and these streamlines
that deviate from the topographic organization of the fiber bundle can also be invalid. The topographic
regularity of the fiber bundles [32, 33, 34] allows us to improve the validity of the tractography. We will
develop a vector-flow-based filer that utilizes the topographic principle of fiber bundles to reduce irregular
streamlines in tractography.
1.3 Contributions
We develop a novel surface-based probabilistic tracking framework for reconstructing fiber connections
in SWM. Our method systematically tackles challenges in existing volume and surface-based methods
for SWM tracking. To track the fibers accurately on the SWM layer, we carefully generate the surface
representation of SWM based on nonlinear registration and deformation of the white matter (WM) surface
5

computed from T1-weighted MRI. After that, we developed a novel method to project the 3D FODs onto the
tangent space of the SWM surface based on an efficient transformation of spherical harmonics (SPHARM)
representation of FODs to the local coordinate system of the tangent space. A probabilistic tracking method
is then developed based on the projected 2D FODs on the SWM surface, where parallel transport is adopted
to propagate the streamlines across triangles and regularize tract smoothness intrinsically. By tracking
based on the intrinsic geometry of the SWM surface, our method avoids the sharp turning angles of SWM
in 3D space and can significantly improve the reconstruction of fibers in the SWM. We demonstrate that
our proposed method outperforms previous tracking algorithms in experiments in terms of providing more
accurate and complete fibers in SWM. Furthermore, our method provides excellent opportunities to study
the variability caused by pathology and neural degeneration in SWM.
We propose a flow-rotation-based interpolation framework based on the local geometry of the fiber
bundles. We decomposed each FOD function with multiple peak lobes into single-peak FODs since they
represent different fiber bundles. Then, we model the consistent single-peak FODs in the neighborhood as
vector field flows fitted by smooth polynomials. Each vector field flow represents the geometry of a fiber
bundle locally and continuously determines the direction of the fibers. This representation enables us to
develop a rotation-based FOD interpolation, which avoids the loss of sharpness of FODs. In experiments,
the proposed interpolation achieves superior performance than the commonly used linear interpolation
approach by providing more realistic interpolated FODs. Furthermore, our interpolation method can be
used to perform upsampling of FOD volume images and significantly improve the performance of FODbased tractography.
We develop a novel geometric framework to process and analyze the FOD functions. This work inherits the key spirits from the flow-rotation-based framework and provides a mathematically rigorous
framework for the FOD functions. We developed a more rigorous and accurate spherical deconvolution
framework to decompose the FOD function into single-peak lobes and an elegant clustering framework to
6

handle the single-peak lobes. We embedded the single-peak lobes into the Wasserstein space, which has
good properties for analyzing and processing the FOD data. We define a new metric for FOD functions,
and interpolation of FOD functions is defined as the Barycenter of the new metric. To compute the interpolation faster, we developed a Fast Approximation that is a generalization of the flow-rotation-based
interpolation. The novel geometric framework can successfully represent the fiber geometry, provide accurate interpolations, and improve the FOD-based tractography.
Besides reducing errors in tracking propagation, we also propose a filtering technique to reduce invalid
streamlines in the post-processing of tractography. We extend the idea of locally representing the fiber
geometry using vector field flow to the entire fiber bundle. We extract an underlying vector field supporting
the whole topographically regular fiber bundle. This work is achieved by formulating and solving a Markov
Random Field (MRF) problem from the peak fields of the FODs. After that, we develop a quantitative metric
called vector flow deviation (VFD) to evaluate how well the fiber streamlines fit as the vector field flows.
Unlike previous filtering algorithms, which need registration streamlines of subjects to a common space
for evaluation or manually defining metrics to measure the distance between streamlines. Our data-driven
FOD filtering utilizes the fiber bundle’s intrinsic geometry to evaluate the streamlines. In experiments, we
demonstrate our method of filtering the tractography fiber bundles with known connectopy, such as the
visual and callosal motor pathways. Our filter successfully removes outliers from the tractography and
significantly improves the regularity of the fiber bundles.
1.4 Organization of the Thesis
Chapter 2 covers the background of diffusion MRI, including the principle of diffusion MR imaging technique, advanced mathematical models in dMRI, and tractography. Necessary mathematical tools are also
introduced. Chapter 3 presents the surface-based probabilistic tractography in superficial white matter.
Chapter 4 presents the flow-rotation-based interpolation for fiber orientation functions. Chapter 5 presents
7

an elegant framework for processing the FOD functions. Chapter 6 covers the vector-flow-based filer that
reduces the invalid streamlines in tractography. Finally, Chapter 7 summarizes the research works and
provides potential research.
1.5 Relevant Publications
Chapter 3 has been published in IEEE Transactions on Medical Imaging [35]. Chapter 4 is published in
The Medical Image Computing and Computer Assisted Intervention Society (MICCAI) 2023 [36]. Chapter
6 has been published on MICCAI 2019 [37].
8

Chapter 2
Background
2.1 Principle of Diffusion MRI
Diffusion is the thermal motion of particles and results in molecular or particle mixing. The diffusion
process of water molecules is the key to the Diffusion MRI. In free water, the squared displacement of
molecules from their starting point over a time t averaged over all the molecules can be characterized by
the Einstein equation [38]:
< r2 >= 6Dt, (2.1)
where < r2 > is the mean squared displacement and the diffusion coefficient is denoted by D. Since we
measure the diffusivity in human tissue instead of free water, the water molecules encounter hindrances,
such as cell membranes and macromolecules. We often refer to the apparent diffusion coefficient (ADC),
which is lower than the diffusion coefficient D.
In MRI, we apply an external magnetic field B0, and then the magnetic dipoles of hydrogen nuclei, also
called spins, will be aligned with the static field B0. After an excitation with a 90-degree radio frequency
(RF) pulse, the spins will process around the B0 at a Lamor frequency ω within the plane perpendicular to
the direction of B0:
9

ω = γB0 (2.2)
where γ is the gyromagnetic ratio. The rotating magnetic field generated by coherently precessing spins
induces a current, which is received by the coils and is the signal to generate the MR image [39]. Then, we
can apply a diffusion weighting gradient field, whose strength changes linearly by position, to the rotated
spins. The molecules’ diffusion process leads to the spins’ dephase since the molecules move to different
locations and process at different Lamor frequencies. The dephasing of the spins causes attenuation of the
received signal. The diffusion signal can be described by:
S(b) = S0 exp(−b.ADC) (2.3)
where S0 is the MR image without any diffusion weighting, ADC is the apparent diffusion coefficient,
and b is "b-factor". For the pulsed gradient spin echo (PGSE), the b-factor is given by the Stejskal–Tanner
expression [40]:
b = γ
2G
2
δ
2
(∆ −
δ
3
), (2.4)
where γ is the gyromagnetic ratio, G is the amplitude of the magnetic field gradient pulses, and the δ and
∆ are their duration and temporal separation, respectively.
2.2 Diffusion Tensor Model
The ADC in Eq. 2.3 depends on the direction of the diffusion weighting gradient field. A single ADC
is inadequate when ordered tissue is in our sample. However, multiple measurements of ADC cannot
10

integrally characterize the ordered tissue. Thus, more advanced mathematical models are needed to characterize diffusion measurements from multiple directions. In [1], the diffusion tensor imaging (DTI) was
proposed to model the diffusion process using a multivariate Gaussian model, and the diffusion signal in
Eq. 2.3 now can be expressed as:
S(g, b) = S0 exp(−bg
T Dg). (2.5)
where g is a unit vector corresponding to the gradient field direction, and the diffusion tensor D is a 3 × 3
symmetric and positive definite matrix:
D =








Dxx Dxy Dxz
Dxy Dyy Dyz
Dxz Dyz Dzz








, (2.6)
where the diagonal elements of this matrix correspond to diffusivities along three Cartesian axes, and the
off-diagonal elements correspond to the correlation between those orthogonal axes.
Six unknown coefficients exist since the matrix D is symmetric. We need at least six diffusion-weighted
images acquired with gradients applied in non-collinear and non-coplanar directions to solve the matrix D.
The tensor model needs low computational cost and can quantify diffusion anisotropy and estimate axon
fibers’ principal direction. However, a fundamental limitation of the tensor model is that it can only represent a single fiber orientation and fails at complicated fiber configurations such as fiber crossings, kissing,
or fanning. More advanced mathematical models are needed to handle the complicated fiber geometry.
11

2.3 Multiple Fiber Models
For better representing the complicated fiber configurations, High Angular Resolution Diffusion Imaging
(HARDI) techniques were developed to acquire more diffusion-weighted images in the q-space. q is referred to as the wavevector, and it depends on the length, strength, and orientation of the gradient pulse
sequence [41]. Advanced mathematical models are proposed using HARDI techniques to represent complicated fiber configurations, and these models attempt to define a probability distribution on the unit sphere,
as each point on the sphere corresponds to a unique orientation. The distribution function either quantifies the fraction of fiber portions or the probability of diffusing water molecules moving in a particular
direction. In the tensor model, the distribution function is parametric, similar to the Gaussian distribution.
The multi-tensor model [42, 43, 44] generalizes the tensor model and replaces the Gaussian model with
several Gaussian densities, each representing different fiber bundles. The ball-and-stick model [45, 46] is a
variation of the multi-sensor model; it differentiates the water molecules into restricted water molecules in
fibers and free molecules that do not interact with fibers. The free molecules are modeled using isotropic
Gaussian, and the restricted molecules in fibers are modeled as a Gaussian model with only one non-zero
eigenvalue so that particles move only in the fiber direction.
Diffusion spectrum imaging (DSI) and QBall imaging (QBall) reconstruct the probability distribution
function called the diffusion orientation distribution function (dODF). DSI [47, 48, 49] attempts to measure the probability of diffusing water molecules moving in a particular direction directly. For a pulse
sequence with infinitesimally short gradient pulses, the probability distribution function p is the Fourier
transform of the measurement A(q), as a function of wavevector q, which is a vector in 3D space that
controlled by varying the gradient strength and orientation. QBall [50] estimates the probability function
as the Funk–Radon transform of the measurement function A(q) on a shell in q-space. QBall acquires the
measurements in q-space using a spherical acquisition instead of the grid in DSI, and the measurements
for QBall imaging are more affordable than DSI.
12

Spherical deconvolution (SD) [51, 7] aims to recover the fiber orientation distribution (FOD) function,
directly quantifying the fraction of fiber portions. The FOD model regards the signal as the sum of measurements from different orientations, with each orientation weighted by the fraction of fibers along that
orientation. Mathematically, the signal A(q) is the convolution between the measurement R(q; x) for a
single fiber population with orientation x and the FOD function f:
A(q) = Z
f(x)R(q, x)dx, (2.7)
where f is the FOD function that needed to be computed and R(q, x) is the response function of a single
fiber bundle. Regularized versions of SD [5, 52] were proposed to overcome the susceptibility of SD to
noise, and these methods also refer to constrained spherical deconvolution (CSD).
Among these mathematical models, the fiber orientation distribution (FOD) model does not require the
pre-specification of the number of underlying fibers, is more accessible to acquire the required data, and
is computationally more efficient [53]. The FOD model directly quantifies the fraction of fiber portions
within a voxel with each orientation. Moreover, the distribution function of FOD provides the sharpest
peaks, which helps the tractography to track the fiber connections more accurately. FOD has become
more popular in the last decade and facilitates the development of FOD-based fiber-tracking techniques
for brain connectivity research. This thesis will focus on the FOD function and the FOD-based techniques
in tractography.
2.4 Spherical Harmonics
FOD function, defined on the unit sphere, is commonly represented and computed as coefficients of the
spherical harmonics up to the order N:
13

f(θ, ϕ) = X
n,m
u
m
n Y
m
n
(θ, ϕ), (2.8)
where θ and ϕ are the polar and azimuth angles of the spherical coordinates in the space of the MRI
scan; Y
m
n
is the mth (−n ≤ m ≤ n) real spherical harmonic at the order n (0 ≤ n ≤ N) and u
m
n
is
the coefficient for the corresponding spherical harmonics. Since the FOD functions are symmetric on the
unit sphere, only even ordered spherical harmonics are used to represent the FODs. The real spherical
harmonics Y
m
n
is defined as:
Y
m
n
(θ, ϕ) =



a
m
n P
|m|
n (cos θ) sin |m|ϕ, if m < 0
b
m
n P
m
n
(cos θ) cos mϕ, if m ≥ 0
(2.9)
where a
m
n
and n
m
n
are normalization constants; P
m
n
is the associated Legendre polynomial at degree n
and order m. The spherical harmonics are the eigenfunctions of the Laplace-Beltrami operator on the unit
sphere, and they form an orthonormal basis for the L
2
space. The spherical harmonics perform similarly
to the Fourier Series on the unit circle, and any function belonging to the unit sphere’s L
2
space can
be approximated by the finite linear combination of the spherical harmonics. Besides the real spherical
harmonics in the Eq. 2.9, the complex spherical harmonics also play a critical role in our work; the mth
complex spherical harmonic at the order n is defined as:
Wm
n
(θ, ϕ) = c
m
n P
m
n
(cos θ)e
imϕ
, (2.10)
where c
m
n
are normalization constants. The relationship between real and complex spherical harmonics is
similar to that between trigonometric and complex exponential functions. The basis transfer from real to
complex of the spherical harmonic is:
14

Wm
n =



(Y
−m
n − iY m
n
)/
√
2, if m < 0
Y
0
n
, if m = 0
(Y
m
n + iY −m
n
)(−1)m/
√
2, if m > 0
(2.11)
This transfer is inevitable, and its inverse is:
Y
m
n =



i(Wm
n − (−1)mW−m
n
)/
√
2, if m < 0
W0
n
, if m = 0
(W−m
n + (−1)mWm
n
)/
√
2, if m > 0
(2.12)
The complex spherical harmonics are efficient in rotating the FOD function, which is the key to efficient
computations in our surface-based tractography in Chap. 3 and interpolations in Chap. 4 and Chap. 5.
We introduce the Funk-Hecke Formula [54], which is commonly used for fast spherical deconvolution
in Eq. 2.7. The spherical convolution between a real spherical harmonic Y
m
n
and a convolution kernel R
can be computed by:
Z
S
2
Y
m
n
(x)R(q · x)dx = Y
m
n
(q)Gn, (2.13)
where R needs to be a function of the inner product q · x and Gn is defined as
Gn = 2π
Z 1
−1
Pn(t)R(t)dt, (2.14)
where Pn is the Legendre polynomial of degree n.
15

2.5 FOD with Compartments Models
The FOD function can more precisely represent the white matter fibers with compartment models, which
characterize diffusion imaging signals into different compartments. Our work uses the FOD reconstruction method [6] with compartment models to compute the FOD function from the diffusion MRI data.
This model has three compartments representing different types of water diffusion: intra-axonal diffusion,
extra-axonal diffusion, and trapped water with negligible diffusion. Mathematically, the diffusion signal
A(b, q) can be expressed as:
A(b, q) = Z
f(x)R(b, q, x)dx + αe−bλiso + γ + n, (2.15)
where the first term corresponds to the intra-axonal compartment, f is the FOD function that needed to be
computed, and R(b, q, x) is the response function. The second term corresponds to the extra-axonal compartment, b is the b-value, λiso is the diffusion coefficient as an unknown parameter, and α is the fraction
of the extra-axonal compartment and is also an unknown parameter. An unknown constant γ represents
the third compartment, corresponding to water molecules trapped inside cell bodies surrounding axons
with negligible diffusion. Finally, n denotes the noise in the signal. The response function R(b, q, x) is
chosen as:
R(b, q, x) = e
−bλa(q·x)
2
, (2.16)
where λa is the diffusion coefficient along the axons and is a pre-selected parameter.
Eq. 2.15 with constraints and regularizations forms an inverse problem, and the FOD function f and
unknown parameters are solved by the iterative application of a quadratic programming step and a gradient
descent step. This FOD reconstruction framework can accurately estimate fiber directions with clean and
16

sharp FOD functions and provide superior performance in tractography. Our work utilized this state-ofthe-art FOD method to develop advanced techniques for improving connectome research.
2.6 Tractography
Tractography is the only available tool for non-invasively studying white matter pathways in the human
brain. These pathways, which consist of axons, transfer information between brain regions and are critical
to our understanding of the function and structure of the brain. The tractography techniques are based on
the fundamental assumption: when many axons align along a common direction, the diffusion of water
molecules will be hindered to a greater extent across this direction than along it [23]. This anisotropy of
diffusion can be measured by diffusion MRI and further represented by advanced fiber orientation models.
Tractography is the technique to characterize and infer the white matter pathways based on diffusion MRI
data and fiber orientation models.
Tractography is the algorithm that integrates voxelwise fiber orientation information into a pathway
that connects brain regions. The pathway usually consists of numerous streamlines; each starts at a seed
point and follows the local fiber orientation information step-by-step, effectively forming a tract connecting brain regions. Depending on whether the algorithm considers the uncertainty in the fiber orientation
model, the tractography techniques can be either deterministic [55, 56] or probabilistic [57, 58, 59]. The
deterministic methods, assuming no uncertainty along the trajectory, integrate along the direction with
the strongest diffusion, for example, the eigenvector direction corresponding to the largest eigenvalue in
the tensor model or the peak value direction in the FOD model. Deterministic algorithms are accessible
and often need low computational costs. However, deterministic tractography is prone to errors caused by
noise and modeling errors and potentially underestimates the fiber connections [60]. In contrast, probabilistic tractography utilizes the uncertainty in the fiber orientation models; for example, the FOD function
quantifies the probability of fiber portions along different directions. Probabilistic tractography algorithms
17

usually sample a direction according to the distribution of the fiber orientation model. Thus, probabilistic tractography is a random walk guided by the local fiber orientation models but with some anatomical
constraints. Probabilistic tractography algorithms often provide more complete fiber connections in the
brain; meanwhile, they generate lots of false positive connections due to their probabilistic nature.
2.7 Opitmal Transport and Wasserstein Space
Optimal transport is an important mathematical problem that studies how to move one mass distribution
to another most efficiently. Let X and Y be two complete and separable metric spaces; P(X) and P(Y )
represent the Borel probability measures on X and Y . T : X → Y is a Borel measurable map, and
µ ∈ P(X) is a measure, the push forward of µ through T is defined by
T#µ(E) = µ(T
−1
(E)), ∀ Borel set E ⊂ Y. (2.17)
For a fixed Borel cost function: c : X → Y . The Monge version [61] of the optimal transport problem is
Problem 1 (Monge optimal transport problem) Let µ ∈ P(X), v ∈ P(Y ). Minimize
T →
Z
X
c(x, T(x))dµ(x). (2.18)
among all transport maps T from µ to v, i.e. all maps T such that T#µ = v.
However, this problem does not always have a solution; for example, if µ is a Dirac delta and v is not,
no admissible T exists. Thus, we often refer to a relaxation problem
Problem 2 (Kantorovich’s formulation of optimal transportation) Minimize
γ →
Z
X×Y
c(x, y)dγ(x, y), (2.19)
18

γ ∈ Γ(µ, v) and Γ(µ, v) is the set of all the couplings of µ to v, i.e., the set of Borel Probability measures on
X × Y such that
γ(A × Y ) = µ(A) ∀ Borel set A ⊂ X, γ(X × B) = v(B) ∀ Borel set B ⊂ Y. (2.20)
The transport plan γ ∈ Γ(µ, v) means that the value of γ(A × B) is the amount of mass initially
in A, which is sent into the set B. Kantorovich’s optimal transportation is finding a transport plan to
minimize the cost. Based on Theorem 1.5 in [61], this problem always has a solution if the cost function
c is continuous. The Kantorovich problem is a relaxation of the Monge problem [62]; namely, for each
transport map T, one can associate a transference plan so that the integral in Eq. 2.18 equals the integral
in 2.19. If c is continuous and µ is nonatomic, it is even possible to show that [61]
inf(Monge) = min(Kantorovich). (2.21)
More importantly, both the objective function Eq. 2.19 and the constraints Eq. 2.20 are linear in γ, so the
Kantorovich problem is actually infinite-dimensional linear programming, and in numerical computation,
the problem becomes finite-dimensional linear programming. In this thesis, we focus on the Kantorovich
problem.
If Y = X in Kantorovich problem, and we choose the cost function as c = d
p
, where d is the metric
of space X, Wasserstein p-distance is the p
th root of min (Kantorovich); namely, Wasserstein p-distance
between Borel probability measures µ and v is
Wp(µ, v) =

 inf
γ∈Γ(µ,v)
Z
X×X
d(x1, x2)
p
dγ(x1, x2)


1
p
. (2.22)
19

Wp is a metric defined on the p-Wasserstein space that is a subspace of the Boreal measure space P(X)
and is defined as
Wp(X) =
µ ∈ P(X) : Z
X
d
p
(x, x0)dµ(x) < ∞ for some, and thus any , x0 ∈ X

. (2.23)
If the space X is bounded, namely d(x1, x2) < M ∀x1, x2 ∈ X with some constant M < ∞,
then P(X) equals Wp(X). The most widely used p-Wasserstein space is W2(X) since it inherits many
important geometric properties of the base space (X, d) and provides excellent properties in statistical
analysis [61, 62].
20

Chapter 3
Surface-based Probabilistic Fiber Tracking in Superficial White Matter
3.1 Introduction
Superficial white matter (SWM) lies directly beneath the cortex and contains the short association fibers,
or U-fibers, connecting neighboring gyri [63, 64, 12]. Over the last decade, significant advances have
been made in connectome imaging techniques with dramatically increased spatial and angular resolution
in diffusion MRI (dMRI) [3]. Current connectome research, however, focuses primarily on long-range
fiber connections in the deep white matter (DWM) with little attention paid to the U-fibers in the SWM
even though the SWM has approximately twice as many fiber connections as the DWM [11] and plays
an essential role in various brain disorders [65, 66, 67, 68]. In this chapter, we will develop a novel fiber
tracking framework based on the surface representation of the SWM and provide dedicated tractography
tools for modeling U-fiber connectivity.
Tractography based on dMRI is the most widely used method for reconstructing fiber connections
in the human brain in vivo. Popular tractography methods typically realize tracking via sampling and
interpolating a fiber trajectory through the 3D image grid of the dMRI based on various mathematical
representations about fiber distributions within voxels [6, 5]. Previous works on fiber tracking in SWM
primarily rely on existing tractography methods. For example, cortical ROIs were used to select fiber tracts
for specific bundles of U-fibers [14, 69, 16, 15, 17]. In addition, whole-brain tractography was applied to
21

cluster U-fibers throughout the cortex [21, 20, 19, 18] and to construct atlases [10]. One known limitation
of existing volume-based tractography methods is that they tend to have many false negatives in the
reconstruction of highly curved fibers, such as the U-fibers [9]. Previous tractography methods can thus
typically only provide a partial representation of the U-fiber connections in the SWM.
There has been increasing interest in advancing the tracking of U-fibers by leveraging the surface representation of the SWM. A surface-based seeding approach was proposed recently to enhance the reconstruction of U-fibers [22], but it still relies on conventional volume-based tractography for fiber tracking
and cannot prevent the reconstructed tracts from deviating away from the sheet-like structure of the SWM.
To ensure the U-fibers follow the SWM beneath the cortex, we proposed a surface-based tractography
method [70] that tracks the U-fibers on the triangular mesh representation of the SWM. This deterministic tracking approach follows the peaks of the interpolated fiber orientation distribution (FOD) [6] on the
SWM mesh to generate the U-fibers. This method outperforms volume-based tractography in reconstructing U-fibers that connect the precentral and postcentral gyrus. While this method performs tracking on
a triangular mesh, it has not fully utilized the intrinsic surface geometry of the SWM to avoid limitations
in existing volume-based tracking. For example, it relies on angular differences from R
3
in quantifying
the changes of directions during the tracking propagation process. Hence, it does not follow the intrinsic
surface geometry to regularize tract smoothness. In addition, the deterministic tracking approach does not
consider the uncertainty in fiber orientations and can affect the completeness of U-fibers’ reconstruction.
In this chapter, we will develop a novel surface-based probabilistic tracking framework for reconstructing U-fibers in SWM. Our method systematically tackles challenges in existing volume- and surface-based
methods for U-fiber tracking. First, we carefully generate the surface representation of SWM based on nonlinear registration and deformation of the white matter (WM) surface computed from T1-weighted MRI.
After that, we develop a novel method to project the 3D FODs onto the tangent space of the SWM surface
based on an efficient transformation of spherical harmonics (SPHARM) representation of FODs to the local
22

Figure 3.1: Accurate placement of the WM cortical surface in the dMRI space. (a) T1-weighted MRI; (b)
the extra-axonal diffusivity calculated from dMRI by our compartment models. After the preprocessing by
HCP-Pipeline, yellow arrows in (c) highlight the residual misalignment of the WM cortical surface (cyan)
and the intra-axonal tissue map calculated by our compartment models. (d) shows the much-improved
alignment of the WM cortex (green) with the same tissue map after nonlinear registration between T1-
weighted MRI in (a) and the extra-axonal diffusivity map in (b).
coordinate system of the tangent space. A probabilistic tracking method is then developed based on the
projected 2D FODs on the SWM surface, where parallel transport is adopted to propagate the streamlines
across triangles and regularize tract smoothness intrinsically. By tracking based on the intrinsic geometry of the SWM surface, our method avoids the sharp turning angles of SWM in the 3D space and can
significantly improve the reconstruction of U-fibers in the SWM. Our experiments demonstrate that our
proposed method performs better than previous tracking algorithms on both in vivo and post-mortem human connectome imaging data. Furthermore, we applied our method to study SWM connectivity changes
in autosomal dominant Alzheimer’s disease (ADAD) patients and demonstrated statistically significant associations between parietal U-fiber connectivity and disease severity. This chapter has a published version
[35] in IEEE Transactions on Medical Imaging.
23

3.2 Methods
3.2.1 Mesh Representation of SWM
The proposed novel method for tracking the U-fibers is based on the surface representation of the SWM,
so the accurate placement of the cortical surfaces in the diffusion MRI (dMRI) space is a prerequisite. As
a first step, the triangular mesh representation of the white matter (WM) and gray matter (GM) cortical
surfaces are reconstructed from the T1-weighted MRI by FreeSurfer [71, 72]. For the dMRI data, we apply
the HCP-Pipeline [73] for distortion correction and its registration with the T1-weighted MRI. Nonlinear
distortions, however, still exist after these processes and affect the accurate alignment of cortical surfaces
with dMRI (Fig. 3.1. (c)). To remove the residual misalignment, we compute the compartment models
from the multi-shell dMRI by our previous work [6], which generates a diffusivity map (Fig. 3.1. (b)) of
the extra-axonal compartment for nonlinear registration with the T1-weighted MRI (Fig. 3.1. (a)) by the
ANTS software [74]. After that, we apply the nonlinear deformation calculated by ANTS to transform the
WM surface from the T1-weighted MRI to the space of the dMRI, which aligns more accurately with the
tissue boundaries, as shown in Fig. 3.1. (d). Finally, the transformed WM mesh is deformed inward toward
the deep white matter (DWM) along the normal direction by a distance of δ (e.g., 0.5mm) to generate the
SWM mesh for surface-based fiber tracking.
3.2.2 FOD Projection onto SWM Surface
The fiber orientation distribution (FOD) represents complicated fiber configurations as a function on a
unit sphere and has been used widely for fiber tracking in the 3D space. Since we track the U-fibers on the
surface rather than R
3
, the whole tracking procedure must be presented in the surface’s local coordinates
and tangent spaces. Thus, we represent the 3D FOD, which we denote as F OD3D, in the surface’s local
coordinates and project the F OD3D onto the tangent space of the SWM surface. We will first construct
24

Figure 3.2: The projection of 3D FODs onto 2D triangles local coordinate systems. (a) The F OD3D function is represented in a local spherical coordinate system defined on the center of a triangle of the SWM
mesh. The reference plane of the system is the triangle plane; meanwhile, the zenith direction z is aligned
with the normal direction. The x−axis is aligned with an edge of the triangle. (b) The projected F OD2D
using the Eq. 3.1. The F OD2D function is parametrized by ϕ in the local coordinate system and can be
normalized to a probability distribution function on the unit circle S
1
.
a local spherical coordinate system at each triangle (tangent space) of the SWM mesh, as shown in Fig.
3.2. (a), where the normal of the triangle is the zenith direction, an edge of the triangle is chosen as the
x − axis whose azimuthal angle ϕ equals 0, and the polar angle θ is the angle with respect to the zenith
direction. Fig. 3.2. (a) shows a F OD3D function at the center of a triangle on the SWM mesh, and it is
projected onto the triangle plane as a F OD2D function in Fig. 3.2. (b). In this local coordinate system,
we parameterize the 3D FOD as F OD3D(θ, ϕ) and project it onto the tangent space as:
F OD2D(ϕ) = Z π
0
F OD3D(θ, ϕ) sin(θ) dθ, (3.1)
The projected FOD2D function is parametrized by the azimuth angle ϕ in the local coordinate system and
defined on a unit circle S
1
in the tangent space of the SWM surface, as shown in Fig. 3.2. (b). The weighted
term sin(θ) in integral reduces the diffusion information perpendicular to the SWM surface. The integral
of the F OD2D function on the circle equals the integral of the F OD3D function on the sphere.
25

Figure 3.3: The local coordinate system Oxyz differs from the physical coordinate system Ox0y
0
z
0 by a
rotation R. (θ0, ϕ0) and (θ, ϕ)represent the same point p in the two systems. The rotation R is decomposed
as three successive rotations around three axes: (a) Ox0y
0
z
0
is rotated by an angle α around Oz0 →
Ox1y
1
z
0
; (b) Ox1y
1
z
0
is rotated by a second angle β around Oy1 → Ox2y
1
z
2
; (c) Ox2y
1
z
2
is rotated by
a third angle γ around Oz2 → Ox3y
3
z
2
, which is Oxyz.
On the other hand, F OD3D at each voxel is commonly parameterized by the spherical coordinates in
the physical space of the MRI scans. As discussed in Sec. 2.4, the F OD3D function is often represented
by the coefficients of real spherical harmonics (SPHARM) up to the order N:
F OD3D(θ0, ϕ0) = X
n,m
u
m
n Y
m
n
(θ0, ϕ0), (3.2)
where θ0 and ϕ0 are the polar and azimuth angles of the spherical coordinates in the space of the MRI
scan; Y
m
n
is the mth (−n ≤ m ≤ n) real spherical harmonic at the order n (0 ≤ n ≤ N) and u
m
n
is the
coefficient for the corresponding spherical harmonics. More details of the SPHARM are introduced in Sec.
2.4.
To track the U-fibers on the surface, we represent and compute functions in Eq. 3.1 based on the
surface’s local coordinates. The critical challenge is to efficiently obtain the SPHARM representation of
F OD3D in the local spherical coordinate system, i.e., we need to solve for the SPHARM coefficients v
m
n
at each triangle:
26

F OD3D(θ, ϕ) = X
n,m
v
m
n Y
m
n
(θ, ϕ) = X
n,m
u
m
n Y
m
n
(θ0, ϕ0), (3.3)
where v
m
n
is the unknown coefficient for the real SPHARM basis in the local system, and u
m
n
is the given
coefficient in the standard physical coordinate system as in Eq. 3.2; (θ0, ϕ0) and (θ, ϕ) represent the same
point in the different coordinates as in Fig. 3.3. With the SPHARM representation of FODs in the local
spherical coordinate system, we can write Eq. 3.1 as:
F OD2D(ϕ) = X
n,m
v
m
n
Z π
0
Y
m
n
(θ, ϕ) sin(θ) dθ, (3.4)
where the integrals are independent of the data and can be precomputed at a high level of numerical
precision, and we only need to compute the coefficients v
m
n
based on the FOD data. Thus, we can obtain
a high-resolution and precise F OD2D function using Eq. 3.4 with minor computation costs rather than
directly expanding and numerically integral of F OD3D with an expensive computation.
The key to computing v
m
n
from u
m
n
is to represent Y
m
n
(θ0, ϕ0) by a linear combination of Y
m
n
(θ, ϕ);
namely, find the transformation of SPHARMs under a coordinate system rotation. Since the rotation transformation of the complex SPHARMs is much easier to drive, we first transfer the real SPHARMs to complex
SPHARMs, shown in Eq. 2.10. Then, we derive the rotation transformation based on complex SPHARMs,
and, finally, go back to real SPHARMs. The transform from real to complex SPHARM is shown in Eq. 2.11,
and the inverse transform is Eq. 2.12.
For a rotation R from the coordinate system Ox0y
0
z
0
(Fig. 3.3) of the MRI physical space to the local
coordinate system Oxyz of the current triangle, we follow Wigner’s work [75] to decompose it as three
successive rotations around three axes: (a) Ox0y
0
z
0
is rotated by an angle α around Oz0 → Ox1y
1
z
0
;
(b) Ox1y
1
z
0
is rotated by a second angle β around Oy1 → Ox2y
1
z
2
; (c) Ox2y
1
z
2
is rotated by a third
angle γ around Oz2 → Ox3y
3
z
2
, which is Oxyz. The angles (α, β, γ) are the Euler angles, where α is
27

computed by finding the angle between Oy and the cross product of Oz and Oz0
, β is the angle between
Oz and Oz0
, and γ is the angle between Ox and Ox2
in step (c) above. In the following, we represent the
successive operations in the order from right to left, and the rotation R operator is now:
R = ZγYβZα, (3.5)
where Zγ is the rotation around the z − axis in the current coordinate system by an angle γ as explained
above, and the other two operators are similarly defined.
Based on a group symmetry argument [75], Wigner has proven that Wigner D-matrices can represent
the transformation of the n
th order complex SPHARMs between two coordinate systems based on the
decomposition in Eq. 3.5:
Wn(θ0, ϕ0) = Wn(θ, ϕ)Dn
(α, β, γ), (3.6)
where Wn is a (2n + 1) complex vector that represents the n
th order complex SPHARMs; Dn
(α, β, γ) is
a (2n + 1) × (2n + 1) matrix:
Dn
m,k(α, β, γ) = e
imγd
n
m,k(β)e
ikα
, (3.7)
where the first and third terms correspond to the rotations Zα and Zγ in Eq. 3.5, and the middle term is
induced from the rotation Yβ. The rotation around the z −axis is trivial since it only changes the azimuth
angle ϕ in the complex SPHARMs. The term d
n
m,k(β), corresponding to a rotation around the y − axis,
is much more challenging. Wigner showed that the rotation Yβ could be decomposed into successive
rotations as:
Yβ = Z−π/2Y−π/2ZβYπ/2Zπ/2
, (3.8)
28

Combining Eq. 3.5 and 3.8, we have:
R = Zγ−π/2Y−π/2ZβYπ/2Zα+π/2
, (3.9)
where the decomposition is much easier to compute since the rotations that depend on angles (α, β, γ)
are around the z − axis; meanwhile, the difficult rotations around the y − axis are fixed, π/2 or −π/2.
The closed-form solution [76, 77] of d
n
m,k(β) is:
d
n
m,k(β) = N
k
m(sin β)
k−m(1 + cos β)
mP
(k−m,k+m)
n−k
(cos β), (3.10)
where P is the Jacobi polynomial and the parameter Nk
m is:
N
k
m =
1
2
k
s
(n − k)!(n + k)!
(n − m)!(n + m)!, 0 ≤ m ≤ k ≤ n, (3.11)
other elements of the (2n + 1) × (2n + 1) matrix d
n
m,k(β) can be induced by the symmetric properties:
d
n
k,m(β) = (−1)k−md
n
m,k(β), (3.12a)
d
n
−m,−k
(β) = (−1)k−md
n
m,k(β), (3.12b)
d
n
m,−k
(β) = (−1)n−md
n
m,k(π − β). (3.12c)
We can now decompose the matrix Dn
(α, β, γ) in Eq. 3.6 and 3.7 as successive matrices:
Dn
m,k(α, β, γ) = E
Z
γ−π/2Dn
1 E
Z
β Dn
2 E
Z
α+π/2
, (3.13)
where Dn
1
and Dn
2
represent the matrix d
n
in Eq. 3.10 at β = −π/2 and β = π/2; EZ
ω
is the diagonal matrix
induced from the rotation by an angle ω around the z − axis. Let vn and un represent the coefficients
29

of the n
th order real SPHARMs under the local coordinate system at the current triangle and standard
coordinate system of the MRI scan; we have the following formula based on Eq. 3.4 and 3.6:
vn = UDn
(α, β, γ)U
−1un, (3.14)
where U is the transformation from real to complex SPHARMs that is defined in Eq. 2.11 and U
−1
is the
inverse transformation defined in Eq. 2.12; vn is a real (2n + 1) vector whose mth element is v
m
n
in Eq.
3.3, and un is the vector of real SPHARM coefficients u
m
n
in Eq. 3.3.
By transforming the real SPHARM coefficients from the common physical space to the local coordinate
system on each triangle, we can represent the F OD3D function using the real SPHARMs in the local
system and compute its 2D projection according to Eq. 3.4. Each F OD2D function is defined on a triangle
plane and represents the diffusion propagation information within this triangle, which allows us to track
the U-fibers on the SWM mesh probabilistically.
3.2.3 Parallel Transport
Streamlining-based tractography methods typically control the angular difference of consecutive tracking
directions for smoothness regularization and reduce false positives; namely, they constrain the new sampling direction according to the track’s previous direction. However, the SWM surface’s geometry structure differs from the R
3
; a direct angle comparison is not meaningful since the track’s previous direction
and new sampling direction lie in different tangent spaces of the surface. To incorporate a regularization mechanism in our surface-based tracking, we develop an intrinsic approach to compare neighboring
tracking directions based on the parallel transport of vectors on surfaces.
An affine connection (or covariant derivative) introduces a geometric structure on a manifold to compare vectors in neighboring tangent spaces [78]. Since the SWM mesh is embedded in R
3
, the Levi-Civita
connection induced by the Riemannian metric can be naturally used. With the Levi-Civita connection,
30

Figure 3.4: For a vector up in the tangent space of point p on a Riemannian manifold, the Levi-Civita
connection parallel transports up along a geodesic γ(t) to a vector uq in the tangent space of point q,
which allows the comparison of vectors in the tangent spaces of nearby points.
a vector in a tangent space can be parallel transported (the covariant derivative along the path is 0) to
another tangent space as in Fig. 3.4; namely, the parallel transported vector uγ(t) along a curve satisfies
the following differential equations:
∇γ
′
(t)
uγ(t) = 0
uγ(0) = up
(3.15)
where ∇ is the Levi-Civita connection, γ(t) is a curve connecting point p and point q; uγ(t)
is the parallel
transported vector along the curve. This differential equation system always has a unique solution based
on the existence and uniqueness theorem of differential equations if the point q is close to the point p.
For a track that leaves a current triangle and enters an adjacent triangle, we will parallel transport the
tracking direction in the current triangle to the new triangle to allow the comparison of angular differences
during the tracking process for smoothness regularization. The parallel transport induced by the LeviCivita connection can be numerically realized for two adjacent triangles on the SWM mesh as shown in
Fig. 3.5; namely, the parallel transport is achieved by unfolding the pair of adjacent triangles, translating
the vector from one triangle to the other along the geodesic (straight line) in the unfolded plane, and finally
folding the triangles back to their original positions together with the transported vector [79]. Since we use
the projected F OD2D function on each triangle for the sampling of fiber tracking directions, we express
the parallel transport based on the local spherical coordinate system of each triangle. For a direction v
31

Figure 3.5: The parallel transport of a vector v in the current triangle T1 to a neighboring triangle T2.
The parallel transport consists of unfolding the triangles to a common plane (a) → (b), translating the red
tangent vector v from T1 to T2 (b), and folding the triangles back to their original positions (b) → (c). The
transported vector vp can then be compared with any sampling direction w in the triangle T2.
parametrized by the azimuth angle ϕ in the local coordinate system of a triangle T1, as in Fig. 3.5. (a), the
vector v is represented as:
v = R
z
(ϕ) ∗ ex, (3.16)
where ex and z are the directional vectors of the x − axis edge and normal in T1 as illustrated in Fig.
3.5, and the matrix R is the operator that rotates by angle ϕ around the normal direction z. The parallel
transport from a triangle T1 into T2 is achieved by rotating the vector v as follows:
vp = R
ec
(θc) ∗ v, (3.17)
where ec is the directional vector of the shared edge of T1 and T2; θc is the angle between the normal
vectors of the two triangles (the two triangles must be oriented consistently). The vector vp lies in the
plane of triangle T2 and can be compared with any newly sampled direction in triangle T2, as in Fig. 3.5.
(c), for the calculation of intrinsic angular differences during the fiber tracking process.
In the extreme case where a track leaves the current triangle at a vertex vi
instead of any edge, we use
the geodesic polar map [80] to flatten all triangles in the one-ring neighborhood of the vertex into a 2D
chart. The geodesic polar map of the vertex vi flattens its one-ring neighborhood into a plane by enforcing
32

Table 3.1: Algorithm 1. Surface-based Probabilistic Tractography
Inputs: FODs, WM mesh, seed ROI, stopping condition, deformation distance δ, angle
threshold θin, F OD2D cut-off threshold F ODmin.
Output: N U-fibers on the SWM mesh.
A. Deform the WM mesh by a distance δ toward DWM to obtain the SWM mesh.
B. Project the FODs onto the SWM mesh to generate the F OD2D at each triangle.
C. Repeat the following steps to track N U-fibers:
1. Start tracking a new U-fiber and set n = 0. Randomly select a seed triangle Tn from the
Seed ROI. Set pn as the center of the seed triangle.
2. Pick a tracking direction v at Tn with its F OD2D > F ODmin via rejection sampling.
If n = 0, save v0 = v. If n > 0, repeat rejection sampling until the angle between v and the
parallel transported direction vp from the previous triangle is less than θin. If the number
of rejection sampling > 50, abort the tracking process and return to step 1.
3. Compute the point p such that p = pn + sv, s > 0, p ∈ ∂Tn. Find the triangle T ∈ NTn
to continue the tracking and parallel transport the current tracking direction v to vp in T
using the method in Sec. 3.2.3.
4. Set pn+1 = p, Tn+1 = T, n = n + 1. Return to step 2 and continue the tracking until a
stopping condition is met.
5. Restart the tracking from the seed triangle in the opposite direction v = −v0 until a
stopping condition is satisfied.
6. Merge the tracking results from the two opposite directions as a complete U-fiber. If the
desired number of tracks is reached, stop the algorithm; otherwise, return to step 1.
the geodesics from vi
(i.e., lines starting at vi
) to remain straight in the plane while rescaling the angles
between geodesics by:
ri =
2π
P
Tjik
θjik
, (3.18)
where θjik is the angle at the vertex vi of all triangle Tjik in the one-ring neighborhood of vi
. In the
2D chart, we can easily find the upcoming triangle to continue the tracking by extending the tracking
line through the vertex. Finally, we map the extended tracking line from the 2D chart to the original
position of the upcoming triangle in the 3D mesh by the inverse geodesic polar map, which will be the
parallel transported vector to compare the angular difference between tracking directions of neighboring
triangles.
33

3.2.4 Tracking Algorithm
Using the triangular mesh representation of SWM and projected F OD2Ds, we will develop a surfacebased probabilistic fiber tracking algorithm of U-fibers as listed in Algorithm 1. Our algorithm starts
from a randomly selected seed triangle and propagates a streamline on the SWM mesh until a user-defined
stopping condition is met. Similar to previous FOD-based tractography algorithms, the stopping condition
can be: (a) the streamline arrives at the boundary of a predefined mask on the SWM mesh; (b) the streamline
meets the conditions of including or excluding ROIs from user-defined protocols; (c) no valid F OD2D data
to continue the tracking. In Algorithm 1, Tn denotes a triangle from the triangular mesh, pn denotes a
point of the streamline, NT denotes the neighboring triangles of a triangle T, and ∂T denotes the boundary
of the triangle T. The main steps of Algorithm 1 are explained as follows.
Sampling: For a point pn and the triangle Tn, the rejection sampling method is used to sample the
direction of the streamline based on the normalized F OD2D function of the triangle Tn. No smoothness
regularization is needed for the sampled direction when pn is the seed point. In other cases, we repeat
sampling the direction until the angle difference between the sampled direction v and the parallel transported vector vp is less than a threshold θin. Moreover, the F OD2D value of the sample direction v is
required to exceed a threshold F ODmin to reduce the propagations in areas with a low probability of valid
fiber connections. If there is no accepted direction until the sampling times reach a threshold (e.g., 50), we
abort the current fiber tracking process and restart it from a newly sampled seed triangle.
Propagation within a triangle For the point pn within the triangle Tn (including the boundary
points), we compute a line segment to propagate the streamline piecewise linearly, which originates from
pn and terminates at the triangle boundary following the sampled direction v. Then, the terminating point
p enters an adjacent triangle T and acts as the starting point of the next line segment in this new triangle.
Propagation across triangles: When our streamline leaves a current triangle and enters one adjacent
triangle, we parallel transport the direction v in the current triangle to vp that lies in the new triangle
34

using the method in Sec. 3.2.3. As discussed in the rejection sampling process, the vector vp is used for
smoothness regularization.
In Algorithm 1, we apply fiber tracking at a seed triangle in symmetric directions. After completing
tracking in both directions, we concatenate the two streamlines to form a complete U-fiber on the SWM
mesh. For practical U-fiber tracking experiments, the user needs to specify the seed ROIs and properly
define the stopping criteria on the cortical surface. We will give more detailed examples of setting up
surface-based tracking protocols in our experimental results.
3.3 Experiments and Results
We will apply our surface-based fiber tracking algorithm to three datasets to demonstrate its effectiveness
in reconstructing U-fibers in the SWM. We first apply the proposed method to the high-resolution dMRI
data from the Human Connectome Project (HCP) [2] and compare it with two other tractography methods.
After that, we applied our method to two post-mortem scans and validated the agreement between the
tracked U-fibers and manually delineated SWM from their high iron contrast on T2-weighted MRI. Finally,
we apply our method to reconstruct the U-fibers in the parietal lobes of 26 subjects with ADAD and show
that a significant association of U-fiber connectivity with disease severity can be successfully detected.
3.3.1 U-fiber Tracking on HCP Subjects
In our first experiment, the T1-weighted MRI and multi-shell dMRI data of 484 subjects from the 500 release of HCP were used. The age range of the 484 subjects is 22 to 36 years old, and there are 195 males and
289 females. T1-weighted MRI of HCP has an isotropic resolution of 0.7mm, and dMRI images have an
isotropic resolution of 1.25mm and were acquired from 270 gradient directions distributed over three different b-values: 1000, 2000, and 3000 s/mm2
. All data were first preprocessed by the HCP-Pipeline, which
reconstructs the cortical surfaces from the T1-weighted MRI and registers the dMRI to the T1-weighted
35

Figure 3.6: Representative HCP examples of U-fibers between the pre-and post-central gyrus reconstructed
by the two surface-based methods (SP: surface-based probabilistic tracking; SD: surface-based deterministic tracking) and the MRtrix software. For better visualization, we downsampled all results to 1000 streamlines. The subject ID is plotted at the top of the results from each HCP subject.
MRI after artifacts correction. After that, FODs and compartmental parameters were computed [6]. In this
experiment, we will reconstruct the U-fibers connecting the pre- and post-central gyrus [14, 16, 20] and
compare the performance with state-of-the-art volume-based fiber tracking from the MRtrix software [59]
and a deterministic surface-based fiber tracking method [70] that our group developed previously.
As the first step for surface-based tracking, we resampled the WM cortical surface from the left hemisphere of all subjects into a mesh with 20,000 triangles and deformed it toward the DWM along the normal
direction by a distance of δ = 0.5mm, which will be validated in our second experiment for the generation of the SWM mesh in surface-based tracking. Following the work in [70], we used the Hamilton-Jacobi
36

Figure 3.7: Quantitative comparisons of the U-fibers reconstructed by the proposed surface-based probabilistic (SP) fiber tracking, surface-based deterministic (SD) fiber tracking, and volume-based fiber tracking
in the MRtrix software on 484 HCP subjects. (a) Box plots of the ‘well-U-connected’ measure of the three
methods. Box plots of the ‘well-distributed’ measure on the precentral and postcentral gyrus are shown
in (b) and (c), respectively. (d) Box plots of the U-ratio measure across subjects. (e) Box plots of the Procrustes distance are used to assess the topographic regularity of the reconstructed U-fibers using the three
methods.
skeleton method [81] to generate the pre- and post-central gyrus’s gyral skeletons according to FreeSurfer
labels. The gyral skeletons from these two gyri will be used as the inclusion ROIs for surface-based tracking. The stopping condition for surface-based tracking is that both inclusion ROIs have been reached.
Furthermore, the seed ROI for surface-based tracking is obtained as triangles within a 2mm geodesic distance to the boundary between the pre-and post-central gyrus. For the MRtrix software, all neighboring
voxels of the triangles in the seed ROI on the SWM mesh are combined to form the seed ROI in volumebased fiber tracking. The gray and white matter masks of the pre-and post-central gyrus from the aseg
results of FreeSurfer were combined to form the mask for fiber tracking with MRtrix, which will stop the
tracking if a tract leaves the mask.
For the deterministic surface-based method in [70], we set the step size to be 0.01mm, the two angle
thresholds as T Hθin = 1◦
and T Hθxing = 10◦
, as used in the same experiment in [70]. We set the
angle threshold θin = 10◦
and the F OD2D cut-off threshold F ODmin = 0.01 for the proposed novel
probabilistic surface-based fiber tracking method. The iFOD1 algorithm is used for MRtrix, and the main
parameters of iFOD1 are set as step size = 0.1mm, angle threshold = 10◦
, and the F OD3D cut-off threshold
used for the MRtrix was 0.05. For all methods, 30,000 seed points are randomly sampled for fiber tracking.
37

Reconstructed U-fibers from five representative HCP subjects are shown in Fig. 3.6. For all subjects, Ufibers reconstructed by the volume-based tractography in MRtrix do not closely follow the cortical folding,
and most streamlines do not touch the gyral skeletons. On the other hand, both surface-based methods can
generate U-shaped streamlines that follow tightly to the cortical surface. Compared to the surface-based
deterministic method, the proposed probabilistic and surface-based fiber tracking method provides more
geometrically regular and topographically organized U-fibers connecting the two neighboring gyri.
To quantitatively evaluate the performance of the fiber tracking methods in the reconstruction of Ufibers in SWM, we use three similar measures proposed in [70]: ‘well-U-connected’, ‘well-distributed’,
‘well-U-shaped’, and an additional metric for topographic regularity. The ‘well-U-connected’ measure
counts the number of reconstructed streamlines connecting two neighboring gyri given a fixed number
of seed points. Considering the cortical thickness of the precentral and postcentral gyrus [82], we count
the number of streamlines with their endpoints within 4mm to the gyral skeletons of the precentral and
postcentral gyrus for each subject. The boxplot of this measure from all 484 HCP subjects is in Fig. 3.7.
(a). The proposed method generates the most valid connections among all three methods, with 48.5%
of seed points generating valid connections. For the surface-based deterministic method, 29.5% of tracts
initiated from the seed points develop into valid connections. For the volume-based tracking from MRtrix,
only 6% of the tracts from the seed points generate valid U-fiber connections even if we choose their seed
points in the sulcal regions to minimize the need for highly curved U-turns during the tracking process.
These results demonstrate the high efficiency of the proposed surface-based probabilistic tractography in
reconstructing the U-fibers. The ‘well-distributed’ measure evaluates how completely the U-fibers connect
the two gyri. Thus, we uniformly partitioned both the precentral and postcentral gyral skeletons into 20
sections and counted how many sections were hit by the streamlines. More specifically, one section of
a gyral skeleton is hit if the endpoint of a streamline is within 4mm of the section. Fig. 3.7. (b) and (c)
show the boxplots of the ‘well-distributed’ measure of all three methods for both the pre-and post-central
38

Figure 3.8: Across the HCP subjects, box plots of the number of touched sections in (a) pre- and (b) postcentral gyrus by MRtrix results with varying numbers of seeds from 3000 to 50,000 are shown. Box plots
of the number of touched sections by results from the proposed surface-based probabilistic (SP) method
(30,000 seeds) are also shown for ease of comparison.
gyrus, which indicates that the fiber tracts reconstructed by the proposed method provide a complete
representation of the U-fibers connecting the two gyri in the SWM. For the ‘well-U-shaped’ measure, we
compute the U-ratio of an individual streamline as the distance between two endpoints divided by the
streamline length. A smaller U-ratio is desirable for the successful reconstruction of valid U-shaped fibers.
For each subject, we compute the mean U-ratio for all the streamlines. Fig. 3.7. (d) shows the box plots
of the mean U-ratio from all 484 HCP subjects. Results from the proposed method have much smaller
U-ratios than the other two methods.
Topographic regularity is an important property widely present in various fiber pathways [34, 32, 33].
In particular, somatotopy is a well-known organizational principle in the motor and sensory system [83,
84]. Inspired by this anatomical property of brain pathways, several quantitative metrics [85, 86, 37] were
proposed to measure the topographic regularity of the reconstructed fiber tracts. In this experiment, we
use the metric proposed in [85] to measure the topographic regularity of the U-fibers. For each subject, we
apply the classical multidimensional scaling (MDS) to project both the start points (precentral gyrus) and
the endpoints (postcentral gyrus) of the streamlines to the 2-D R
2 plane, respectively. Then, we compute
the Procrustes distance between the projected start and end points to measure their shape difference in
39

the R
2 plane. For the reconstructed tracts from all three methods, the box plots of the Procrustes distance
of the U-fibers from all subjects are shown in Fig. 3.7. (e). The lowest distance of the U-fibers from
the proposed method suggests they are topographically the most organized. On the other hand, tracts
from volume-based tractography have much higher Procrustes distances. These results demonstrate more
topographically regular U-fibers can be generated by constraining the fiber tracking process along the
SWM surface.
Compared to our novel surface-based tracking methods, we can see from Fig.3.7 that volume-based
MRtrix tracking has lower performances in all measures. In these previous experiments, we fixed the
number of seeds for all methods to 30,000 for each subject. Among all the measures, the ‘well-U-connected’
(ratio of valid connections), ‘well-U-shaped’, and topographic regularity are relatively stable with respect to
tract numbers, but the ‘well-distributed’ measure in Fig. 3.7. (b) and (c) can vary with respect to the number
of valid tracts in MRtrix results. To further demonstrate that our proposed method provides more complete
U-connections even if MRtrix dramatically increases the number of seeds; we ran the same experiment
using MRtrix, varying the number of seeds from 3000 to 500,000. Boxplots of the number of sections in
the two gyri touched by the valid U-fibers in MRtrix results are shown for each seed number in Fig. 3.8.
For ease of comparison, the boxplot of results from our method (SP) obtained previously from only 30,000
seeds was also plotted in Fig. 3.8. For MRtrix results, we can see the number of touched sections in both
gyri increases as the number of seeds grows. However, the number of touched sections stabilizes once the
number of seeds arrives at around 200,000. Even if we generated tracts from 500,000 seeds using MRtrix,
the number of touched sections on both gyri remains significantly lower than our method. This result
shows that it is hard for traditional volume-based tractography to reconstruct complete U-connections
even with excessive computations, mainly because they will keep generating connections with repeated
patterns. On the contrary, the proposed surface-based tracking method can efficiently generate a more
complete representation of SWM connections.
40

Figure 3.9: T2-weighted 7T MRIs of two post-mortem brains are shown in (a) Subject 1 (91-year-old male)
and (b) Subject 2 (97-year-old female). SWM with higher iron concentrations is highlighted with blue
arrows. Manually delineated SWM masks in the lateral frontal lobe for subjects 1 and 2 are shown in (c)
and (d), respectively.
3.3.2 Validation with 7T Post-mortem MRI
In this experiment, we will apply our surface-based fiber tracking algorithms to high-resolution dMRI of
two post-mortem brains and compare the reconstructed fiber tracts with manually segmented SWM masks
on T2-weighted MRIs, where the higher iron concentration in SWM [87, 88] provides excellent contrast
for manual segmentation.
Subjects were recruited from the Death Verification Service of the Capital in São Paulo (SVOC), and
after obtaining informed consent from the family in situ MRI was performed right before autopsy at the
PISA facility (Imaging Platform at the Autopsy Room) located at the Medical School of the University of
São Paulo of Brazil. The cadavers were scanned in cranio on a 7T SIEMENS Magnetom scanner with a 32
41

Figure 3.10: Overlay of reconstructed U-fibers with the SWM on T2-weighted post-mortem MRIs. Top row:
subject 1; Bottom row: subject 2. (a) and (d) show a zoomed view of the overlay of the U-fibers with the
MRI in the red box highlighted in (b) and (e), respectively. (c) and (f) show a zoomed view of the overlay
of the U-fibers with the MRI in the blue box highlighted in (b) and (e), respectively.
receiver channel head coil (Nova Medical, USA) before brain procurement at room temperature. Following
the HCP imaging protocol, we acquired T1-weighted MRI (0.75mm isotropic), multi-shell dMRI (1.3mm
isotropic, 76 gradient directions over three b-values 1000 s/mm2
, 2000 s/mm2, and 3000 s/mm2), and
high-resolution T2-weighted MRI (0.45x0.45x2mm) that provides excellent contrast for SWM (Fig. 3.9.
(a) and (b)). In this experiment, we manually delineated the SWM with the ITK-SNAP software [89] on
the frontal lobe of two post-mortem subjects (Fig. 3.9. (c) and (d)) based on the T2-weighted MRI. These
masks will be used to evaluate the reconstructed U-fibers quantitatively. Subject 1 is a 91-year-old male,
and Subject 2 is a 97-year-old female.
The same preprocessing steps used in the first experiment were applied to the post-mortem MRIs to
generate the WM cortical surface from T1-weighted MRI and FODs from the dMRI. We tested different
choices of the deformation distance parameter δ from 0.3mm to 0.7mm to generate the SWM mesh and
evaluate its impact on the reconstructed U-fibers compared to the manually segmented SWM masks. For
42

Figure 3.11: (a) The mean distance from the U-fibers reconstructed by the proposed algorithm to the SWM
mask as a function of the deformation distance parameter δ. Boxplots of the distance of the reconstructed
U-fibers to the SWM masks in (b) subject 1 and (c) subject 2. MRtrix: iFOD1 algorithm in the MRtrix
software. SP: proposed surface-based probabilistic tracking method with δ = 0.5mm.
the proposed surface-based probabilistic fiber tracking on the SWM meshes, we use the same tractography
parameters as in the first experiment: θin = 10◦
and F ODmin = 0.01. On the lateral frontal lobe of each
subject, we applied the same Hamilton-Jacobi skeleton analysis methods as in the first experiment and
extracted the gyral skeletons as the stopping condition and sulcal basins as the seed ROI.
For each SWM mesh generated from a given deformation distance δ, we randomly sampled 105
seed
points to compute the U-fibers in the lateral frontal lobe with our algorithm. For the results with δ =
0.5mm, we overlay the reconstructed U-fibers of the proposed method with the high-contrast T2-weighted
MRI of the two post-mortem brains and plotted the results in Fig. 3.10, where we can see U-fibers generated
by our method agree very well with the SWM in these images. To quantitatively measure how accurately
the reconstructed U-fibers follow the SWM, we computed the distance between the U-fibers and the SWM
mask. For each point q from the U-fibers, its distance to the mask M is defined as follows:
d(q, M) = min
m∈M
||q − m||2, (3.19)
43

For accurate distance computation, we interpolated the SWM mask to an isotropic resolution of 0.1 mm.
For each SWM mesh from a given δ, the mean distance of all the points on the reconstructed U-fibers is
calculated to assess the impact of the deformation distance parameter. As shown in Fig. 3.11. (a), the best
deformation distance is around 0.5mm for both subjects. For this parameter choice (δ = 0.5mm), boxplots
of the distance from the U-fibers to the SWM mask on the lateral frontal lobe are shown in Fig. 3.11. (b) and
(c) for the two subjects. We also compared the volume-based tractography algorithm in MRtrix with the
same setup as in Sec. 3.3.1 for seed region and mask generation and parameter choices. Like our surfacebased tracking, 105
seed points were randomly sampled for the volume-based tractography method. In
addition, we removed invalid streamlines from volume-based tractography that do not enter the gyral
cortex. Boxplots of the distance from the fiber tracts generated by MRtrix to the SWM masks are plotted
in Fig. 3.11. (b) and (c). For both subjects, the reconstructed U-fibers from the proposed surface-based
tracking method are distributed much closer to the SWM mask than the fiber tracts generated by volumebased tractography.
3.3.3 U-Fiber Connectivity Changes in ADAD
In the third experiment, we will apply our surface-based fiber tracking method to a cohort of 26 ADAD
patients. There are 13 males and 13 females in this cohort, and the age range is 21 to 47 years old. Unlike
late-onset AD (LOAD), these ADAD patients all carry the A431E mutation in the PSEN1 gene and are
predetermined to develop dementia at a similar age range [90]. In addition, these patients tend to present
atypically high tau pathology in parietal cortices [91] instead of medial temporal cortices as in typical
LOAD. It is thus interesting to examine how the SWM connectivity changes in the parietal cortices of
these ADAD patients.
44

Figure 3.12: The reconstructed U-fibers in the parietal lobe of two ADAD patients are superimposed over
their SWM meshes. The left and right hemispheres of subject 1 are shown in (a) and (b), respectively. The
left and right hemispheres of subject 2 are shown in (c) and (d), respectively.
Following the HCP protocol, T1-weighted MRI and multi-shell dMRI were acquired from the 26 ADAD
patients on a 3T Siemens Prisma scanner at the University of Southern California (USC). All subject recruitment and data acquisition were approved by the Institutional Review Board (IRB) at USC. For each
subject, we use the Sum of Boxes (SOB) score of the Clinical Dementia Rating Scale (CDR) assessment [92]
as the severity of dementia of Alzheimer’s type.
The same preprocessing steps, as in Sec. 3.3.1, were first conducted to generate the SWM mesh and
FODs for surface-based fiber tracking. The same parameters were used: δ = 0.5mm, θin = 10◦
, and
F ODmin = 0.01 for surface-based tracking with our novel algorithm in the parietal lobes of each subject, which include the postcentral, superior parietal, inferior parietal, and supramarginal gyrus based on
45

FreeSurfer labels. For all subjects, 30000 seed points were sampled from the sulcal regions for surface-based
fiber tracking. In Fig. 3.12, we show the reconstructed U-fibers for two representative subjects: subject 1
is a 43-year-old female with a CDR-SOB score of 5.5, and subject 2 is a 37-year-old male with a CDR-SOB
score of 18.
Next, we quantitatively evaluated the association between the connectivity of the reconstructed Ufibers and disease severity. Motivated by the apparent fiber density based on FOD magnitude in 3D [93],
we computed the peak value of the F OD2D function along the U-fibers. We obtained the average peak
values of F OD2D over the parietal lobe for each hemisphere. After that, we first fitted a linear regression
model between the mean peak value and the CDR-SOB score of the 26 subjects in each hemisphere. The
results are shown in Fig. 3.13. (a) and (b), where the p-values of the linear coefficient for the left and right
hemispheres are 1.34e-3 and 1.62e-2, respectively. For 25 of the 26 ADAD patients, tau PET scans based on
the 18F-AV-1451 tracer were acquired [94]. Following the processing steps of PET-Surfer in FreeSurfer [95],
standard uptake value ratios (SUVRs) of tau PET signals were calculated in cortical areas as their severity
measure of tau pathology. Then, a linear regression model was fitted between the mean F OD2D peak
value along the U-fibers and the tau SUVR of the parietal lobe in each hemisphere. The results are shown
in Fig. 3.13. (c) and (d), where the p-values of the linear coefficient for the left and right hemispheres are
7.67e-5 and 2.25e-4, respectively. These statistically significant results suggest that the U-fiber connectivity
decreases with the increase of clinically assessed disease severity (CDR-SOB) and tau pathology in cortical
areas with known AD pathology for ADAD.
46

Figure 3.13: The linear regression results between the CDR-SOB score, tau SUVR, and the mean F OD2D
peak value on the U-fibers of the parietal cortices of the ADAD patients. (a), (c) are results from the left
hemisphere and (b), (d) are results from the right hemisphere.
3.4 Discussion
3.4.1 Parameters
We deform the WM inward by a distance δ, chosen as 0.5mm, to generate the SWM mesh across all
experiments. In experiment 3.3.2, we found that the SWM lies around 0.5mm below the WM mesh based
on manually delineated SWM mask in post-mortem data. Our measurement is consistent with previous
SWM imaging research [88], which confirms the SWM lies around 0.5mm below the WM surface based on
a high iron contrast and ultrahigh resolution MR imaging technique. For future work, it will be interesting
to vary the deformation distance for different cortical regions and examine whether this would enhance
the performance in U-fiber reconstruction and modeling SWM connectivity changes in brain disorders.
To obtain the F OD2D on the SWM mesh, we interpolate the F OD3D at each triangle from a neighborhood of voxels. At the current dMRI resolution of HCP data, the width of the neighborhood is greater
than the typical SWM thickness, which provides a sufficient window to incorporate related information of
surrounding SWM at each surface point. For future work, if we have dMRI data with much higher spatial
resolution than typical SWM thickness, it will be essential to adapt the interpolation techniques to include
more neighboring voxels to ensure the interpolated F OD2Ds have sufficient information related to SWM.
47

For fiber tracking on the SWM mesh, the angular threshold θin is chosen to regularize the smoothness
of the fiber tracts. In all experiments, we set θin to 10◦
to balance the streamlines’ regularity and the Ufibers’ completeness. For example, a lower θin for the experiment in Sec. 3.3 will generate smoother tracts,
but the overall U-fiber bundle will be less ‘well-connected’ between the pre-and post-central gyrus for a
given number of seeds. This phenomenon is unsurprising since a strict regularization in fiber tracking
often reduces the number of successfully tracked streamlines in tractography.
Another parameter related to fiber tracking is the minimum accepted F OD2D value F ODmin in the
rejection sampling procedure, which we set as 0.01 as the default value for all experiments in this work.
This parameter is similar to the FOD cut-off thresholds in MRtrix and other FOD-based fiber tracking
methods. However, as we project the F OD3D onto the tangent space of the SWM surface using Eq.
3.1, the F OD2D value is not directly comparable to the F OD3D value. Indeed, the F OD3D function
distributed sparsely (concentrated within small areas) on the unit sphere, an integral in Eq. 3.1, which
performs similarly to the weighted average, would reduce its amplitude. The default F ODmin is around
1/10 of the mean F OD2D peak value in HCP data and is similar to the often-used cut-off value in volumebased tractography. Increasing F ODmin will decrease the number of successful streamlines but improve
regularity and consistency.
In experiments, we resampled each WM cortical surface into a mesh with 20,000 triangles. After the
resampling, the triangles’ dimensions are comparable to the dMRI resolution, and the areas of the triangles
are more uniform. Moreover, after resampling, the triangles’ areas between different subjects are similar,
and we can use fixed tracking parameters for all subjects.
3.4.2 More General SWM Connectivity
Following the most widely used U-fiber reconstruction protocols in previous research [16, 21, 20, 19, 18,
10], we adopted the gyral crown as the termination point for surface-based fiber tracking in this work
48

so that U-fibers will connect adjacent gyri. However, attributing to the surface-based representation and
parallel-transport-based regularization, our method can choose ROIs and termination conditions on any
cortical regions if anatomically needed. A post-mortem dissection study [12] reported intra-gyral SWM
fibers connecting different parts of the same gyri. Moreover, the intra-gyral and inter-gyri U-fibers can
converge into junctional areas and form pyramid-shaped structures. Another histological and dMRI-based
macaque brain research [96] found longitudinal fibers propagating along the gyral crown and forming a
sheet structure. For future work, we will apply our method to explore these complex SWM fiber structures
with more general anatomical protocols informed by histology findings.
3.4.3 Code and Computation Speed
We implemented the proposed method in MATLAB and C++ and distributed it publicly to the research
community (https://github.com/Xinyu-Nie/SurfTracker). The C++ code can be compiled on Windows,
Mac, and Linux. Our code automatically uses multiple threads, and the C++ version is efficient for fast
computation. For example, our method needs less than two minutes to generate about 30, 000 streamlines
in experiment 3.3.1 on a Laptop with a 2.61GHz eight-core Intel CPU. The computational cost is primarily from the fiber tracking processes. The computation of F OD2Ds is speedy (within several seconds)
because of our closed-form computation of F OD2Ds using Eq. 3.4 and Eq. 3.14.
3.4.4 Clinical Applications
We applied the proposed method to an ADAD dataset in experiment 3.3.3 and found that the F OD2D peak
amplitudes were significantly associated with the CDR-SOB clinical scores and tau pathology deposition.
The reduction of the F OD2D amplitudes implies possible degeneration of fiber connectivity in the SWM
of the parietal lobe of ADAD patients with the progression of disease severity. In future work, we will apply
our method to more brain imaging datasets and examine the change of SWM connectivity due to aging
49

and various pathological changes, including amyloid-beta deposition and hippocampal atrophy. Moreover,
we will also apply our method for constructing U-fiber bundle atlases across cortical areas using HCP data
and distributing related tools in the same software.
50

Chapter 4
Flow-based Geometric Interpolation of Fiber Orientation Distribution
Functions
4.1 Introduction
Diffusion MRI (dMRI) is the most widely used technique for studying human brain structural connectivity in vivo [97]. Significant improvements in imaging techniques dramatically increased the spatial and
angular resolution of dMRI [3] and provided opportunities for advanced models such as fiber orientation
distribution (FOD) [5], which facilitates the development of FOD-based fiber tracking for brain connectivity research. However, well-known challenges in current tractography methods generate large amounts
of false positives and negatives [98]. While there have been considerable efforts in developing novel fiber
tracking methods [99], a critical step in tractography, FOD interpolation, has received rare attention.
In popular FOD-based tractography, linear interpolation is commonly adopted for numerical efficiency.
Still, it often generates artificial directions and ignores rotations between neighboring FODs, as shown in
Fig. 4.1. (B), which can lead to false positive streamlines. To enhance FOD interpolation, a Riemannian
framework was proposed in [13]; under the square root reparameterization, the space of FOD functions
can form the positive orthant of the unit Hilbert sphere. However, this framework is computationally
expensive and sometimes fails to provide anatomically meaningful interpolations [24]. A rotation group
51

Figure 4.1: An illustrative example of FOD interpolations. (A) FODs of four neighboring voxels from a
bending fiber bundle were highlighted by red circles in (B) and (C). The blue arrow shows the direction
of the fiber bundle. (B) Interpolated FODs are generated by linear interpolation, where artificial peaks are
generated. (C) The interpolated result by our proposed method correctly accounts for rotation and follows
the bending geometry of the fiber bundle.
action-based framework [24] was proposed that simultaneously averages the shape and rotation of FODs.
A later work [25] proposed a rotation-induced Riemannian metric for FODs and introduced a weighted
mean for FOD interpolation. However, since only one rotation is used for the whole FOD, these methods
cannot handle more general situations where individual FOD peaks experience different rotations. More
importantly, these methods have not been adopted in a tractography framework to advance fiber tracking
performance due to their numerical complexity.
In this chapter, we develop a novel framework to perform geometrically consistent interpolation of
FODs and demonstrate its effectiveness in enhancing the performance of fiber tracking. We decompose
each FOD function with multiple peak lobes into components, each with only one peak lobe. Then, we
locally model neighboring voxels’ single-peak components, consistent in direction, as a vector field flow
fitted by polynomials. Each vector field locally represents the geometry of an underlying fiber bundle and
continuously determines the direction of single-peak components within the support. Then, a closed-form
solution is developed to account for rotations of FODs represented as spherical harmonics and realize the
geometrically consistent interpolation of each FOD component, as shown in Fig. 4.1. (C). The interpolation
52

Figure 4.2: A FOD function with two peak lobes is decomposed into two SPHARM-based FOD functions,
each with only a salient peak lobe.
of a complete FOD function with multiple peak lobes is obtained by merging the single-peak interpolations from all the covering vector fields. In our experiments, we use HCP data to quantify the accuracy
of the proposed FOD interpolation algorithm and show that it achieves superior performance than the
commonly used linear interpolation approach. Furthermore, we apply our interpolation method to perform upsampling of FOD fields and significantly improve the performance of FOD-based tractography
both qualitatively and quantitatively. This chapter has a published version [36].
4.2 Method
4.2.1 FOD Decomposition
The fiber orientation distribution (FOD), which is introduced in Sec. 2.3 and Sec. 2.5, can represent the
complicated crossing fiber’s geometry [6]. However, the multiple peak lobes of the FOD function pose a
challenge for image processing. Our solution is to decompose the FOD function into several independent
components, each containing only one peak lobe. A FOD function is conventionally represented under the
real spherical harmonics (SPHARMs) basis, which is introduced in Sec. 2.4, up to the order N:
F OD(θ, ϕ) = X
n,m
u
m
n Y
m
n
(θ, ϕ) = U
T Y (θ, ϕ), (4.1)
53

Figure 4.3: A vector field within a tube centering at a voxel P0. (A) Tube T
k
p0
(yellow box) centering at p0
along the peak direction vp0 of the k
th FOD component F
k
p0
(green). (B) The FODs within the tube and the
{F
k
pt
} are in green. (C) The vectors {vpt} are picked at voxels within the tube, and these vectors are used
to solve Eq. 4.3 and compute the supporting vector field V
k
p0
.
where Y
m
n
is the mth (−n ≤ m ≤ n) real SPHARM basis at the order n (0 ≤ n ≤ N) and u
m
n
is the
coefficient for the corresponding basis, U is the vector that represents all the coefficients u
m
n
, and θ and ϕ
are the polar and azimuth angles of the spherical coordinates in R
3
. For any FOD function, we expand it
using Eq. 4.1 on a unit sphere represented by a triangular mesh, search the peaks on the mesh, and accept
the peaks whose value is higher than a threshold T HD (e.g., 0.1). For a FOD function with K peak lobes,
we solve the following optimization problem for its decomposition:
arg min
U1,··· ,UK
||X
K
k=1
Uk − U||2
2 + λ1
X
K
k=1
||AkUk − AkU||2
2 + λ2
X
K
k=1
||BkUk||2
2
, (4.2)
where Uk are the coefficients for the decomposed single-peak FOD components, and Ak is the matrix
that represents the values of SPHARMs at neighboring directions around the k
th peak (vertices within
two rings of each peak). The first term enforces the sum of the decomposed single-peak components to
equal the original FOD; the second term enforces each component to equal the original FOD near the
corresponding peak; the third term suppresses each component around other peaks. We show an example
of a FOD function decomposition in Fig. 4.2, where a FOD function is decomposed into two single-peak
components.
54

4.2.2 Modeling Single Peak FOD Components as Flow of Vector Fields
For each single-peak FOD component, we model it with the flow of a smooth vector field, which supports
geometrically consistent interpolations of FOD components.
We represent the k
th single-peak component of the FOD function at a voxel p0 as F
k
p0
. We choose the
peak direction of F
k
p0
as the seeding vector vp0 of the local supporting vector field. Then we compute a tube
T
k
p0
, centering at p0, along the direction vp0 with a radius r and a height h (Fig. 4.3. (A)). For each voxel pt
within the tube T
k
p0
, we choose the single-peak component F
k
pt
(Fig. 4.3. (B)) whose peak direction vpt
is
closest to vp0
, and the peak direction vpt
is a vector at pt
(Fig. 4.3. (C)). We do not pick any vector for voxels
without a valid peak direction whose angular difference is less than a threshold θ to the seeding vector
vp0
. These peak vectors {vpt} form a vector field within this tube, and we use a second-order polynomial
to fit each component of this vector field:
arg min
a0,a1,··· ,a9
card(T
k
p0 X
)
t=0
||v
d
pt − a0 −
X
3
i=1
aixi −
X
l+h=2
anx
l
ix
h
j
||2
2 + λ3
X
9
n=4
||an||2
2
, (4.3)
where v
d
pt
(1 ≤ d ≤ 3) represent d
th component of the vector at voxel pt
. The second term regulates the
second-order coefficients for smoothness. The polynomials are used to model the vector field V
k
p0 within
the tube T
k
p0
that represents the k
th underlying fiber bundle locally around the voxel p0.
4.2.3 Rotation Calculation for SPHARM-based FODs
For a target point q where we perform the interpolation, we choose the nearest voxel p0, which has been
augmented with a set of tubes {T
k
p0
} and vector fields {V
k
p0
} through the computation of Sec. 4.2.2. For
each vector field V
k
p0
, we compute the vector vq at point q using its polynomial representation. Each of
the corresponding k
th single-peak FOD components F
k
pt
from voxels within one voxel distance to q are
55

Figure 4.4: In the flowchart, the original FODs are highlighted by red circles at the four corners. The FOD
functions are decomposed into single-peak components in the first row and locally fitted by vector field
flows. In the second row, we rotate the single-peak components at four corners to align with the vectors
at each target point, and then we compute their weighted mean at each target point. Finally, we combine
the single-peak components into complete FODs
used for interpolation. First, we rotate each single-peak component F
k
pt
to align its peak direction with the
vector vq.
An easy way to compute the rotation is R = exp([r]
×), where r is a vector with its direction determined by the crossing product between the peak vectors vpt
and vq, and its length is the angle between
vpt
and vq; [r]
× is the cross-product matrix of a vector r:
[r]
× =








0 −rz ry
rz 0 −rx
−ry rx 0








. (4.4)
Since the rotated single-peak FOD components are now aligned in direction, we can compute the weighted
mean of SPHARM coefficients, which is the interpolated FOD component corresponding to the k
th peak
around voxel p0. The weights can be inverse distance or linear interpolation weights. After interpolating
56

all the FOD single-peak components independently, we combine them into the complete interpolated FOD
function at point q. The flowchart of the method is shown in Fig. 4.4. Our framework independently
handles the single-peak components of different FODs and successfully obtains geometrically consistent
interpolation of complicated crossing fiber geometry.
An essential step for the interpolation above is to transform the FOD function by a rotation R. A
straightforward numerical way is to expand the FOD function on a spherical triangular mesh, rotate the
mesh to rotate the function, and compute the inner products between the rotated FOD function and each
of the SPHARM bases to obtain the coefficients. However, the numerical method is computationally expensive. Instead, use the closed-form solution proposed in Sec. 3.2.2 to derive a matrix from the rotation
R that can be applied to the coefficients of the SPHARMs. Let FODR represent the FOD function after
applying the rotation R. We have the following relation:
F ODR(θ, ϕ) = F OD(θr, ϕr), (4.5)
where (θr, ϕr) is the coordinate acquired by rotating the coordinate (θ, ϕ) with the inverse rotation R−1
.
We represent Eq. 4.5 using the SPHARMs:
X
n,m
v
m
n Y
m
n
(θ, ϕ) = X
n,m
u
m
n Y
m
n
(θr, ϕr), (4.6)
where v
m
n
and u
m
n
are coefficients for FODR and FOD, respectively. To solve v
m
n
, we follow the work in
Sec. 3.2.2; namely, we represent the SPHARM function Y
m
n
(θr, ϕr) by a linear combination of Y
m
n
(θ, ϕ).
For rotation R−1
, we follow Eq. 3.5 to decompose it as three successive rotations around three axes:
R
−1 = ZγYβZα, (4.7)
57

Figure 4.5: One interpolated slice of the FODs. For FODs within the red box, a zoomed view of the proposed
interpolation, the linear interpolation, and the ground truth are plotted.
where Zγ and Zα are the rotations around the current z − axis by angles γ and α, respectively; Yβ is the
rotation around the current y − axis by an angle β. Details of this decomposition are introduced in Sec.
3.2.2. Now we can compute the n
th-order coefficients of the rotated FODR in Eq. 4.6 using Eq. 3.14 :
vn = UDn
(α, β, γ)U
−1un, (4.8)
where U is the transformation from real to complex SPHARMs that is defined in Eq. 2.10 and U
−1
is the
inverse transformation defined in Eq. 2.11; vn is a real (2n + 1) vector whose mth element is v
m
n
in Eq.
4.6, and un is the vector of real SPHARM coefficients u
m
n
in Eq. 4.6. Dn
(α, β, γ) is the Wigner D-matrix
that can be computed using Eq. 3.13.
The computation achieved by the closed-form representation is efficient because it only involves smallsize matrix operations. For example, the largest matrix in coefficients computation for a FOD function
represented by up to 16th-order SPHARMs is 33 × 33.
4.2.4 Evaluation Methods
We compare the proposed FOD interpolation method with the linear interpolation of SPHARM coefficients, the most used method in FOD-based tractography. We measure the quality of the interpolated FOD
58

functions based on down-sampling; we down-simple a ground truth FOD volume data to half the resolution, interpolate the down-sampled data to the original resolution, and measure the interpolated FOD
functions against the ground truth data based on two metrics. The first metric measures the sharpness of
the interpolated FOD functions, which indicates the specificity and accuracy of the FOD function. Inspired
by the full width at half maximum (FWHM) in signal processing, we define the full area at half maximum
(FAHM) of an FOD function f as:
F AHM(f) = area({x : f(x) > max(
f
2
)})/4π. (4.9)
FAHM is more sensitive to boating effects than entropy [24] and generalized fractional anisotropy
[25]. Another metric is the relative error between the interpolated FOD function and the ground truth
FOD function. The relative error is the L
2 distance of two FOD functions divided by the L
2 norm of the
ground truth FOD function.
We also evaluate the effectiveness of the proposed method on tractography. We up-sample the FOD
volume images to super-resolution images, including the cortical spinal tract (CST) area that connects the
cortical surface to the internal capsule. Then, we run the popular tractography from the MRtrix [59] on the
original and super-resolution data. We use Topographic Regularity, an essential property widely presented
in motor and sensory pathways [85, 86, 37] and introduced in Sec. 3.3.1, to show the improvements of the
tractography on up-sampled FOD data. We measure the topographic regularity using the intuitive metric
discussed in 3.3.1, where the classical multidimensional scaling (MDS) is used to project both the beginning
(cortical surface) and ending points (internal capsule) of the streamlines of each CST bundle to R
2
. Then,
the Procrustes distance between the projected beginning and ending points is computed to characterize
how well topographic regularity is preserved during fiber tracking.
59

Figure 4.6: Boxplots from quantitative comparisons using data from 40 HCP subjects. (A) Mean FAHM of
the ground-truth FODs, upsampled FODs from the proposed and linear interpolation. (B) The relative error
between up-sampled FODs and ground truth. (C) Procrustes distance of CST bundles from fiber tracking
based on linear interpolation (original) and our up-sampled FODs.
4.3 Experiment Results
We evaluated the FOD interpolation using 40 HCP subjects [2], including 20 females and 20 males. We
reconstructed 16th-order SPHARM-based FODs [6] from the HCP data with an isotropic resolution of
1.25mm. For parameters in our method, we set λ1, λ2, and λ3 Eq. 4.2 and 4.3 to be 1; the radius r, height
h, and θ of the tubes to be three times of voxel size, five times of voxel size, and 10◦
.
The HCP FOD data is used as the ground truth for down-sampling-based evaluation. We show the
FODs from one interpolated slice of a subject in Fig. 4.5 and highlight the FODs from an ROI (red box)
where several fiber bundles cross. Contrasting to the proposed method, FODs from linear interpolation
tend to lose their sharpness. For each subject, we computed the FAHM and relative L
2
error for each
interpolated FOD, which was then used to compute the mean FAHM and mean relative L
2
error among
all interpolated FODs. We show the boxplots of these measures from the 40 HCP subjects in Fig. 4.6. (A)
and (B). The FAHM measurement shows our approach avoids the bloating effects and preserves a similar
level of sharpness to the ground truth FODs; the lower mean relative L
2
error to the ground truth from
our method further shows the proposed interpolation achieved more accurate interpolation.
60

Figure 4.7: Tractography of CST on original HCP data and upsampled data for three subjects, and for each
bundle, we downsample the number of streamlines to 1000 for visualization.
We up-sampled the 40 HCP FOD volume images around the CST region to super-resolution images
with an isotropic resolution of 0.25mm. Then, we ran FOD-based probabilistic tractography on the original
and up-sampled FOD data using the iFOD1 algorithm of the MRtrix software [59]. In each run, 104
seed
points are randomly selected, and the main parameters of iFOD1 are set as: step size = 0.02mm, which is
around 0.1 times the voxel size of the up-sampled image, and angle threshold = 7◦
. The same parameters
were used for the original HCP dataset to avoid the bias of parameters. Three representative examples
of the reconstructed CST bundle from the motor cortex to the internal capsule are shown in Fig. 4.7,
where we can easily see that the tracts from the super-resolution FOD by our proposed interpolation
method are much smoother and better reflect the somatotopic organizational principles of the CST from
neuroanatomy than the baseline tracking results. In Fig. 4.6. (C), boxplots of the results from Procrustes
analyses further confirm this observation and demonstrate that our geometric FOD interpolation algorithm
can significantly enhance the anatomical consistency of fiber tracking results.
61

Chapter 5
A Geometric Framework for Processing Fiber Orientation Distribution
Functions
5.1 Introduction
Diffusion MRI (dMRI) can quantify the diffusion of water molecules in biological tissue, which is restricted
by fibrous structures such as membranes. In the last decades, advanced improvements in image acquisition
techniques significantly increased the spatial and angular resolution of dMRI [2, 3, 4] and provided preconditions for advanced mathematical models such as diffusion orientation distribution function (ODF) [50],
diffusion spectrum imaging (DSI) [49], fiber orientation distribution (FOD) [5, 6, 7]. Compared with other
mathematical models, the FOD model does not require the pre-specification of the number of underlying
fibers, is more accessible to acquire the required data, and is computationally more efficient [53]. FOD has
become more popular in the last decade and facilitates the development of FOD-based tractography for
brain connectivity research.
However, fewer works have been proposed for analyzing and processing the FOD data. In particular,
interpolation, the essential operation in image processing [13], has received rare attention. The crucial
transformation, such as convolution, translation, rotation, scaling, or warping, requires accurate interpolation. We introduced in Chap. 4, the FOD interpolation also plays a critical role in the tractography
62

algorithms and highly affects the tractography’s performance. In FOD-based tractography, linear and
nearest-neighbor interpolations are commonly used. The nearest neighbor interpolation is the fastest in
computation and is often used in asymmetric-FOD-based tractography [100, 101]. Linear interpolation
is more prevalent in state-of-the-art tractography [59, 8, 99] since it is adopted for efficiency and continuity. Still, the nearest neighbor interpolation induces severe discontinuity. Linear interpolation suffers
from bloating effects [13], which degenerates the specificity of the high angular resolution data. Moreover,
it is discussed in the Chap. 4, linear interpolation generates artificial directions, which can lead to false
positive streamlines. In Chap. 4, we developed a vector-flow-based framework that performs the geometrically consistent interpolation of FOD functions to overcome the difficulties of previous interpolations.
We decomposed FOD functions with multiple peak lobes into single-peak components handled separately.
The proposed framework in Chap. 4 outperforms widely used linear interpolation and provides more
geometrically proper interpolation, but it still has some limitations that need to be improved. First, to
construct the vector field flow, we need a heuristic process to construct the tube, which supports the
vector field. This heuristic process reduces the consistency of the algorithm since the structure of the tube
highly depends on the regions of the brain. Second, we need a minimum number of FOD functions from
neighboring voxels for vector field fitting. These limit the usage and accuracy of the interpolations in
small or boundary regions where constructing a tube and vector field flow is difficult. More importantly,
a unified framework is still needed to analyze and process the FOD functions. We expect to embed the
FOD functions into a proper metric space so that we can quantitatively analyze and compare them, and
the operators, such as interpolation, can naturally be defined within this metric space.
This chapter inherits the spirit of the Chap. 4 and develops a proper space for FOD functions; namely,
we handle the peak lobes of a FOD function separately and embed them into the Wasserstein space, which
has excellent properties for analyzing and processing the peak lobes. Besides utilizing the Wasserstein
63

Figure 5.1: A FOD function with three peak lobes plotted in spherical radiation pattern (a) and on the unit
sphere (b).
space, we improved and modified the critical procedures in the vector-flow-based framework with advanced techniques. Firstly, we more accurately detect the peak lobes of each FOD function based on a
spherical deconvolution framework. Secondly, we improve the decomposition of the FOD functions using
an optimization framework with adaptive constraints. Thirdly, the decomposed single-peak FOD functions are clustered based on an advanced framework rather than the heuristic tube construction. Then,
single-peak FOD functions are embedded and processed in the Wasserstein space, and the interpolation
is computed based on the Wasserstein Barycenter. For efficient computation, we propose a Fast Approximation, which is a generalization of the vector-flow-based framework, for the Wasserstein Barycenter.
Instead of interpolating the shift information using polynomial-based vector field flow fitting, we compute the shift interpolation using the Karcher Mean of the directional vectors on the unit sphere, and
this new method leads to more flexible and easier usage. In experiments, we used synthetic datasets to
demonstrate the proposed framework can handle peak lobes with large angle differences and represent the
highly curved fiber geometry. We also use HCP data to quantify the accuracy of the proposed interpolation
and show that it achieves superior performance than the linear interpolation. Furthermore, we show our
interpolation method can significantly improve the performance of FOD-based tractography.
64

5.2 Method
This section will introduce the properties of FOD functions and the difficulties of processing them. Then,
we will develop new frameworks to handle these difficulties.
5.2.1 Challenges of FOD Functions Processing
The fiber orientation distribution (FOD), which is introduced in Sec. 2.3 and Sec. 2.5, represents the complicated cross or bifurcate fiber geometry and quantifying the fraction of fiber portions [5, 6]. The FOD
function is commonly represented as a function defined on the unit sphere S
2
at each voxel, as shown in
Fig. 5.1. For reconstruction computation and data storage, a FOD function f is conventionally represented
under the basis of real spherical harmonics (SPHARMs), which is introduced in Sec 2.4:
f(q) = X
n,m
u
m
n Y
m
n
(q) = Y (q)
TU, (5.1)
where Y
m
n
is the mth (−n ≤ m ≤ n) real SPHARM basis at the order n (0 ≤ n ≤ N) and u
m
n
is the
coefficient of the corresponding basis, U is the vector that represents all the coefficients {u
m
n }, and q is the
vector representation of a point on the unit sphere S
2
.
A FOD function is often sparse on the S
2
and comprises multiple sharp peak lobes, each representing
a different fiber bundle, as shown in Fig. 5.1. The FOD function’s sparsity and mixture peak structure
pose significant challenges for image processing and analysis. Indeed, to qualitatively analyze and process
the data, we commonly represent the data in metric or norm space. In this way, people can compare
and analyze the data, and general operators on data, such as interpolation, smoothing, and convolution,
can be defined. However, the FOD functions do not form a well-defined metric or norm space unless
we parameterize the FOD function using a mixture of known probability distributions such as the von
Mises-Fisher distribution used in [102]. FOD functions are commonly assumed to be nonparametric and
65

Figure 5.2: Probability distribution functions for the illustration of different distances. (a) shows three
disjoint distributions with single sharp peak lobes. (b) shows two single peak lobe distributions and a
distribution with two peak lobes.
are often embedded in a larger space. For example, FOD functions are commonly considered as a subset
of the L
2
(S
2
) since they are represented under the SPHARM basis. Linear interpolation is the trilinear
weighted mean based on the L
2 norm. By linearity, any weighted interpolation of SPHARM coefficients
is equivalent to the pointwise interpolation of the FOD function on the S
2
:
f(q) = X
i
w
i
f
i
(q) = Y (q)
T
(
X
i
w
iU
i
), (5.2)
where {w
i} are the non-negative weights summing to 1; {U
i} and {f
i} are coefficeints and FOD functions
from neighboring voxels, respectively; f is the interpolated FOD function. Besides interpolation, pointwise
smoothing and convolution can also be defined. However, the L
2 norm is improper for the functions with
multiple sharp peak lobes. For example, the pointwise interpolation only roughly cuts off each peak and
sums them up, as shown in Fig. 5.3.
66

Table 5.1: Distances between Probability distributions
Pairs (p1, p2) (p1, p3) (p2, p3) (p4, p5) (p4, p6) (p5, p6)
L
2
5.31 5.31 5.31 3.75 5.31 4.61
Arcdist 1.57 1.57 1.57 0.98 1.53 1.49
W1 0.2 0.7 0.5 0.47 0.1 0.42
W2 0.2 0.7 0.5 0.56 0.1 0.48
The Arcdist between probability distributions pi and pj is arcdist(
√pi
,
√pj ),
which is defined in (5.4).
Generally, the normalized FOD functions are regarded as the probability functions on S
2
. In the nonparametric Riemannian framework [13], the square root reparameterized probability functions distributed
on a Hilbert sphere:
Ψ(M) = {ψ : M → R≥0 |
Z
M
ψ
2
(x)dx = 1}, (5.3)
where M is a manifold. The normalized FOD function p
f(x)/||f||1 lies on Ψ(S
2
). The distance between
two functions on Ψ(M) is the inverse cosine of their inner product:
arcdist(ψ1, ψ2) = cos−1
(
Z
M
ψ1(x)ψ2(x)dx). (5.4)
However, we can easily show this arc distance is topologically equivalent to the L
2 distance except for
a square root reparameterization. Thus, we cannot expect this distance to work significantly differently
from the L
2 distance.
Kullback–Leibler divergence (KLD) [103] and Jensen–Shannon divergence (JSD) [104] are used to compare normalized FOD functions. However, the KLD and JSD are improper; these distances achieve the maximum for two distributions with disjoint supports, regardless of how close or how far their peak lobes are
from each other [105]. Compared with other distances, the Wasserstein p-distance [106], which considers
both the shape and shift differences, enables us to differentiate and characterize FOD functions better. The
Wasserstein p-distance comes from an important mathematical problem called optimal transport problem;
67

for more about this problem, please refer to Sec. 2.7. Let (M, d) be a manifold equipped with a metric, the
Wasserstein p-distance is a distance defined between two Borel probability measures on M:
Wp(p1, p2) =

 inf
γ∈Γ(p1,p2)
Z
M×M
d(x, y)
p
dγ(x, y)


1
p
, (5.5)
where p1 and p2 are Borel probability measures on M, d(x, y) is the distance between x and y, and Γ is
the set of all the couplings of p1 and p2, as defined in Eq. 2.20. A coupling γ is a joint probability measure
on M × M whose marginals are p1 and p2 on the first and second factors, respectively. If p1 = δ(x − η1)
and p2 = δ(y − η2) are Dirac delta functions, the only possible coupling is γ = δ(x − η1, y − η2) so
that Wp(p1, p2) equals the distance between η1 and η2 on the manifold M. Intuitively, the metric Wp is
the minimum distance of turning one probability distribution into another. W1 and W2 are commonly
used in various appilcations; W1 is also called earth mover’s distance, and W2 has excellent structures for
statistical analysis [106, 62]. The weighted Fréchet functional of several probability distributions {ps} is
defined as:
F(µ) = X
S
s
w
sWp(µ, ps)
p
, (5.6)
where {w
s} are the non-negative weights summing to 1. A Wasserstein Barycenter of {ps} is a minimizer
µ
∗ of F in P(M) (if it exists), where P(M) is the space of Borel probability measures defined on M. The
interpolation is naturally defined as the Wasserstein Barycenter.
We show the distances between several probability distributions based on different metrics in Fig. 5.2
and Table5.1. We show the probability distribution functions defined on interval M = [0, 1] for ease
of visualization. Probability distributions in Fig. 5.2. (a) show both L
2 norm and Arcdist achieve the
maximum for two sharp peak lobes with disjoint supports, regardless of their shift distance. On the other
hand, Both W1 and W2 equal the shift distance between peak lobes. Fig.5.2. (b) shows both L
2 norm and
Arcdist cannot identify the probability distributions with different numbers of peak lobes; we expect p4
68

Figure 5.3: Interpolations between Probability distribution functions based on different distances. Two
single peak lobe distributions are shown in (a). (b) to (e) show the interpolations between them based on
different metrics. One distribution with two peak lobes and another with three are shown in (f). (g) to (j)
show their interpolations.
and p6 to be closer since they represent single fiber bundles with close directions, while p5 contains a
crossing fiber pattern, which should be far from p4 and p6.
Moreover, we demonstrate the mean interpolations between two probability distributions. We design
the two probability functions under different situations: (a) Two probability functions contain only one
sharp peak lobes, represented as p1 and p2, respectively. (b) One function p3 contains two peak lobes, and
another function p4 contains three peaks. All these situations are common in the interpolation of FOD
functions. Results of interpolations are shown in Fig. 5.3. In the first row (a), the mean of L
2 norm (linear
interpolation) and the mean of Arcdist roughly cut off each peak and sum them up, and the interpolated
functions are bloating. The W2 barycenter is a shift of the short peak lobe to the middle of p3 and p4, and
this is the interpolation we expect for single-peak FOD functions. The W1 barycenter does not work well
because it is not unique; for example, both the L
2 norm mean and the W2 barycenter are W1 barycenters.
None of the interpolations in the second row (b) are satisfying. For the W2 barycenter, we expect the most
right peak lobe of p4 to be reduced in magnitude instead of moved to a new location. The mass-preserving
69

Figure 5.4: (a) The FOD function f with peak lobes mixed. (b) Deconluluted function fs from the FOD
function.
property of W2 distance can potentially mix information from different fiber bundles. Situations will be
more tricky for real FOD functions, which can contain over four peak lobes.
5.2.2 FOD Function Decomposition
Based on the analysis in the last section, the Wasserstein 2-distance is the best for processing the singlepeak normalized FOD functions. However, it does not work for multiple peak lobe FOD functions. To
analyze and handle the multiple-peak structure of the FOD function, we decompose the FOD function into
several independent components, each containing only one peak lobe. Since the angular resolution of the
diffusion MRI is limited, the peak lobes with close directions can be mixed, as in Fig. 5.4. (a). The local
maximum detection cannot identify these mixed FOD peak lobes. To finely differentiate the FOD peak
lobes, we developed a spherical deconvolution (SD) framework for peak detection. Besides reconstructing
70

the FOD function directly, the SD framework was applied to sharpen the orientation distribution function
(ODF) to acquire the FOD function [58, 107]. Similarly, we form the spherical convolution as:
f(q) = Z
S
2
fs(x)R
k
(q, x)dx, (5.7)
where Rk
(q, x) is the kernel function, and fs is an approximation of the mixed Dirac delta functions that
represent the peaks. The kernel function Rk
(q, x) is chosen as:
R
k
(q, x) = Ek

1 + q · x
2
k
, (5.8)
where Ek is a scaling constant to ensure the integral of Rk
(q, x) on S
2
is 1 for any q ∈ S
2
. Rk
(q, x) is an
approximation of the Dirac delta function δ(1 − q · x):
fs(q) = lim
k→∞ Z
S
2
fs(x)R
k
(q, x)dx, ∀q ∈ S
2
, (5.9)
where the convergence is in both uniform and L
2 norm sense for any function fs ∈ L
2
(S
2
) ∩ C(S
2
)
based on Theorem 2.31 and 2.34 in [54]. The deconvoluted function fs in (5.7) will become bloating and
finally equals the FOD functions as k goes to infinity. We choose k so that the kernel Rk
(q, x) performs
the same sharpness as the anisotropic FOD functions; namely, we ensure the kernel has the same Full
Area at Half Maximum (FAHM), which is defined in Eq. 4.9, to the FOD functions from the body of corpus
callosum. The deconvoluted fs (Fig. 5.4. (b)) is much sharper than the FOD function. The fast computation
of deconvolution (5.7) depends on the real SPHARMs representation of fs and the Funk-Hecke Formula
introduced in Eq. 2.13:
Z
S
2
X
n,m
v
m
n Y
m
n Ek

1 + q · x
2
k
dx =
X
n,m
v
m
n Y
m
n
(q)Gn, (5.10)
71

where {v
m
n } are the unknown SPHARM coefficients of function fs and Gn is defined as:
Gn = 2πEk
Z 1
−1
Pn(t)

1 + t
2
k
dt, (5.11)
where Pn is the Legendre polynomial of degree n. Eq. 5.7 becomes:
f(q) = X
n,m
v
m
n Y
m
n
(q)Gn. (5.12)
We solve Eq. 5.12 as follows:
arg min
v
||Av − s||2
2 + λ0||Bv||2
2
,
(5.13)
where v is the vector representation of unknown coefficients {v
m
n }; s represents the FOD function values
at locations Q = {q1, q2, . . . , qH}, where the values are larger than the commonly used FOD cutoff
threshold (e.g., 0.1); Matrix A comes from {Y
m
n
(qh)} and {Gn} in Eq. 5.12, and B suppresses the value of fs
uniformly at complementary set of Q. The deconvolution is a least square problem. Since the deconvoluted
function fs is much sharper than the FOD function, we use higher-order SPHARMs in Eq. 5.12 (e.g., up
to 24th order). If the maximum SPHARM order of fs is the same as the FOD function, we can derive a
closed-form solution for Eq. 5.12, but the function fs will be more bloating and noisy. After solving Eq.
5.13, the peaks of the FOD function are detected based on local maximums of the sharp function fs.
For the FOD function detected with K peak lobes, we solve the following optimization problem for its
decomposition:
arg min
U1,··· ,UK
||X
K
k=1
Uk − U||2
2 + λ1
X
K
k=1
||AkUk − AkU||2
2 + λ2
X
K
k=1
||BkUk||2
2
, (5.14)
72

Figure 5.5: The FOD function with multiple peak lobes is decomposed into five single-peak component
FOD functions, and the decomposed components are rescaled for better visualization.
where {Uk} are the coefficients for the decomposed single-peak FOD components; Ak is the matrix representing the values of SPHARMs at neighboring directions of k
th peak where the FOD values are higher
than half of the peak value. Bk is the matrix representing the values of SPHARMs except for neighboring directions where the FOD values are higher than 30% of k
th the peak value. The first term enforces
the sum of the decomposed single-peak components to equal the original FOD; the second term enforces
each component to equal the original FOD near the corresponding peak; the third term suppresses each
component except the neighboring of the corresponding peak. This equation is an improvement of Eq.
4.2 since both Ak and Bk are adaptive constraints depending on the sharpness of the k
th peak lobe, and
the decomposition can be applied to FOD functions represented under different orders of SPHARMs. We
show an example of the FOD function decomposition in Fig. 5.5, where the FOD function is decomposed
into five single-peak components.
5.2.3 New Metric and Interpolation for FOD functions
We choose the Wasserstein 2-distance as the underlying metric. (S
2
, dS
2 ) is a complete, separable and
compact manifold, where dS
2 is the great-circle distance on S
2
. The metric space (P(S
2
), W2), the probability measure space on S
2
equipped with W2, is complete, separable and compact according to Theorem
2.7 in [61].
73

Since W2 distance ignores the magnitude, we regard each normalized FOD function as a function in
P(S
2
). We define a new metric dF between two FOD function f1 and f2 as:
df (f1, f2) = q
W2(p1, p2)
2 +

||f1||1 − ||f2||1
2
,
(5.15)
where p1 and p2 are normalized probability functions from f1 and f2, respectively. It is defined on P(S
2
)×
R+, and we consider a FOD function represented as f = (f /||f||1, ||f||1). We can check that df is a welldefined metric that satisfies the triangle inequality and measures differences in shape, shift, and magnitude.
As discussed in Sec. 5.2.1, the Wasserstein 2 Barycenter is improper as the interpolation for FOD
functions with multiple peak lobes, but it works fine for interpolation between peak lobes. Interpolation
is generally used among local FOD functions, and there are underlying geometrical fiber structures of
the neighboring FOD functions; namely, the FOD functions from neighboring voxels belong to several
different fiber bundles. We assume the nonzero peak lobes {f
n
s } from neighboring {Fs} can be regrouped
based on distance in Eq. 5.15. Assume we regroup {f
n
s } into M new groups, the Barycenter of the mth
group {f
l
m} with L elements is computed by:
arg min
hm∈P(S
2)×R+
X
L
l
w
l
d
2
f
(f
l
m, hm), (5.16)
where {w
l} is a rearrange of {w
s} according to the regroup. The Barycenters {h
∗
m} composed the peak
lobes of the interpolated FOD function. The minimizer of Eq. 5.16 is the Wasserstein Barycenter between
L normalized functions {f
l
m/||f
l
m||1} multiplied by the magnitude PL
l w
l
||f
l
m||1.
The metric space (P(S
2
), W2) is complete, separable, and compact, and the Fréchet functional in Eq.
5.6 with p = 2 is continuous in this space. Thus, the minimizer (Wasserstein Barycenter) of the Fréchet
functional always exists due to compactness. The Wasserstein Barycenter is not always unique; for example, the Barycenter of two Dirac delta measures at antipodal points is not unique. However, the normalized
7

Figure 5.6: The flowchart of the geometric interpolation framework. The blue, purple, green, and red
arrows mean hierarchical clustering, peak lobe field generation, interpolation, and summation. The FOD
functions from neighboring voxels are decomposed and clustered at the first row. Each peak lobe is used as
the generator for a peak lobe field within the same cluster. The same color border highlights the peak lobe
field and its generator peak lobe. Interpolation is computed within each peak lobe field from the second
to the third row. Finally, we combine interpolations from all peak lobe fields to get the complete FOD
functions.
FOD peak lobes work as absolutely continuous measures since they are density functions. The strict convexity of Fréchet functional guarantees the uniqueness of the Wasserstein Barycenter between absolutely
continuous measures. The strict convexity of Fréchet functional is defined and proved for (P(R
n
), W2) in
Proposition 3.1.6 of [61]. But, this proposition is also true for (P(S
2
), W2). Indeed, Theorem 1.34 in [61]
shows there exists only one transport plan (Eq. 2.19), which is included by a transport map (Eq. 2.18), between absolutely continuous measures on S
2
. The proof in Proposition 3.1.6 of [62] works for (P(S
2
), W2)
based on the existence and uniqueness of the transport map.
75

5.2.4 Numerical Computation
Using a spherical triangular mesh, we discretize S
2
for numerical computation. In the discrete case, the
Wasserstein 2-distance between two probability measures p1 and p2 is computed by solving a linear programming problem:
W2(p1, p2)
2 = min
γ
X
i,j
d(xi
, yj )
2
γ(xi
, yj )



P
i
γ(xi
, yj ) = p2(yj )
P
j
γ(xi
, yj ) = p1(xi)
γ(xi
, yj ) ≥ 0,
(5.17)
where {xi} and {yi} are supports of p1 and p2. The equality constraints ’ dimensions are low since the
normalized peak lobes’ supports are sparse. The Wasserstein Barycenter between probability measures
{ps} is also computed by linear programming:
min
u,v,γ1,··· ,γS
X
S
s
w
sX
i,j
d(x
s
i
, yj )
2
γs(x
s
i
, yj )



P
j
γs(x
s
i
, yj ) = ps(x
s
i
)
P
i
γs(x
s
i
, yj ) = u(yj )
γs(x
s
i
, yj ) ≥ 0
−ϵ ≤ u − Cv ≤ ϵ,
(5.18)
where {x
s
i
} are supports of {ps}, C is the matrix representing the values of SPHARMs at location {yj};
the optimal u
∗
and v
∗
are the Barycenter measure and its SPHARM coefficients, respectively; ϵ controls the
precision of the SPHARM approximation of u. The last constraint is critical since it provides regularity to
76

Table 5.2: Algorithm: Geometric Interpolation
Inputs: FOD functions {F1, · · · , FS} from neighboring voxels {1, . . . , S} and target points
{pt}.
Output: Interpolated FOD functions {Ft} at target points.
A. Deconvolute and decompose each Fs into peak lobes {f
1
s
, · · · , f Ns
s } by solving Eq. 5.13 and
Eq. 5.14.
B. Computing the spatial weights {w
s
t
, 1 ≤ s ≤ S} at each pt
(trilinear weights).
C. Regroup all peak lobes with the following steps:
1. Compute the distance Eq. 5.15 between pairs of all peak lobes {f
n
s
, 1 ≤ n ≤ Ns, 1 ≤ s ≤ S}
by solving Eq. 5.17.
2. Cluster all peak lobes {f
n
s } into K clusters {f
(k,h)
s , 1 ≤ h ≤ Hk, 1 ≤ k ≤ K} using
hierarchical clustering with single-linkage.
3. Viewing each f
(k,h)
s as a generator and selecting the closest peak lobe (to f
(k,h)
s ) at each other
voxel, these peak lobes form a peak lobe field P
(k,h)
s .
4. Divide the magnitude of f
(k,h)
s by N
(k,h)
s , which is the number of peak lobe fields contain
f
(k,h)
s .
D. Interpolating at each of {pt}:
1. Computing the interpolation of peak lobes within peak lobe field P
(k,h)
s by solving Eq. 5.18
with spatial weights {w
s
t }.
2. Combine the interpolations from all the peak lobe fields to obtain the complete interpolation
at pt
.
the solution and enables direct computation of SPHARM coefficients. The size of the linear programming
Eq. 5.18 is higher than Eq. 5.17 and is still affordable due to the sparse supports.
In this work, we regroup the nonzero peak lobes {f
n
s } from neighboring voxels {1, . . . , S} with corresponding FOD functions {Fs} based on distance defined in Eq. 5.15. The flowchart of the regroup and
interpolation is shown in Fig. 5.6. To reduce computations, we first cluster {f
n
s } into K clusters using
hierarchical clustering with single-linkage, as shown in the first row in Fig. 5.6. The clustering threshold
T hC is approximately the greatest angle at which we allow the fibers belonging to the same fiber bundle
to bifurcate within one voxel. Within a cluster, each peak lobe is used as a generator; we pick up the
closest peak lobe (to the generator) at any other voxel, and these peak lobes form a peak lobe field (second
row in Fig. 5.6). Since each peak lobe can belong to multiple peak lobe fields, to ensure that each peak
lobe contributes properly, we multiply its magnitude by 1/N, where N is the number of peak lobe fields
77

containing the peak lobe. For example, in the third row of Fig. 5.6, the magnitude weight of the red peak
lobe from the left-up corner is multiplied by 1
4
since it belongs to 4 peak lobe fields. Then, interpolation
for each peak lobe field is computed based on Eq. 5.16 and Eq. 5.18. Finally, we combine all the peak
lobe field interpolations to obtain the complete interpolation. The interpolation algorithm, including the
preprocessing steps, is demonstrated in Table5.2.
5.2.5 Fast Approximation
When we need real-time interpolation from the neighboring FOD functions. The proposed interpolation
in Eq. 5.16, which needs solving linear programming Eq. 5.18, is not practical in big data processing. The
Wasserstein Barycenter simultaneously interpolates the shift and shape deformation information. Instead,
we interpolate them independently for fast computation.
Given {f
l
m} and weights {w
l}, since each of {f
l
m} is antipodal symmetrical, we only consider their
half components in the same hemisphere. We first compute the weighted mean among {f
l
m}, which is
achieved by computing the Karcher mean of the peaks on S
2
:
arg min
y∈S
2
X
L
l
w
l
d
2
S
2 (x
l
m, y), (5.19)
where {x
l
m} are peak locations of peak lobes {f
l
m}. The Karcher mean is unique when {x
l
m} are located
within the same hemisphere, and it is computed using iterative Algorithm 1 in [13]. Since {x
l
m} are close
to each other, we chose the normalized Euclidean center as the initial for fast convergence. Let y
∗ be the
Karcher mean, and the Dirac delta function at y
∗
is the Wasserstein Barycenter of delta functions at {x
l
m}.
Compared with the method in Chap. 4, which utilizes vector field flow to interpolate the shift information,
the Karcher mean has no minimum requirement of FOD functions’ number and can be applied to more
general interpolations. Next, we rotate each f
l
m to align its peak direction with the Karcher mean y
∗
. The
fast computation of each peak lobe f
l
m rotation is computed following the works in Eqs. 4.5, 3.8 and 4.8.
78

Figure 5.7: The interpolated FOD functions between the first peak lobe (first column in each figure) and the
second peak lobe (last column in each figure) using different methods. The second peak lobes are rotated
to perform different angle differences, shown on the left, to the first peak lobes. ’Linear’ represents the
linear interpolation, ’RF’ is the Riemannian Framework, ’WB’ is the Wasserstein Barycenter, and ’FA’ is
the Fast Approximation.
5.3 Experiments
We demonstrate the proposed geometry framework in synthetic and Human Connectome Project (HCP)
data [2]. We also apply the proposed interpolations to the tractography.
5.3.1 Synthetic Data
The diffusion MR images of HCP were acquired from 270 gradient directions distributed over three different
b-values: 1000, 2000, and 3000 s/mm2
. The resolution of the diffusion MR images is 1.25mm isotropic.
The HCP-Pipeline was used for preprocessing. The FOD functions with 8
th order and 16th order SPHARMs
are reconstructed using the method [6].
79

Figure 5.8: The Full Area at Half Maximum (FAHM) of interpolated FOD peak lobes in the middle columns
of each figure in Fig. 5.7. ’Linear’ represents the linear interpolation, ’RF’ is the Riemannian Framework,
’WB’ is the Wasserstein Barycenter, and ’FA’ is the Fast Approximation, ’Baseline’ represents the baseline
peak lobe. (a) and (b) are the 8
th and 16th FOD functions, respectively.
The first experiment shows the interpolations when we change the angles between two FOD peak
lobe functions. We select the baseline single-peak FOD function from the corpus callosum body of subject 100307 in HCP data. We first fix the baseline peak lobe to align with the y-axis and then rotate the
baseline to perform different angles (3
◦
to 60◦
) to the y − axis with a step size of 3
◦
. We test the linear
interpolation, the interpolation in Riemannian Framework (RF) [13], and the Wasserstein Barycenter and
its Fast Approximation proposed in Sec. 5.2.5. The results of the interpolations are shown in Fig. 5.7,
where three interpolations with weights [0.75 0.25], [0.5 0.5], and [0.25 0.75] between the first and second
peak lobes (first and last in each row) were performed. The linear and RF interpolations simply cut off two
FOD functions and sum them up. The Wasserstein Barycenter and Fast Approximation provide reasonable
interpolations. We use the metric Full Area at Half Maximum (FAMH) defined in Eq. 4.9, which measures
the sharpness, to measure the bloating effect quantitatively. We measure FAMH of the middle interpolation in each row of each figure in Fig. 5.7 and plot the degrees vs FAMH in Fig. 5.8. We can observe that
linear interpolation and RF interpolation dramatically lose sharpness when the angle difference exceeds
9
◦
, particularly for 16th order FOD functions. The interpolations of Wasserstein Barycenter show a little
80

Figure 5.9: The interpolated FOD functions of the synthetic half-cycle data are based on different methods.
Figure 5.10: The boxplots of the distances between the interpolated FOD functions and the ground truth
FOD functions in Fig. 5.9. ’Linear’ represents the linear interpolation, ’RF’ is the Riemannian Framework
interpolation, ’WB’ is the Wasserstein Barycenter, and ’FA’ is the Fast Approximation. (a) and (b) are the
8
th and 16th FOD functions, respectively.
bloating effect when the angle difference exceeds 50◦
in 8
th order FOD functions. The Fast Approximation
avoids the bloating effect since we interpolate the shift and shape information independently.
Next, we use synthetic data to show that the Wasserstein-metric-based framework and its fast approximation can successfully model the bending geometry of the fiber bundles. We designed a synthetic
half-circle FOD image (ground truth in Fig. 5.9) with an inner radius of 1 mm and an outer radius of 6
mm. This dataset contains different fiber bending levels (radius of curvature) in practical tractography,
particularly the highly curved U-fibers. The resolution of this synthetic image is 1mm isotropic, and we
Table 5.3: Mean Distances between the interpolated FOD functions and the ground truth
Methods Linear RF WB FA
8
th 0.093 (5.3
◦
) 0.096 (5.5
◦
) 0.058 (3.3
◦
) 0.056 (3.2
◦
)
16th 0.126 (7.22◦
) 0.127 (7.28◦
) 0.041 (2.34◦
) 0.048 (2.76◦
)
81

Figure 5.11: One interpolated slice of the FODs. For FODs within the red box, a zoomed view of the Ground
truth, linear interpolation, Wasserstein Barycenter, and Fast Approximation.
downsample it to 2mm (downsampled in Fig. 5.9). We interpolate the downsampled data to 1mm by different methods and compare them with the ground truth data. The FOD functions in this synthetic image
also come from the corpus callosum body of subject 100307 in HCP data. In Fig. 5.9, we can observe that
the linear and the RF interpolations generate the FOD function with crossing peak lobes, particularly in
the highly bending regions; the Wasserstein Barycenter and its approximation provide FOD functions that
follow the bending geometry. Moreover, we measure the proposed distance df , which is defined in Eq.
5.15, between the interpolated FOD functions and the ground truth FOD functions for quantitative analysis. The boxplots of the measurements are shown in Fig. 5.10, and the mean distances to the ground truth
are presented in Table5.3. We can observe that the linear and RF interpolations generate more significant
errors in 16th order FOD functions that are sharper. At the same time, both Wasstersterin Barycenter and
its Fast Approximation outperform the other two methods, particularly for high-order FOD functions.
5.3.2 HCP Data
We evaluated the FOD interpolation using 40 HCP subjects, including 20 females and 20 males. For parameters in our method, we set k in Eq. 5.7 to be 100 to perform the same FAMH to the FOD functions from
the body of the corpus callosum. λ1 and λ2 in Eq. 5.14 to be 1. The clustering threshold T hc is set to be
0.8 (around 45◦
); namely, assume the fibers from the same fiber bundle can not bifurcate over 45◦ within
one voxel size.
82

Figure 5.12: Boxplots of the FAHM of the interpolated FOD functions in (a) and the distance between the
interpolated FOD functions to the ground truth FOD functions in (b). ’Linear’ represents the linear interpolation, ’GR’ is the ground truth, ’WB’ is the Wasserstein Barycenter, and ’FA’ is the Fast Approximation.
We choose the Region of Interest (ROI) within the Corpus Callosum for affordable computations, particularly the part connecting the two lower limb regions, two upper limb regions, and two trunk regions. We
down-sampled the ground-truth HCP FOD volumes to half their resolution, interpolated the down-sampled
data to their original resolution, and measured the interpolated FOD functions against the ground-truth
HCP volume. Since the RF interpolation performs almost the same as the linear interpolation, as shown in
the last section, we only compare linear interpolation, Wasserstein Barycenter, and Fast Approximation in
this experiment. We show the FODs from one interpolated slice of the HCP subject 100307 in Fig. 5.11 and
highlight the FODs from an ROI (red box) near the cortex. Contrasting to the Wasserstein Barycenter and
its approximation, FODs from linear interpolation tend to lose their sharpness. For each subject, we computed the mean FAHM and mean distance df (5.15) to ground truth for each interpolated FOD. We show
the boxplots of these measures from the 40 HCP subjects in Fig. 5.12. The FAHM measurement shows the
Wasserstein Barycenter and the Fast Approximation avoids the bloating effects and preserves a sharpness
similar to the ground truth FODs. At the same time, the linear interpolations generate bloating effects.
Moreover, the lower distances to the ground truth show that the proposed two interpolations achieved
interpolations that are closer to the ground truth data.
83

Figure 5.13: The tractography results on the synthetic half-circle data based on different interpolations.
5.3.3 Application in Tractography
In this experiment, we apply the proposed interpolations to tractography. We develop the basic probabilistic walk-based tractography algorithm; namely, at each step of a track, we sample the direction based
on the interpolated FOD from neighboring FOD functions. We use rejection sampling to sample the direction based on the distribution of the interpolated FOD function. We compare linear interpolation and Fast
Approximation for efficiency, which performs almost the same as Wasserstein Barycenter.
We first apply the tractography algorithm to the downsampled synthetic half-cycle data in Fig. 5.9. The
step size of the tracking algorithm is set to be 0.025 times the voxel size. Since the largest turn angle of FOD
functions between adjacent voxels in this data is 45◦
, we set the angle threshold of the tracking algorithm to
45◦
. We uniformly select the seeds from the left end of the half circle. The streamlines of the tractography
are shown in Fig. 5.13. The streamlines based on the Fast Approximation interpolation smoothly follow the
circle geometry with less oscillation. In contrast, the streamlines of linear interpolation tend to perform
84

Figure 5.14: Tractography of CST on original HCP data and upsampled data for three subjects, and for
each bundle, we downsample the number of streamlines to 200 for visualization.
weird sharp turns near the circle center, and streamlines near the outside perform much better. This
experiment shows that linear interpolation makes reconstructing highly curved fibers in the brain hard.
Next, we up-sample the FOD volume images to super-resolution images using the Fast Approximation
interpolation. We up-sampled the 40 HCP FOD volume images around the cortical spinal tract (CST)
region to super-resolution images with an isotropic resolution of 0.25mm. Since we assume the peak
lobes belonging to the same fiber bundle can not bifurcate over 45◦ within 1.25mm (resolution of HCP
data), we expect the peak lobes belonging to the same fiber bundle in the up-sampled data don’t bifurcate
over 9
◦
. As shown in Fig. 5.9 and Fig. 5.10, when the angle difference between the 16th order peak lobes is
less than 9
◦
, the linear interpolation performs very close to the Fast Approximation and the bloating effect
is weak. Thus, we apply the developed probabilistic tractography with linear interpolation instead of Fast
Approximation to the up-sampled data for computational efficiency. The parameters are set as step size =
0.025mm, which is 0.1 times the voxel size of the up-sampled image, and angle threshold = 7
◦
. We also
directly apply the developed probabilistic tractography algorithm with linear interpolation to the original
HCP FOD data to compare with the tractography results on up-sampled data. The tractography based on
the original HCP data used the same parameters to avoid parameter bias. Three representative examples of
the reconstructed CST bundle from the motor cortex to the internal capsule are shown in Fig. 5.14, where
we can easily see that the streamlines from the up-sampled image are smoother and better reflect the
85

Figure 5.15: Boxplots of the Procrustes distances for CSTs of 40 subjects, HCP represents the reconstructed
CST on the original HCP images, and Upsampled presents the reconstructed CST on the Upsampled superresolution images.
somatotopic organizational principles of the CST from neuroanatomy than the direct linear interpolationbased results. We use Topographic regularity to quantitatively measure the regularity of the streamlines,
which is an important property present in the motor and sensory system [83, 84] and introduced in Sec.
3.3.1. We use the intuitive metric introduced in Sec. 3.3.1, where the classical multidimensional scaling
(MDS) is used to project both the beginning (cortical surface) and ending points (internal capsule) of the
streamlines of each CST bundle to R
2
. Then, the Procrustes distance between the projected beginning
and ending points is computed to characterize how well topographic regularity is preserved during fiber
tracking. The boxplots of the Procrustes distances are shown in Fig. 5.15. The results demonstrate that
our geometric FOD interpolation algorithm can significantly enhance the anatomical consistency of fiber
tracking results.
86

Chapter 6
Topographic Filtering of Tractograms as Vector Field Flows
6.1 Introduction
Diffusion MRI (dMRI) is currently the most widely used in vivo imaging tool for studying human brain
connectivity. While there has been tremendous progress in improving both the spatial and angular resolution of dMRI, especially the multiband and multi-shell imaging popularized by the Human Connectome
Project (HCP) [2], there is still much debate about how to relate these complicated connectome imaging
data to the underlying brain circuits [108], especially the validity of the commonly used tractography
techniques [9]. For several significant neurocircuitries, a striking property is their fiber pathways follow
a topography-preserving relation when they travel between anatomical regions [33, 32]. The wide presence of such topography-preserving connectivity, i.e., connectopy [32], provides a great opportunity to
examine and improve the biological validity of tractography techniques [85, 86]. Using the fiber orientation distributions (FODs) reconstructed from connectome imaging data, we develop a novel method to
construct a vector field flow representation of topographically regular fiber pathways here. We then apply
this novel representation to develop a tractogram filtering method to improve the topographic regularity
of fiber bundles with known connectopy.
The wide presence of connectopy in various brain regions was first identified in post-mortem studies, such as the retinotopy of the visual pathways [109, 110]. The in vivo characterization of anatomical
87

topography in the early visual system was first achieved with the popularization of electrophysiological studies [111]. With task-based functional MRI (tfMRI), the topographically ordered mapping between
anatomy and function in various brain regions has been revealed in vivo, including the retinotopic organization of axons in the visual system [112], somatotopic organization of the motor and sensory pathways
[83], the tonotopic organization of the auditory pathways [113], and the anatomy of cortico-striatal pathways [114]. In contrast to the rich literature on functional connectopy, very limited in vivo research has
been conducted on the structural connectopy of tractography techniques. Recently, a novel FOD-based
tractography method was developed [85] by incorporating topographic regularity and achieved improved
regularity for the reconstruction of visual pathways and corticospinal tracts. In addition, A quantitative
measure called group graph spectral distance (G2
SD) was proposed [86] to evaluate the topographical
regularity of fiber tracts from any tractography techniques and then applied for the filtering of irregular
tracts.
Instead of modeling the regularity among fiber tracts, which can be affected by tractography techniques
and various parameters, the novel method we develop in this work directly considers how the connectopy
of the anatomical pathways are represented in the dMRI data. Using the FODs (Fig. 6.1. (B)) reconstructed
from the connectome imaging data, our main idea is to extract an underlying vector field that supports
the topographically regular fiber bundle (Fig. 6.1. (A)). This work is achieved by formulating and solving a
Markov Random Field (MRF) problem [115] from the peak fields of the FODs (Fig. 6.1. (C)). After that, we
develop a quantitative measure to evaluate how well the fiber tracts fit as the vector field flows (Fig. 6.1.
(D)). Using HCP data, we demonstrate our method in filtering the tractograms of two fiber bundles with
known connectopy: the visual pathway [112] and the callosal motor pathways [116]. We also quantify the
effectiveness of our method with the independent (G2
SD) measure from previous work on topographic
tract filtering [86]. This chapter has a published version [37].
88

Figure 6.1: A vector field flow representation for the visual pathway bundle. (A) Fiber bundle from FODbased tractography. The red box is the ROI used for visualization in (B)-(D). (B) The FODs in the ROI. (C)
The multiple peaks of the FODs. (D) The vector field extracted from the FOD peaks. Note that we use the
QIT software ( http://cabeen.io/qitwiki) to visualize all FOD peaks and vector fields as sticks in this work.
6.2 Method
6.2.1 Vector Field Flows for Topographically Regular Tracts
In neuroscience, topographic regularity is commonly understood as the preservation of spatial relationships of axonal fibers as they traverse the space and connect distant brain regions [33]. Following this idea,
the tracts of a topographically regular fiber bundle should not intersect, and the points on the tracts need
to preserve their relative topological relation along the bundle. This characteristic is consistent with the
flows induced by a smooth vector field. More specifically, for a continuously differentiable vector field v(x)
89

in R
3
, its flows are integral curves defined as the solution of the following ordinary differential equation
(ODE):
dx
dt
= v(x), x ∈ Ω, (6.1)
where Ω is a bounded open set in R
3
.
From the theory of the ODE and differential manifolds [117], the flows of a vector field intersect (asymptotically) only at the stationary points (v(x) = 0). Moreover, the local flows are diffeomorphisms that
preserve the topological structure (details can be found in [117]). For example, the local flows map a
connected open set to another connected open set, and the map is one-to-one. This property means the
flows automatically satisfy the neuroscientific understanding of topographic regularity. Interestingly, this
definition was naturally adopted in the original tractography algorithm by Basser et al. [55], where the
streamlines were defined as the flows of the field of principal eigenvectors from the tensor model. The
challenge, however, is that the field of eigenvectors from the tensor model is not smooth vector fields accurately depicting the underlying white matter anatomy. The flows from the eigenvector fields are thus
inappropriate for studying fiber pathways passing through regions with crossings. Modern tractography
techniques use advanced models such as fiber orientation distribution (FOD) to resolve better crossing
fibers [5, 6]. With the ability to more faithfully represent multiple fiber directions, the FOD representation
provides the opportunity to truly analyze fiber tracts as flows of vector fields and study their topo-graphic
regularity. Given any fiber bundle from FOD-based tractography, the main technical hurdle is thus to extract the associated vector field from the FOD fields. Once this vector field is computed, we can quantify
the degree of the tracts as vector field flows.
90

Figure 6.2: Main steps in our topographic tract filtering algorithm. (A) Input tracts of the callosal motor
(CM) pathway from FOD-based tractography. (B) The field of MVD from FOD peaks for the ROI highlighted
by the yellow box in (A). (C) The PVF for the CM pathways is computed by solving the MRF model. (D)
The preserved fiber tracts after removing 30% of tracts with high VFD value. (E) For the ROI highlighted
by the white box in (A), the 10% of tracts with the highest VFD values are plotted as white tubes together
with the original input tracts. This helps visualize the irregularity of the filtered tracts by our algorithm.
6.2.2 Vector Field Computation from FODs
For the FOD represented as spherical harmonics at each voxel, we extract its peaks as local maxima on a
spherical mesh and denote them as a set of normalized vectors v1, v2, · · · ,, which is called a multi-valued
directional (MVD) in graphics literature [118]. Because the FOD is symmetric with respect to the sphere’s
origin, we only take one direction from each pair of symmetric peaks. For each direction, we associate a
magnitude. For example, they can be the FOD values in corresponding directions. The collection of the
directionals over the image volume is called an MVD field (Fig. 6.2. (B)). Typically, the number of peaks
in FODs is chosen to be less than four, which matches the sparsity assumption in the FOD reconstruction
algorithm [6].
Given a tractogram of a fiber bundle, we will compute a principal vector field (PVF) to perform vector
field flow analysis. To extract the PVF directly from an MVD field, we formulate and solve a Markov
91

Random Field (MRF) problem. Especially once the MVD field is computed, we only need to assign an
index to the set of vectors (directionals) at each voxel to define a vector field over the brain volume. For
a voxel p, we denote lp as the index into the directional at this voxel. The vector field decomposition,
then, is essentially an integer programming problem. To extract a smooth vector field, we will introduce
regularization between neighboring voxels by measuring the angular differences of vectors. Overall, we
maximize the log-likelihood of the posterior distribution of MRF [119] as defined:
E =
X
p∈Ω
θp(lp) + X
(p,q)∈N
θpq(lp, lq), (6.2)
where Ω is the image domain under consideration, N denotes the neighborhood of voxels, and lp is the
index of the directional at the p
th voxel. Two-factor terms are the unary factor θp and the pairwise factor
θpq.
At each voxel, the decomposed vector should be close to the tangent vector of a tract passing through
this voxel. We use the average tangent vectors if more than one tract passes through this voxel. Moreover,
the vector field should follow a relatively high probability of direction (strong diffusion) whenever possible.
The unary factor θp is thus defined as follows:
θp(lp) = λ1F(lp) + λ2(p)| < vlp
, up > |, (6.3)
where F is the associated FOD magnitude of the directional vlp
, and up is the normalized average tangent
vector of the tracts that pass voxel p. For a voxel p, where there is no tract passing, up is set to 0. When
the directional vector is aligned to the tangent vector, the absolute value of the inner product is large, and
a higher probability directional contributes more to the unary factor. The regularization parameter λ2 is
set to be in proportion to the density of fiber tracts in the voxel, i.e., λ2(p) = k ∗ d(p), where d(p) is the
number of tracts passing through voxel p, and k is a positive number.
92

For smoothness regularization, the pairwise factor θpq is defined as follows:
θpq(lp, lq) = λ3 | < vlp
, vlq > |, (6.4)
The pairwise factor measures the similarity of the directional vectors of each pair in the neighborhood
system N. In our implementation, we use the 6-connected neighborhood system.
Numerically, we use the max-sum Belief Propagation (BP) algorithm [115] to solve the MRF problem for
efficient computation, which generates a smooth PVF that supports the tractogram under consideration.
6.2.3 Measure for Topographic Regularity
The smooth principle vector field (PVF) generated in Sec. 6.2.2 is consistent with the trend of the fiber
bundle. To remove irregular fiber tracts that do not follow the PVF, we define a measure called Vector
Flow Deviation (V F D) to quantify the consistency between a tract T(t) and the PVF:
V F D(T) =
(
R
(vT(t) − ut)
2dt)
1/2
L(T)
, (6.5)
where the integral is defined over the tract T, ut
is the normalized tangent vector of T at point T(t), vT(t)
is the normalized PVF at point T(t), and L(t) is the length of the tract. The V F D measure computes the
mean square error between the tangent vector of the tract and the vector field. Based on the uniqueness
theorem from ODE, V F D(T) = 0 (tangent vector of T is perfectly aligned along the vector field pointwise)
if and only if T is an integral curve of the vector field.
Given a set of precomputed input tracts and the corresponding FOD data, we can develop a topographic
tract filtering method based on the V F D measure. As summarized in Fig. 6.2, we first compute the MVD
(Fig. 6.2. (B)) from the FODs and then obtain the PVF (Fig. 6.2. (C)) by solving the MRF optimization
problem. Finally, we compute the V F D to measure the topographic regularity of each tract and generate
93

Figure 6.3: Representative tract filtering results from 3 HCP subjects for the visual (left) and CM (left)
pathway. For each case, regions with outlier tracts were highlighted on the input tracts.
the tracts to be preserved (Fig. 6.2. (D)) and removed (Fig. 6.2. (E)). Our algorithm is independent of the
tractography algorithm used to generate the fiber tracts. We only need to adjust the threshold value for
different fiber bundles for the V F D measure. Solving the MRF contributes most to the computational cost
of our algorithm. The algorithm usually takes 5-8 minutes for 10000 tracts on a 2.8G HZ 2-core Laptop.
6.3 Experimental Results
This section presents experimental results to demonstrate our topographic tract filtering algorithm based
on vector field flows. The preprocessed MRI data of 30 subjects released from the HCP were used in our
experiments. From the multi-shell diffusion MRI data, FODs were first computed using the method in [6]
before the vector field analysis. We applied our method to filter two fiber bundles: optic radiation of the
visual pathway and the colossal motor (CM) pathway, which are known to have topographic regularity.
The inputs to our method are tractograms of around 104
tracts reconstructed by FOD-based tractography
in MRtrix [59] and related anatomical ROIs. For the optic radiation, the lateral geniculate nucleus (LGN)
region detected based on the method in [120] was used as the seed ROI, and the V1 cortex region from
94

Figure 6.4: Box plots of the average G2
SD measures of the removed and preserved bundle of 30 HCP
subjects. (A) Visual pathway; (B) CM pathway.
Free-Surfer was used as the inclusion ROI. For the CM pathways, the corpus callosum ROI from FreeSurfer
segmentation was used as the seed ROI, and the precentral gyrus from both hemispheres was used as the
inclusion ROIs. For our method, the parameters used in the MRF model for vector field extraction are set
as λ1 = 1, λ3 = 10, and k = 0.1 for all experiments. The Belief Propagation (BP) algorithm typically
converges within 4-6 iterations. The filtering threshold for the CM pathway is set to remove 30% of the
tracts, and the visual pathway to remove 40% of the tracts for all subjects.
Three representative examples of the filtering of the visual and CM pathway are shown in Fig. 6.3.
Our method can clearly remove the outlier tracts and generate much cleaner fiber bundles based on the
parameters discussed above. Some regions have been highlighted for visual inspection with the dashed
white ellipsoids on the original input bundles.
Besides visual assessment, we also conducted quantitative evaluations using the group graph spectral
distance (G2
SD) proposed in [86]. The G2
SD distance is low for topographically organized tracts and high
for irregular tracts. For each bundle, we calculated the G2
SD measure for the removed tracts after filtering
and the tracts that were preserved. Based on our filtering parameters, 40% of the tracts were removed,
and 60% of the tracts were preserved for the visual pathway; 30% of the tracts were removed, and 70% of
the tracts were preserved for the CM pathway. The G2
SD was calculated for each tract, and a mean value
was calculated for the removed and pre-served bundle for each subject. The boxplots of the mean G2
SD
95

values of the removed and preserved bundle are shown in Fig. 6.4. We can see that our method correctly
removed the tracts with higher G2
SD measures for both bundles, which shows the high consistency of our
method with previous work on topographic filtering [86].
96

Chapter 7
Conclusions
In this thesis, we developed four main works that improve dMRI data analysis and tractography performance in connectome research. In Chap. 3, we developed a novel surface-based probabilistic tractography
for reconstructing the fiber connections in the superficial white matter (SWM) layer. This new tractography outperforms traditional tractography and provides us with opportunities to study the connectivity
in SWM. In Chap. 4, we develop a vector field flow framework to interpolate the fiber orientation distribution (FOD) functions. The new interpolation overcomes critical difficulties in traditional interpolations
and alleviates the propagation errors in tractography algorithms. In Chap. 5, we proposed a geometric
framework that is based on the Wasserstein space to process and analyze the FOD functions. This geometric framework provides a novel interpolation to generalize the vector field flow framework. In Chap.
6, to reduce the invalid connections in tractography, we develop a flow-based filter algorithm that aims to
improve the consistency between the tractography and the fiber geometry. We present some promising
works to extend our current works.
The most straightforward work is to combine the proposed interpolations and the surface-based tractography. We can apply the proposed advanced interpolations to provide more accurate FOD functions
around the SWM layer. This will improve the connectivity research in the SWM and provide us with a
more accurate pathologic analysis of SWM.
97

Current direct FOD-based research is rare, and most works use the FOD functions for tractography
and analyze the tractography results. However, tractography is error-prone and algorithm-dependent; it
is not reliable for quantitative statistical analysis. The critical limitation of direct FOD-based research is
that no excellent tools exist to analyze the FOD functions. Kinds of literature [62, 61] have shown many
excellent properties in the Wasserstein 2 space for data analysis. These works inspired us to develop
tools for analyzing the FOD functions. We expect to develop statistical tools based on FOD functions for
longitudinal and pathological analysis.
Integrating the proposed interpolations with existing state-of-the-art tractography algorithms [99, 59]
is appealing and can potentially improve tractography performance in various fiber bundles.
98

Bibliography
[1] Peter J Basser, James Mattiello, and Denis LeBihan. “MR diffusion tensor spectroscopy and
imaging”. In: Biophysical journal 66.1 (1994), pp. 259–267.
[2] David C Van Essen, Kamil Ugurbil, Edward Auerbach, Deanna Barch, Timothy EJ Behrens,
Richard Bucholz, Acer Chang, Liyong Chen, Maurizio Corbetta, Sandra W Curtiss, et al. “The
Human Connectome Project: a data acquisition perspective”. In: Neuroimage 62.4 (2012),
pp. 2222–2231.
[3] Kamil Uğurbil, Junqian Xu, Edward J Auerbach, Steen Moeller, An T Vu,
Julio M Duarte-Carvajalino, Christophe Lenglet, Xiaoping Wu, Sebastian Schmitter,
Pierre Francois Van de Moortele, et al. “Pushing spatial and temporal resolution for functional
and diffusion MRI in the Human Connectome Project”. In: Neuroimage 80 (2013), pp. 80–104.
[4] Stamatios N Sotiropoulos, Saad Jbabdi, Junqian Xu, Jesper L Andersson, Steen Moeller,
Edward J Auerbach, Matthew F Glasser, Moises Hernandez, Guillermo Sapiro, Mark Jenkinson,
et al. “Advances in diffusion MRI acquisition and processing in the Human Connectome Project”.
In: Neuroimage 80 (2013), pp. 125–143.
[5] J-Donald Tournier, Fernando Calamante, and Alan Connelly. “Robust determination of the fibre
orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical
deconvolution”. In: Neuroimage 35.4 (2007), pp. 1459–1472.
[6] Giang Tran and Yonggang Shi. “Fiber orientation and compartment parameter estimation from
multi-shell diffusion imaging”. In: IEEE transactions on medical imaging 34.11 (2015),
pp. 2320–2332.
[7] 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”. In: Neuroimage 23.3 (2004), pp. 1176–1185.
[8] Dogu Baran Aydogan, Russell Jacobs, Stephanie Dulawa, Summer L Thompson,
Maite Christi Francois, Arthur W Toga, Hongwei Dong, James A Knowles, and Yonggang Shi.
“When tractography meets tracer injections: a systematic study of trends and variation sources of
diffusion-based connectivity”. In: Brain Structure and Function 223 (2018), pp. 2841–2858.
99

[9] KH Maier-Hein, PF Neher, JC Houde, MA Côté, E Garyfallidis, J Zhong, M Chamberland, FC Yeh,
YC Lin, Q Ji, et al. The challenge of mapping the human connectome based on diffusion tractography.
Nat Commun 8: 1349. 2017.
[10] Miguel Guevara, Pamela Guevara, Claudio Román, and Jean-François Mangin. “Superficial white
matter: A review on the dMRI analysis methods and applications”. In: Neuroimage 212 (2020),
p. 116673.
[11] Almut Schüz and Valentino Braitenberg. “The human cortical white matter: quantitative aspects
of cortico-cortical long-range connectivity”. In: Cortical Areas. CRC Press, 2002, pp. 389–398.
[12] Harumichi Shinohara, Xiaoliang Liu, Riho Nakajima, Masashi Kinoshita, Noriyuki Ozaki,
Osamu Hori, and Mitsutoshi Nakada. “Pyramid-shape crossings and intercrossing fibers are key
elements for construction of the neural network in the superficial white matter of the human
cerebrum”. In: Cerebral Cortex 30.10 (2020), pp. 5218–5228.
[13] Alvina Goh, Christophe Lenglet, Paul M Thompson, and René Vidal. “A nonparametric
Riemannian framework for processing high angular resolution diffusion images and its
applications to ODF-based morphometry”. In: NeuroImage 56.3 (2011), pp. 1181–1201.
[14] Marco Catani, Flavio Dell’Acqua, Francesco Vergani, Farah Malik, Harry Hodge, Prasun Roy,
Romain Valabregue, and Michel Thiebaut De Schotten. “Short frontal lobe connections of the
human brain”. In: cortex 48.2 (2012), pp. 273–291.
[15] Joshua D Burks, Lillian B Boettcher, Andrew K Conner, Chad A Glenn, Phillip A Bonney,
Cordell M Baker, Robert G Briggs, Nathan A Pittman, Daniel L O’Donoghue, Dee H Wu, et al.
“White matter connections of the inferior parietal lobule: A study of surgical anatomy”. In: Brain
and behavior 7.4 (2017), e00640.
[16] Alexandre Pron, Christine Deruelle, and Olivier Coulon. “U-shape short-range extrinsic
connectivity organisation around the human central sulcus”. In: Brain Structure and Function
226.1 (2021), pp. 179–193.
[17] Francesco Vergani, Luis Lacerda, Juan Martino, Johannes Attems, Christopher Morris,
Patrick Mitchell, Michel Thiebaut De Schotten, and Flavio Dell’Acqua. “White matter connections
of the supplementary motor area in humans”. In: Journal of Neurology, Neurosurgery & Psychiatry
85.12 (2014), pp. 1377–1385.
[18] Miguel Guevara, Claudio Román, Josselin Houenou, Delphine Duclap, Cyril Poupon,
Jean François Mangin, and Pamela Guevara. “Reproducibility of superficial white matter tracts
using diffusion-weighted imaging tractography”. In: Neuroimage 147 (2017), pp. 703–725.
[19] Minhui Ouyang, Tina Jeon, Virendra Mishra, Haixiao Du, Yu Wang, Yun Peng, and Hao Huang.
“Global and regional cortical connectivity maturation index (CCMI) of developmental human
brain with quantification of short-range association tracts”. In: Medical Imaging 2016: Biomedical
Applications in Molecular, Structural, and Functional Imaging. Vol. 9788. SPIE. 2016, pp. 328–334.
100

[20] Claudio Román, Miguel Guevara, Ronald Valenzuela, Miguel Figueroa, Josselin Houenou,
Delphine Duclap, Cyril Poupon, Jean-François Mangin, and Pamela Guevara. “Clustering of
whole-brain white matter short association bundles using HARDI data”. In: Frontiers in
neuroinformatics 11 (2017), p. 73.
[21] Fan Zhang, Ye Wu, Isaiah Norton, Laura Rigolo, Yogesh Rathi, Nikos Makris, and
Lauren J O’Donnell. “An anatomically curated fiber clustering white matter atlas for consistent
white matter tract parcellation across the lifespan”. In: NeuroImage 179 (2018), pp. 429–447.
[22] Dmitri Shastin, Sila Genc, Greg D Parker, Kristin Koller, Chantal MW Tax, John Evans,
Khalid Hamandi, William P Gray, Derek K Jones, and Maxime Chamberland. “Surface-based
tracking for short association fibre tractography”. In: NeuroImage 260 (2022), p. 119423.
[23] Timothy EJ Behrens, Stamatios N Sotiropoulos, and Saad Jbabdi. “MR diffusion tractography”. In:
Diffusion MRI. Elsevier, 2014, pp. 429–451.
[24] H Ertan Çetingül, Bijan Afsari, Margaret J Wright, Paul M Thompson, and René Vidal. “Group
action induced averaging for HARDI processing”. In: 2012 9th IEEE International Symposium on
Biomedical Imaging (ISBI). IEEE. 2012, pp. 1389–1392.
[25] Junning Li, Yonggang Shi, and Arthur W Toga. “Diffusion of fiber orientation distribution
functions with a rotation-induced Riemannian metric”. In: Medical Image Computing and
Computer-Assisted Intervention–MICCAI 2014: 17th International Conference, Boston, MA, USA,
September 14-18, 2014, Proceedings, Part III 17. Springer. 2014, pp. 249–256.
[26] Kesshi M Jordan, Bagrat Amirbekian, Anisha Keshavan, and Roland G Henry. “Cluster confidence
index: A streamline-wise pathway reproducibility metric for diffusion-weighted MRI
tractography”. In: Journal of Neuroimaging 28.1 (2018), pp. 64–69.
[27] Eleftherios Garyfallidis, Matthew Brett, Marta Morgado Correia, Guy B Williams, and
Ian Nimmo-Smith. “Quickbundles, a method for tractography simplification”. In: Frontiers in
neuroscience 6 (2012), p. 175.
[28] Eleftherios Garyfallidis, Marc-Alexandre Côté, Francois Rheault, Jasmeen Sidhu, Janice Hau,
Laurent Petit, David Fortin, Stephen Cunanne, and Maxime Descoteaux. “Recognition of white
matter bundles using local and global streamline-based registration and clustering”. In:
NeuroImage 170 (2018), pp. 283–295.
[29] Lauren J O’Donnell and Carl-Fredrik Westin. “Automatic tractography segmentation using a
high-dimensional white matter atlas”. In: IEEE transactions on medical imaging 26.11 (2007),
pp. 1562–1575.
[30] Robert E Smith, Jacques-Donald Tournier, Fernando Calamante, and Alan Connelly. “SIFT:
Spherical-deconvolution informed filtering of tractograms”. In: Neuroimage 67 (2013),
pp. 298–312.
[31] Jon Haitz Legarreta, Laurent Petit, François Rheault, Guillaume Theaud, Carl Lemaire,
Maxime Descoteaux, and Pierre-Marc Jodoin. “Filtering in tractography using autoencoders
(FINTA)”. In: Medical Image Analysis 72 (2021), p. 102126.
101

[32] Saad Jbabdi, Stamatios N Sotiropoulos, and Timothy E Behrens. “The topographic connectome”.
In: Current opinion in neurobiology 23.2 (2013), pp. 207–215.
[33] Gaurav H Patel, David M Kaplan, and Lawrence H Snyder. “Topographic organization in the
brain: searching for general principles”. In: Trends in cognitive sciences 18.7 (2014), pp. 351–363.
[34] Jean-Philippe Thivierge and Gary F Marcus. “The topographic brain: from neural connectivity to
cognition”. In: Trends in neurosciences 30.6 (2007), pp. 251–259.
[35] Xinyu Nie, Jialiang Ruan, Maria Concepción Garcia Otaduy, Lea Tenenholz Grinberg,
John Ringman, and Yonggang Shi. “Surface-Based Probabilistic Fiber Tracking in Superficial
White Matter”. In: IEEE Transactions on Medical Imaging 43.3 (2024), pp. 1113–1124. doi:
10.1109/TMI.2023.3329451.
[36] Xinyu Nie and Yonggang Shi. “Flow-Based Geometric Interpolation of Fiber Orientation
Distribution Functions”. In: International Conference on Medical Image Computing and
Computer-Assisted Intervention. Springer. 2023, pp. 46–55.
[37] Xinyu Nie and Yonggang Shi. “Topographic filtering of Tractograms as vector field flows”. In:
Medical Image Computing and Computer Assisted Intervention–MICCAI 2019: 22nd International
Conference, Shenzhen, China, October 13–17, 2019, Proceedings, Part III 22. Springer. 2019,
pp. 564–572.
[38] Derek K Jones. “Gaussian modeling of the diffusion signal”. In: diffusion MRI. Elsevier, 2009,
pp. 37–54.
[39] Peter J Basser and Evren Özarslan. “Introduction to diffusion MR”. In: Diffusion MRI (2009),
pp. 2–10.
[40] Edward O Stejskal and John E Tanner. “Spin diffusion measurements: spin echoes in the presence
of a time-dependent field gradient”. In: The journal of chemical physics 42.1 (1965), pp. 288–292.
[41] Kiran K Seunarine and Daniel C Alexander. “Multiple fibers: beyond the diffusion tensor”. In:
Diffusion Mri. Elsevier, 2014, pp. 105–123.
[42] David S Tuch, Timothy G Reese, Mette R Wiegell, Nikos Makris, John W Belliveau, and
Van J Wedeen. “High angular resolution diffusion imaging reveals intravoxel white matter fiber
heterogeneity”. In: Magnetic Resonance in Medicine: An Official Journal of the International Society
for Magnetic Resonance in Medicine 48.4 (2002), pp. 577–582.
[43] Yunmei Chen, Weihong Guo, Qingguo Zeng, Guojun He, Baba Vemuri, and Yijun Liu. “Recovery
of intra-voxel structure from HARD DWI”. In: 2004 2nd IEEE International Symposium on
Biomedical Imaging: Nano to Macro (IEEE Cat No. 04EX821). IEEE. 2004, pp. 1028–1031.
[44] Daniel C Alexander and Gareth J Barker. “Optimal imaging parameters for fiber-orientation
estimation in diffusion MRI”. In: Neuroimage 27.2 (2005), pp. 357–367.
102

[45] Timothy EJ Behrens, Mark W Woolrich, Mark Jenkinson, Heidi Johansen-Berg, Rita G Nunes,
Stuart Clare, Paul M Matthews, J Michael Brady, and Stephen M Smith. “Characterization and
propagation of uncertainty in diffusion-weighted MR imaging”. In: Magnetic Resonance in
Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 50.5
(2003), pp. 1077–1088.
[46] Tim Hosey, Guy Williams, and Richard Ansorge. “Inference of multiple fiber orientations in high
angular resolution diffusion imaging”. In: Magnetic Resonance in Medicine 54.6 (2005),
pp. 1480–1489.
[47] David Solomon Tuch et al. “Diffusion MRI of complex tissue structure”. PhD thesis.
Massachusetts Institute of Technology, 2002.
[48] VJ Wedeen, TG Reese, DS Tuch, MR Weigel, JG Dou, RM Weiskoff, and D Chessler. “Mapping
fiber orientation spectra in cerebral white matter with Fourier-transform diffusion MRI”. In:
Proceedings of the 8th Annual Meeting of ISMRM, Denver. 2000, p. 82.
[49] Van J Wedeen, Patric Hagmann, Wen-Yih Isaac Tseng, Timothy G Reese, and Robert M Weisskoff.
“Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging”. In:
Magnetic resonance in medicine 54.6 (2005), pp. 1377–1386.
[50] David S Tuch. “Q-ball imaging”. In: Magnetic Resonance in Medicine: An Official Journal of the
International Society for Magnetic Resonance in Medicine 52.6 (2004), pp. 1358–1372.
[51] A Anderson and Z Ding. “Sub-voxel measurement of fiber orientation using high angular
resolution diffusion tensor imaging”. In: Book of abstracts: Tenth Annual Meeting of the
International Society for Magnetic Resonance in Medicine. Berkeley, CA: ISMRM. Vol. 10. 2002, p. 440.
[52] Daniel C Alexander. “Maximum entropy spherical deconvolution for diffusion MRI”. In: Biennial
International Conference on Information Processing in Medical Imaging. Springer. 2005, pp. 76–87.
[53] Bing Jian and Baba C Vemuri. “A unified computational framework for deconvolution to
reconstruct multiple fibers from diffusion weighted MRI”. In: IEEE transactions on medical
imaging 26.11 (2007), pp. 1464–1471.
[54] Kendall Atkinson and Weimin Han. Spherical harmonics and approximations on the unit sphere: an
introduction. Vol. 2044. Springer Science & Business Media, 2012.
[55] Peter J Basser, Sinisa Pajevic, Carlo Pierpaoli, Jeffrey Duda, and Akram Aldroubi. “In vivo fiber
tractography using DT-MRI data”. In: Magnetic resonance in medicine 44.4 (2000), pp. 625–632.
[56] Susumu Mori, Barbara J Crain, Vadappuram P Chacko, and Peter CM Van Zijl.
“Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging”.
In: Annals of Neurology: Official Journal of the American Neurological Association and the Child
Neurology Society 45.2 (1999), pp. 265–269.
[57] Timothy EJ Behrens, H Johansen Berg, Saad Jbabdi, Matthew FS Rushworth, and
Mark W Woolrich. “Probabilistic diffusion tractography with multiple fibre orientations: What
can we gain?” In: neuroimage 34.1 (2007), pp. 144–155.
103

[58] Maxime Descoteaux, Rachid Deriche, Thomas R Knosche, and Alfred Anwander. “Deterministic
and probabilistic tractography based on complex fibre orientation distributions”. In: IEEE
transactions on medical imaging 28.2 (2008), pp. 269–286.
[59] J-Donald Tournier, Robert Smith, David Raffelt, Rami Tabbara, Thijs Dhollander,
Maximilian Pietsch, Daan Christiaens, Ben Jeurissen, Chun-Hung Yeh, and Alan Connelly.
“MRtrix3: A fast, flexible and open software framework for medical image processing and
visualisation”. In: Neuroimage 202 (2019), p. 116137.
[60] Manabu Kinoshita, Kei Yamada, Naoya Hashimoto, Amami Kato, Shuichi Izumoto, Takahito Baba,
Motohiko Maruno, Tsunehiko Nishimura, and Toshiki Yoshimine. “Fiber-tracking does not
accurately estimate size of fiber bundle in pathological condition: initial neurosurgical experience
using neuronavigation and subcortical white matter stimulation”. In: Neuroimage 25.2 (2005),
pp. 424–429.
[61] Luigi Ambrosio, Alberto Bressan, Dirk Helbing, Axel Klar, Enrique Zuazua, Luigi Ambrosio, and
Nicola Gigli. “A user’s guide to optimal transport”. In: Modelling and Optimisation of Flows on
Networks: Cetraro, Italy 2009, Editors: Benedetto Piccoli, Michel Rascle (2013), pp. 1–155.
[62] Victor M Panaretos and Yoav Zemel. An invitation to statistics in Wasserstein space. Springer
Nature, 2020.
[63] Francesco Vergani, Sajedha Mahmood, Cristopher M Morris, Patrick Mitchell, and
Stephanie J Forkel. “Intralobar fibres of the occipital lobe: a post mortem dissection study”. In:
Cortex 56 (2014), pp. 145–156.
[64] Marco Catani, Naianna Robertsson, Ahmad Beyh, Vincent Huynh, Francisco de Santiago Requejo,
Henrietta Howells, Rachel LC Barrett, Marco Aiello, Carlo Cavaliere, Tim B Dyrby, et al. “Short
parietal lobe connections of the human and monkey brain”. In: cortex 97 (2017), pp. 339–357.
[65] Owen R Phillips, Shantanu H Joshi, Fabrizio Piras, Maria Donata Orfei, Mariangela Iorio,
Katherine L Narr, David W Shattuck, Carlo Caltagirone, Gianfranco Spalletta, and
Margherita Di Paola. “The superficial white matter in Alzheimer’s disease”. In: Human brain
mapping 37.4 (2016), pp. 1321–1334.
[66] Thomas Veale, Ian B Malone, Teresa Poole, Thomas D Parker, Catherine F Slattery,
Ross W Paterson, Alexander JM Foulkes, David L Thomas, Jonathan M Schott, Hui Zhang, et al.
“Loss and dispersion of superficial white matter in Alzheimer’s disease: a diffusion MRI study”.
In: Brain communications 3.4 (2021), fcab272.
[67] Shuyue Wang, Fan Zhang, Peiyu Huang, Hui Hong, Yeerfan Jiaerken, Xinfeng Yu, Ruiting Zhang,
Qingze Zeng, Yao Zhang, Ron Kikinis, et al. “Superficial white matter microstructure affects
processing speed in cerebral small vessel disease”. In: Human Brain Mapping 43.17 (2022),
pp. 5310–5325.
[68] Arash Nazeri, M Mallar Chakravarty, Tarek K Rajji, Daniel Felsky, David J Rotenberg,
Mikko Mason, Li N Xu, Nancy J Lobaugh, Benoit H Mulsant, and Aristotle N Voineskos.
“Superficial white matter as a novel substrate of age-related cognitive decline”. In: Neurobiology
of aging 36.6 (2015), pp. 2094–2106.
104

[69] Yupeng Wu, Dandan Sun, Yong Wang, Yunjie Wang, and Yibao Wang. “Tracing short connections
of the temporo-parieto-occipital region in the human brain using diffusion spectrum imaging and
fiber dissection”. In: Brain research 1646 (2016), pp. 152–159.
[70] Jin Kyu Gahm and Yonggang Shi. “Surface-based Tracking of U-fibers in the Superficial White
Matter”. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2019: 22nd
International Conference, Shenzhen, China, October 13–17, 2019, Proceedings, Part III 22. Springer.
2019, pp. 538–546.
[71] Bruce Fischl, André Van Der Kouwe, Christophe Destrieux, Eric Halgren, Florent Ségonne,
David H Salat, Evelina Busa, Larry J Seidman, Jill Goldstein, David Kennedy, et al. “Automatically
parcellating the human cerebral cortex”. In: Cerebral cortex 14.1 (2004), pp. 11–22.
[72] Anders M Dale, Bruce Fischl, and Martin I Sereno. “Cortical surface-based analysis: I.
Segmentation and surface reconstruction”. In: Neuroimage 9.2 (1999), pp. 179–194.
[73] Matthew F Glasser, Stamatios N Sotiropoulos, J Anthony Wilson, Timothy S Coalson,
Bruce Fischl, Jesper L Andersson, Junqian Xu, Saad Jbabdi, Matthew Webster,
Jonathan R Polimeni, et al. “The minimal preprocessing pipelines for the Human Connectome
Project”. In: Neuroimage 80 (2013), pp. 105–124.
[74] Nicholas J Tustison, Philip A Cook, Arno Klein, Gang Song, Sandhitsu R Das, Jeffrey T Duda,
Benjamin M Kandel, Niels van Strien, James R Stone, James C Gee, et al. “Large-scale evaluation
of ANTs and FreeSurfer cortical thickness measurements”. In: Neuroimage 99 (2014), pp. 166–179.
[75] Eugene Paul Wigner and U Fano. “Group theory and its application to the quantum mechanics of
atomic spectra”. In: American Journal of Physics 28.4 (1960), pp. 408–409.
[76] Shan-Tao Lai, P Palting, and Ying-Nan Chiu. “On the closed form of Wigner rotation matrix
elements”. In: Journal of Mathematical Chemistry 19 (1996), pp. 131–145.
[77] G Aubert. “An alternative to Wigner d-matrices for rotating real spherical harmonics”. In: AIP
Advances 3.6 (2013).
[78] Stephanie Alexander. “Michael Spivak, A comprehensive introduction to differential geometry”.
In: (1978).
[79] Fernando De Goes, Mathieu Desbrun, and Yiying Tong. “Vector field processing on triangle
meshes”. In: ACM SIGGRAPH 2016 Courses. 2016, pp. 1–49.
[80] Eugene Zhang, Konstantin Mischaikow, and Greg Turk. “Vector field design on surfaces”. In:
ACM Transactions on Graphics (ToG) 25.4 (2006), pp. 1294–1326.
[81] Yonggang Shi, Paul M Thompson, Ivo Dinov, and Arthur W Toga. “Hamilton–Jacobi skeleton on
cortical surfaces”. In: IEEE transactions on medical imaging 27.5 (2008), pp. 664–673.
[82] TJ Biega, RR Lonser, and JA Butman. “Differential cortical thickness across the central sulcus: a
method for identifying the central sulcus in the presence of mass effect and vasogenic edema”. In:
American journal of neuroradiology 27.7 (2006), pp. 1450–1453.
105

[83] Tarek A Yousry, Urs D Schmid, Andre G Jassoy, Dorothea Schmidt, Wilhelm E Eisner,
Hanns-Jurgen Reulen, Maximilian F Reiser, and Joseph Lissner. “Topography of the cortical
motor hand area: prospective study with functional MR imaging and direct motor mapping at
surgery.” In: Radiology 195.1 (1995), pp. 23–29.
[84] Wolfgang Grodd, Ernst Hülsmann, Martin Lotze, Dirk Wildgruber, and Michael Erb.
“Sensorimotor mapping of the human cerebellum: fMRI evidence of somatotopic organization”.
In: Human brain mapping 13.2 (2001), pp. 55–73.
[85] Dogu Baran Aydogan and Yonggang Shi. “Tracking and validation techniques for topographically
organized tractography”. In: NeuroImage 181 (2018), pp. 64–84.
[86] Junyan Wang, Dogu Baran Aydogan, Rohit Varma, Arthur W Toga, and Yonggang Shi. “Modeling
topographic regularity in structural brain connectivity with application to tractogram filtering”.
In: NeuroImage 183 (2018), pp. 87–98.
[87] BPRSRG Drayer, P Burger, R Darwin, S Riederer, R Herfkens, and GA Johnson. “MRI of brain
iron”. In: American Journal of Roentgenology 147.1 (1986), pp. 103–110.
[88] Evgeniya Kirilina, Saskia Helbling, Markus Morawski, Kerrin Pine, Katja Reimann,
Steffen Jankuhn, Juliane Dinse, Andreas Deistung, Jürgen R Reichenbach, Robert Trampel, et al.
“Superficial white matter imaging: Contrast mechanisms and whole-brain in vivo mapping”. In:
Science Advances 6.41 (2020), eaaz9281.
[89] Paul A Yushkevich, Joseph Piven, Heather Cody Hazlett, Rachel Gimpel Smith, Sean Ho,
James C Gee, and Guido Gerig. “User-guided 3D active contour segmentation of anatomical
structures: significantly improved efficiency and reliability”. In: Neuroimage 31.3 (2006),
pp. 1116–1128.
[90] Jill Murrell, Bernardino Ghetti, Elizabeth Cochran, Miguel Angel Macias-Islas, Luis Medina,
Arousiak Varpetian, Jeffrey L Cummings, Mario F Mendez, Claudia Kawas, Helena Chui, et al.
“The A431E mutation in PSEN1 causing familial Alzheimer’s disease originating in Jalisco State,
Mexico: an additional fifteen families”. In: Neurogenetics 7 (2006), pp. 277–279.
[91] Brian A Gordon, Tyler M Blazey, Jon Christensen, Aylin Dincer, Shaney Flores, Sarah Keefe,
Charles Chen, Yi Su, Eric M McDade, Guoqiao Wang, et al. “Tau PET in autosomal dominant
Alzheimer’s disease: relationship with cognition, dementia and other biomarkers”. In: Brain 142.4
(2019), pp. 1063–1076.
[92] Charles P Hughes, Leonard Berg, Warren Danziger, Lawrence A Coben, and Ronald L Martin. “A
new clinical scale for the staging of dementia”. In: The British journal of psychiatry 140.6 (1982),
pp. 566–572.
[93] David Raffelt, J-Donald Tournier, Stephen Rose, Gerard R Ridgway, Robert Henderson,
Stuart Crozier, Olivier Salvado, and Alan Connelly. “Apparent fibre density: a novel measure for
the analysis of diffusion-weighted magnetic resonance images”. In: Neuroimage 59.4 (2012),
pp. 3976–3994.
106

[94] John M Ringman, Naghmeh Dorrani, Sara Gutiérrez Fernández, Rebecca Signer,
Julian Martinez-Agosto, Hane Lee, Emilie D Douine, Yuchuan Qiao, Yonggang Shi, Lina D’Orazio,
et al. “Characterization of spastic paraplegia in a family with a novel PSEN1 mutation”. In: Brain
Communications 5.2 (2023), fcad030.
[95] Douglas N Greve, Claus Svarer, Patrick M Fisher, Ling Feng, Adam E Hansen, William Baare,
Bruce Rosen, Bruce Fischl, and Gitte M Knudsen. “Cortical surface-based analysis reduces bias
and variance in kinetic modeling of brain PET data”. In: Neuroimage 92 (2014), pp. 225–236.
[96] Colin Reveley, Anil K Seth, Carlo Pierpaoli, Afonso C Silva, David Yu, Richard C Saunders,
David A Leopold, and Frank Q Ye. “Superficial white matter fiber systems impede detection of
long-range cortical connections in diffusion MR tractography”. In: Proceedings of the National
Academy of Sciences 112.21 (2015), E2820–E2828.
[97] Brian A Wandell. “Clarifying human white matter”. In: Annual review of neuroscience 39 (2016),
pp. 103–128.
[98] Cibu Thomas, Frank Q Ye, M Okan Irfanoglu, Pooja Modi, Kadharbatcha S Saleem,
David A Leopold, and Carlo Pierpaoli. “Anatomical accuracy of brain connections derived from
diffusion MRI tractography is inherently limited”. In: Proceedings of the National Academy of
Sciences 111.46 (2014), pp. 16574–16579.
[99] Dogu Baran Aydogan and Yonggang Shi. “Parallel transport tractography”. In: IEEE transactions
on medical imaging 40.2 (2020), pp. 635–647.
[100] Ye Wu, Yoonmi Hong, Yuanjing Feng, Dinggang Shen, and Pew-Thian Yap. “Mitigating gyral bias
in cortical tractography via asymmetric fiber orientation distributions”. In: Medical image analysis
59 (2020), p. 101543.
[101] Matteo Bastiani, Michiel Cottaar, Krikor Dikranian, Aurobrata Ghosh, Hui Zhang,
Daniel C Alexander, Timothy E Behrens, Saad Jbabdi, and Stamatios N Sotiropoulos. “Improved
tractography using asymmetric fibre orientation distributions”. In: Neuroimage 158 (2017),
pp. 205–218.
[102] Tim McGraw, Baba Vemuri, Robert Yezierski, and Thomas Mareci. “Segmentation of high angular
resolution diffusion MRI modeled as a field of von Mises-Fisher mixtures”. In: Computer
Vision–ECCV 2006: 9th European Conference on Computer Vision, Graz, Austria, May 7-13, 2006,
Proceedings, Part III 9. Springer. 2006, pp. 463–475.
[103] Maxime Descoteaux and Rachid Deriche. “High angular resolution diffusion MRI segmentation
using region-based statistical surface evolution”. In: Journal of Mathematical Imaging and Vision
33.2 (2009), pp. 239–252.
[104] Ming-Chang Chiang, Marina Barysheva, Agatha D Lee, Sarah Madsen, Andrea D Klunder,
Arthur W Toga, Katie L McMahon, Greig I De Zubicaray, Matthew Meredith, Margaret J Wright,
et al. “Brain fiber architecture, genetics, and intelligence: a high angular resolution diffusion
imaging (HARDI) study”. In: Medical Image Computing and Computer-Assisted
Intervention–MICCAI 2008: 11th International Conference, New York, NY, USA, September 6-10, 2008,
Proceedings, Part I 11. Springer. 2008, pp. 1060–1067.
107

[105] Charles Zheng, Franco Pestilli, and Ariel Rokem. “Quantifying error in estimates of human brain
fiber directions using Earth Mover’s Distance”. In: arXiv preprint arXiv:1411.5271 (2014).
[106] Gabriel Peyré, Marco Cuturi, et al. “Computational optimal transport: With applications to data
science”. In: Foundations and Trends® in Machine Learning 11.5-6 (2019), pp. 355–607.
[107] Fang-Cheng Yeh, Van Jay Wedeen, and Wen-Yih Isaac Tseng. “Estimation of fiber orientation and
spin density distribution by diffusion deconvolution”. In: Neuroimage 55.3 (2011), pp. 1054–1062.
[108] Saad Jbabdi, Stamatios N Sotiropoulos, Suzanne N Haber, David C Van Essen, and
Timothy E Behrens. “Measuring macroscopic brain connections in vivo”. In: Nature neuroscience
18.11 (2015), pp. 1546–1555.
[109] U Ebeling and H -J Reulen. “Neurosurgical topography of the optic radiation in the temporal
lobe”. In: Acta neurochirurgica 92 (1988), pp. 29–36.
[110] Gordon Holmes and WT Lister. “Disturbances of vision from cerebral lesions, with special
reference to the cortical representation of the macula”. In: Brain 39.1-2 (1916), pp. 34–73.
[111] Roger B Tootell, Eugene Switkes, Martin S Silverman, and Susan L Hamilton. “Functional
anatomy of macaque striate cortex. II. Retinotopic organization”. In: Journal of neuroscience 8.5
(1988), pp. 1531–1568.
[112] Stephen A Engel, David E Rumelhart, Brian A Wandell, Adrian T Lee, Gary H Glover,
Eduardo-Jose Chichilnisky, and Michael N Shadlen. “fMRI of human visual cortex.” In: Nature
(1994).
[113] Deniz Bilecen, Klaus Scheffler, Nena Schmid, Kurt Tschopp, and Joachim Seelig. “Tonotopic
organization of the human auditory cortex as detected by BOLD-FMRI”. In: Hearing research
126.1-2 (1998), pp. 19–27.
[114] Suzanne N Haber and Brian Knutson. “The reward circuit: linking primate anatomy and human
imaging”. In: Neuropsychopharmacology 35.1 (2010), pp. 4–26.
[115] Judea Pearl. Probabilistic reasoning in intelligent systems: networks of plausible inference. Morgan
kaufmann, 1988.
[116] Mathias Wahl, Birgit Lauterbach-Soon, Elke Hattingen, Patrick Jung, Oliver Singer, Steffen Volz,
Johannes C Klein, Helmuth Steinmetz, and Ulf Ziemann. “Human motor corpus callosum:
topography, somatotopy, and link between microstructure and function”. In: Journal of
neuroscience 27.45 (2007), pp. 12132–12138.
[117] L Conlon. “Differentiable Manifolds Birkhäuser Advanced Texts”. In: Basler Lehrbacher (1st ed.
1994 or 2nd ().
[118] Amir Vaxman, Marcel Campen, Olga Diamanti, Daniele Panozzo, David Bommes,
Klaus Hildebrandt, and Mirela Ben-Chen. “Directional field synthesis, design, and processing”. In:
Computer graphics forum. Vol. 35. 2. Wiley Online Library. 2016, pp. 545–572.
108

[119] Richard Szeliski, Ramin Zabih, Daniel Scharstein, Olga Veksler, Vladimir Kolmogorov,
Aseem Agarwala, Marshall Tappen, and Carsten Rother. “A comparative study of energy
minimization methods for markov random fields with smoothness-based priors”. In: IEEE
transactions on pattern analysis and machine intelligence 30.6 (2008), pp. 1068–1080.
[120] Alexandra Kammen, Meng Law, Bosco S Tjan, Arthur W Toga, and Yonggang Shi. “Automated
retinofugal visual pathway reconstruction with multi-shell HARDI and FOD-based analysis”. In:
Neuroimage 125 (2016), pp. 767–779.
109

Abstract (if available)

Abstract Diffusion MRI (dMRI) is an essential non-invasive imaging tool that allows us to study the microstructure and connectivity of the human brain {\it in vivo}. Tractography, a critical technique for reconstructing fiber connections based on dMRI, has become a critical tool in brain connectome research, but it is error-prone and suffers from various limitations. In this thesis, we will develop novel computational techniques to overcome several limitations in current tractography techniques and advance the state-of-the-art in human connectome research.

A critical challenge for tractography in connectome research is reconstructing fiber connections in the superficial white matter (SWM) distributed beneath the cortical surface. Indeed, complex structures and the high curvature of the fibers in SWM lead to a false-negative reconstruction in the conventional tractography algorithms. Another weakness of tractography is that it generates enormous anatomically invalid connections that reduce the accuracy of the connectivity analysis, even with advanced mathematical models such as fiber orientation distribution (FOD). In this thesis, we propose a novel surface-based tractography algorithm that successfully reconstructs the fiber connectivity in the SWM. On the other hand, we develop geometric frameworks to process and analyze the FOD functions and provide novel interpolations of the FOD that help reduce the implausible propagations in tractography. Finally, we proposed a post-processing vector-flow-based filter to reduce the unrealistic streamlines in tractography.
To overcome the difficulties of traditional volume-based tractography in SWM, we propose a novel surface-based framework for the probabilistic tracking of fibers on the triangular mesh representation of the SWM. By deriving a closed-form solution to transform the spherical harmonics (SPHARM) coefficients of 3D fiber orientation distributions (FODs) to local coordinate systems on each triangle, we develop a novel approach to project the FODs onto the tangent space of the SWM. After that, we utilize parallel transport to realize the intrinsic propagation of streamlines on SWM following probabilistically sampled fiber directions. Our intrinsic and surface-based method eliminates the need to perform the necessary but challenging sharp turns in 3D compared with conventional volume-based tractography methods. We performed quantitative comparisons to demonstrate the proposed algorithm can more effectively reconstruct the fiber connections in SWM than previous methods.

Interpolation plays a vital role in the tractography algorithms since the step size of the tracking is much smaller than the resolution of the dMRI data. Previous interpolation methods of FODs often generate false artificial information, leading to anatomically incorrect propagations in tractography. We propose a flow-rotation-based interpolation framework based on the local geometry of the fiber bundles. Our approach decomposes each FOD function with multiple peak lobes into single-peak components, and neighboring single-peak components consistent in direction are modeled as vector field flows. The vector field flow representation enables geometrically consistent interpolation of FOD functions in a rotation-based framework. Experiments demonstrate that our method produces anatomically meaningful FOD interpolation and improves state-of-the-art tractography.

We developed an advanced geometric framework to improve the flow-rotation-based framework and put the essential concepts in more efficient and mathematically rigorous frameworks. The decomposed single-peak FOD functions are embedded and processed in Wasserstein space, a well-defined metric space with excellent properties. We define a new metric to compare and analyze the FOD functions, and the interpolation is computed as the Barycenter of the new metric. For efficient computation, we propose a Fast Approximation, which is a generalization of the flow-rotation-based interpolation, for the Barycenter. In experiments, we use synthetic and real datasets to demonstrate that the new geometric framework can handle FOD functions representing complicated fiber geometry and provide accurate interpolations that improve tractography performance.

To reduce the invalid streamlines in the post-processing of tractography, we proposed a filter based on the vector flow representation, which is used to model the local fiber geometry in flow-rotation-based interpolation. The vector flow representation is consistent with the topographic principle in brain connections, and it can be naturally extended to the global representation of the topographically regular fiber bundles. Thus, we developed a framework that models the streamlines of fiber bundles as vector field flows. We compute a fiber bundle's principal vector field (PVF) by solving a Markov Random Field problem based on the FOD peaks. We propose a metric named Vector Flow Deviation (VFD) for tractography by quantifying the consistency between the PVF and the tractography streamlines and applying the VFD to filter out topographically irregular streamlines. In the experiments, we visually and quantitatively show that our method successfully removes the fiber bundles' outliers while significantly improving the topographical regularity.

In this thesis, we developed four novel computational techniques for connectome analysis: a surface-based tractography for reconstructing SWM connections, a flow-rotation-based FOD interpolation that overcomes the difficulties in commonly used interpolations, a geometric framework to process and analyze the FOD functions efficiently, and a filter that reduces the invalid connections in tractography. All four techniques overcome critical difficulties in tractography and connectome analysis and outperform state-of-the-art techniques in visual and quantitative comparisons. We develop software tools for the proposed algorithms and will distribute them publicly to the research community.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Methods for improving reliability and consistency in diffusion MRI analysis

PDF

Correction, coregistration and connectivity analysis of multi-contrast brain MRI

PDF

Improved brain dynamic contrast enhanced MRI using model-based reconstruction

PDF

Novel theoretical characterization and optimization of experimental efficiency for diffusion MRI

PDF

High-performance distributed computing techniques for wireless IoT and connected vehicle systems

PDF

New methods for carotid MRI

PDF

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

PDF

Shape, pose, and connectivity in subcortical networks across the human lifespan

PDF

Functional connectivity analysis and network identification in the human brain

PDF

Signal processing methods for brain connectivity

PDF

3D vessel mapping techniques for retina and brain as an early imaging biomarker for small vessel diseases

PDF

Advanced techniques for stereoscopic image rectification and quality assessment

PDF

Graph-based models and transforms for signal/data processing with applications to video coding

PDF

Machine learning techniques for outdoor and indoor layout estimation

PDF

Contributions to structural and functional retinal imaging via Fourier domain optical coherence tomography

PDF

Functional real-time MRI of the upper airway

PDF

Architecture design and algorithmic optimizations for accelerating graph analytics on FPGA

PDF

Detection and decoding of cognitive states from neural activity to enable a performance-improving brain-computer interface

PDF

Advanced machine learning techniques for video, social and biomedical data analytics

PDF

Neighborhood and graph constructions using non-negative kernel regression (NNK)

Asset Metadata

Creator Nie, Xinyu (author)

Core Title Novel computational techniques for connectome analysis based on diffusion MRI

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Degree Conferral Date 2024-05

Publication Date 05/17/2024

Defense Date 04/23/2024

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

Tag brain connectivity,diffusion MRI,fiber orientation distribution,OAI-PMH Harvest,spherical harmonic,superficial white matter,tractography,tractography filtering,Wasserstein distance

Format theses (aat)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Shi, Yonggang (committee chair)

Creator Email 1993nxy@gmail.com,xnie@usc.edu

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

Unique identifier UC113939761

Identifier etd-NieXinyu-12914.pdf (filename)

Legacy Identifier etd-NieXinyu-12914

Document Type Dissertation

Format theses (aat)

Rights Nie, Xinyu

Internet Media Type application/pdf

Type texts

Source 20240517-usctheses-batch-1151 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu