# Manifold-aware Synthesis of High-resolution Diffusion from Structural
Imaging

Benoit Anctil-Robitaille  Antoine Théberge  Pierre-Marc Jodoin  Maxime
Descoteaux  Christian Desrosiers  Hervé Lombaert  Department of Computer and
Software Engineering, ETS Montreal, Canada  Department of Computer Science,
Sherbrooke University, Canada

###### Abstract

The physical and clinical constraints surrounding diffusion-weighted imaging
(DWI) often limit the spatial resolution of the produced images to voxels up
to 8 times larger than those of T1w images. Thus, the detailed information
contained in T1w images could help in the synthesis of diffusion images in
higher resolution. However, the non-Euclidean nature of diffusion imaging
hinders current deep generative models from synthesizing physically plausible
images. In this work, we propose the first Riemannian network architecture for
the direct generation of diffusion tensors (DT) and diffusion orientation
distribution functions (dODFs) from high-resolution T1w images. Our
integration of the Log-Euclidean Metric into a learning objective guarantees,
unlike standard Euclidean networks, the mathematically-valid synthesis of
diffusion. Furthermore, our approach improves the fractional anisotropy mean
squared error (FA MSE) between the synthesized diffusion and the ground-truth
by more than 23% and the cosine similarity between principal directions by
almost 5% when compared to our baselines. We validate our generated diffusion
by comparing the resulting tractograms to our expected real data. We observe
similar fiber bundles with streamlines having less than 3% difference in
length, less than 1% difference in volume, and a visually close shape. While
our method is able to generate high-resolution diffusion images from
structural inputs in less than 15 seconds, we acknowledge and discuss the
limits of diffusion inference solely relying on T1w images. Our results
nonetheless suggest a relationship between the high-level geometry of the
brain and the overall white matter architecture.

###### keywords:

Diffusion synthesis , Manifold-valued data learning, Riemannian geometry

##  1 Introduction

Diffusion MRI is of crucial importance in multiple challenging tasks,
including the diagnosis of complex cognitive disorders [31, 30, 39], the study
of neurodegenerative diseases [24, 19] and neurosurgical planning [10].
Nonetheless, diffusion-weighted imaging (DWI) suffers from a low signal-to-
noise ratio (SNR) and a poor spatial resolution arising from physical and
clinical limitations such as the use of echo-planar imaging (EPI) and limited
patient scanning time. Indeed, the induced trade-off between image resolution,
SNR and imaging time [41] in the acquisition of DWI often results in images
with voxel size up to 8 times larger than other common modalities such as
structural T1w images, e.g. 1 mm iso. for DWI vs 2 mm iso. for T1w [45].
Hence, it has been shown that the increased presence of partial volume effect
(PVE) in DW images, due to their low spatial resolution, impairs their
subsequent analysis [40, 2].

These limitations have triggered the development of post-processing methods
that aim to improve the spatial resolution of low-resolution diffusion
volumes. Tackling the estimation of high-resolution (HR) diffusion from low-
resolution (LR) images was first explored with interpolation-based methods [5,
54, 13]. They resample existing images to a higher-resolution grid and is
nowadays a default step in diffusion MRI processing tools such as in
TractoFlow [49] and MRtrix3 [51]. Although fine anatomical details can be
enhanced with such approach, interpolations will always be limited by the
inherent coarseness of the original diffusion data as they exclusively rely on
intra-image information.

Machine learning offers an effective way to leverage the rich information
contained in HR images for the synthesis of diffusion imaging in the same
resolution, thus going beyond interpolation. In [3], a fully supervised image
quality transfer (IQT) framework using random forests is proposed to learn a
non-linear mapping between paired low-quality and high-quality diffusion data.
Similarly, in [14], a supervised 2D SRCNN is used for the same objective. The
authors demonstrate that learning a mapping from an LR input to its HR version
not only helps recovering anatomical details better than interpolation, but
can also help in downstream tasks such as tractography. However, such
approaches rely on limited high-resolution diffusion data which are costly and
challenging to acquire. Moreover, the methods in [3, 14] have only been tested
on small datasets comprising a maximum of 23 subjects and 3 subjects
respectively.

In parallel, deep neural networks offer unsupervised learning techniques that
only require few paired training samples to train specific synthesis tasks.
More particularly, Generative Adversarial Networks (GANs) [21] have been
successfully used for the synthesis of missing modalities [11], image-to-image
translation [58, 35] and image super-resolution [43], just to name a few.
Therefore, deep generative models could be a key solution for the synthesis of
high-resolution diffusion from unpaired images expressing a higher level of
structural details such as T1w images.

While the synthesis of raw DWI signals with a proper angular resolution is an
extremely resource-intensive task, current deep learning architectures
struggle to generate plausible diffusion reconstruction schemes, notably
Diffusion Tensor (DT) or Orientation Distribution Functions (ODFs) because of
their non-Euclidean nature [26]. Indeed, each voxel of a DT image lies on a
Riemannian manifold of symmetric positive definite (SPD) 3$\times$3 matrices
[5], and ODFs can be represented as points on an n-Sphere manifold
$\mathbb{S}^{n}$ [8]. In the context of image synthesis, the inability of
networks to capture the underlying non-linear Riemannian manifold geometry of
the data results in the generation of implausible images that miss the
important mathematical properties of diffusion imaging [26]. Consequently, the
limitations of current deep neural networks (DNN) have impeded the development
of generative models in diffusion imaging, which have been mostly restricted
to the synthesis of DT scalar maps such as Fractional Anisotropy (FA) and Mean
Diffusivity (MD).

Figure 1: An example of the input and output of our network. The HR T1w input
image (first column), the real LR diffusion (second column), the generated HR
diffusion (third column) and the recovered HR T1w image (last column) of a
test subject.

###  1.1 Related Work

##### Structural-to-Diffusion Synthesis

Amidst the literature, [22] studies the generation of diffusion-derived scalar
maps from downsampled structural images. To do so, the authors use a CycleGAN
to learn the intermodal relationships between T1w images and FA/MD maps and
successfully translate one to another. They demonstrate that structural images
share sufficient information with the diffusion anisotropy of tissues to
synthesize plausible 2D FA and MD slices. Similarly, in [33], a Self-attention
Conditional GAN (SC-GAN) is used to generate FA and MD maps from different
input modalities including structural T1w images. Their results indicate that
both the 3D contextual information and the adversarial objective are important
building blocks for the synthesis of diffusion data. In [57], dual GANs with a
Markovian discriminator [36] are employed for the harmonization of inter-site
DT-derived metrics. Finally, in [46] functional MRI in combination with
structural T1w inputs are fed to a CNN network to generate DT. This body of
work demonstrates the potential of generative models for the structural-to-
diffusion synthesis of imaging data. Nevertheless, they remain limited [22,
33, 46] in not exploiting the high-resolution information contained in the
structural images to their full extent by either considering downsampled
version of the T1w inputs or 2D slices with limited context. In addition, even
though diffusion scalar maps are clinically useful, they mostly ignore fiber
orientations and are of limited interest for tasks such as tractography or
connectome visualization. With regards to the generated tensors in [46], the
authors provide no guarantee on their mathematical validity, such as symmetric
positive definiteness, nor on their usability in a downstream task such as
tractography.

##### Manifold-Valued Data Learning

Deep learning models are well suited to model data lying in an Euclidean
vector space. However, the Euclidean operations from which they are built
upon, e.g., convolutions or pooling, are not well defined on curved manifolds.
Moreover, the application of Euclidean geometry to manifold-valued data, such
as DT and ODF, has well-documented side effects [5]. Consequently, studies
that use neural networks for the accurate processing of data on Riemannian
manifolds have started to emerge [7, 6]. However, these works require
substantial modifications of known deep learning models and call for further
investigation in a broader set of scenarios.

Another avenue for the processing of manifold-valued data resides in the
design of computationally efficient Riemannian metrics. To that purpose, [5]
proposed a Log-Euclidean metric to process data lying on the
$\mathcal{S}_{++}^{*}$ manifold with applications to diffusion tensors. With
the help of the $\log$ and $\exp$ maps defined in [5], one can process tensors
using Euclidean operations and guarantee that the processed tensors keep their
SPD properties. Likewise, a Log-Euclidean framework has also been proposed in
[8] for the computation of orientation distribution function and applied to
diffusion ODF. These two frameworks, combined with the matrix backpropagation
of spectral layers presented in [27], constitutes the fundamentals of the
following manifold-valued data learning approaches. For instance, in [25], the
authors have integrated the Log-Euclidean metric into their deep learning
model called SPDNet to learn compact and discriminative SPD matrices. Although
SPDNet offers a way to learn data on $\mathcal{S}_{++}^{*}$, it has not been
designed for spatially organized and volumetric SPD matrices learning as in
DT.

More recently, [26] proposed a Wasserstein GAN (WGAN) [4] leveraging the Log-
Euclidean metric to synthesize plausible DT, among other manifold-valued data
type. To ensure the validity of the generated data, the authors project the
output of their generator network to a vector space using the $\log$ map from
[5] prior to the discriminator assessment. The $\exp$ operation is then used
on the synthesized output to recover valid DT. Despite its ability to generate
mathematically valid diffusion, the model in [26] outputs images that are not
conditioned by any real subject specific information (e.g., a T1w image) and,
thus, are less clinically valuable. In addition, this prior work only focuses
on the generation of DT in 2D, which once again limits the value of the
generated data.

###  1.2 Contributions

This work proposes a novel deep learning architecture that leverages the
detailed information of high-resolution structural images to guide the
synthesis of DT and ODF in the same high-resolution space. Based on the
CycleGAN architecture [58], our solution exploits the inherent cross-modality
representations of structural and diffusion images to learn functions that map
one to the other in a weakly supervised manner. To do so, we address the
current limitations of deep learning models built upon Euclidean operations by
integrating a Riemannian framework, namely the Log-Euclidean framework, for
statistical computations on DT and ODF directly into the model. Such
Riemannian framework within the learning procedure of the network enforces a
valid synthesis of diffusion data that lies on a desired Riemannian manifold.
By constraining the generated diffusion to lie on the Riemannian manifold of
$3\times 3$ symmetric positive definite matrices ($\mathcal{S}_{++}^{3}$) for
DT and on the n-Sphere $\mathbb{S}^{\mathbf{n}}$ manifold for ODF, we
guarantee the mathematical coherence of the solution. As opposed to existing
deep generative models that focused on the generation of diffusion derived
scalar maps [36, 22, 33], our model outputs complete diffusion schemes, here
the $3\times 3$ DT and the spherical harmonic coefficients of ODF. This
important difference allows one to perform tractography, tractogram
visualization and fiber bundle segmentation in addition to scalar maps
computation directly from our network output while only requiring a T1w image
as input.

Specifically, our contributions are as follows:

  * 1.

The first Riemannian network for the cycle-consistent mapping between real-
valued images and data lying on the $\mathcal{S}_{++}^{3}$ and the
$\mathbb{S}^{\mathbf{n}}$ manifolds;

  * 2.

The first deep learning model for the guided super-resolution of DT and ODF
from unpaired high-resolution structural images and limited priors;

  * 3.

A comprehensive analysis of synthesized diffusion imaging including the
evaluation of full-valued diffusion data, scalar maps and tractography.

##  2 Method

