# Analytical and Fast Fiber Orientation Distribution Reconstruction in
3D-Polarized Light Imaging

Abib Alimi
corresponding author
Athena Project-Team
Inria Sophia Antipolis-Méditerranée
Université Côte d’Azur, France
<EMAIL_ADDRESS>&Samuel Deslauriers-Gauthier
Athena Project-Team
Inria Sophia Antipolis-Méditerranée
Université Côte d’Azur, France &Felix Matuschke
Institute of Neuroscience and Medicine (INM-1)
Research Center Jülich, Germany &Andreas Müller
Simulation Lab Neuroscience
Jülich Supercomputing Centre
Institute for Advanced Simulation
JARA
Research Center Jülich, Germany &Sascha E. A. Muenzing
Institute of Neuroscience and Medicine (INM-1)
Research Center Jülich, Germany &Markus Axer
Institute of Neuroscience and Medicine (INM-1)
Research Center Jülich, Germany &Rachid Deriche
Athena Project-Team
Inria Sophia Antipolis-Méditerranée
Université Côte d’Azur, France

###### Abstract

Three dimensional Polarized Light Imaging (3D-PLI) is an optical technique
which allows mapping the spatial fiber architecture of fibrous postmortem
tissues, at sub-millimeter resolutions. Here, we propose an analytical and
fast approach to compute the fiber orientation distribution (FOD) from high-
resolution vector data provided by 3D-PLI. The FOD is modeled as a sum of $K$
orientations/Diracs on the unit sphere, described on a spherical harmonics
basis and analytically computed using the spherical Fourier transform.
Experiments are performed on rich synthetic data which simulate the geometry
of the neuronal fibers and on human brain data. Results indicate the
analytical FOD is computationally efficient and very fast, and has high
angular precision and angular resolution. Furthermore, investigations on the
right occipital lobe illustrate that our strategy of FOD computation enables
the bridging of spatial scales from microscopic 3D-PLI information to macro-
or mesoscopic dimensions of diffusion Magnetic Resonance Imaging (MRI), while
being a means to evaluate prospective resolution limits for diffusion MRI to
reconstruct region-specific white matter tracts. These results demonstrate the
interest and great potential of our analytical approach.

_K_ eywords fiber orientation distribution $\cdot$ 3D-PLI $\cdot$ polarized
light imaging $\cdot$ diffusion MRI

##  1 Introduction

To improve our understanding of the human brain structure and function as well
as its normal development and complex disorders, it is compulsory to develop
and combine different imaging techniques at multiple spatial scales. With a
millimeter resolution, diffusion Magnetic Resonance Imaging (MRI) (35, 11, 55,
28, 23, 50, 36) is the only in vivo and non invasive imaging technique able to
assess the structural architecture of human fibrous organs such as the brain
and the heart and to infer their respective connectivity (41, 18, 27, 26).
However, in vivo diffusion MRI suffers from a lack of ground truth (38) and a
relatively low spatial resolution. Therefore, 3D-Polarized Light Imaging
(3D-PLI) (8, 29, 34, 44, 4) has been investigated and, thanks to its high
spatial resolution, it is presented as a complementary and potential technique
for the validation and guidance of diffusion MRI fiber orientation estimates
and tractography.

3D-PLI is an optical imaging technique that utilizes the birefringence
property of brain tissue to map its spatial organization at micrometer
resolution (8, 29). The term ‘3D’ in this context is first of all related to
the within-section 3D fiber orientation estimation, and only secondly in the
sense of 3D volume reconstruction of serial sections. The technique, thus,
provides high-resolution 3D microscopic fiber orientation measurements
recovered from each unstained histological tissue section. Complementing
diffusion MRI with 3D-PLI measurements, therefore, requires to bridge the
spatial scale information from micro to millimeter dimensions. Thus, the fiber
orientation distribution (FOD), which is widely used in diffusion community
(52, 1, 51, 14, 15), is a suitable tool to compare diffusion MRI and 3D-PLI
fiber orientation estimates in 3D. For validation purpose, 3D-PLI has been
used to reconstruct brain white matter (WM) fiber pathways and presented as an
adequate means to validate diffusion MRI fiber tractography (3). Moreover,
(13, 60, 26) used high-resolution 3D-PLI, or more specifically the generated
Fiber Orientation Map (FOM) which reflects the determined 3D fiber courses
best, in order to identify and to follow the tiny fiber tracts to complement
diffusion MRI results from the same region, rather than addressing long-
distance tractography in both modalities. Concerning the brain’s fiber
architecture, histology-based studies compared FODs reconstructed from
diffusion MRI to those estimated either in 2D from stained tissue sections
(37, 12), or more recently, in 3D via optical coherence tomography (57) and
confocal microscopy (31, 48). However, to our knowledge, 3D-PLI has not yet
been used to evaluate 3D fiber orientation distributions obtained from
diffusion MRI, and building PLI-based FODs paves the way for it.

On this account, (10) proposed an estimate of the FOD derived from high-
resolution 3D-PLI data, called pliODF. This technique enables the bridging of
spatial scales from micro to macroscopic dimensions, for instance, to the
millimeter resolution of diffusion MRI. pliODF benefits from the important
concept of super-voxel (SV) which allows downsampling of the high-resolution
data, and then estimating the FOD on the unit sphere by expanding a normalized
directional histogram with a series of spherical harmonics (SH) functions. By
changing the SV size, the FOD can be computed at different spatial scales (10,
5, 2). (10) refer to a tissue voxel containing a single high-resolution fiber
orientation as ‘native voxel’ and the super-voxel regroups $K=n_{x}\times
n_{y}\times n_{z}$ native voxels. The pliODF, though, is derived from an
empirically parameterized directional histogram, which makes it a discrete
method. It relies on a histogram binning procedure which can introduce angular
errors in the recovered FODs (5), and since it is represented on a SH basis,
it can also limit the order of its SH expansion (6). In order to overcome
these issues, (6) proposed to regularize the pliODF using the Laplace-Beltrami
operator and therefore improved its results.

In this work, we present an analytical procedure to determine the FOD from
3D-PLI data in order to assess the spatial distribution of fiber orientations
within a super-voxel. Our method is based on defining each fiber orientation
in a native voxel as a 2D Dirac delta function on the unit sphere and the FOD
as a sum of these Diracs within each super-voxel. The FOD is continuously
described on a spherical harmonics basis and analytically computed via the
spherical Fourier transform and by means of the Diracs. We remind that the FOD
is the spatial distribution of the fiber orientations within an image voxel as
known in diffusion MRI (52, 15). In the context of 3D-PLI, (10) introduced a
discrete approximate of the FOD called the ‘pliODF’ and, here, we propose an
analytical FOD denoted ‘aFOD’.

This paper reviews and extends our previous works in which the analytical
method was introduced (5) and its performance evaluated (2). Here, we give (1)
a full description of our method to analytically compute the aFOD from high-
resolution 3D-PLI vector data, then (2) the used datasets and the different
evaluation criteria are presented for a further/extended comparison with
pliODF. Results demonstrate (3) the high angular resolution and precision of
the aFOD and the dependence of pliODF on its empirical histogram binning to
perform similarly, and (4) that the aFOD is computationally very efficient and
very fast. We finally (5) compute aFODs on human brain data at different
spatial scales while preserving the original high-resolution organization of
WM fiber bundles and (6) present it as a means for diffusion MRI to assess
potential resolution limits to track specific fascicles in the human brain.

##  2 Theory

The fiber orientation distribution provides a comprehensive statistical
description of integrated 3D-PLI vector data (10, 5, 2) and positions itself
as a suitable tool for multimodal analysis and comparisons (48). As a function
defined on the unit sphere $\mathbb{S}^{2}$, the FOD can be expanded as a
linear combination of spherical harmonics. In this section, we give a brief
background information on spherical harmonics before presenting the pliODF
estimate by (10) and the analytical approach –aFOD– we propose.

###  2.1 Spherical harmonics