Figure 2: The forward cycle of our manifold-aware CycleGAN. $G_{Y}$ generates
high resolution DT on the $\mathcal{S}_{++}^{3}$ manifold using the
$\log_{\mathrm{Id}}$ and $\exp_{\mathrm{Id}}$ mapping. $D_{Y}$ assesses the
generated images quality and provides feedback to $G_{Y}$. $G_{X}$ tries to
reconstruct the original T1w images from
$\log_{\mathrm{Id}}(\exp_{\mathrm{Id}}(G_{Y}(\mathbf{x})))$.

Let $X$ be the real-valued domain of structural images (e.g., T1w) and $Y$ be
the manifold-valued domain of diffusion images (e.g., DT or ODF). We aim at
learning mapping functions $G_{Y}:X_{\mathrm{HR}}\mapsto Y_{\mathrm{HR}}$ and
$G_{X}:Y_{\mathrm{HR}}\mapsto X_{\mathrm{HR}}$ that translate high-resolution
T1w images to high-resolution diffusion images and the other way around.
Nevertheless, we mainly have access to unpaired HR T1w images and LR diffusion
data that differ in terms of subject, resolution and nature. To address this
problem, we propose a Manifold-Aware CycleGAN (MA-CycleGAN) architecture (see
Figure 2) that inherently handles both the domain translation and the super-
resolution of diffusion, while accounting for the Riemannian geometry of the
data.

We train our network with unpaired training samples
$\\{\mathbf{x}_{i}\\}^{N}_{i=1}$ where $\mathbf{x}_{i}\in X_{\mathrm{HR}}$ is
a 3D structural image, and $\\{\mathbf{y}_{j}\\}^{M}_{j=1}$ where
$\mathbf{y}_{j}\in Y_{\mathrm{LR}}$ is a diffusion image in a lower spatial
resolution. The mathematical validity of the generated DT and ODF is ensured
by projecting every synthesized voxel on the tangent plane at their isotropic
counterpart (i.e., the $3\times 3$ identity matrix for DT and the uniform
distribution for ODF). To do so, the exponential and logarithm maps borrowed
from their underlying Riemannian manifolds geometry are used. Two
discriminators $D_{X}$ and $D_{Y}$ evaluate the synthesized HR T1w images
$G_{X}(\mathbf{y})$ and downsampled HR diffusion images, using a learned
residual function $\mathcal{F}$ as follows $\mathcal{F}(G_{Y}(\mathbf{x}))$,
with regards to their real data distribution
$G_{X}(\mathbf{y})\sim{\mathbb{P}_{X_{\mathrm{HR}}}}$ and
$\mathcal{F}(G_{Y}(\mathbf{x}))\sim{\mathbb{P}_{\log(Y_{\mathrm{LR}})}}$. By
combining pixel-wise reconstruction losses and higher-level adversarial
feedback in a single objective, our MA-CycleGAN is able to exploit the local
tissue information and global geometry of high-resolution structural images to
produce realistic and usable sharp diffusion data.

In the following sections, we detail the Riemannian frameworks embedded in our
architecture for DT and ODF learning. Moreover, we frame our cycle-consistent
and adversarial objectives incorporating both the $\exp$ and $\log$ maps of
the aforementioned Riemannian frameworks and an up-and-down sampling strategy.
Then, we present our anisotropy-based attention mechanism that helps the
network to focus on meaningful fiber tracts information.

###  2.1 Riemannian Framework for Diffusion Tensors Learning

Diffusion tensors are $3\times 3$ symmetric positive matrices $\mathbf{M}$
that can be decomposed in three real and positive eigenvalues
$\boldsymbol{\Upsigma}$ and three corresponding eigenvectors $\mathbf{U}$
using eigendecomposition such that
$\mathbf{M}=\mathbf{U}\boldsymbol{\Upsigma}{\mathbf{U}}^{\top}$. The
eigenvector $\mathbf{u_{1}}$ associated with the largest eigenvalue of
$\mathbf{M}$ represents the principal direction of diffusion and aligns with
the underlying fibers population. DT-derived metrics, describing the shape of
the tensor, are computed from the positive eigenvalues
$\lambda_{1}>\lambda_{2}>\lambda_{3}\in\boldsymbol{\Upsigma}$. One of the most
important DT-derived metric, fractional anisotropy, measures how far the shape
of the diffusion tensor is from a sphere (i.e., how anisotropic the diffusion
is). This metric is computed as follows:

$\begin{split}\mathrm{FA}=\sqrt{\frac{1}{2}}\frac{\sqrt{(\lambda_{1}-\lambda_{2})^{2}+(\lambda_{2}-\lambda_{3})^{2}+(\lambda_{1}-\lambda_{3})^{2}}}{\sqrt{(\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2})}}.\end{split}$
(1)

Because of their SPD properties, tensors lie on a non-linear manifold denoted
as

$\mathcal{S}_{++}^{3}=\big{\\{}\mathbf{M}\in{}^{3\times
3},\,\mathbf{M}={\mathbf{M}}^{\top},\,{\mathbf{x}}^{\top}\mathbf{M}\mathbf{x}>0,\,\forall\mathbf{x}\in{}^{3},\,\left\|\mathbf{x}\right\|_{2}>0\big{\\}}\\\
$ (2)

$\mathcal{S}_{++}^{3}$ is not a vector space (i.e., the linear combination of
two elements in $\mathcal{S}_{++}^{3}$ may lie outside this space), thus using
standard Euclidean operations to process statistics on diffusion tensors can
lead to undesirable effects like the well-documented swelling effect in [5].
Moreover, not considering the $\mathcal{S}_{++}^{3}$ manifold while
synthesizing DT with deep neural network can lead to the generation of non-SPD
tensors as can be seen in Table 3. Such tensors are physically incorrect and
must be avoided. To accurately process DT, [5] proposed a Log-Euclidean metric
that will be presented in the following sections.

Figure 3: We use the Log-Euclidean metric to project the manifold-valued
diffusion data to a tangent plane before processing them with Euclidean
operations. On the left, the projection of diffusion tensors from the
$\mathcal{S}_{++}^{3}$ manifold to the tangent plane at the isotropic tensor
$\mathcal{T}_{\mathrm{Id}}$. On the right, the projection of square-root re-
parameterized ODFs from the $\mathbb{S}^{n}$ manifold to the tangent plane at
the uniform distribution $\mathcal{T}_{\mathbf{u}}$.

####  2.1.1 Log-Euclidean Metric

Diffusion tensor matrices are well defined in the Log-Euclidean metric, where
a matrix logarithm and exponential can be conveniently processed in a metric
and can always be mapped back to a valid symmetric diffusion tensor [5]. Let
$\mathbf{M}=\mathbf{U}\boldsymbol{\Upsigma}{\mathbf{U}}^{\top}$ be the
eigendecomposition of a symmetric matrix $\mathbf{M}$. The computation of the
logarithm and the exponential of a tensor noted as $\log_{\mathrm{Id}}$ and
$\exp_{\mathrm{Id}}$ are defined as follows:

$\displaystyle\forall\mathbf{P}\in\mathcal{S}_{++}^{3},\
\log_{\mathrm{Id}}(\mathbf{P})$ $\displaystyle\ =\
\mathbf{U}\,\log(\boldsymbol{\Upsigma})\,{\mathbf{U}}^{\top}\,\in\,\mathcal{T}_{\mathrm{Id}}$
(3) $\displaystyle\forall\mathbf{S}\,\in\,\mathcal{T}_{\mathrm{Id}},\
\exp_{\mathrm{Id}}(\mathbf{S})$ $\displaystyle\ =\
\mathbf{U}\,\exp(\boldsymbol{\Upsigma})\,{\mathbf{U}}^{\top}\,\in\,\mathcal{S}_{++}^{3}$
(4)

Using these definitions, the geodesic distance between two points on the
$\mathcal{S}_{++}^{3}$ manifold, $\mathbf{P}_{1}$ and $\mathbf{P}_{2}$, can
then be expressed as

$\displaystyle
dist(\mathbf{P}_{1},\mathbf{P}_{2})\,=\,\big{\|}\log_{\mathrm{Id}}(\mathbf{P}_{1})-\log_{\mathrm{Id}}(\mathbf{P}_{2})\big{\|}_{2}.$
(5)

We provide in A the definition of the matrix backpropagation for diffusion
tensors learning used to train our network and the link to our open-source
implementation. For more details, the reader is referred to [27].

###  2.2 Riemannian Framework for ODF Learning

Orientation distribution functions, here represented as
$\mathbf{p}(\mathbf{s})$, are probability density functions modelling the
diffusion of water molecules at any angle $\mathbf{s}$ on the 2-sphere
$\mathbb{S}^{2}$. The space $\mathcal{P}$ of such PDFs forms the set

$\mathcal{P}=\big{\\{}\mathbf{p}:\mathbb{S}^{2}\to\mathbb{R}^{+}\,|\,\forall\mathbf{s}\in\mathbb{S}^{2},\mathbf{p}(\mathbf{s})\geq
0;\,\int_{\mathbf{s}\in\mathbb{S}^{2}}\mathbf{p}(\mathbf{s})d\mathbf{s}=1\big{\\}}$
(6)

The constrained function space $\mathcal{P}$ is not a vector space but a
nonlinear differentiable manifold that, just like the aforementioned
$\mathcal{S}_{++}^{3}$ manifold, needs to be equipped with an efficient
Riemannian metric to accurately process statistics on it [48]. Fortunately,
PDFs can be re-parameterized in multiple ways leading to known manifolds with
closed-form and computationally-efficient Riemannian operations. The square-
root re-parameterization of PDFs is a particularly convenient one as it
results in a unit Hilbert sphere manifold with an $\mathbb{L}^{2}$ metric
[48].

####  2.2.1 Square Root Re-Parameterization of ODF

With the help of the square-root re-parameterization of ODF
$\psi(\mathbf{s})=\sqrt{\mathbf{p}(\mathbf{s})},\forall\mathbf{s}\in\mathbb{S}^{2}$,
the space $\boldsymbol{\psi}$ can be viewed as the positive orthant of a unit
Hilbert sphere:

$\boldsymbol{\psi}=\big{\\{}\psi:\mathbb{S}^{2}\to\mathbb{R}^{+}\,|\,\forall\mathbf{s}\in\mathbb{S}^{2},\psi(\mathbf{s})\geq
0;\,\int_{\mathbf{s}\in\mathbb{S}^{2}}\psi^{2}(\mathbf{s})d\mathbf{s}=1\big{\\}}$
(7)

where the geodesic, exponential and logarithm maps are defined in a closed
form. In practice, ODFs are represented as histograms of $N$ bins where each
bin represents an orientation on a discretized 2-sphere. In the context of
learning spatially organized ODFs in 3D, working directly with the re-
parameterized ODFs $\psi(\mathbf{s})$ would require considerable resources.
Indeed, $N$ is typically in the order of few hundreds. To alleviate this
burden, $\psi(\mathbf{s})$ can be represented in a more compact form using a
spherical harmonic basis [12, 8] as follows:

$\psi(\mathbf{s})=\sum_{i=1}^{K}c_{i}\mathbf{B}_{i}(\mathbf{s})$ (8)

Here $K\ll N$ is the number of orthonormal basis functions used to represent
$\psi(\mathbf{s})$ and $\\{\mathbf{B}_{i}\\}_{i\in K}$ is the set of spherical
harmonic basis functions as in [12]. From this parametric representation, an
efficient Log-Euclidean framework has been proposed in [8] and is presented in
Section 2.2.2. Similarly to the FA of the diffusion tensor model, the
generalized fractional anisotropy (GFA) [52] of the ODF can be computed from
the spherical harmonic coefficients as in Eq. (9) below:

$\begin{split}\mathrm{GFA}=\sqrt{1-\frac{(c_{0}^{0})}{\sum_{k=0}^{L}\sum_{m=-k}^{k}(c_{k}^{m})^{2}}}.\end{split}$
(9)

Here, the GFA measure how far is the ODF from the uniform distribution.

####  2.2.2 Log-Euclidean Metric

Given the parametric representation of $\psi(\mathbf{s})$ in Eq. (8), the
square root of any ODF can be expressed by its Riemannian coordinate
$\mathbf{c}=(c_{1},c_{2},\ldots,c_{K})^{\top}$ and gives the probability
family $PF_{K}$:

$\footnotesize
PF_{K}=\Big{\\{}p(\mathbf{s}\,|\,\mathbf{c})=\Big{(}\sum_{i=1}^{K}c_{i}\mathbf{B}_{i}(\mathbf{s})\Big{)}^{2}:\int_{\mathbf{s}}p(\mathbf{s}\,|\,\mathbf{c})d\mathbf{s}=\sum_{i=1}^{K}c_{i}^{2}=1,\sum_{i=1}^{K}c_{i}\mathbf{B}_{i}(\mathbf{s})\geq
0,\forall\mathbf{s}\in\mathbb{S}^{2}\Big{\\}}.$ (10)

Following Eq. (10), the parameter space $PS_{K}$ can be defined as

$PS_{\\!K}=\big{\\{}\mathbf{c}\
|\,\left\|\mathbf{c}\right\|=\sum_{i=1}^{K}c_{i}^{2}=1,\
\sum_{i=1}^{K}c_{i}\mathbf{B}_{i}(\mathbf{s})\geq
0,\forall\mathbf{s}\in\mathbb{S}^{2}\big{\\}}$ (11)

which is also a subset of the sphere manifold $\mathbb{S}^{K-1}$. The sphere,
being a simple and well-studied manifold, makes the Log-Euclidean framework
for ODFs computation straight-forward and efficient as seen in Eq. (12) and
Eq. (13) below:

$\displaystyle\forall\mathbf{c}\in PS_{\\!K}\subset
S^{k-1},\,\log_{\mathbf{u}}(\mathbf{c})=\frac{\mathbf{c}-\mathbf{u}\,\mathrm{cos\Psi}}{\,\|\mathbf{c}-\mathbf{u}\,\mathrm{cos\Psi}\|_{2}}\Psi,$
$\displaystyle\qquad\text{ where
}\Psi=\mathrm{arcos}(\left\langle\mathbf{u}|\mathbf{c}\right\rangle)$ (12)
$\displaystyle\forall\mathbf{v_{c}}\in\mathcal{T}_{\mathbf{u}},\,\exp_{\mathbf{u}}(\mathbf{v_{c}})=\mathbf{u}\,\mathrm{cos\,\Psi}+\frac{\mathbf{v_{c}}}{\,\|\mathbf{v_{c}}\|_{2}}\,\mathrm{sin}\,\Psi,\text{
where }\Psi=\|\mathbf{v_{c}}\|_{2}.$ (13)

Here, $\mathbf{u}$ is the uniform orientation distribution function defined as
$\mathbf{u}=(1,0,\ldots,0)$. We use the maps in Eq. (12) and Eq. (13) to
accurately learn ODFs and ensure their validity throughout the training
process. Furthermore, the Log-Euclidean framework offers a simple geodesic
estimation between two parameterized ODFs $p(\,\cdot\,|\,\mathbf{c})$ and
$p(\,\cdot\,|\,\mathbf{c}^{\prime})$ as follows:

$dist\big{(}p(\,\cdot\,|\,\mathbf{c}),p(\,\cdot\,|\,\mathbf{c}^{\prime})\big{)}\,=\,\left\|\log_{\mathbf{u}}(\mathbf{c})-\log_{\mathbf{u}}(\mathbf{c}^{\prime})\right\|_{2}.$
(14)

###  2.3 Adversarial Training

Figure 4: On the left, our forward cycle: A T1w image is translated to a high-
resolution DT on the $\mathcal{S}_{++}^{3}$ manifold and back to the T1w
domain where the cycle consistent loss is computed. On the right, our backward
cycle: An upsampled DT on $\mathcal{T}_{\mathrm{Id}}$ is translated to the T1w
domain and back to $\mathcal{T}_{\mathrm{Id}}$ where the cycle consistent loss
is computed.

In a standard GAN setup [21], a generator network $G$ tries to generate
samples so close to the true data distribution that a discriminator network
$D$ is unable to distinguish between real and fake examples. Following the
CycleGAN architecture, our method uses two generators and two discriminators
denoted as $G_{X},G_{Y},D_{X}$ and $D_{Y}$. Here, $G_{X}$ takes a batch of
upsampled diffusion volumes as input and tries to fool $D_{X}$ by generating
realistic high-resolution T1w volumes. Similarly, $G_{Y}$ takes HR T1w volumes
as input and tries to generate plausible diffusion volumes in the same HR
space.

Because we only have access to LR diffusion data, the synthesized HR diffusion
is downsampled using a learned residual function $\mathcal{F}$ prior to
$D_{Y}$’s assessment as shown in Figure 2. Our adversarial objectives follow
the LSGAN formulation in [38] and are expressed as follows:

$\displaystyle\operatorname*{min}_{D_{Y}}\,\mathcal{L}_{\mathrm{LSGAN}}(G_{Y},D_{Y},X_{\mathrm{HR}},Y_{\mathrm{LR}})\
=$
$\displaystyle\qquad\tfrac{1}{2}\,\mathbb{E}_{\mathbf{y}\sim\mathbb{P}_{Y_{\mathrm{LR}}}}\big{[}(D_{Y}(\log(\mathbf{y}))-1)^{2}\big{]}+\tfrac{1}{2}\,\mathbb{E}_{\mathbf{x}\sim\mathbb{P}_{X_{\mathrm{HR}}}}\big{[}(D_{Y}(\mathcal{F}(G_{Y}(\mathbf{x}))))^{2}\big{]}$
(15)
$\displaystyle\operatorname*{min}_{D_{X}}\,\mathcal{L}_{\mathrm{LSGAN}}(G_{X},D_{X},Y_{\mathrm{LR}},X_{\mathrm{HR}})\
=$
$\displaystyle\qquad\tfrac{1}{2}\,\mathbb{E}_{\mathbf{x}\sim\mathbb{P}_{X_{\mathrm{HR}}}}\big{[}(D_{X}(\mathbf{x})-1)^{2}\big{]}+\tfrac{1}{2}\,\mathbb{E}_{\mathbf{y}\sim\mathbb{P}_{Y_{\mathrm{HR}}}}\big{[}(D_{X}(G_{X}(\\!\uparrow\\!\log(\mathbf{y}))))^{2}\big{]}$
(16)

where $\uparrow$ represents trilinear upsampling and $\log$ is the logarithm
map defined in either Eq. (3) or Eq. (12). It should be noted that both the
upsampling of the real data and the discriminator evaluation of the generated
diffusion information are performed in the Log-Euclidean domain to account for
the underlying data manifold.

###  2.4 Cycle-consistency Loss

The adversarial losses alone are not sufficient to drive the generation of HR
diffusion data. Indeed, $D_{Y}$ only evaluates downsampled data and, thus,
cannot help $G_{Y}$ improving beyond a certain precision level. Therefore, the
cycle-consistency loss denoted in Eq. (2.4) not only helps ensuring the
structural coherence of the synthesized images across modalities, but also
provides important high-resolution gradients to train $G_{Y}$.

Our cycle-consistency loss is threefold: 1) the error between the original HR
structural volume $\mathbf{x}$ and the reconstructed volume
$G_{X}(G_{Y}(\mathbf{x}))$, 2) the error between the upsampled diffusion
$\\!\uparrow\\!\\!\log(\mathbf{y})$ and its HR reconstruction
$G_{Y}\big{(}G_{X}(\\!\uparrow\\!\log(\mathbf{y}))\big{)}$, and 3) the error
between the original LR diffusion $\mathbf{y}$ and the downsampled recovered
volume
$\mathcal{F}\big{(}G_{Y}\big{(}G_{X}(\uparrow\\!\log(\mathbf{y}))\big{)}\big{)}$.
Combining these in a single loss gives

$\displaystyle\mathcal{L}_{\mathrm{cyc}}(G_{Y},G_{X})\ =\
\lambda_{\mathrm{cyc}_{X}}\underbrace{\mathbb{E}_{\mathbf{x}\sim\mathbb{P}_{X_{\mathrm{HR}}}}\big{[}\left\|G_{X}(G_{Y}(\mathbf{x}))\,-\,\mathbf{x}\right\|_{1}\big{]}}_{\text{Forward
Cycle HR}}$ $\displaystyle\qquad\quad\ +\
\tfrac{1}{2}\lambda_{\mathrm{cyc}_{Y}}\underbrace{\mathbb{E}_{\mathbf{y}\sim\mathbb{P}_{Y_{\mathrm{LR}}}}\big{[}\left\|G_{Y}\big{(}G_{X}(\uparrow\\!\log(\mathbf{y}))\big{)}\,-\,\\!\uparrow\\!\log(\mathbf{y})\right\|_{1}\big{]}}_{\text{Backward
Cycle HR}}$ $\displaystyle\qquad\qquad\quad\ +\
\tfrac{1}{2}\lambda_{\mathrm{cyc}_{Y}}\underbrace{\mathbb{E}_{\mathbf{y}\sim\mathbb{P}_{Y_{\mathrm{LR}}}}\big{[}\left\|\mathcal{F}\big{(}G_{Y}\big{(}G_{X}(\uparrow\\!\log(\mathbf{y}))\big{)}\big{)}\,-\,\log(\mathbf{y})\right\|_{1}\big{]}}_{\text{Backward
Cycle LR}}$ (17)

We employ the $\ell_{1}$ norm in Eq. (2.4) to measure both the forward and
backward cycle reconstruction errors, as it is less sensitive to large errors
than the $\ell_{2}$ norm [56]. Again, the log map is used to project the
generated and the real manifold-valued data onto a tangent plane before
computing the cycle-consistency loss. Furthermore, two parameters
$\lambda_{\mathrm{cyc}_{X}}$ and $\lambda_{\mathrm{cyc}_{Y}}$ control the
contribution of both cycles and have been empirically tuned. Full cycles,
including the manifold mappings and the loss computation in the tangent plane,
are illustrated in Figure 4.

###  2.5 Image Prior Regularization

By using a cycle-consistency loss in both directions, the CycleGAN model is
able to learn a bijective mapping between two domains using unpaired examples
[58]. However, for many cross-domain translation problems, the solution space
is extremely large and the model does not necessarily converge to a solution
that satisfies important domain-specific properties [37]. This is problematic,
especially in the case of medical images synthesis where the generated images
must not only be realistic from the discriminator’s point of view, but also be
faithful to expected results of the downstream tasks and known anatomical
properties. Thus, to ensure the model’s convergence towards plausible
solutions, we introduce a prior loss as follows:

$\displaystyle\mathcal{L}_{\mathrm{prior}}(G_{Y},G_{X})\ =\
\lambda_{\mathrm{prior}_{X}}\mathbb{E}_{\mathbf{x}\sim\mathbb{P}_{X_{\mathrm{HR}}},\mathbf{y}\sim\mathbb{P}_{Y_{\mathrm{LR}}}}\big{[}\left\|G_{Y}(\mathbf{x}_{i})\,-\,\uparrow\\!\log(\mathbf{y}_{i})\right\|_{1}\big{]}$
$\displaystyle\qquad\quad\ +\
\lambda_{\mathrm{prior}_{Y}}\mathbb{E}_{\mathbf{x}\sim\mathbb{P}_{X_{\mathrm{HR}}},\mathbf{y}\sim\mathbb{P}_{Y_{\mathrm{LR}}}}\big{[}\left\|G_{X}(\uparrow\\!\log(\mathbf{y}_{i}))\,-\,\mathbf{x}_{i}\right\|_{1}\big{]}$
(18)

where $\mathbf{x}_{i}$ and $\mathbf{y}_{i}$ are paired volumes taken from a
limited number of subjects. With this loss, the super-resolved diffusion stays
close to the real upsampled diffusion while integrating high-frequency
elements from the HR structural images.

###  2.6 Diffusion Anisotropy Weighted Loss

Figure 5: Comparison of the real FA (blue) and the generated FA (orange)
distribution density of the same subject with and without diffusion anisotropy
weighted loss.

Voxels expressing fiber tracts information typically have higher FA values
than those representing other tissues like grey matter (GM) or cerebrospinal
fluid (CSF). Therefore, we would like our diffusion synthesis method to be
particularly accurate in regions with higher fractional anisotropy. However,
as seen in Figure 5, voxels with high FA are underrepresented compared to
those with lower values. Consequently, this imbalance problem drives the
network’s generation towards diffusion with FA close to the mean. To alleviate
this issue, we weight the diffusion error in $\mathcal{L}_{\mathrm{cyc}}$ and
$\mathcal{L}_{\mathrm{prior}}$ by the FA/GFA of the target volume at every
voxel. The benefit of such mechanism can be observed on the density plots at
the bottom of Figure 5 which clearly exhibit a more faithful FA distribution
when using the proposed diffusion anisotropy weighting scheme.

###  2.7 Full Objective

Combining all loss terms, our full objective function is given by

$\displaystyle\mathcal{L}\ $ $\displaystyle=\
-\,\mathcal{L}_{\mathrm{LSGAN}}(G_{X},D_{X},Y_{\mathrm{LR}},X_{\mathrm{HR}})\,-\mathcal{L}_{\mathrm{LSGAN}}(G_{Y},D_{Y},X_{\mathrm{HR}},Y_{\mathrm{LR}})$
$\displaystyle\qquad\qquad+\,\mathcal{L}_{\mathrm{cyc}}(G_{Y},G_{X})\,+\,\mathcal{L}_{\mathrm{prior}}(G_{Y},G_{X})$
(19)

As in standard adversarial learning approaches, we train the generators and
discriminators concurrently by solving a mini-max problem:

$\begin{split}G_{X}^{*},G_{Y}^{*}\,=\,\arg\operatorname*{min}_{G_{X},G_{Y}}\operatorname*{max}_{D_{X},D_{Y}}\mathcal{L}(G_{X},G_{Y},D_{X},D_{Y})\end{split}$
(20)

Hence, Eq. (19) combines the error feedback from both the HR structural images
and the LR diffusion in an adversarial and voxel-wise manner.

##  3 Results

###  3.1 Data and Pre-Processing

We employ the T1w and diffusion MRI data of 1,065 subjects from the HCP1200
release of the Human Connectome Project [53] to evaluate our manifold-aware
CycleGAN. The T1w (0.7mm3 voxels, FOV=224mm, matrix=320, 256 sagittal slices
in a single slab) and diffusion (1.25mm3 voxels, sequence=Spin-echo EPI,
repetition time (TR)=5520 ms, echo time (TE)=89.5ms) data were acquired with a
Siemens Skyra 3T scanner [47] and minimally processed following [20].

Diffusion tensors were fitted using the DSI Studio toolbox software [29] and
the dODFs estimated using the constant solid angle (CSA) method [1] from the
DIPY library [17]. Diffusion ODFs were then re-parameterized following Section
2.2.1 and further estimated using 4th order spherical harmonics [12]. Both the
DT and the ODF volumes were transformed to the Log-Euclidean domain using
their respective Riemannian framework described in Section 2. Once in the Log-
Euclidean domain, these volumes were upsampled to the T1w spatial resolution
using trilinear interpolation and aligned to their corresponding HR T1w
images. In experiments, we consider these upsampled and aligned diffusion
volumes as the “ground truth” diffusion. The T1w images have been rescaled to
the [0,1] range by min-max normalization. Finally, both the structural and the
diffusion volumes, in high and low-resolution, were decomposed in overlapping
patches of $32^{3}$ and $18^{3}$ voxels respectively. Volumes are processed
patch-wise by our model for two important reasons: 1) limiting the memory
required by the model to compute network activations and outputs, and 2)
increasing the amount of training examples.

Figure 6: The architecture of our residual downsampling function $\mathcal{F}$
used to reduce the resolution of our synthesized HR diffusion before $D_{Y}$
assessment.

###  3.2 Implementation Details

Our two generator networks, $G_{X}$ and $G_{Y}$, follow the U-Net
implementation in [9] where the last activation layer has been replaced to fit
the different output data scale. In all the experiments, we used a sigmoid in
$G_{X}$ to generate T1w images within the [0,1] range. For the generation of
diffusion, the last activation function of $G_{Y}$ varies depending on the
generated diffusion model. For the generation of DT, we used a hard hyperbolic
tangent function and for the generation of ODF, a tanh activation.

Moreover, we set the amount of input and output channels of our networks
according to our input and output data shape as described in Table 1. The DT
inputs are of shape $9\times 32\times 32\times 32$, where the 9 channels
represent the flattened $3\times 3$ diffusion tensors at every voxel. The ODF
inputs are of shape $C\times 32\times 32\times 32$ where C depends on the
order of the spherical harmonic basis that is used to represent them. Due to
limitations on computational resources, we estimated the diffusion ODFs using
4th order spherical harmonics, yielding a total of $C=15$ coefficients.

Our discriminator networks $D_{X}$ and $D_{Y}$ follow the SRGAN discriminator
architecture in [34] where we replaced the 2D convolutions by 3D ones.
Furthermore, we reduced the number of feature maps in convolution layers to
32, 64, 128 and 256 as suggested in [43] for volumetric data. In all
scenarios, $D_{X}$ assesses T1w volumes of shape $1\times 32\times 32\times
32$ while $D_{Y}$ takes as inputs volumes of shape $9\times 18\times 18\times
18$ for DT and $15\times 18\times 18\times 18$ for ODFs.

###  3.3 Training Setup

To train our network, we randomly selected 70% (746 subjects) of the 1,065
subjects for training, 20% (213 subjects) for validation and 10% (106
subjects) for testing. From the 745 training subjects, we kept aside between 0
and 75 subjects as paired priors in Eq. (2.5). To form our training,
validation and test sets, we randomly selected 50,000 unpaired HR T1w and LR
diffusion patches from remaining training subjects, 10,000 patches from
validation subjects and 5,000 patches from test subjects. The same number of
patches (50,000) was selected from our subjects kept as paired priors, i.e.
aligned HR T1w and upsampled diffusion, to form a paired training set.
Although we kept the number of patches in the paired training set constant, we
experimentally adjusted the number of subjects used to randomly extract them.
To do so, we trained our network with an increasing number of paired subjects,
as reported in Table 2, and kept the best setup for our following experiments.

Table 1: Implementation details of our generator networks $G_{X}$ and $G_{Y}$ and discriminator networks $D_{X}$ and $D_{Y}$. Network |  Input Shape |  Output Shape |  Last Activation
---|---|---|---
$G_{Y}$ (T1w $\mapsto$ DT) |  $\phantom{0}1\times 32\times 32\times 32$ |  $\phantom{0}9\times 32\times 32\times 32$ |  Hardtanh
$G_{X}$ (DT $\mapsto$ T1w) |  $\phantom{0}9\times 32\times 32\times 32$ |  $\phantom{0}1\times 32\times 32\times 32$ |  Sigmoid
$G_{Y}$ (T1w $\mapsto$ ODF) |  $\phantom{0}1\times 32\times 32\times 32$ |  $15\times 32\times 32\times 32$ |  Tanh
$G_{X}$ (ODF $\mapsto$ T1w) |  $15\times 32\times 32\times 32$ |  $\phantom{0}1\times 32\times 32\times 32$ |  Sigmoid
$D_{X}$ (T1w) |  $\phantom{0}1\times 32\times 32\times 32$ |  1 |  Linear
$D_{Y}$ (DT) |  $\phantom{0}9\times 18\times 18\times 18$ |  1 |  Linear
$D_{Y}$ (ODFs) |  $15\times 18\times 18\times 18$ |  1 |  Linear

We trained the networks using an Adam optimizer [32] with a learning rate of
$10^{-4}$ and beta1, beta2 values of 0.5 and 0.999. Hyper-parameters
$\lambda_{\mathrm{prior}_{X}}$, $\lambda_{\mathrm{prior}_{Y}}$,
$\lambda_{\mathrm{cyc}_{X}}$ and $\lambda_{\mathrm{cyc}_{Y}}$ were
experimentally set to 10, 0.5, 5 and 0.25. Furthermore, we used a reduce on
plateau learning rate scheduler with a patience of 10 epochs and a factor of
10. Batches of 8 patches were used and the models were trained for 35 epochs
($\sim 220$k steps) on an NVIDIA TITAN XP GPU with 12 GB of VRAM. All
experiments were repeated three times with a different initialization seed.

###  3.4 Baselines

As mentioned before, deep learning models for the synthesis of manifold-valued
data are just starting to emerge. Consequently, the number of baselines
requiring minimal adaptation for the evaluation of our model is limited.
Nevertheless, we compare our model to the three approaches described below.

##### Manifold-Aware WGAN

Our first baseline is an adaptation of the manifold-aware WGAN presented in
[26] for the conditional generation of diffusion from structural T1w images.
This method denoted as “MA-WGAN” in our results, can generate plausible
manifold-valued images by incorporating the Log-Euclidean maps within the
network and therefore provides a natural point of comparison for our method.
For this baseline, we use the same generator $G_{Y}$ and discriminator $D_{Y}$
as in our proposed model. The manifold mapping used in the network is changed
according to the generated diffusion scheme following Section 2.

##### Manifold-Aware U-Net

We also compare our method to a supervised U-Net model. For this baseline, we
use the same generator as in our own architecture, i.e. $G_{Y}$ [9], but train
it in a supervised manner with paired HR T1w and upsampled diffusion volumes
in the Log-Euclidean domain. Similar to the _MA-WGAN_ baseline, we change the
manifold mapping of the network according to the generated diffusion
reconstruction scheme (DT or ODF). This baseline, denoted as “MA-U-Net” in
results, helps us measure the effect of our adversarial and cycle-consistent
losses.

##### U-Net

Our last baseline is a standard supervised U-Net [9] trained with paired HR
T1w images and upsampled diffusion without manifold-awareness. With this
method, we aim at measuring the performance gain of our method induced by both
the manifold-mapping and the use of additional unpaired samples. In addition,
we validate that manifold-awareness is necessary to synthesize realistic
samples strictly lying on the data manifold.