The spherical harmonics of order l and phase factor m are defined as

$Y_{l}^{m}(\theta,\phi)=N_{l}^{m}P_{l}^{m}(\cos\theta)e^{jm\phi}$ (1)

where $N_{l}^{m}$ is a normalization coefficient, $P_{l}^{m}$ are the
associated Legendre polynomials, the angles $\theta\in[0,\pi]$ and
$\phi\in[0,2\pi)$ are the colatitude and azimuth, respectively. The spherical
harmonics $\\{Y_{l}^{m}:-l\leqslant m\leqslant l,l=0,1,...\\}$ form an
orthonormal basis over $L^{2}(\mathbb{S}^{2})$ and any square integrable
function $f(\theta,\phi)\in L^{2}(\mathbb{S}^{2})$ can be expressed as

$f(\theta,\phi)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}c_{lm}Y_{l}^{m}(\theta,\phi)$
(2)

where $c_{lm}$ are the SH coefficients of $f$. These coefficients are obtained
by computing the spherical Fourier transform (25) of $f$ defined as

$c_{lm}=\int_{\mathbb{S}^{2}}f({\bf w})\overline{Y_{l}^{m}}({\bf
w})d\mathbf{w}$ (3)

for ${\bf w}\in\mathbb{S}^{2}$ and the over-bar denotes the conjugation. The
unit vector ${\bf w}$ is related to $(\theta,\phi)$ by

${\bf
w}(\theta,\phi)=[\sin\theta\cos\phi\quad\sin\theta\sin\phi\quad\cos\theta]^{T}$
(4)

where $\cdot^{T}$ indicates transposition. In this paper, both functions
$f(\theta,\phi)$ and $f(\bf w)$ are assumed to be equal and written
interchangeably to simplify notations.

###  2.2 pliODF, a discrete solution

(10) proposed an estimate of the FOD from high-resolution 3D-PLI vector data,
which they called pliODF. For this, the collection of all native voxels, which
represents the fiber orientation map (FOM) (8, 29), is downsampled into cubic
compartments or super-voxels (SV) containing $n_{x}\times n_{y}\times n_{z}$
native voxels each. This allows to determine the spatial distribution of 3D
fiber orientations in a given SV. Thus, starting from high-resolution FOMs, as
illustrated in Fig. (3) in (10), their approach is implemented as follows:

  1. 1.

Define super-voxels of the FOMs measured from unstained histological tissue
sections.

  2. 2.

Create a normalized directional histogram on the unit sphere from discretized
distribution of fiber orientation vectors in each super-voxel.

  3. 3.

Approximate the directional histogram using spherical harmonics, see Eq. (2).

In practice, the SH approximation is truncated at a maximum order or bandlimit
$L_{max}$. Moreover, in the original implementation of the pliODF technique,
the SH coefficients are determined with a linear least square method by
minimizing $||\mathbf{h}-\mathbf{B}\mathbf{c}||^{2}$, yielding

$\mathbf{c}=(\mathbf{B}^{T}\mathbf{B})^{-1}\mathbf{B}^{T}\mathbf{h}$ (5)

where, in matrix format, $\mathbf{B}$ is the basis matrix of the truncated SH,
$\mathbf{h}$ the vector entries of the directional histogram, and $\mathbf{c}$
the coefficients vector of the expansion up to order $L_{max}$.

This technique benefits not only from the 3D nature of 3D-PLI dataset but also
from the super-voxel which allows the estimation of the FODs at different
spatial scales. Nevertheless, the reconstructed pliODF is an empirical and
discrete approximation of the FOD. Indeed, it is a directional histogram
defined by its bin centers and the number of 3D-PLI fiber orientation vectors
falling in each bin is dependent on the bin size or corresponding dihedral
angle. Therefore, the $L_{max}$ of the SH expansion of the pliODF can be
limited as shown in (6).

###  2.3 Analytical FOD in 3D-PLI

We demonstrate, in this section, that without any discretization which leads
to loss of orientation information, the aFOD can be analytically reconstructed
from high-resolution 3D-PLI orientation measurements. We show how we model our
aFOD before we obtain the exact SH coefficients that uniquely characterize it.

####  2.3.1 aFOD as sum of k Diracs

In each voxel of the fiber orientation map (visualizing the dominant fiber
orientation), we propose to describe the unit fiber orientation vector as a
two dimensional Dirac delta function $\delta$ on the sphere. Therefore, in a
super-voxel containing $K$ orientations $(\theta_{k},\phi_{k})$, the fiber
orientation distribution function $f$ can be modeled as the sum of $K$ Diracs
(20, 19), that is

$f(\theta,\phi)=\frac{1}{K}\sum_{k=1}^{K}\delta(\cos\theta-\cos\theta_{k})\delta(\phi-\phi_{k})$
(6)

Note that $f(\theta,\phi)$ is completely defined by parameters
$(\theta_{k},\phi_{k})$ of the $K$ orientations located in the SV. In the next
step, the coefficients of the SH expansion of $f$ are determined.

####  2.3.2 Computation of the SH coefficients

Unlike (10) who used a linear least square method to get the coefficients of
the SH expansion, we directly recover the coefficients that uniquely define
our aFOD $f$ by computing the spherical Fourier transform through Eq. (3)
since we know the parameters $(\theta_{k},\phi_{k})$ of all fiber orientation
vectors determined from 3D-PLI analysis (8, 29, 4) , in a given super-voxel.
That is,

$\begin{split}c_{lm}&=\int_{\mathbb{S}^{2}}f({\bf w})\overline{Y_{l}^{m}}({\bf
w})d\mathbf{w}\\\
&=\int_{0}^{2\pi}\int_{0}^{\pi}f(\theta,\phi)\overline{Y_{l}^{m}}(\theta,\phi)\sin\theta
d\theta d\phi\end{split}$

then $f$ is replaced by the collection of Diracs in Eq. (6) and $Y_{l}^{m}$ by
its expression in Eq. (1) to get

$\begin{split}c_{lm}&=\frac{N_{l}^{m}}{K}\sum_{k=1}^{K}\int_{-1}^{1}\delta(\cos\theta-\cos\theta_{k})P_{l}^{m}(\cos\theta)d\cos\theta\\\
&\times\int_{0}^{2\pi}\delta(\phi-\phi_{k})e^{-jm\phi}d\phi\end{split}$

and finally applying the sifting property of the Diracs gives

$\displaystyle c_{lm}$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}N_{l}^{m}P_{l}^{m}(\cos\theta_{k})e^{-jm\phi_{k}}$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\overline{Y_{l}^{m}}(\theta_{k},\phi_{k})$
(7)

This solution for the $c_{lm}$ is the analytical and continuous equivalent of
the discrete one approximated for pliODF. Eq. (7) yields, with respect to
$L_{max}$, the exact solution which requires no directional histogram and,
therefore, all the information from the high-resolution fiber orientation
vectors are taken into account in the FOD computation. Furthermore, it can be
evaluated at any $L_{max}$.

Using Eq. (2), consequently, the analytical FOD $\mathbf{f}$ is simply
computed through matrix-vector multiplication:

$\mathbf{f}=\mathbf{B}\mathbf{c}$ (8)

with $\mathbf{B}$ the SH basis matrix and $\mathbf{c}$ the vector of the
analytically recovered SH coefficients.

Note that although our reconstruction approach uses the SH to represent the
aFOD, it is independent of the SH basis used and here the real and symmetric
basis defined in (17) is chosen. In this case, the number of terms of the SH
basis of order $l$ is $L=(1/2)(l+1)(l+2)$, where only the even orders are
considered.

##  3 Methods

Figure 1: (A) 3D visualization of the simulated fiber bundles crossing at
angle $X~{}=~{}90^{\circ}$ by the FAConstructor (45). (B)-(D) The resulting
Fiber Orientation Maps of 3D-PLI simulations through $60\text{ }\mu m$
sections with a crossing at $X=30^{\circ},60^{\circ}$ and $90^{\circ}$,
respectively. The 3D orientations of entire fiber structures (A) and
determined local fiber orientations (B)-(D) are color coded as indicated by
the color sphere at the bottom right.