###  3.5 DT and ODF Synthesis

We first test our network and baselines for the task of DT and diffusion ODF
synthesis. In this setup, we train our network with unpaired HR T1w patches
and LR diffusion patches in the Log-Euclidean domain. Moreover, we use a
paired training set of 50 000 patches from 50 randomly chosen subjects. The
paired training set is used as prior for our method and as the training set
for our baselines. Hence, both our method and baselines are trained with the
same amount of paired information.

##### Evaluation Metrics

Three metrics were considered to quantitatively evaluate the generated
diffusion. First, we use the cosine similarity to compare the principal fiber
orientation of every synthesized tensor and ODF to their expected real
orientation:

$\begin{split}\mathrm{similarity}(\mathbf{a},\mathbf{b})\,=\,\left\lvert\frac{\mathbf{a}\cdot\mathbf{b}}{\left\|\mathbf{a}\right\|\left\|\mathbf{b}\right\|}\right\rvert\end{split}$
(21)

To retrieve the main orientation of the ODF, we first calculate the spherical
coordinates at which the value of the ODF is maximum using a discretized
sphere of 724 vertices. We then convert these spherical coordinates to
Euclidean coordinates to obtain their principal orientation vector. For DT,
the main orientation is given by the eigenvector associated with the largest
eigenvalue of the DT as described in Section 2.1.

Figure 7: Metrics on the sagittal, coronal and axial slices between the
generated HR diffusion of a random test subject and its interpolated ground-
truth. (Top row) FA, (middle row) geodesic distance and (bottom row) cosine
similarity.

In addition, we compare the generated and real diffusion with the mean square
error (MSE) between their FA/GFA. The FA and GFA are computed following Eq.
(1) and Eq. (9). This help evaluating the shape of the synthesized DT/ODF
independently of their orientation.

Finally, we measure the mean geodesic distance using Eq. (5) for DT and Eq.
(14) for ODF between our generated data and the real diffusion. This latter
metric encode both the orientation and the shape error in a single measure.

###  3.6 Tractography

To further assess the integrity of the synthesized diffusion volumes by the
proposed method, we performed whole-brain tractography on both the real and
the generated data, and segmented the resulting tractograms into bundles. We
then posed the tractograms generated on real data as ground truth and
extracted quantitative measures from whole-brain tractograms. Likewise, we
segmented bundles to measure how much is tractography impacted. In the
following subsections, we describe each step of the analysis.

Figure 8: Mean streamline length (in mm) for the segmented bundles from DT
tractography (top row) and ODF tractography (bottom row).

##### Streamlines Generation and Bundling

Tracking was performed using the EuDX algorithm [16] with a step-size of 0.5
mm. A maximum angle of 60 degrees was used between steps, using the principal
direction of the diffusion tensor and maxima of ODF. Maxima were extracted
from the ODF using scilpy 111https://github.com/scilus/scilpy. Seeding was
done at 2 seeds per voxel on the whole white-matter mask, which was computed
from the ground-truth T1w image using Dipy [17]. Streamlines with a length
below 10 mm or above 300 mm were discarded. Whole brain tractograms were then
segmented using RecobundlesX [18], using 80 bundles from Yeh et al. [55] as
reference. To allow for a more robust comparison, and because initial
streamline points placement depends on randomness which may have an impact on
the reconstructed streamlines, tracking was performed five times with
different random seeds on each volume.

##### Streamlines Assessment

Figure 9: Visualization of three segmented bundles from real and synthesized
tractograms. (On the left), left Corticospinal Tract (CST), Arcuate Fasciculus
(AF) and Inferior Longitudinal Fasciculus (ILF) segmented from real and
generated DT data. (On the right), left CST, AF and ILF segmented from real
and synthesized ODF data.

The similarity between the streamlines reconstructed from the real and
generated diffusion volumes is assessed using three measures: streamline
length, bundle volume and bundle shape. First, we compare the streamline
length (in mm) between all real and synthesized whole-brain tractograms, as
well as each segmented bundle. We also compare the volume occupied by the
whole-brain tractograms and bundles in a voxel-wise matter. Furthermore, the
shape similarity of reconstructed tractograms is measured by their voxel-wise
Dice, overlap (OL) and overreach (OR). The OL, defined as

$\begin{split}\mathrm{OL}\,=\,\frac{|B\cap A|}{|A|},\end{split}$ (22)

where $A$, $B$ are binary bundle masks, quantifies how much the volume of
bundle $A$ is reconstructed by bundle $B$. The OR, expressed as

$\begin{split}\mathrm{OR}\,=\,\frac{|B\cup A|-|B\cap A|}{|A|},\end{split}$
(23)

evaluates how much of bundle $B$ goes over bundle $A$. Segmented bundles were
merged to allow for pairwise comparison between real and generated data.
Merged bundles with fewer than 100 streamlines were discarded from the
analysis. To evaluate the overall reconstruction quality, we also report Dice,
OL and OR between the reconstructed tractograms and the atlas used for bundle
segmentation.

Figure 9 presents some of the reconstructed bundles from real and synthesized
data.

###  3.7 Diffusion Synthesis Analysis

We report in Table 2 the performance of our model with an increasing number of
paired subjects used as prior in Eq. (2.5). We observe that for both DT and
ODF, the mean cosine similarity increases with the number of subjects and
reaches its peak when extracting the paired patches from 50 subjects. As we
keep the number of extracted patches constant (50,000 patches), 50 subjects
represent an average of 1,000 patches per subject which seems a good trade-off
between sampling diversity and sparsity. In this setup, our model yields a
mean cosine similarity of 0.8648 $\pm\,$5.5\text{\times}{10}^{-3}$$ and 0.8846
$\pm\,$4.4\text{\times}{10}^{-3}$$ in voxels with $\mathrm{FA}\geq 0.2$ for DT
and ODF respectively. In voxels with $\mathrm{FA}\geq 0.5$, a mean cosine
similarity of 0.9167 $\pm\,$3.2\text{\times}{10}^{-3}$$ is reached for DT and
0.9425 $\pm\,$1.4\text{\times}{10}^{-4}$$ for ODF. As reported in Table 4,
this corresponds to FA MSE values of 0.0089 $\pm\,$1.5\text{\times}{10}^{-3}$$
and 0.0159 $\pm\,$1.4\text{\times}{10}^{-4}$$ for DT and GFA MSE values of
0.0229 $\pm\,$2.8\text{\times}{10}^{-4}$$ and 0.0614
$\pm\,$1.3\text{\times}{10}^{-3}$$ for ODF.

Table 2: The cosine similarity between the principal orientation obtained by our method with different amount paired subjects used as prior in Eq. (2.5). Since the metrics are more relevant in regions that typically encode fibers information, we report them at increasing FA thresholds of 0.2 and 0.5. Paired [2pt] Subjects | Cosine Sim (DT) | Cosine Sim (ODF)
---|---|---
FA $\geq$ 0.2 |  FA $\geq$ 0.5 |  GFA $\geq$ 0.2 |  GFA $\geq$ 0.5
0 | 0.7315 | 0.7955 | 0.7029 | 0.7380
10 | 0.7795 | 0.8599 | 0.8664 | 0.9058
25 | 0.8152 | 0.8863 | 0.8745 | 0.9081
50 | 0.8648 | 0.9167 | 0.8846 | 0.9425
75 | 0.8317 | 0.8993 | 0.8674 | 0.9051

As a point of comparison, we give in Table 4 the performance of our baselines
when they are trained with the same 50 paired subjects. As can be seen, the
proposed model obtains the highest performance for all metrics when trained on
DT images, as well as better geodesic and fiber orientation estimation than
baselines when trained on ODF.

Compared to the fully-supervised U-Net, which does not enforce manifold
consistency on the output, our model improves FA MSE by 23.14%, mean geodesic
distance by 81.11% and mean cosine similarity by 4.23% for regions with
$\mathrm{FA}\geq 0.5$ for DT. This shows the benefit of imposing manifold-
awareness constraints on the network’s output. Without these constraints,
U-Net generates an average of 3,843.5 non-SPD tensors and 1,559.7 non-PDF ODF
as noted in Table 3. Our model also provides statistically better performance,
for both DT and ODF data, on the mean cosine similarity and mean geodesic
distance compared to MA-U-Net (paired t-test p $<$ 0.05), that does not
include cycle-consistency and is only trained with paired data.

Table 3: The amount of non-SPD tensors and non-PDF dODFs generated by the compared methods for five randomly selected test subjects. The manifold-awareness ensures that the generated diffusion schemes lie on their respective data manifold. Method |  Manifold-Awareness |  Non-SPD Tensors |  Non-PDF ODF
---|---|---|---
U-Net |  ✗ |  $3843.5\pm 481.22$ |  $1559.7\pm 122.3$
MA-U-Net |  ✓ |  0 |  0
MA-WGAN |  ✓ |  0 |  0
Ours |  ✓ |  0 |  0

Improvements are particularly important for voxels with $\mathrm{FA}\geq 0.5$,
where our method obtains a 3.43% higher mean cosine similarity, 11.19% lower
geodesic and 13.58% better FA MSE for DT and 4.18% higher mean cosine
similarity for ODF. This demonstrates the impact of our anisotropy-weighted
loss described in Section 2.6, which gives more importance to voxels with
higher anisotropy values.

Figure 10: Comparison of the color FA of the generated diffusion of compared methods on a sagittal, axial and coronal slice of a random test subject. Figure 11: Qualitative comparison of the generated ODF by the compared methods. First row, the color encoded GFA of a coronal slice of a test subject. Second row, the generated ODF and GFA in a region with crossing fibers. Third row, a close-up look of the generated ODF and GFA. Fourth row, a close-up look at the ODF peaks. Table 4: The fractional anisotropy mean squared error (FA MSE), geodesic, and cosine similarity between the principal orientations obtained by the different compared methods when trained with 50 paired subjects. Since the metrics are more relevant in regions that typically encode fibers information, we report them at increasing FA thresholds of 0.2 and 0.5. | Method | FA$\\!{}^{*}$ MSE | Geodesic | Cosine Similarity
---|---|---|---|---
| FA$\\!{}^{*}$ $\geq$ 0.2 | FA$\\!{}^{*}$ $\geq$ 0.5 | FA$\\!{}^{*}$ $\geq$ 0.2 | FA$\\!{}^{*}$ $\geq$ 0.5 | FA$\\!{}^{*}$ $\geq$ 0.2 | FA$\\!{}^{*}$ $\geq$ 0.5
DT | U-Net | 0.0100 | 0.0207 | 2.2559 | 2.5457 | 0.8266 | 0.8795
MA-WGAN | 0.0294 | 0.0781 | 0.7791 | 1.046 | 0.5724 | 0.6246
MA-U-Net | 0.0088 | 0.0184 | 0.3947 | 0.5415 | 0.8288 | 0.8863
Ours | 0.0089 | 0.0159 | 0.3531 | 0.4809 | 0.8648 | 0.9167
ODF | U-Net | 0.0220 | 0.0559 | 0.3639 | 0.6299 | 0.8702 | 0.9045
MA-WGAN | 0.0443 | 0.1173 | 0.5351 | 0.9121 | 0.6823 | 0.7278
MA-U-Net | 0.0219 | 0.0563 | 0.3611 | 0.6249 | 0.8710 | 0.9047
Ours | 0.0229 | 0.06145 | 0.3467 | 0.6238 | 0.8846 | 0.9425
_* : FA is used for DT and GFA for ODF_.