###  3.1 Synthetic data generation

The simulated 3D-PLI dataset consists of two 3D fiber bundles intersecting at
an angle $X$ in the XY-plane. A 3D collision solving algorithm, described in
(39), is performed to ensure that the individual fibers do not overlap. For
this purpose, the fibers are divided into short segments so that two segments
can be checked for a collision. In the event of a collision, the two segments
are slowly pushed apart. This happens simultaneously for all colliding objects
and is repeated until there is no more collisions. Fig 1A shows an example of
the resulting fiber configuration for the $X=90^{\circ}$ crossing case. Note
that $X$ is varied from $0^{\circ}$ to $90^{\circ}$ with $10^{\circ}$
increments to give $10$ different crossing situations. More details on the
fiber generation can be found in A.

The 3D models generated in this way are then used for the simulation using
SimPLI tool (22, 40) to produce synthetic 3D-PLI datasets. The noise of the
CCD camera is modelled as in (49) i.e. with a negative binomial distribution
in which the variance $\sigma$ depends on the expected light intensity $\mu$
by $\sigma=3\mu$. From the generated 3D-PLI signal, high-resolution
orientations are reconstructed based on the approaches in (8, 49). The
generated data consists of 10 slices of $60\text{ }\mu m$ thick, with a pixel
resolution of $64\times 64\text{ }\mu m^{2}.$

###  3.2 Experimental brain data

The post mortem human brain tissue used for this study was acquired in
accordance with the local ethic committee of our partner university at the
Heinrich Heine University Düsseldorf. As confirmed by the ethic committee,
postmortem human brain studies do not need any additional approval, if a
written informed consent of the subject is available. For the research carried
out here, this consent is available.

The brain was removed $24\text{ }h$ after the donor’s death (male, age: 54
years, suffered from multiple system atrophy). The right hemisphere was fixed
in $4\%$ buffered formaldehyde solution for 4 years. After immersion in a
$20\%$ solution of glycerin with Dimethyl sulfoxide (DMSO) for cryoprotection,
the brain was deep frozen and cryo-sectioned at $70\text{ }\mu m$ thickness
(with Polycut CM 3500, Leica, Germany) along the coronal plane, orthogonal to
the line connecting the anterior and posterior commissures determined from a
post mortem turbo spin-echo MRI scan pre-sectioning. Before each sectioning
step, en face images (referred to as blockface images) were taken from the
frozen surface of the brain block to generate reference images used for serial
section image alignment. Each section was finally mounted unstained on a glass
slide, immersed in a $20\%$ solution of glycerin to avoid dehydration, and
cover-slipped.

50 consecutive sections were selected for the present study and measured with
the Large-area Polarimeter (LAP) setup with tilting stage, as described in
(9). Pixel size of all images was set to $64\times 64\text{ }\mu m^{2}$. The
experimental protocol provided one planar measurement and four tilting
position measurements (in north, west, east, and south direction with $8\text{
}^{\circ}$ tilting angle), and each measurement was composed of a set of 18
images taken at $10\text{ }^{\circ}$ equidistant polarization filter rotation
angles. Image analysis included calibration, segmentation, and tilting
analysis to determine the local fiber structure orientations, utilizing the
analysis method described in (49).

Three-dimensional reconstruction of the brain slab volume consisting of the 50
sections was achieved by a multi-step registration approach: Firstly, the
3D-PLI images were aligned to their corresponding blockface images via affine
transformations utilizing in-house developed software based on the software
packages ITK (59) and elastix (32). Secondly, non-linear registration using
the ANTs toolkit (7) was performed. Finally, all 3D-PLI modalities were
spatially transformed using the obtained deformation fields and orientations
were corrected accordingly.

One tissue section was measured with the Polarizing Microscope (PM) setup
without tilting stage, as described in (44). Pixel size of was set to
$1.3\times 1.3\text{ }\mu m^{2}$. The experimental protocol provided just the
planar measurement with each measurement being composed of a set of 18 images
taken at $10\text{ }^{\circ}$ equidistant polarization filter rotation angles.
Image analysis included calibration, segmentation, and flat image analysis to
determine the local fiber structure orientations, utilizing the analysis
method described in (44).

Table 1: Imaging spatial resolutions resulting from super-voxels. Super-voxel size | Resolution [$mm^{3}$]
---|---
1 $\times$ 1 $\times$ 1 | .064 $\times$ .064 $\times$ .07
2 $\times$ 2 $\times$ 2 | .128 $\times$ .128 $\times$ .14
5 $\times$ 5 $\times$ 5 | .32 $\times$ .32 $\times$ .35
10 $\times$ 10 $\times$ 10 | .64 $\times$ .64 $\times$ .7
25 $\times$ 25 $\times$ 25 | 1.6 $\times$ 1.6 $\times$ 1.75
50 $\times$ 50 $\times$ 50 | 3.2 $\times$ 3.2 $\times$ 3.5

###  3.3 Evaluations

To assess the performance of aFOD and pliODF reconstruction approaches and see
the relationship between them, we do tests on synthetic and on experimental
human brain 3D-PLI data.

####  3.3.1 Angular resolution

In diffusion MRI, it is generally assumed that the local maxima or peaks of
the FOD coincide with the fiber directions (52, 54, 17, 14). Here, the angular
resolution corresponds to the critical angle $\tau$ under which the peaks of
the two crossing fiber bundles are confounded, that is, they can not be
detected separately. To find the peaks of the FOD, the
sh2peaks111https://www.mrtrix.org algorithm implemented in Mrtrix3 (53) is
used, with the minimum peak amplitude set to 0.5 and the other parameters to
default.

####  3.3.2 Angular precision

Let ‘single fiber’ denote a bundle of coherently organized/oriented fibers.
The angular precision $\epsilon$ is the deviation or angular error of the peak
of the reconstructed FOD from the simulated ground-truth. Here we investigate
the $\epsilon$ of both aFOD and pliODF, as well as the effect of the
directional histogram bin width on $\epsilon$. The angular precision
$\epsilon$ is evaluated in degrees as follows,

$\epsilon=\frac{180}{\pi M}\sum_{m=1}^{M}\arccos({\bf v}_{m}^{T}{\bf w}_{m})$
(9)

which is the average distance between ${\bf v}_{m}$ the $m^{th}$ ground truth
orientation and ${\bf w}_{m}$ the closest of the recovered peak orientations
of the FOD and $M$ is the number of fibers in the super-voxel.

####  3.3.3 Computational runtimes

In order to ensure a fair comparison and due to the large size of 3D-PLI human
brain data, both reconstruction algorithms are performed on the JURECA (30)
supercomputer at Forschungszentrum Jülich, in Germany. The supercomputer
consists of 1872 compute nodes and each node is equipped with two Intel Xeon
E5-2680 v3 Haswell $24$-core processors operating at $2.5$ GHz. The considered
compute nodes contain $128$ GB of main memory each, and their number is set to
$5$, for a total of $120$ cores in order to keep the compute time in an
acceptable range. The super-voxel size and the spherical harmonics bandlimit
$L_{max}$ are taken as parameters, plus the histogram bin size for the pliODF
generation.

##  4 Results

###  4.1 Single Fiber Orientation

We recall that single fiber refers to a single bundle of fibers. Fig. 2
depicts, using this simulated single fiber dataset, the FOD angular precision
$\epsilon$ from both methods as a function of the directional histogram bin
size in degrees. The red dotted curve represents the $\epsilon$ of pliODF
while the green horizontal line is the constant $\epsilon$ of aFOD. The
$\epsilon$ of discrete pliODF converges to the $\epsilon$ of aFOD as the bin
size decreases, regardless of the considered spherical harmonics $L_{max}$,
see (A) to (C). Moreover, in (D) is zoomed in the blue square in (C) to show
that with $L_{max}=12$ and small bin sizes ($<5^{\circ}$), both pliODF and
aFOD have angular precision smaller than $1^{\circ}$ and when the bin size
gets smaller, $\epsilon$ of pliODF decreases as well.

Figure 2: Angular precision $\epsilon$ of ‘single fiber’ FODs as a function of
the directional histogram bin size, in degrees. (A)–(C) For decreasing bin
size, the $\epsilon$ of pliODF (red dotted line) converges to the constant
$\epsilon<1^{\circ}$ of aFOD, in green horizontal line. At $L_{max}=12$, when
bin size is very small, for both FODs $\epsilon<1^{\circ}$, this is
illustrated in (D) where the blue square in (C) is zoomed in.

Therefore, in the rest of the experiments, in order to get results closer to
our analytical approach, bin widths of $3.6^{\circ}$ (consistent with (10))
and $1.8^{\circ}$ will the considered for the estimation of the pliDOFs.

The previous observations –and choices of bin width– are corroborated in Table
2 which quantitatively compares the $\epsilon$ of recovered aFOD and pliODF
for the single fiber population with varying SH $L_{max}$. The angular
precision of the aFOD is always $\epsilon<1^{\circ}$ regardless of $L_{max}$.
This applies to the $\epsilon$ of pliODF, except for $L_{max}=4$ when bin size
is $1.8^{\circ}$, and $L_{max}\leq 6$ for a bin width of $3.6^{\circ}$. Note,
though, that the pliODF produces more erroneous results (e.g. up to
$\epsilon>30^{\circ}$) when larger bin width are used (Fig. 2).

Angular Error $\epsilon[^{\circ}]$
---
$L_{max}$ |  |  aFOD | pliODF $(1.8^{\circ})$ | pliODF $(3.6^{\circ})$
4 |  |  0.55∘ | 1.33∘ | 2.12∘
6 |  |  0.57∘ | 0.78∘ | 1.24∘
8 |  |  0.59∘ | 0.34∘ | 0.52∘
10 |  |  0.62∘ | 0.03∘ | 0.09∘
12 |  |  0.64∘ | 0.33∘ | 0.60∘
Table 2: Single fiber angular precision $\epsilon$ in degrees [∘] as a
function of $L_{max}$. For aFODs $\epsilon$ is always $<1^{\circ}$, and the
difference with $\epsilon$ of pliODFs is not very significant with regard to
the considered bin sizes.

###  4.2 Crossing Fiber Orientations

####  4.2.1 Angular Resolution

Fig. 3 displays the angular resolution of both aFOD and pliODF relative to
$L_{max}$. It can be observed that $90^{\circ}$-crossings can be separated at
any SH $L_{max}\geq 4$, and at least $L_{max}=6$ is necessary to discriminate
$60^{\circ}$-crossings and $L_{max}=10$ for $40^{\circ}$-crossings. These
observations are valid for the aFOD as well as the pliODF with the selected
bin widths, though only results for bin size $3.6^{\circ}$ are shown for
illustration. Moreover, one can also notice shape differences between the aFOD
and the pliODF at $40$\- and $60^{\circ}$-crossings.

Figure 3: Effects of SH $L_{max}$ on the angular resolution of aFOD and pliODF
(with bin size $3.6^{\circ}$) estimated at different crossing configurations.

####  4.2.2 Angular Precision

Fig. 4 displays the angular precision $\epsilon$ of the FODs reconstructed
from crossing fiber bundles using both methods at different spherical
harmonics $L_{max}$. The green curve represents the $\epsilon$ of aFOD and the
red and black dotted curves represent the $\epsilon$ of pliODF estimated with
bin sizes of $3.6^{\circ}$ and $1.8^{\circ}$, respectively. The following
observations can be made. First, globally for all algorithms, $\epsilon$
decreases as the crossing angle $X$ increases: it goes from around $7^{\circ}$
at $X=40^{\circ}$ and $L_{max}=10$ down to less than $1.5^{\circ}$ at
$X=90^{\circ}$ and $L_{max}=6$. Second, in all combinations of $X$ and
$L_{max}$, the aFODs have a smaller angular error than pliODFs have, except
for crossing $X=40^{\circ}$ and $L_{max}=10$ as can be seen in the plot at the
right of Fig. 4. Third, as for pliODF, the smaller the bin size, the lower the
angular error. Finally, as seen qualitatively in Fig. 3 for both methods, the
fiber crossing at $X=90^{\circ}$ can be correctly reconstructed at any SH
order with a maximum error $\epsilon=1.3^{\circ}$, the $X=60^{\circ}$ is
resolved when $L_{max}\geq 6$ with maximum $\epsilon<6^{\circ}$ and, the
$X=40^{\circ}$ only when $L_{max}\geq 10$ with maximum error $\epsilon\approx
7^{\circ}$.

Figure 4: Angular precision of the FODs reconstructed at varying SH $L_{max}$
from different crossing fiber configurations $X$.

###  4.3 Computational runtime

The FODs’ computational runtime measurements on the system JURECA (30) are
presented in Table 3 in minutes:seconds. The test is performed on human brain
data consisting of $50$ coronal slices of $1950\times 1350$ native voxels
222The dataset is actually masked and only $38\%$ of the voxels are
considered. These voxels contain the brain tissue. each (see Sect. 3.2).
Concerning pliODF, the directional histogram bin size is set to $1.8^{\circ}$,
equivalent to $20000$ bins in order to get closest results to the analytical
solution. The following observations are made. First, the execution time
increases with the SH maximum order $L_{max}$ for both reconstruction
techniques as previously reported in (10, 2). Second, at high resolution, the
analytical method is faster than the pliODF technique by a factor between $24$
and $73$ at SV size of $5\times 5\times 5$ corresponding to $320\times
320\times 350\text{ }\mu m^{3}$ resolution (in Table 1). Third, when the SV
size is $>25$ isotropic native voxels, it is $3$ to $6$ times quicker to
generate pliODFs than aFODs, and both reconstruction approaches are completed
within seconds. In between these SV sizes, the analytical solution runs
faster. This could be accounted for the continuous setting of the analytical
approach, which makes it compute the SH coefficients from every native voxel
located in the super-voxel, regardless of its size. However, the discrete
pliODF approach estimates the SH coefficients from the defined directional
histogram in the considered super-voxel (10, 24). Therefore, pliODF runs
moderately faster when the SV size is large, at least $50$ isotropic native
voxels since the number of super-voxels in the whole dataset decreases i.e.
there is a few directional histograms to estimate the SH coefficients from.
Conversely, with a small SV size, the number of super-voxels as well as
histograms is large, thus the analytical technique is faster.

Table 3: A comparison of runtime measurements in minutes:seconds. All the
aFODs are run in the order of seconds regardless of the super-voxel size and
the SH $L_{max}$. pliODF takes minutes to run in small sized SV and a few
seconds when the resolution is coarser. The execution time generally increases
with $L_{max}$. Super-voxel size [in voxels]
---
| $L_{max}$ | 50 $\times$ 50 $\times$ 50 | 25 $\times$ 25 $\times$ 25 | 10 $\times$ 10 $\times$ 10 | 5 $\times$ 5 $\times$ 5
pliODF | 4 | 0:01.23 | 0:01.57 | 0:17.29 | 2:07.52
| 6 | 0:01.90 | 0:02.89 | 0:33.50 | 4:08.71
| 8 | 0:02.92 | 0:04.98 | 0:58.72 | 7:16.63
| 10 | 0:04.26 | 0:07.71 | 1:30.81 | 11:15.95
| 12 | 0:06.00 | 0:11.10 | 2:12.14 | 16:23.47
aFOD | 4 | 0:04.80 | 0:01.27 | 0:01.51 | 0:05.25
| 6 | 0:09.03 | 0:02.41 | 0:02.34 | 0:07.69
| 8 | 0:15.66 | 0:04.13 | 0:03.73 | 0:09.66
| 10 | 0:24.97 | 0:06.45 | 0:05.82 | 0:11.38
| 12 | 0:37.47 | 0:09.70 | 0:07.87 | 0:13.37