Our results can be further appreciated in Figure 7, where we report the
metrics yielded by our method on a sagittal, axial and coronal slice of a
randomly chosen test subject. We see that our network is able to recover most
of the fibers orientation, especially in regions with typically higher FA/GFA
like the corpus callosum. Moreover, the estimated FA/GFA is generally faithful
to the real data except in the corticospinal tract where the error is higher.

Finally, as a qualitative evaluation, we compare in Figure 10 the generated
color encoded FA of compared methods. From this figure, we can see how our
cycle-consistent and prior losses help our model converging towards plausible
solutions better than baselines, especially compared to the MA-WGAN that only
relies on an adversarial objective. We also observe a visually more faithful
orientations estimation by our method in the splenium of the corpus callosum
and in the pons. In Figure 12, we compare the real GFA map of a test subject
to the generated maps by the tested methods. We also compare the real and
synthesized GFA maps to the associated HR T1w image of this said subject. It
can be seen from Figure 12 that our method generates GFA maps that are more
consistent with the real HR T1w image of the subject while visually reducing
boundary artifacts. Our method also seems to recover fine anatomical details
that are present in the HR T1w but not in the real diffusion. Indeed, the GFA
maps generated by our method exhibit sharper edges and better tissue
delineation particularly at the bottom of the coronal slice where the real
diffusion lacks details. We provide in C the color encoded FA and FA maps of
five additional test subjects to visually assess the generated images by our
method on different brain geometries. In Figure 11, we compare the ODFs
produced by our method and baselines in a region with crossing fibers. One can
see that our method generates ODFs that smoothly transition between
orientations and plausible crossing estimation.

###  3.8 Streamlines Length and Volume

We report in Figure 13 the mean streamline length (in mm) and volume (in
voxels) for the reconstructed wholebrain tractograms of real and generated
diffusion. We further detail in Figure 8 the mean streamline length and, in
Figure 14, the mean bundle volume for each segmented bundles.

We can observe from the results that reconstructed whole brain tractograms are
very similar in size, but that individual bundles segmented from synthesized
data tend to be shorter. Indeed, looking at the mean bundle lengths and
volumes reported in Table 5, it can be seen that bundles segmented from the
generated DT/ODF tractograms are respectively 14.98% and 2.32% shorter than
real data. Nonetheless, generated ODFs tend to produce bundles that are
slightly more voluminous with mean volume of $\sim$24,408 voxels compared to
$\sim$24,210. While some bundles were only segmented on the real data or the
generated data (IFOF vs. TPT, for example), we can observe that the bundles
recovered in both cases exhibit similar statistics.

Table 5: Measures on the real and generated segmented bundles from tractography. |  | Real | Generated
---|---|---|---
|  | (mean ± std) | (mean ± std)
DT | Length (mm) | 99.99 ± 54.37 | 85.01 ± 50.02
Volume (voxels) | 15963.62 ± 24199.55 | 14119.52 ± 24274.11
ODF | Length (mm) | 120.94 ± 72.29 | 118.13 ± 71.57
Volume (voxels) | 24210.36 ± 42082.48 | 24408.78 ± 44255.32

###  3.9 Streamlines Shape

We report in Figure 13 the Dice, OL and OR between the wholebrain tractograms
of real and generated diffusion. The same measures are then detailed in Figure
15 for all segmented bundles. We can see from these two figures that, despite
their similarity in length and volume, the space occupied by segmented bundles
from real and generated data may vary. This disparity can be observed in
Figure 9 where we compare three real and generated segmented bundles.
Furthermore, we observe a high inter-bundle metrics performance variability.
For instance, bundles such as the Frontal Aslant Tract (FAT), Frontopontine
(FPT) and Vertical Occipital Fasciculus (VOF) reach a high agreement (i.e.
Dice, OL and OR close to 1) whereas the Superior Cerebellar Peduncle (SCP) and
Occipitopontine Tract (OPT) are hardly matched. This high inter-bundle
variability can be further appreciated in Table 6 where we note the mean Dice,
OL and OR and their standard deviation for all segmented bundles.

Table 6: Average metrics from the bundle shape comparison between real and synthesized data. |  Dice |  OL |  OR
---|---|---|---
|  (mean ± std) |  (mean ± std) |  (mean ± std)
DT |  0.46 ± 0.24 |  0.38 ± 0.21 |  1.36 ± 0.86
ODF |  0.54 ± 0.24 |  0.42 ± 0.19 |  1.23 ± 0.58

From Table 6, we can also see that the segmented bundles from the generated
ODF are more faithful to their real counterpart than DT with a mean Dice of
0.54 ± 0.24, a mean OL of 0.42 ± 0.19 and a mean OR of 1.23 ± 0.58.

##  4 Discussion

Figure 12: GFA of the real low-resolution diffusion (first column), the real
upsampled diffusion (second column) and the generated diffusion of the
compared methods. One can notice the ability of our network to recover fine
anatomical details that are clearly visible in the HR T1w image but not in the
original diffusion.

In this work, we propose a novel Riemannian deep learning architecture for the
synthesis of 3D manifold-valued data and have tested its performance on two
tasks: 1) the generation of diffusion tensors (DTs) and 2) the generation of
diffusion orientation diffusion functions (ODFs). Specifically, we have
explored the feasibility of generating high-resolution DT and ODF from high-
resolution structural T1w images and unpaired LR diffusion. We show in Table 3
that a standard model relying on Euclidean operations fails to capture the
geometry of the diffusion data manifold which leads to the estimation of
physically incorrect diffusion. To alleviate this issue, we have built a
framework on top of recent advances in manifold-valued data processing and
Riemannian geometry [5, 8, 26] to ensure the validity of the generated
diffusion. We have evaluated the generated volumes properties using mean
squared errors of FA/GFA maps, geodesic distances and cosine similarities
between real and predicted principal fiber orientation. To further evaluate
the integrity of the synthesized diffusion in a typical diffusion application,
we have performed tractography and assessed the lengths, volumes and shapes of
resulting tractrograms.

###  4.1 Diffusion Synthesis Performance

The generation of DT/ODF solely relying on a T1w image is an ill-posed problem
for which a single T1w intensity can correspond to several fiber arrangements.
However, by providing the contextual information required for the network to
localize the structural input, we observe that strong fiber patterns can
successfully be recovered by our method and baselines. Hence, we believe that,
to a certain extent, the high-level geometry of the brain globally impacts its
underlying fibers organization. This seems to be particularly true in regions
of higher anisotropy where fiber tracts are strongly organized. As a result,
we observed a better estimation of the principal fibers orientation in regions
with higher FA/GFA and a generally poor estimation of the principal
orientations in regions with high inter-subject variance such as the
ventricles.

Using HR structural images to drive the synthesize of diffusion helps
recovering fine anatomical details and sharp edges better than interpolation
based methods. By leveraging the detailed information contained in HR
structural images, our network is not limited by the coarseness of the low-
resolution input diffusion signal such as in interpolation.

This transfer of information from HR to LR images is enforced by our cycle-
consistency loss that preserves a high structural coherence between the HR
structural inputs and the generated diffusion. In addition, our adversarial
and prior loss help recovering plausible fiber patterns better than our
baselines by leveraging unpaired examples of real diffusion. Furthermore, the
manifold-awareness yields on-par or better metrics when compared with
equivalent Euclidean architectures. The manifold-awareness, however, comes
with the benefit of ensuring the mathematical properties of the synthesized
diffusion regardless of the amount of network training.

###  4.2 Tractography Performance

####  4.2.1 DT vs. ODF

We observe from Table 5 that bundles segmented from ODF tractography tend to
be longer and more voluminous which is an expected behavior. While a thorough
comparison between DT and ODF tractography is out-of-scope for this work (c.f.
[15, 50, 28]), we can nevertheless appreciate that the synthesized DT and ODF
behave in a manner similar to their real counterparts in the context of
tractography.

Moreover, we observed that the bundles segmented from ODF tractography have a
higher mean Dice, higher overlap and lower overreach than bundles segmented
from DT tractography. Since the same algorithm was used to perform
tractography on both datasets, the main difference is that DT tractography
makes use of a single direction while ODF tractography may use multiple local
maxima to propagate streamlines. As such, the lower discrepancies in ODF
bundle shapes can be explained by the multiple directions used in tractography
compensating for local errors. At the opposite, DT tractography is known to be
sensitive to local estimation error [23]. This sensitivity, which often lead
to the early termination or to the switch to a wrong adjacent tract of the
tracking algorithm [28], can greatly affect the final tractograms shape.

####  4.2.2 Bundle shape analysis

Since DT/ODF generation studies are still few, it is hard to provide a
definitive conclusion on the quality of the streamlines generated on
synthesized data. However, we observe from Figure 13 that whole brain
tractograms have a similar streamline length, occupy the same number of voxels
and have the same shape. Segmented bundles, if extracted from both real and
generated data, also tend to exhibit the same average length and volume.

While the reported bundle shape metrics might seem to indicate a poor
correspondence between bundles generated from real and synthesized data, it
should be noted that bundle segmentation is a highly variable operation. For
example, Rheault et al. [42], which analyzed the reproducibility of the
segmentation of a single bundle between human experts and non-experts, reports
a significant difference between the volume of segmented bundles between
experts and non-experts as well as median Dice scores around 0.77 for intra-
rater reproducibility, 0.65 for inter-rater reproducibility and 0.8 for
reproducibility with a gold standard. In a similar manner, Schilling et al.
[44] analyzed the variability in the segmentation of 14 bundles between 42
groups using both manual and automatic segmentation. While few actual metrics
are reported, figures indicate a generally low Dice score, as well as a high
variability in bundle volume and streamline length for inter-protocol
responsibility. Analysis for specific pathways report Dice scores between 0.4
and 0.6 for inter-protocol and inter-subject reproducibility and Dice scores
between 0.6 and 0.8 for intra-protocol reproducibility. As such, we can
theorize that the reported Dice scores in the present work could be impacted
by the inherent variability in the segmentation process.

##  5 Conclusion

This work presents a novel Riemannian network architecture for the cycle-
consistent synthesis of diffusion tensors and diffusion ODF in high-resolution
structural T1w space. The results have demonstrated that our Riemannian
architecture can synthesize valid diffusion images with a 5% improvement in
principal fibers orientation and a 23% improvement in FA MSE with respect to
our baselines. The better performance of our approach over compared methods
shows the benefit of using both paired and unpaired samples in a single
objective. Furthermore, as opposed to standard Euclidean deep learning models,
which generate an average of 3,844 invalid tensors and 1,560 invalid ODFs per
volume, our method guarantees the mathematical coherence of the synthesized
diffusion schemes, free of invalid tensors or ODFs.

Moreover, we have evaluated qualitatively our generated diffusion volumes by
comparing their tractograms with their real counterparts. It was observed that
our generated T1w-driven diffusion shares similarities with the real diffusion
in terms of streamline length, volume and fiber bundles shape. We have also
shown the ability of our network to transfer fine anatomical details from the
high-resolution T1w images to diffusion images. This transfer of information
allows the generation of images with sharper edges and a higher level of
details that could not be achieved with image interpolation.