###  4.4 On human brain data

Figure 5: FOD map of U-shaped fibers connecting adjacent gyri at varying
resolutions from (B) to (D) with a high-resolution retardation map displayed
in (A). The global fiber organization is preserved through the scales. FODs
overlay a recomputed retardation map and are colour-coded according to their
local orientation (red: left-right; green: anterior-posterior; blue: inferior-
superior).

Figure 6: Stratum sagittale (SS) in deep WM described by FODs at different
super-voxel sizes i.e. at different spatial resolutions. (A) pictures a high-
resolution retardation map while (B) to (D) show the SS at SV sizes of $5$ and
$50$ isotropic native voxels. The global fiber organization is preserved
through the scales. FODs are colour-coded according to their local orientation
(red: left-right; green: anterior-posterior; blue: inferior-superior). The
green square RIO in (D) is zoomed in (E) to present a voxel containing a
crossing FOD computed at different SH $L_{max}$ and the white rectangular RIO
is zoomed in on Fig. 7.

The coronal brain slices with a pixel size of $64\text{ }\mu m$ isotropic
presented in Sect. 3.2 are used to perform the previous computational running
time test, and here, to illustrate the ability of the analytical method to
reconstruct the FODs from real data at different imaging scales while bridging
the resolution gap. The corresponding FOM contains $1350\times 1950\times 50$
native voxels 333Again, the dataset is masked..

####  4.4.1 U-shaped fibers

aFODs of U-shaped fibers also known as short arcuate or integrals (in the
occipital lobe) are pictured in Fig. 5. The subplot (A) displays the
retardation map (8) of a high-resolution section and the green square is the
region of interest (RIO). (B) to (D) overlaying recomputed retardation maps,
show the U-shaped fiber connection of neighbouring gyri at different spatial
resolutions, from $SV=5\times 5\times 5$ to $SV=25\times 25\times 25$ (see
Table 1), while preserving their original high-resolution structure. From
these same maps, one can also appreciate the fibers fanning into the gyri.

####  4.4.2 Stratum sagittale

In Fig. 6 is shown a FOD map of the stratum sagittale (SS), yellow RIO in (A).
The SS consists of longitudinal fibers (going out of the plane), lateral to
the ventricles and located in the deep WM region of the occipital lobe. It can
be reported that: 1) Going from smaller $SV=5\times 5\times 5$ (B) to larger
$SV=50\times 50\times 50$ (C,D), the spatial resolution gets ($10$ times)
coarser ($3.2\times 3.2\times 3.5\text{ }mm^{3}$, see Table 1), however, the
global WM tissue structure of the SS is conserved by both the pliODF (C) and
the aFOD (D) with crossings. 2) Overall, the shapes of the aFOD and pliODF are
very similar, although the latter sometimes presents some spurious lobes at
the center of the crossings. 3) The aFOD is more sharply defined as the
spherical harmonics $L_{max}$ increases as shown in (E). This observation is
also valid for pliODF (but not displayed, with eventual small spurious lobes).
A zoom in on the white RIO is displayed in Fig. 7.

####  4.4.3 Diffusion MRI resolution limit

Figure 7: Zoom in on the white rectangle of Fig. 6: details about the layers
of the stratum sagittale. The spatial resolution is $128\times 128\times
140\text{ }\mu m^{3}$ in (A), $320\times 320\times 350\text{ }\mu m^{3}$ in
(B) and $640\times 640\times 700\text{ }\mu m^{3}$ in (C). The blue and the
red arrows point the thinner and the more pronounced layers of the SS,
respectively. By reducing the scale, the thinner layer form crossing aFODs
with a vertical fiber bundle.

Fig. 7 displays a zoom in on the aFOD map of the stratum sagittale. At high
resolution (A), the two layers of the SS are distinguishable as shown by (56)
in their postmortem dissection of the occipital lobe, even though they
mentioned the difficulty to appreciate them on coronal sections. From our FOD
maps build from high-resolution 3D-PLI vector data, one can separately
identify them up to a coarser resolution of $640\times 640\times 700\text{
}\mu m^{3}$ in (C) where the thinner layer intersect another vertical fiber
bundle to form crossing FODs. We can see loss of WM fibers at this scale,
suggesting, thus, that higher resolutions of $128\times 128\times 140\text{
}\mu m^{3}$ and $320\times 320\times 350\text{ }\mu m^{3}$ could be considered
as suitable spatial scales for diffusion MRI to reach and apply bundle-
specific tractography (46) of the SS. Indeed, using a SH $L_{max}=6$ at these
particular scales, the specific WM tracts of the SS can be tracked, maybe more
preferably at SV size of $2\times 2\times 2$ equivalent to $128\times
128\times 140\text{ }\mu m^{3}$ resolution as observed in our study.

Figure 8: Fiber orientation map of a coronal brain section in the region of
the stratum sagittale (A) observed with 3D-PLI at $1.3\times 1.3\times\mu
m^{2}$ resolution. Details of crossings (radial and horizontal fibers) in a
transition zone between white and gray matter is shown in (B). (C) highlights
the crossing of three pathway crossings at the bottom right (dark spot).

##  5 Discussion

In this paper, we propose an analytical and fast approach to resolve the FOD
reconstruction problem in 3D-Polarized Light Imaging. We quantitatively
compare the performances of our approach and the discrete pliODF method on
both highly rich synthetic and human brain 3D-PLI datasets.

##### On Angular Precision, Angular Resolution

As pointed out by (42), “characterisation of any error in the resolution of
the single fiber orientation will provide an easy way of identifying any
underlying systematic problems inherent in certain (FOD reconstruction)
techniques…”. Thus, both reconstruction techniques are performed on a single
fiber synthetic model. Results show that, unlike our aFOD, pliODF is dependent
on the number bins of its discretized directional histogram. This demonstrates
that pliODF requires high number of bins to reach results similar to the
analytical solution. Indeed, in this experiment (see Sect. 4.1), the angular
precision $\epsilon$ of aFOD is always smaller than $1^{\circ}$ while the
$\epsilon$ of pliODF converges to the aFOD’s by decreasing its bin width.

Our aFOD is analytically defined, therefore, the produced angular errors in
the single and crossing fiber experiments could be explained by the truncation
of the spherical harmonics basis as reported in the previous chapter (and by
potential biases from the FOD peaks detection method). In the multiple fiber
context, except for $X=60^{\circ}$-crossing where $\epsilon$ increases when
$L_{max}$ goes from 6 to 8, the angular precision improves when the SH
truncation order increases, independent of the crossing configuration. In the
particular cases of crossing at $X=0^{\circ}\text{ and }X=90^{\circ}$, the
deviation is always very low, i.e. less than $1^{\circ}$. Moreover, using high
SH orders helps sharpen the aFOD (and pliODF with proper bin width) and gives
more defined peaks very useful in fiber-tracking applications (50).

As for pliODF, the presence of angular error can be attributed to the SH basis
truncation as well, however, it could be mainly due to its discrete nature.
pliODF by definition uses an arbitrary discretized directional histogram and
only provides an approximate of the FOD (10), which is the main limitation of
the approach. Indeed, the directional histogram is defined by its bin centers
and the number of native orientation vectors falling in each bin is dependent
on the bin size or corresponding dihedral angle. This could explain the fact
that, in the single fiber experiment when the histogram bin size is
$<4^{\circ}$, pliODF has lower $\epsilon$ than aFOD does. In fact, the closer
this bin center to the orientation of single fiber ground truth, the smaller
the angular error of pliODF. Conversely, the further the bin center, the
larger the $\epsilon$ as shown in Fig. 2 and Table 2. In the crossing fiber
experiment, results in Fig. 4 show again that the aFOD presents less angular
deviation than the pliODF.

The use of the discretized directional histogram has also a significant
influence on the angular resolution of the pliODF in the crossing fiber
experiment, but the bin sizes we have selected ($3.6^{\circ}\text{ and
}1.8^{\circ}$) are ‘empirically’ chosen so as to get pliODF results close to
the aFOD solution. As said earlier, when the bin width gets smaller, i.e. the
number of histogram bins gets larger, pliODF is expected to behave as well as
our the aFOD (Fig. 3, Fig. 4). Nevertheless, one should be cautious about
possible over-binning of the directional histogram which could eventually
cause some undesirable effects such as the appearance of spurious lobes (Fig.
3, Fig. 6) or eventual extended computation runtimes.

As aforementioned, the value of the spherical harmonics order $L_{max}$
determines the maximum achievable FOD angular resolution. By increasing
$L_{max}$, more acute crossing fiber populations could be resolved as shown in
Sect. 4.2. Here again we set the directional histogram bin size and both
approaches give very similar angular resolutions $\tau$ but aFOD in general
demonstrates better angular precision $\epsilon$ (Fig. 4).

As for imaging noise influence, in a previous study (2), we tested our
analytical method on a synthetic dataset and it was revealed very little
influence of noise on both the angular resolution and angular precision of the
aFOD. We expect pliODF to perform as well, if the good number of radial
histogram bins is used. Indeed, the noise in 3D-PLI data we deal with comes
mainly from the CCD camera and its removal is performed during the estimation
of the high-resolution fiber orientation vector (8, 9, 33, 58, 4, 49).

##### On computational runtimes

The analytical approach is very fast and computationally efficient. This could
be due to the choice of defining each fiber orientation as a $2D$ Dirac
function on the unit sphere, which simplify and accelerate the computation of
the SH coefficients, analytically. Note that in our study, only a subset of
the human brain was analysed. An entire human brain of about $1200\text{
}cm^{3}$ has, however, about $100$ times more voxels when scanned at $64\times
64\times 70\text{ }\mu m^{3}$ resolution and it would theoretically last about
$100$ times longer to calculate the FODs. Moreover, at a (ultra) higher
resolution of $1.3\times 1.3\times 70\text{ }\mu m^{3}$ the dataset gets about
$2\times 10^{5}$ times larger and the computations will require more time. We
recall that to reach millimeter resolutions, the SV size should be very large,
therefore, pliODF runs faster as seen in Sect. 4.3 and might still be
applicable, given a suitable bin size is used. However, since the calculation
time at this scale is already quite fast (seconds for the current dataset),
the use of the aFOD is still advocated, since it gives the exact solution and
has no bias of the histogram binning.

##### Limitation

One limitation of representing the fiber orientation distribution in a
truncated SH basis is the presence of Gibbs ringing effects as discussed in a
recent review by (15). However, this issue can be resolved or drastically
mitigated for the aFOD by introducing an apodized delta function in its
definition as recently proposed by (43) and (21) in diffusion MRI.

##### On human brain data

Finally, our analytical algorithm preserves the global integrity of the WM
fiber architecture in the human brain. Indeed, experiments in deep WM region
(Fig. 6) as well as in the cortex (Fig. 5) where the fiber structures are well
known (47, 56), show the conservation of the general anatomical organization
of the stratum sagittale (SS) and the U-fibers at different spatial
resolutions. It is noticed though that pliODF sometimes presents some spurious
lobes, which could be due to the empirical (over- or under-) binning of the
directional histogram. Conversely, high SH orders improve the sharpness of the
FOD as well as its angular resolution and precision as demonstrated in the
simulations (Fig. 3, Fig. 4), and commonly known in diffusion MRI (51, 16).

Further analysis of the SS illustrates the loss of its thinner layer or simply
its intersection with a differently oriented fiber population to form more
complex fiber configurations. This is visible at a coarse resolution of about
$640\times 640\times 700\text{ }\mu m^{3}$ in Fig. 7C. Whereas at higher
resolutions of around $128\times 128\times 140\text{ }\mu m^{3}$, both layers
can be clearly appreciated. This observation would suggest a potential spatial
resolution limit, for diffusion MRI to reach, at which the WM fibers of the SS
could be reconstructed and possibly improve their tracking.

Note, that the state-of-the-art 3D-PLI signal analysis only provides a single
orientation estimate, i.e. one orientation vector per native voxel.
Consequently, an initial image resolution of $64\times 64\text{ }\mu m^{2}$
leads to fiber orientation descriptions (including their fiber tract
crossings) at larger spatial scale. This does not mean, that there is no
crossing present in human brain data shown in Figs. (5,6,7). As demonstrated
in Fig. 8, the studied region of interest in the sagittal stratum, for
example, is full of fiber crossings which are revealed by 3D-PLI at axonal
resolution, i.e. at $1.3\times 1.3\text{ }\mu m^{2}$ pixel size. Future
research will certainly address aFOD determination in serial 3D-PLI
measurements at the one micrometer scale.

##  6 Conclusion

We have presented an analytical and fast method to define and reconstruct the
fiber orientation distribution in 3D-PLI. To achieve it, we model each high-
resolution fiber orientation vector contained in the super-voxel as 2D Diracs
delta function on the unit sphere and the aFOD as their sum. The aFOD is
analytically described on a spherical harmonics basis and efficiently computed
via the spherical Fourier transform and by means of the Diracs. This makes the
aFOD not only computationally very fast but also with high angular resolution
and precision and the ability to close the resolution gap with diffusion MRI.
This work serves as a important step towards a major goal which is the
validation of diffusion MRI fiber orientation estimates and histology guided
tractography via 3D-Polarized Light Imaging technique, more specifically, it
leads to multimodal brain imaging analysis across the scales, through the
fiber orientation distribution.

## Acknowledgements

This work was partly supported by ANR "MOSIFAH" under ANR-13-MONU-0009-01, the
ERC under the European Union’s Horizon 2020 Research and Innovation Program
(ERC Advanced Grant agreement No. 694665:CoBCoM). M.S, F.M and S.M were partly
funded by the European Union’s Horizon 2020 Research and Innovation Program
under Grant Agreement No. 785907 (HBP SGA2). The authors gratefully
acknowledge the computing time granted through JARA-HPC on the supercomputer
JURECA (30) at Forschungszentrum Jülich.

##  Appendix A Synthetic Data Generation

The nerve fiber model imitates two fiber bundles in the XY-plan crossing at an
angle $X$ which varies from $0^{\circ}$ to $90^{\circ}$, with $10^{\circ}$
increments. To build this configuration, the following steps were performed:

  1. 1.

Two inital fibers bundles with a length of $3840\text{ }\mu m$ and a radius of
$640\text{ }\mu m$ were placed at an crossing angle $X$.

  2. 2.

These initial fiber bundles were then filled with smaller fiber bundles with a
radius of $64\text{ }\mu m$, on which the collision algorithm (39) was
performed for $25$ steps, to allow a ‘macroscopic’ weave pattern.

  3. 3.

The resulting fiber bundles were then filled one last time with fibers with a
radius of $10\text{ }\mu m$. In this state, each fiber is spatially divided
into $20\text{ }\mu m$ segments. To limit the curvature of each fiber, the
minimum radius of a bending fiber is limited to $40\text{ }\mu m$.

This configuration runs with the collision solving algorithm until no
collision are detected.

The resulting fiber configuration, for instance as pictured in Fig. 1A where
$X=90^{\circ}$, is then used for the 3D-PLI simulation to generate the
resulting FOM. For the simulation, the model is discretized on a volume of
$3840\times 3840\times 600\text{ }\mu m^{3}$ to allow the calculations
performed by simPLI (22, 40). Each fiber is modelled as a filled cylindrical
object with an absorption coefficient $\mu=1/{mm}^{3}$ and a birefringence
$\Delta n=-0.001$. The deriving microscopic images have a in-plane resolution
of $60\times 60\text{ }\mu m^{2}$. The signal noise was calculated with a gain
factor of $g=3$ (see (49)) and an initial light intensity of $26000\text{
}a.u$. The signal from the 3D-PLI pipeline is then analyzed to calculate the
FOM.