Our results suggest that the high-level geometry of the brain, encoded in
structural T1w images, can be used to predict its global fiber bundles
organization. Leveraging this principle, our method could enable the fast
synthesis of DT and ODF in situations where the acquisition of diffusion
imaging is not available. More generally, it offers the basis of a framework
targeting any real-to-manifold-valued image translation tasks. For instance,
our method could be used for missing modalities synthesis and datasets
completion, manifold-valued image inpainting or manifold-valued population
based statistics that rely on non-Euclidean metrics.

## Acknowledgments

This work was supported financially by the Canada Research Chair on Shape
Analysis in Medical Imaging, the Research Council of Canada (NSERC), the Fonds
de Recherche du Québec (FQRNT), the Réseau de Bio-Imagerie du Québec (RBIQ),
and ETS Montreal.

## References

  * Aganj et al. [2010] Aganj, I., Lenglet, C., Sapiro, G., Yacoub, E., Ugurbil, K., Harel, N., 2010\.  Reconstruction of the orientation distribution function in single- and multiple-shell q-ball imaging within constant solid angle.  Magnetic Resonance in Medicine .
  * Alexander et al. [2001] Alexander, A.L., Hasan, K.M., Lazar, M., Tsuruda, J.S., Parker, D.L., 2001.  Analysis of partial volume effects in diffusion-tensor mri.  Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 45, 770–780.
  * Alexander et al. [2017] Alexander, D.C., Zikic, D., Ghosh, A., Tanno, R., Wottschel, V., Zhang, J., Kaden, E., Dyrby, T.B., Sotiropoulos, S.N., Zhang, H., Criminisi, A., 2017.  Image quality transfer and applications in diffusion MRI.  NeuroImage .
  * Arjovsky et al. [2017] Arjovsky, M., Chintala, S., Bottou, L., 2017.  Wasserstein generative adversarial networks, in: Precup, D., Teh, Y.W. (Eds.), Proceedings of the 34th International Conference on Machine Learning, PMLR. pp. 214–223.
  * Arsigny et al. [2006] Arsigny, V., Fillard, P., Pennec, X., Ayache, N., 2006\.  Log-Euclidean metrics for fast and simple calculus on diffusion tensors.  Magnetic Resonance in Medicine .
  * Brooks et al. [2019] Brooks, D., Schwander, O., Barbaresco, F., Schneider, J.Y., Cord, M., 2019.  Riemannian batch normalization for spd neural networks, in: Advances in Neural Information Processing Systems, pp. 15489–15500.
  * Chakraborty et al. [2019] Chakraborty, R., Bouza, J., Manton, J., Vemuri, B.C., 2019\.  A deep neural network for manifold-valued data with applications to neuroimaging, in: Information Processing in Medical Imaging, Springer International Publishing. pp. 112–124.
  * Cheng et al. [2009] Cheng, J., Ghosh, A., Jiang, T., Deriche, R., 2009\.  A Riemannian framework for orientation distribution function computing, in: Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), pp. 911–918.
  * Çiçek et al. [2016] Çiçek, Ö., Abdulkadir, A., Lienkamp, S.S., Brox, T., Ronneberger, O., 2016.  3D U-Net: Learning Dense Volumetric Segmentation from Sparse Annotation.  Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) .
  * Costabile et al. [2019] Costabile, J.D., Alaswad, E., D’Souza, S., Thompson, J.A., Ormond, D.R., 2019.  Current applications of diffusion tensor imaging and tractography in intracranial tumor resection.  Frontiers in Oncology .
  * Dar et al. [2019] Dar, S.U., Yurt, M., Karacan, L., Erdem, A., Erdem, E., Cukur, T., 2019\.  Image Synthesis in Multi-Contrast MRI with Conditional Generative Adversarial Networks.  IEEE Transactions on Medical Imaging 38, 2375–2388.
  * Descoteaux et al. [2007] Descoteaux, M., Angelino, E., Fitzgibbons, S., Deriche, R., 2007\.  Regularized, fast, and robust analytical Q-ball imaging.  Magnetic Resonance in Medicine .
  * Dyrby et al. [2014] Dyrby, T.B., Lundell, H., Burke, M.W., Reislev, N.L., Paulson, O.B., Ptito, M., Siebner, H.R., 2014.  Interpolation of diffusion weighted imaging datasets.  NeuroImage .
  * Elsaid and Wu [2019] Elsaid, N.M., Wu, Y.C., 2019\.  Super-Resolution Diffusion Tensor Imaging using SRCNN: A Feasibility Study, in: Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBS, Institute of Electrical and Electronics Engineers Inc.. pp. 2830–2834.
  * Farquharson et al. [2013] Farquharson, S., Tournier, J.D., Calamante, F., Fabinyi, G., Schneider-Kolsky, M., Jackson, G.D., Connelly, A., 2013.  White matter fiber tractography: why we need to move beyond dti.  Journal of neurosurgery 118, 1367–1377.
  * Garyfallidis [2013] Garyfallidis, E., 2013.  Towards an accurate brain tractography.  Ph.D. thesis. University of Cambridge.
  * Garyfallidis et al. [2014] Garyfallidis, E., Brett, M., Amirbekian, B., Rokem, A., Van Der Walt, S., Descoteaux, M., Nimmo-Smith, I., 2014.  Dipy, a library for the analysis of diffusion mri data.  Frontiers in Neuroinformatics .
  * Garyfallidis et al. [2018] Garyfallidis, E., Côté, M.A., Rheault, F., Sidhu, J., Hau, J., Petit, L., Fortin, D., Cunanne, S., Descoteaux, M., 2018.  Recognition of white matter bundles using local and global streamline-based registration and clustering.  NeuroImage 170, 283–295.  Segmenting the Brain.
  * Gattellaro et al. [2009] Gattellaro, G., Minati, L., Grisoli, M., Mariani, C., Carella, F., Osio, M., Ciceri, E., Albanese, A., Bruzzone, M., 2009.  White matter involvement in idiopathic parkinson disease: a diffusion tensor imaging study.  American Journal of Neuroradiology 30, 1222–1226.
  * Glasser et al. [2013] Glasser, M.F., Sotiropoulos, S.N., Wilson, J.A., Coalson, T.S., Fischl, B., Andersson, J.L., Xu, J., Jbabdi, S., Webster, M., Polimeni, J.R., Van Essen, D.C., Jenkinson, M., 2013\.  The minimal preprocessing pipelines for the Human Connectome Project.  NeuroImage .
  * Goodfellow et al. [2014] Goodfellow, I.J., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., Bengio, Y., 2014\.  Generative adversarial nets, in: Advances in Neural Information Processing Systems, Neural information processing systems foundation. pp. 2672–2680.
  * Gu et al. [2019] Gu, X., Knutsson, H., Nilsson, M., Eklund, A., 2019\.  Generating diffusion mri scalar maps from t1 weighted images using generative adversarial networks, in: Felsberg, M., Forssén, P.E., Sintorn, I.M., Unger, J. (Eds.), Image Analysis, Springer International Publishing. pp. 489–498.
  * Huang et al. [2004] Huang, H., Zhang, J., Van Zijl, P.C., Mori, S., 2004.  Analysis of noise effects on DTI-based tractography using the brute-force and multi-ROI approach.  Magnetic Resonance in Medicine 52, 559–565.
  * Huang et al. [2007] Huang, J., Friedland, R., Auchus, A., 2007.  Diffusion tensor imaging of normal-appearing white matter in mild cognitive impairment and early alzheimer disease: preliminary evidence of axonal degeneration in the temporal lobe.  American journal of neuroradiology 28, 1943–1948.
  * Huang and Van Gool [2017] Huang, Z., Van Gool, L., 2017\.  A riemannian network for SPD matrix learning.  31st AAAI Conference on Artificial Intelligence, AAAI 2017 .
  * Huang et al. [2019] Huang, Z., Wu, J., Van Gool, L., 2019.  Manifold-Valued Image Generation with Wasserstein Generative Adversarial Nets.  Proceedings of the AAAI Conference on Artificial Intelligence .
  * Ionescu et al. [2015] Ionescu, C., Vantzos, O., Sminchisescu, C., 2015.  Matrix backpropagation for deep networks with structured layers, in: Proceedings of the IEEE International Conference on Computer Vision (ICCV), pp. 2965–2973.
  * Jeurissen et al. [2019] Jeurissen, B., Descoteaux, M., Mori, S., Leemans, A., 2019\.  Diffusion mri fiber tractography of the brain.  NMR in Biomedicine 32, e3785.
  * Jiang et al. [2006] Jiang, H., van Zijl, P.C., Kim, J., Pearlson, G.D., Mori, S., 2006.  Dtistudio: Resource program for diffusion tensor computation and fiber bundle tracking.  Computer Methods and Programs in Biomedicine 81, 106–116.
  * Kantarci et al. [2017] Kantarci, K., Murray, M.E., Schwarz, C.G., Reid, R.I., Przybelski, S.A., Lesnick, T., Zuk, S.M., Raman, M.R., Senjem, M.L., Gunter, J.L., Boeve, B.F., Knopman, D.S., Parisi, J.E., Petersen, R.C., Jack, C.R., Dickson, D.W., 2017\.  White-matter integrity on DTI and the pathologic staging of Alzheimer’s disease.  Neurobiology of Aging .
  * Kelly et al. [2018] Kelly, S., et al., 2018.  Widespread white matter microstructural differences in schizophrenia across 4322 individuals: Results from the ENIGMA Schizophrenia DTI Working Group.  Molecular Psychiatry .
  * Kingma and Ba [2015] Kingma, D.P., Ba, J.L., 2015\.  Adam: A method for stochastic optimization, in: 3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings, International Conference on Learning Representations, ICLR.
  * Lan et al. [2020] Lan, H., Toga, A., Sepehrband, F., 2020.  SC-GAN: 3D self-attention conditional GAN with spectral normalization for multi-modal neuroimaging synthesis.  bioRxiv .
  * Ledig et al. [2017] Ledig, C., Theis, L., Huszár, F., Caballero, J., Cunningham, A., Acosta, A., Aitken, A., Tejani, A., Totz, J., Wang, Z., Shi, W., 2017.  Photo-realistic single image super-resolution using a generative adversarial network, in: Proceedings - 30th IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017.
  * Lei et al. [2019] Lei, Y., Harms, J., Wang, T., Liu, Y., Shu, H.K., Jani, A.B., Curran, W.J., Mao, H., Liu, T., Yang, X., 2019.  MRI-only based synthetic CT generation using dense cycle consistent generative adversarial networks.  Medical Physics 46, 3565–3581.
  * Li and Wand [2016] Li, C., Wand, M., 2016.  Precomputed real-time texture synthesis with markovian generative adversarial networks, in: Leibe, B., Matas, J., Sebe, N., Welling, M. (Eds.), Computer Vision – ECCV 2016, Springer International Publishing. pp. 702–716.
  * Lu et al. [2019] Lu, G., Zhou, Z., Song, Y., Ren, K., Yu, Y., 2019\.  Guiding the one-to-one mapping in cyclegan via optimal transport.  Proceedings of the AAAI Conference on Artificial Intelligence .
  * Mao et al. [2017] Mao, X., Li, Q., Xie, H., Lau, R.Y., Wang, Z., Paul Smolley, S., 2017.  Least squares generative adversarial networks, in: Proceedings of the IEEE International Conference on Computer Vision (ICCV), pp. 2794–2802.
  * Neuner et al. [2010] Neuner, I., Kupriyanova, Y., Stöcker, T., Huang, R., Posnansky, O., Schneider, F., Tittgemeyer, M., Shah, N.J., 2010\.  White-matter abnormalities in Tourette syndrome extend beyond motor pathways.  NeuroImage .
  * Oouchi et al. [2007] Oouchi, H., Yamada, K., Sakai, K., Kizu, O., Kubota, T., Ito, H., Nishimura, T., 2007.  Diffusion Anisotropy Measurement of Brain White Matter Is Affected by Voxel Size: Underestimation Occurs in Areas with Crossing Fibers.  American Journal of Neuroradiology .
  * Poot et al. [2013] Poot, D.H., Jeurissen, B., Bastiaensen, Y., Veraart, J., Van Hecke, W., Parizel, P.M., Sijbers, J., 2013.  Super-resolution for multislice diffusion tensor imaging.  Magnetic Resonance in Medicine .
  * Rheault et al. [2020] Rheault, F., et al., 2020.  Tractostorm: The what, why, and how of tractography dissection reproducibility.  Human Brain Mapping 41, 1859–1874.
  * Sánchez and Vilaplana [2018] Sánchez, I., Vilaplana, V., 2018\.  Brain mri super-resolution using 3d generative adversarial networks.  arXiv preprint arXiv:1812.11440 .
  * Schilling et al. [2020] Schilling, K.G., et al., 2020.  Tractography dissection variability: what happens when 42 groups dissect 14 white matter bundles on the same dataset?  bioRxiv .
  * Shi et al. [2016] Shi, F., Cheng, J., Wang, L., Yap, P.T., Shen, D., 2016\.  Super-resolution reconstruction of diffusion-weighted images using 4D low-rank and total variation, in: Mathematics and Visualization, Springer Heidelberg. pp. 15–25.
  * Son et al. [2019] Son, S.J., yong Park, B., Byeon, K., Park, H., 2019\.  Synthesizing diffusion tensor imaging from functional MRI using fully convolutional networks.  Computers in Biology and Medicine 115, 103528.
  * Sotiropoulos et al. [2013] Sotiropoulos, S.N., et al., 2013.  Advances in diffusion MRI acquisition and processing in the Human Connectome Project.  NeuroImage .
  * Srivastava et al. [2007] Srivastava, A., Jermyn, I., Joshi, S., 2007.  Riemannian analysis of probability density functions with applications in vision, in: Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, pp. 1–8.
  * Theaud et al. [2020] Theaud, G., Houde, J.C., Boré, A., Rheault, F., Morency, F., Descoteaux, M., 2020\.  Tractoflow: A robust, efficient and reproducible diffusion mri pipeline leveraging nextflow and singularity.  NeuroImage 218, 116889\.
  * Thomas et al. [2014] Thomas, C., Frank, Q.Y., Irfanoglu, M.O., Modi, P., Saleem, K.S., Leopold, D.A., Pierpaoli, C., 2014.  Anatomical accuracy of brain connections derived from diffusion mri tractography is inherently limited.  Proceedings of the National Academy of Sciences 111, 16574–16579.
  * Tournier et al. [2019] Tournier, J., Smith, R., Raffelt, D., Tabbara, R., Dhollander, T., Pietsch, M., Christiaens, D., Jeurissen, B., Yeh, C.H., Connelly, A., 2019\.  Mrtrix3: A fast, flexible and open software framework for medical image processing and visualisation.  bioRxiv .
  * Tuch [2004] Tuch, D.S., 2004.  Q-ball imaging.  Magnetic Resonance in Medicine .
  * Van Essen et al. [2013] Van Essen, D.C., Smith, S.M., Barch, D.M., Behrens, T.E., Yacoub, E., Ugurbil, K., 2013\.  The WU-Minn Human Connectome Project: An overview.  NeuroImage .
  * Yap et al. [2014] Yap, P.T., An, H., Chen, Y., Shen, D., 2014.  Fiber-driven resolution enhancement of diffusion-weighted images.  NeuroImage .
  * Yeh et al. [2018] Yeh, F.C., Panesar, S., Fernandes, D., Meola, A., Yoshino, M., Fernandez-Miranda, J.C., Vettel, J.M., Verstynen, T., 2018\.  Population-averaged atlas of the macroscale human structural connectome and its network topology.  NeuroImage 178, 57–68.
  * Zhao et al. [2016] Zhao, H., Gallo, O., Frosio, I., Kautz, J., 2016\.  Loss Functions for Image Restoration With Neural Networks.  IEEE Transactions on Computational Imaging .
  * Zhong et al. [2020] Zhong, J., Wang, Y., Li, J., Xue, X., Liu, S., Wang, M., Gao, X., Wang, Q., Yang, J., Li, X., 2020\.  Inter-site harmonization based on dual generative adversarial networks for diffusion tensor imaging: Application to neonatal white matter development.  BioMedical Engineering Online 19\.
  * Zhu et al. [2017] Zhu, J.Y., Park, T., Isola, P., Efros, A.A., 2017\.  Unpaired Image-to-Image Translation Using Cycle-Consistent Adversarial Networks.  Proceedings of the IEEE International Conference on Computer Vision .