## References

  * (1) Daniel C Alexander.  Multiple-fiber reconstruction algorithms for diffusion mri.  White Matter in Cognitive Neuroscience: Advances in Diffusion Tensor Imaging and Its Applications, 1064:113–133, 2005.
  * (2) Abib Alimi, S Deslauriers-Gauthier, F Matuschke, D Schmitz, M Axer, and R Deriche.  Analytical fiber odf reconstruction in 3d polarized light imaging: Performance assessment.  In ISBI 2019-IEEE International Symposium on Biomedical Imaging, 2019.
  * (3) Abib Alimi, Samuel Deslauriers-Gauthier, and Rachid Deriche.  Towards validation of diffusion mri tractography: bridging the resolution gap with 3d polarized light imaging.  In ISMRM 2019 - International Society for Magnetic Resonance in Medicine, 2019.
  * (4) Abib Alimi, Marco Pizzolato, Rutger HJ Fick, and Rachid Deriche.  Solving the inclination sign ambiguity in three dimensional polarized light imaging with a pde-based method.  In 2017 IEEE 14th International Symposium on Biomedical Imaging (ISBI 2017), pages 737–740. IEEE, 2017.
  * (5) Abib Alimi, Yves Ussou, Pierre-Simon Jouk, Gabrielle Michalowicz, and Rachid Deriche.  An analytical fiber odf reconstruction in 3d polarized light imaging.  In 2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018), pages 1276–1279. IEEE, 2018.
  * (6) Abib OY Alimi, Samuel Deslauriers-Gauthier, and Rachid Deriche.  Regularizing the ODF estimate with the Laplace-Beltrami operator in 3D Polarized Light Imaging.  CoBCoM 2017 - Computational Brain Connectivity Mapping Winter School Workshop, November 2017.
  * (7) Brian B Avants, Nicholas J Tustison, Gang Song, Philip A Cook, Arno Klein, and James C Gee.  A reproducible evaluation of ants similarity metric performance in brain image registration.  Neuroimage, 54(3):2033–2044, 2011.
  * (8) Markus Axer, Katrin Amunts, David Grässel, Christoph Palm, Jürgen Dammers, Hubertus Axer, Uwe Pietrzyk, and Karl Zilles.  A novel approach to the human connectome: ultra-high resolution mapping of fiber tracts in the brain.  Neuroimage, 54(2):1091–1101, 2011.
  * (9) Markus Axer, David Grässel, Melanie Kleiner, Jürgen Dammers, Timo Dickscheid, Julia Reckfort, Tim Hütz, Björn Eiben, Uwe Pietrzyk, Karl Zilles, et al.  High-resolution fiber tract reconstruction in the human brain by means of three-dimensional polarized light imaging.  Frontiers in neuroinformatics, 5:34, 2011.
  * (10) Markus Axer, Sven Strohmer, David Gräßel, Oliver Bücker, Melanie Dohmen, Julia Reckfort, Karl Zilles, and Katrin Amunts.  Estimating fiber orientation distribution functions in 3d-polarized light imaging.  Frontiers in neuroanatomy, 10:40, 2016.
  * (11) Peter J Basser, James Mattiello, and Denis LeBihan.  Mr diffusion tensor spectroscopy and imaging.  Biophysical journal, 66(1):259–267, 1994.
  * (12) Matthew D Budde and Joseph A Frank.  Examining brain microstructure using structure tensor analysis of histological sections.  Neuroimage, 63(1):1–10, 2012.
  * (13) Svenja Caspers, Markus Axer, Julian Caspers, Christiane Jockwitz, Kerstin Juetten, Julia Reckfort, David Graessel, Katrin Amunts, and Karl Zilles.  Target sites for transcallosal fibers in human visual cortex–a combined diffusion and polarized light imaging study.  Cortex, 72:40–53, 2015.
  * (14) Alessandro Daducci, Erick Jorge Canales-Rodrı, Maxime Descoteaux, Eleftherios Garyfallidis, Yaniv Gur, Ying-Chia Lin, Merry Mani, Sylvain Merlet, Michael Paquette, Alonso Ramirez-Manzanares, et al.  Quantitative comparison of reconstruction methods for intra-voxel fiber recovery from diffusion mri.  IEEE transactions on medical imaging, 33(2):384–399, 2014.
  * (15) Flavio Dell’Acqua and J-Donald Tournier.  Modelling white matter with spherical deconvolution: How and why?  NMR in Biomedicine, 32(4):e3945, 2019.
  * (16) Maxime Descoteaux.  High angular resolution diffusion MRI: from local estimation to segmentation and tractography.  PhD thesis, 2008.
  * (17) Maxime Descoteaux, Elaine Angelino, Shaun Fitzgibbons, and Rachid Deriche.  Regularized, fast, and robust analytical q-ball imaging.  Magnetic resonance in medicine, 58(3):497–510, 2007.
  * (18) Maxime Descoteaux, Rachid Deriche, Thomas R Knosche, and Alfred Anwander.  Deterministic and probabilistic tractography based on complex fibre orientation distributions.  IEEE transactions on medical imaging, 28(2):269–286, 2008.
  * (19) Samuel Deslauriers-Gauthier and Pina Marziliano.  Sampling signals with a finite rate of innovation on the sphere.  IEEE Transactions on Signal Processing, 61(18):4552–4561, 2013\.
  * (20) Samuel Deslauriers-Gauthier, Pina Marziliano, Michael Paquette, and Maxime Descoteaux.  The application of a new sampling theorem for non-bandlimited signals on the sphere: Improving the recovery of crossing fibers for low b-value acquisitions.  Medical image analysis, 30:46–59, 2016.
  * (21) Thijs Dhollander, Louise Emsell, Wim Van Hecke, Frederik Maes, Stefan Sunaert, and Paul Suetens.  Track orientation density imaging (todi) and track orientation distribution (tod) based tractography.  NeuroImage, 94:312–336, 2014.
  * (22) Melanie Dohmen, Miriam Menzel, Hendrik Wiese, Julia Reckfort, Frederike Hanke, Uwe Pietrzyk, Karl Zilles, Katrin Amunts, and Markus Axer.  Understanding fiber mixture by simulation in 3d polarized light imaging.  NeuroImage, 111:464–475, 2015.
  * (23) Patric Hagmann, Lisa Jonasson, Philippe Maeder, Jean-Philippe Thiran, Van J Wedeen, and Reto Meuli.  Understanding diffusion mr imaging techniques: from scalar diffusion-weighted imaging to diffusion tensor imaging and beyond.  Radiographics, 26(suppl_1):S205–S223, 2006.
  * (24) Claudia Hänel, Ali C Demiralp, Markus Axer, David Gräßel, Bernd Hentschel, and Torsten W Kuhlen.  Interactive level-of-detail visualization of 3d-polarized light imaging data using spherical harmonics.  2017\.
  * (25) Dennis M Healy Jr, Harrie Hendriks, and Peter T Kim.  Spherical deconvolution.  Journal of Multivariate Analysis, 67(1):1–22, 1998.
  * (26) Dylan JHA Henssen, Jeroen Mollink, Erkan Kurt, Robert van Dongen, Ronald HMA Bartels, David Grä$\beta$el, Tamas Kozicz, Markus Axer, and Anne-Marie Van Cappellen van Walsum.  Ex vivo visualization of the trigeminal pathways in the human brainstem using 11.7 t diffusion mri combined with microscopy polarized light imaging.  Brain Structure and Function, 224(1):159–170, 2019.
  * (27) Ben Jeurissen, Maxime Descoteaux, Susumu Mori, and Alexander Leemans.  Diffusion mri fiber tractography of the brain.  NMR in Biomedicine, 32(4):e3785, 2019.
  * (28) Derek K Jones.  Diffusion mri.  Oxford University Press, 2010.
  * (29) P-S Jouk, Yves Usson, Gabrielle Michalowicz, and Laurence Grossi.  Three-dimensional cartography of the pattern of the myofibres in the second trimester fetal human heart.  Anatomy and embryology, 202(2):103–118, 2000.
  * (30) Jülich Supercomputing Centre.  JURECA: General-purpose supercomputer at Jülich Supercomputing Centre.  Journal of large-scale research facilities, 2, A62, 2016.
  * (31) Ahmad Raza Khan, Anda Cornea, Lindsey A Leigland, Steven G Kohama, Sune Nørhøj Jespersen, and Christopher D Kroenke.  3d structure tensor analysis of light microscopy data for validating diffusion mri.  Neuroimage, 111:192–203, 2015.
  * (32) Stefan Klein, Marius Staring, Keelin Murphy, Max A Viergever, and Josien PW Pluim.  Elastix: a toolbox for intensity-based medical image registration.  IEEE transactions on medical imaging, 29(1):196–205, 2009.
  * (33) Melanie Kleiner, Markus Axer, David Gräßel, Julia Reckfort, Uwe Pietrzyk, Katrin Amunts, and Timo Dickscheid.  Classification of ambiguous nerve fiber orientations in 3d polarized light imaging.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 206–213. Springer, 2012.
  * (34) Luiza Larsen, Lewis D Griffin, David GRäßel, Otto W Witte, and Hubertus Axer.  Polarized light imaging of white matter architecture.  Microscopy research and technique, 70(10):851–863, 2007.
  * (35) Denis Le Bihan, Eric Breton, Denis Lallemand, Philippe Grenier, Emmanuel Cabanis, and Maurice Laval-Jeantet.  Mr imaging of intravoxel incoherent motions: application to diffusion and perfusion in neurologic disorders.  Radiology, 161(2):401–407, 1986.
  * (36) Alexander Leemans.  Diffusion mri of the brain: The naked truth.  NMR in Biomedicine, 32(4):e4084, 2019.
  * (37) Trygve B Leergaard, Nathan S White, Alex De Crespigny, Ingeborg Bolstad, Helen D’Arceuil, Jan G Bjaalie, and Anders M Dale.  Quantitative histological validation of diffusion mri fiber orientation distributions in the rat brain.  PloS one, 5(1):e8595, 2010.
  * (38) Klaus H Maier-Hein, Peter F Neher, Jean-Christophe Houde, Marc-Alexandre Côté, Eleftherios Garyfallidis, Jidan Zhong, Maxime Chamberland, Fang-Cheng Yeh, Ying-Chia Lin, Qing Ji, et al.  The challenge of mapping the human connectome based on diffusion tractography.  Nature communications, 8(1):1349, 2017.
  * (39) Felix Matuschke, Kévin Ginsburger, Cyril Poupon, Katrin Amunts, and Markus Axer.  Dense fiber modeling for 3d-polarized light imaging simulations, 2019\.
  * (40) Miriam Menzel, Kristel Michielsen, Hans De Raedt, Julia Reckfort, Katrin Amunts, and Markus Axer.  A jones matrix formalism for simulating three-dimensional polarized light imaging of brain tissue.  Journal of the Royal Society Interface, 12(111):20150734, 2015.
  * (41) Susumu Mori and Peter CM Van Zijl.  Fiber tracking: principles and strategies–a technical review.  NMR in Biomedicine: An International Journal Devoted to the Development and Application of Magnetic Resonance In Vivo, 15(7-8):468–480, 2002\.
  * (42) Greg D Parker, David Marshall, Paul L Rosin, Nicholas Drage, Stephen Richmond, and Derek Kenton Jones.  A pitfall in the reconstruction of fibre odfs using spherical deconvolution of diffusion mri data.  Neuroimage, 65:433–448, 2013.
  * (43) David Raffelt, J-Donald Tournier, Stuart Crozier, Alan Connelly, and Olivier Salvado.  Reorientation of fiber orientation distributions using apodized point spread functions.  Magnetic Resonance in Medicine, 67(3):844–855, 2012.
  * (44) Julia Reckfort, Hendrik Wiese, Uwe Pietrzyk, Karl Zilles, Katrin Amunts, and Markus Axer.  A multiscale approach for the reconstruction of the fiber architecture of the human brain based on 3d-PLI.  Frontiers in Neuroanatomy, 9, September 2015.
  * (45) Jan André Reuter, Felix Matuschke, Miriam Menzel, Nicole Schubert, Kévin Ginsburger, Cyril Poupon, Katrin Amunts, and Markus Axer.  Faconstructor: an interactive tool for geometric modeling of nerve fiber architectures in the brain.  International journal of computer assisted radiology and surgery, pages 1–9, 2019.
  * (46) Francois Rheault, Etienne St-Onge, Jasmeen Sidhu, Klaus Maier-Hein, Nathalie Tzourio-Mazoyer, Laurent Petit, and Maxime Descoteaux.  Bundle-specific tractography with incorporated anatomical and orientational priors.  NeuroImage, 186:382–398, 2019.
  * (47) Heinrich Sachs.  Das Hemisphärenmark des menschlichen Grosshirns: Der Hinterhauptlappen. Mit einem Vorwort von C. Wernicke.  G. Thieme, 1892.
  * (48) Kurt Schilling, Vaibhav Janve, Yurui Gao, Iwona Stepniewska, Bennett A Landman, and Adam W Anderson.  Comparison of 3d orientation distribution functions measured with confocal microscopy and diffusion mri.  Neuroimage, 129:185–197, 2016.
  * (49) Daniel Schmitz, Sascha EA Muenzing, Martin Schober, Nicole Schubert, Thomas Lippert, Katrin Amunts, and Markus Axer.  Derivation of fiber orientations from oblique views through human brain sections in 3d-polarized light imaging.  Frontiers in neuroanatomy, 12:75, 2018.
  * (50) J-Donald Tournier.  Diffusion mri in the brain–theory and concepts.  Progress in Nuclear Magnetic Resonance Spectroscopy, 2019.
  * (51) 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.  Neuroimage, 35(4):1459–1472, 2007.
  * (52) J-Donald Tournier, Fernando Calamante, David G Gadian, and Alan Connelly.  Direct estimation of the fiber orientation density function from diffusion-weighted mri data using spherical deconvolution.  NeuroImage, 23(3):1176–1185, 2004.
  * (53) 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.  NeuroImage, page 116137, 2019.
  * (54) David S Tuch.  Q-ball imaging.  Magnetic resonance in medicine, 52(6):1358–1372, 2004.
  * (55) 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.  Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine, 48(4):577–582, 2002\.
  * (56) Francesco Vergani, Sajedha Mahmood, Cristopher M Morris, Patrick Mitchell, and Stephanie J Forkel.  Intralobar fibres of the occipital lobe: a post mortem dissection study.  Cortex, 56:145–156, 2014.
  * (57) Hui Wang, Christophe Lenglet, and Taner Akkin.  Structure tensor analysis of serial optical coherence scanner images for mapping fiber orientations and tractography in the brain.  Journal of biomedical optics, 20(3):036003, 2015.
  * (58) Hendrik Wiese, David Gräßel, Uwe Pietrzyk, Katrin Amunts, and Markus Axer.  Polarized light imaging of the human brain: a new approach to the data analysis of tilted sections.  In Polarization: Measurement, Analysis, and Remote Sensing XI, volume 9099, page 90990U. International Society for Optics and Photonics, 2014\.
  * (59) Terry S Yoo, Michael J Ackerman, William E Lorensen, Will Schroeder, Vikram Chalana, Stephen Aylward, Dimitris Metaxas, and Ross Whitaker.  Engineering and algorithm design for an image processing api: a technical report on itk-the insight toolkit.  Studies in health technology and informatics, pages 586–592, 2002\.
  * (60) Michael M Zeineh, Nicola Palomero-Gallagher, Markus Axer, David Grä$\beta$el, Maged Goubran, Andreas Wree, Roger Woods, Katrin Amunts, and Karl Zilles.  Direct visualization and mapping of the spatial course of fiber tracts at microscopic resolution in the human hippocampus.  Cerebral cortex, 27(3):1779–1794, 2016.