##  Appendix A Backpropagation for Diffusion Tensors Learning

The $\log_{\mathrm{Id}}(\,\cdot\,)$ and $\exp_{\mathrm{Id}}(\,\cdot\,)$ maps
ensure that our generator network synthesizes valid diffusion tensors. Both
maps involves the EIG operator on $3\times 3$ symmetric matrices whose
gradients must be defined to train the network with the standard
backpropagation algorithm. Employing the matrix generalization of
backpropagation [Ionescu et al., 2015, Huang and Van Gool, 2017], we define
the partial derivatives of the objective function at the $k^{th}$ layer with
respect to a generated diffusion tensor $\mathbf{M}$ as follows:

$\footnotesize\forall\mathbf{M}_{k-1}\in\mathbf{X_{k-1}},\frac{\partial{L^{(k)}}}{\partial{\mathbf{M}_{k-1}}}=2\mathbf{U}\bigg{(}{\mathbf{K}}^{\top}\circ{\bigg{(}{\mathbf{U}}^{\top}\frac{\partial{L^{(k^{{}^{\prime}})}}}{\partial{\mathbf{U}}}\bigg{)}}_{sym}\bigg{)}{\mathbf{U}}^{\top}\bigg{)}+\mathbf{U}{\bigg{(}\frac{\partial{L^{(k^{{}^{\prime}})}}}{\partial{\mathbf{\Sigma}}}\bigg{)}}_{diag}{\mathbf{U}}^{\top},$
(24)

with

$\mathbf{M}_{k-1}=\mathbf{U}\boldsymbol{\Upsigma}{\mathbf{U}}^{\top},$ (25)

and

$\mathbf{K}(i,j)=\left\\{\begin{matrix}\frac{1}{\sigma_{i}-\sigma_{j}}&i\neq
j\\\ 0&i=j,\end{matrix}\right.$ (26)

where ${\mathbf{A}}_{sym}=\frac{1}{2}(\mathbf{A}+{\mathbf{A}}^{\top})$ and
${\mathbf{A}}_{diag}$ is $\mathbf{A}$ with nonzero elements only in its
diagonal. Furthermore, $\mathbf{X_{k-1}}\in{}^{B\times 9\times D\times H\times
W}$ is a batch of 3D DT patches with size $D\times H\times W$ and $\mathbf{K}$
is built upon the eigenvalues $\sigma$ in $\mathbf{\Sigma}$ of $\mathbf{M}$.

The variations $\frac{\partial{L^{(k^{{}^{\prime}})}}}{\partial{\mathbf{U}}}$
and $\frac{\partial{L^{(k^{{}^{\prime}})}}}{\partial{\mathbf{\Sigma}}}$ in Eq.
(24) are defined for every tensor $\mathbf{M}$ for both maps
$\log_{\mathrm{Id}}(\cdot)$ and $\exp_{\mathrm{Id}}(\cdot)$ respectively as:

$\frac{\partial{L^{(k^{{}^{\prime}})}}}{\partial{\mathbf{U}}}=2{\bigg{(}\frac{\partial{L^{(k+1)}}}{\partial{\mathbf{M_{k}}}}\bigg{)}}_{sym}\mathbf{U}\log(\boldsymbol{\Upsigma}),$
(27)
$\frac{\partial{L^{(k^{{}^{\prime}})}}}{\partial{\mathbf{\Sigma}}}=\boldsymbol{\Upsigma}^{-1}{\mathbf{U}}^{\top}{\bigg{(}\frac{\partial{L^{(k+1)}}}{\partial{\mathbf{M_{k}}}}\bigg{)}}_{sym}\mathbf{U},$
(28)

and

$\frac{\partial{L^{(k^{{}^{\prime}})}}}{\partial{\mathbf{U}}}=2{\bigg{(}\frac{\partial{L^{(k+1)}}}{\partial{\mathbf{M_{k}}}}\bigg{)}}_{sym}\mathbf{U}\exp(\boldsymbol{\Upsigma}),$
(29)
$\frac{\partial{L^{(k^{{}^{\prime}})}}}{\partial{\mathbf{\Sigma}}}=\exp(\boldsymbol{\Upsigma}){\mathbf{U}}^{\top}{\bigg{(}\frac{\partial{L^{(k+1)}}}{\partial{\mathbf{M_{k}}}}\bigg{)}}_{sym}\mathbf{U}.$
(30)

With the help of Eq. (24) and the successive application of the
$\exp_{\mathrm{Id}}(\cdot)$ and $\log_{\mathrm{Id}}(\cdot)$ maps, we make sure
that the generated tensors are and stay SPD throughout the training. An open-
source implementation of these differentiable functions is available on our
Github 222https://github.com/banctilrobitaille/torch-
vectorized333https://torch-vectorized.readthedocs.io/en/latest/.

##  Appendix B Tractography Assessment

This appendix presents additional figures, referenced in Section 3 and in
Section 4, that compare the generated and the expected real tractograms.
First, we compare in Figure 13 the mean volume, streamline length, Dice, OL
and OR when considering the whole brain tractograms.

Figure 13: In A), the mean streamline length (in mm). In B), the mean volume
(in voxels) and in C), the mean Dice, OL and OR.

We then compare the mean volume of each recovered bundle in Figure 14 below:

Figure 14: Top row: Mean bundle volumes (in voxels) from real and generated
DT. Bottom row: Mean bundle volumes (in voxels) from real and generated ODF.

Finally, Figure 15 presents the Dice, OL and OR yielded by each recovered
bundles.

Figure 15: Top row: Dice, OL and OR between bundles and tractograms from real
and generated DT. Bottom row: Dice, OL and OR between bundles and tractograms
from real and generated ODF. The lines indicate the mean metric values for all
bundles.

##  Appendix C Visual Assessment of Synthesized Volumes

We provide, in the following appendix, the FA and color FA maps of five
additional test subjects. We first compare in Figure 16 the generated and real
FA of a sagittal, axial and coronal slice. We then compare in Figure 17 the
generated and real color FA of the same five subjects. One can appreciate how
our method adapts to different brain geometries and generates plausible FA and
principal orientations despite the differences in the anatomy of the subjects.

Figure 16: The real and synthesized FA maps of five test subjects in the
sagittal, axial and coronal planes. Figure 17: The real and synthesized color
FA maps of five test subjects in the sagittal, axial and coronal planes.