# A Riemannian framework for incorporating white matter bundle priors in ODF-based tractography algorithms.

Thomas Durantel, Gabriel Girard, Emmanuel Caruyer, Olivier Commowick,

Julie Coloigner

► To cite this version:

Thomas Durantel, Gabriel Girard, Emmanuel Caruyer, Olivier Commowick, Julie Coloigner. A Riemannian framework for incorporating white matter bundle priors in ODF-based tractography algorithms. PLoS ONE, In press, pp.1-10. hal-04246380

HAL Id: hal-04246380

https://hal.science/hal-04246380v1

Submitted on 17 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A Riemannian framework for incorporating white matter bundle priors in ODF-based tractography algorithms.

Thomas Durantel, Gabriel Girard, Emmanuel Caruyer, Olivier Commowick and Julie Coloigner

**Abstract**— Diffusion magnetic resonance imaging (dMRI) tractography is a powerful approach to study brain structural connectivity. However, its reliability in a clinical context is still highly debated. Recent studies have shown that most classical algorithms achieve to recover the majority of the existing true bundles. However, the generated tractograms contain many invalid bundles. This is due to the crossing fibers and bottleneck problems which increase the number of false positives fibers. In this work, we proposed to overpass this limitation with a novel method to guide the algorithms in those challenging regions with prior knowledge of the anatomy. We developed a method of creating and combination of anatomical prior applicable to any orientation distribution function (ODF)-based tractography algorithms. The proposed method captures the track orientation distribution (TOD) from an atlas of segmented fiber bundles and incorporates it during the tracking process, using a Riemannian framework. We tested the prior incorporation method on two ODF-based state-of-the-art algorithms, iFOD2 and Trekker PTT, on the diffusion-simulated connectivity (DiSCo) dataset and on the Human Connectome project (HCP) data. We showed that our method improves the overall spatial coverage and connectivity of a tractogram on the two datasets, especially in crossing fiber regions. Moreover, the fiber reconstruction may be improved on clinical data, informed by prior extracted on high quality data, and therefore could help in the study of brain anatomy and function.

**Index Terms**— Anatomical prior, Diffusion MRI, Track Orientation Distribution, Tractography, Tractometry

## I. INTRODUCTION

DIFFUSION magnetic resonance imaging (dMRI) is an MRI modality that allows to measure the thermal agitation of the water molecule in the brain [1]. This agitation being constrained by the tissues micro-structures, typically the nervous system axons, dMRI enables the voxelwise estimation

T. D. is with the University of Rennes, INRIA, CNRS, INSERM, IRISA UMR 6074, Empenn ERL U-1228, F-35000, Rennes, France.

G. G. was with the Signal Processing Laboratory (LTS5), Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland and is with the Department of Computer Science, Université de Sherbrooke, Québec, Canada.

E. C. is with the University of Rennes, INRIA, CNRS, INSERM, IRISA UMR 6074, Empenn ERL U-1228, F-35000, Rennes, France.

O. C., deceased, was with the University of Rennes, INRIA, CNRS, INSERM, IRISA UMR 6074, Empenn ERL U-1228, F-35000, Rennes, France.

J. C. is with the University of Rennes, INRIA, CNRS, INSERM, IRISA UMR 6074, Empenn ERL U-1228, F-35000, Rennes, France (e-mail: julie.coloigner@irisa.fr).

of the orientations of the white matter fibers [2], [3]. By randomly choosing seeds from within the brain white matter, then following, from one voxel to the next, the local fiber orientations, one can achieve to estimate the brain connectivity and characterize the physical connections that mediate information transfer between cortical regions [3]. This process is called fibers tracking, or fiber tractography, and the resulting set of white matter trajectories is called tractograms [4]. The ability of those approaches to delineate the white matter fiber pathways offers unprecedented insight into the structural connections of the human brain and enormous potential for the study of human brain anatomy, development and function [5]. Moreover, tractography has been proven particularly useful to neurosurgeons for the planning of surgery, especially to preserve important white matter pathways during resections [6].

A multitude of new tracking algorithms has been proposed to improve the quality of the tractograms [7], [8]. The first methods were deterministic methods, where only the principal diffusion orientation of the diffusion tensor was used in each voxel [9]. However, studies have shown the difficulties that these methods have to represent complex brain regions, like crossing or fanning fibers [10]. In response, new diffusion model, like the orientation distribution function (ODF) which encodes continuously all the diffusion orientations within a voxel and thus can characterize the complex fiber structure, have been developed [11], [12]. Modern methods have then sought to estimate the fiber dispersion estimation using probabilistic methods [5], [13], [14]. Among them, existing local tractography approaches estimates the local fiber orientation at each voxel independently by simultaneously fitting the local model and propagating in the most consistent direction, such as Kalman filtering method [15], particle filtering method [16], [17], random walks methods [11] and graph theory method [18]. Other approaches have proposed more global approach which compute multiple fibers pathways and select the best ones based on the diffusion characteristics, like using Bayesian model [19] or the Hough transform [20] and machine learning more recently [21]. However, despite the improvement of the dMRI acquisition, the modeling and the tracking, white matter fibers tractography still can't overcome some limitations [22], [23] which prevents those methods to be usable in the clinical context. Maier-Hein et al (2017). have shown that most modern algorithms can recover up to 90% of the ground truth white matter bundles, but these are poorly recovered in terms of
spatial extent [23]. Moreover, those tractograms also contained large amounts of invalid bundles [23]. Overall, the iFOD2 [24] and PTT [25] algorithms had good performances in recent international challenges [23], [26].

To overpass the poor spatial extent of white matter bundles, one proposed solution is to increase the numbers of the total generated streamlines leading to an over representation of easy tracks and an under-representation of difficult tracks, inducing a density bias [27]. For example, the local tractography methods aim to recover a global connectivity by inferring from local directional information that causes the algorithm to follow, in complex regions, the easiest path available and sometimes the non-realistic decisions [23]. Various studies have suggested that advanced diffusion microstructure modeling [28], [29], streamline filtering techniques [30] or advances in machine-learning-driven tractography [31] could make improvement on the false positive rate.

Other approaches have proposed to compensate the lack of global information by adding them in the form of anatomical priors whose purpose would be to guide the algorithm in complex regions. As for segmentation and label fusion approach, recent methods are based on the use of an anatomical atlas such as the TRACULA method [32], with promising results combining deterministic tractography and anatomical prior. Then, diffusion priors during the tracking process were proposed in [33]–[35], improving the delineation of white matter bundles. In [35], a bundle-specific method incorporates anatomical and orientational priors based on a template, to improve the reconstruction of long fibers and increase of reproducibility, sensitivity and specificity. On the other hand, in [36], another well-known approach based on a machine learning method was proposed to automatically segment, with high precision, the overall shape of a bundle, being a potentially great source of prior.

Based on these promising results, we developed a method of anatomical prior creation and combination, which is usable with any tractography algorithm based on orientation distribution function (ODF) [37]. Our method use pre-segmented fiber bundles by [36], to agglomerate global information from several different brains and capture the orientational variability in complex brain region. In this work, prior are computed on this anatomical atlas and expressed in the form of voxel-wise TOD [38] and then combined with the orientation distribution function (ODF) [37] data using a Riemannian framework [39]. We decided to incorporate those priors in two state-of-the-art algorithms using traditional ODF data, MRtrix iFOD2 [24] and Trekker PTT [25]. We then evaluated those methods on the diffusion-simulated connectivity (DiSCo) dataset [40] and on the Human Connectome project (HCP)<fnref n="1" /> data to show the increase quality of the tractogram and more specifically on the spatial extent of the reconstructed bundles.

## II. METHODS

Our method can be separated into 3 distinct parts: i) the construction of an anatomical atlas from segmented fibers, considered as gold standard, ii) the extraction and estimation

of the TOD anatomical priors, from the atlas and iii) the combination between the priors and the subject data. The entire framework is illustrated in Figure 1.

### A. Atlas construction

The first step of our framework is the creation of a reference fiber atlas from a set of high-resolution diffusion images (see Figure 1.a). In this paper, we chose high-resolution diffusion images from the Human Connectome Project (HCP) young adult study [41], which were acquired at high spatial resolution and b-values with 90 directions and thus contains high quality information on the fibers orientations.

The proposed atlasing method follows a modified version of Guimond et al (2000). [42] which was adapted for diffusion data, to compute an atlas of Tensor images from a set of control subjects. This approach is based on a procedure which compute iteratively the atlas by registering the tensor images of HCP images onto a current reference. After each iteration, it performs an average of all the previous unbiased atlases to compute the next one at the following iteration. The main modifications from the original Guimond et al. (2000) method is to use diffeomorphisms encoded as Stable Vector Fields (SVF) and the log-Euclidean framework [43] to compute the average transformations and the approach was adapted to diffusion data. The entire procedure is detailed in [44].

For each HCP data, 72 bundles were obtained by Wasserthal et al. (2018), using the TractSeg algorithm [36] which perform fiber tracking and fiber filtering, first by regions of interest (ROIs), then by an expert. The overall process used to segment these bundles is detailed in [36] and the data is available online [45].

The transformation field obtained by the atlasing method from each individual space to the atlas space are then applied to 72 segmented fiber bundles of interest of each subject in order to all align them to the fiber atlas.

### B. Priors Extraction

Once aligned to a common space, the bundles of interest are combined in order to represent the general shape of bundles. In order to capture complex orientation, a local fibre orientation prior was estimated using the track orientation distribution (TOD) [38], in each voxel. The tractogram is represented as a voxel-wise probability distribution function (PDF) in the image domain, rather than a set of individual tracks (samples from this distribution). The TOD, thus, captures the expected fibre directions [38].

Our TOD imaging method operates as follows : first, all the fibers directions of the voxel are extracted, then clustered using a k-means algorithm to define the main directions within the voxel and to correct for the density bias. This way, up to 4 principal directions can be extracted, allowing us to correct for the density bias as well as to characterize complex crossing regions. Then, the TOD in each voxel is represented using a set of modified spherical harmonics (SH) basis functions (see [37]) and constructed by projecting one spherical point spread function (PSF) per extracted main direction along the z direction.

<footnote n="1">¹https://www.humanconnectome.org/</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Overall method pipeline. (a) Atlas creation. (b) TOD estimation. (c) Data and prior combination.</figcaption></figure>

The PSF, along the z direction at position $\epsilon$, can be easily obtained in SH basis:

$$ \delta_{z(\epsilon)}(\theta, \phi) = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} c_l^m Y_l^m(\theta, \phi), \quad (1) $$

where $0 \le \theta \le \pi$ and $-\pi \le \phi \le \pi$ are the spherical coordinates, $Y_l^m$ the spherical harmonic of degree $l$ and order $m$ and $c_l^m$, the coefficients, given as:

$$ c_l^m = \int_{-\pi}^{\pi} \int_{0}^{\pi} \delta_{z(\epsilon)}(\theta, \phi) Y_l^m(\theta, \phi) d\theta d\phi \quad (2) $$

We found that resolving the equation 2 with the PSF pointing towards the direction $\theta = \frac{\pi}{2}$ and then rotating them to match the direction extracted from the atlas greatly simplifies the calculations. In [46], authors showed a simple method to rotate functions expressed in SH basis. Since 4 different directions can be represented in a voxel, the final step of the prior extraction is to average the projected PSF. The averaging of distributions defined on the sphere is performed through a Riemannian framework [39]. To do that, we define $\psi(\theta, \phi) = \sqrt{\delta_{z(\epsilon)}(\theta, \phi)}$, the square-root density function of the PDF $\delta(\theta, \phi)$. The square-root is used to ensure that the logarithm maps are available in closed form. Thus, $\delta(\theta, \phi)$ has to be strictly positive. For this reason, in order to find the right PSF, several different distributions were tested (i.e. a Dirac distribution, a Watson distribution and a Gaussian distribution) and then were expressed on the SH basis. We found out that a 2D Gaussian distribution, defined on the unit sphere, is the only one whose number of negative values do not diverge with the degree $l$.

Then, in order to calculate the average distribution, we used the weighted Karcher mean, $\bar{\psi}$, of a set of $n$ points $\{\psi_i\}_{i=0}^n$ in a Riemannian manifold defined by:

$$ \bar{\psi} = \operatorname{argmin} \frac{1}{2} \sum_{i=0}^{n} \omega_i \operatorname{dist}(\psi, \psi_i)^2 \quad (3) $$

with $\omega_i \le 0$ and $\sum_{i=0}^n \omega_i = 1$

As described in [39], the unique solution is $\bar{\psi}$, such as:

$$ \sum_{i=0}^{n} \omega_i \log \bar{\psi}(\psi_i) = 0 \quad (4) $$

where $\log \bar{\psi}$ is the logarithm map from $\bar{\psi}$ to $\psi_i$, given by:

$$ \log_{\bar{\psi}}(\psi_i) = \frac{\psi_i - \langle \psi_i, \bar{\psi} \rangle \bar{\psi}}{\sqrt{1 - \langle \psi_i, \bar{\psi} \rangle}} \cos^{-1} \langle \psi_i, \bar{\psi} \rangle \quad (5) $$

with $\langle \cdot, \cdot \rangle$, the normal dot product.

At the end of this step, using the equation 5 and $\omega_i = 1/n, \forall i \in \{1, ..., n\}$, we obtained $\bar{\psi}$, which is the square-root of the prior for each voxel. An example of the prior with the SH coefficients, truncated at a degree of $l = 8$, of the PSF described previously are shown in Figure 4. The prior creation process in order to average all extracted directions has to be calculated only once.

### C. Data and prior combination

As displayed in Figure 1, the ODF image of the individual diffusion dataset is calculated to estimate the next direction. During tractography, in order to inform the ODF with the anatomical prior, as previously explained in the previous section, a weighted Karcher mean is computed between the ODF and the prior using equation 4, to obtain the enhanced-odf (EODF). However, this method involves calculating, at each step of the tracking, the dot product between distribution on the sphere, that implies an integral over the sphere which would complicate the process. To simplify the computations, Goh et al. (2009) proposed to discretize the PDF and to work with the square root of histograms [39]. In this case, the dot product become summations.

Since the overall purpose of our method is only to guide the tractography algorithms and not to reflect the structure of the priors, the weighting factor $\omega_i$ must be well studied. We want to be able to distinguish between the simple linear region, where not much guidance is needed, especially in easier-to-track bundles with only one fiber direction, and the complex regions with crossing fibers, where the use of the prior is more needed. In order to address this specification, we choose to use 2 measures to weight the prior: the generalised fractional anisotropy (GFA), for the prior and the Akaike information criterion (AIC) [47] for the data.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Left: DiSCo dataset pathways. Right: Prior weighting map example for the DiSCo dataset. Higher value denotes higher usage of the prior and thus more complex regions.</figcaption></figure>

The GFA, is given by:

$$GFA = \sqrt{\frac{1 - c_0^2}{\sum_{i=0}^{N} c_i^2}} \quad (6)$$

where the $c_i$ are the SH coefficients. The GFA is a measure of the anisotropy within the considered voxel, that is, in the case of the prior, this value can reflect the complexity of a region.

The AIC, is described by:

$$AIC = 2k - 2 \ln(L) \quad (7)$$

where k is the numbers of parameters in the model and L the maximized likelihood. The AIC is a measurement of the estimation quality of a model.

Empirically, to well differentiate linear and complex regions, we defined the prior weight given by:

$$\omega_{prior} = \alpha(1 - GFA) + \beta\exp\left(\frac{AIC_{min} - AIC}{2}\right) \quad (8)$$

where $\alpha$ represents the amount of information extracted from the priors, via the GFA, and $\beta$ the amount of information extracted from the data, via the AIC. Those parameters need to be adapted to the data. For the scenario presented in this work, we choose $\alpha = 0.35$ and $\beta = 0.65$ (see Figure 3).

<figure><img src="image_4.png" /><figcaption>Fig. 3. Evolution of the Dice score versus $\alpha$ (from the Equation 8) for the tracking of the CST and for the iFOD and PTT algorithms. The method performs most effectively with an $\alpha$ value of 0.35.</figcaption></figure>

An example of a weighting map for the DiSCo dataset can be seen on Figure 2 with higher $\omega_{prior}$ in crossing fiber areas.

The result of this averaging process is then an EODF, expressed in SH basis, usable in any ODF-based tractography algorithm. Figure 4 shows examples of data and priors combination, as well as the effect of weighting on the results. The first row represents a voxel where the ODF and the prior

<figure><img src="image_5.png" /><figcaption>Fig. 4. Examples of data and prior combination. On the first row: prior and data are aligned. The EODF main orientation are unchanged. On the second and third rows: Prior and data are crossing. The corresponding EODF peak is more or less attenuated using the weight of the prior.</figcaption></figure>

have the same direction, corresponding to an easy-to-track regions with one-way crossing. On the other hand, in the second and third rows, the orientation of the ODF and the prior are different. In this case, according to the value of $\omega_{prior}$, the EODF is almost equal to the ODF or a mixture of the prior and ODF.

### III. EXPERIMENTS

In order to validate our method using priors, we conducted our first experiments on the DiSCo challenge dataset [26]. Then, we tested our algorithms on HCP data. For each experiment, we implemented our prior approach on two state-of-the-art algorithms, MRtrix iFOD2 [24] and Trekker PTT [25]. Thus, on the two dataset, 4 tractography reconstruction were performed, both with and without using priors. For iFOD2, we used a step size of 0.2 voxel and a maximum angle between successive steps of 20 degrees. For PTT, we used the default parameters.

#### A. DiSCo challenge dataset

The aim of this experiment is to test the ability of our novel method to improve the connectivity estimation and the effect of noise. In this context, we do not have access to several datasets to build an anatomical atlas, instead, we computed the TOD: (i) using the ground-truth fiber pathways as anatomical priors and (ii) using fiber tracked on high quality data as anatomical priors.

1) Data: We performed two different variants of this experiment: in the first one, priors are estimated on ground-truth fibers, not available on in vivo data; and the second experiment is used to demonstrate the benefit of our method, without having a ground-truth but priors build on high resolution data.

(i) Fiber pathways ground-truth: The DiSCo dMRI images have a grid of 40 × 40 × 40 with voxel size of 1 mm² and are acquired with 4 different b-shell of 1000, 2000, 3000 and 1300 s/mm², 90 directions per b-shell and 4 b0. In order to obtain synthetic data with the same quality as those
<figure><img src="image_6.png" /><figcaption>Fig. 5. Visual comparison of the resulting tractograms, with and without prior. From left to right, in column: iFOD2, PTT, iFOD2 with priors and PTT with priors. Last column on the right: Prior weighting map, where brighter pixel denotes a higher usage of the prior. From top to bottom: The Arcuate fascicle (AF), the Cingulum (CG), the Corticospinal tract (CST), the Optic radiation (OR) and the Superior longitudinal fascicle I (SLF I).</figcaption></figure>

in clinical context, we used degraded data with a grid of 20×20×20 voxels and only using the 3000 s/mm² b-shell with 60 directions. The prior is computed with the original data. In addition, to study the effect of noise, we added Rician noise on the dMRI. The resulting images have an average SNR of 0dB, 10dB, 20dB, 30dB, 40dB and 50dB. The individual ODF were computed using the method described in [48] and represented in a spherical harmonic basis truncated at the 8th order. Then, EODF were calculated using the method described in this paper.

(ii) *multi-shell ODF*: In this second set of experiment, we computed the priors on fibers tracked with the multi-shell multi-tissue constrained spherical deconvolution (CSD) fiber ODFs. The idea is to simulate a clinical context by computing the ODF on clinical dMRI data and combined with anatomical priors calculated on high resolution data. The anatomical prior are calculated on 40 × 40 × 40 voxels DiSCo dMRI at SNR 20. From those images, fiber ODF were computing using the MRtrix implementation of the multi-shell multi-

tissue constrained spherical (msmt-CSD) method described in [49] using default parameters and using 8th order SH basis. Then, fibers were tracked using the MRtrix iFOD2 algorithm with a step size of 0.2 voxel and default parameters. Only the streamlines that reach the endings ROIs were kept. Finally, the TODs were extracted from those fibers using the method detailed in the previous section. After prior estimation on high resolution data, we estimated the ODF of a subject data with 20 × 20 × 20 grid and SNR= 20. Then, EODFs were obtained by combining the TOD with the clinical data fiber ODF.

2) *Metrics*: To quantify the results, we computed the Pearson correlation coefficient (r) between the ground truth connectivity matrices of the three DiSCo dataset and the resulting tractograms connectivity matrix. These matrices are computed by counting the numbers of streamlines that reach both start and end ROIs for each DiSCo fibers bundles.
## B. HCP data

1) **Data:** To test our method on in vivo high quality data, we used the HCP young adult data, acquired with 90 gradients distributed on 3 shells of b=1000, 2000, and 3000 s/mm² with 6 b= 0 acquisitions and a final resolution of 1.25mm³. Among them, the 105 pre-segmented HCP images were used, as described in II-A, to build fiber atlases of the Arcuate fascicle (AF), the Cingulum left (CG), the Corticospinal tract (CST), the Optic radiation (OR) and the Superior longitudinal fascicle I (SLF I) to study the influence of the addition of prior on bundles of different degrees of complexity and in different region of the brain. A cross validation was performed on the 105 subjects in which 100 subjects are used to build the atlas and the 5 remaining to perform the tractography algorithm. This process is repeated 8 times, giving us a total of 40 subjects to test the method on. For each image belonging to training set, ODFs were computed on the shell b= 3000 s/mm² using the method described in [48] and represented in a spherical harmonic basis truncated at the 8th order. Then, the TOD images for each bundle were registered on the subject space. Those priors were incorporating during the tracking process for each subject. The tractography algorithms were tested on the 5 remaining subjects.

In order to compare our method with state of the art anatomical priors methods, we used the BST algorithm ([35]), using the default parameters, to compute enhanced ODF for the CST and the AF bundles. We then tracked those bundles with the iFOD2 and PTT algorithms using the same parameters as before.

2) **Metrics:** To quantify the overall shape quality of the tractograms, we computed the generalized Dice score [50] between the segmentation of fibers obtained with tractography algorithm and the reference fibers. In order to measure the improvement on the numbers of streamlines reaching both ROI endpoints, we computed the percentage of streamlines that correctly connects the beginning and the end regions of bundles over the total number of streamlines in the reference bundle (noted valid streamlines VS). It is worth noting that this measure does not mean anything in absolute, being biased by the numbers of streamlines in the reference tracks. However, since the same numbers of streamlines is generated for each bundle, it allows us to compare the overall bundles shape and quality with and without the addition of priors.

## IV. RESULTS

### A. DiSCo challenge dataset

Figure 6 displays the correlation with the ground-truth connectivity of the two experiments, for the three DiSCo dataset. First, on Figure 6.a, we observed that lower the noise is, better are the correlation improvement for the experiment using the fiber pathways as ground-truth. Next, both algorithms performed better with the addition of priors for all level of noise and the three datasets. We can also notice that the improvement decreases when the quality of the data increases. In average, the correlation increases of 0.23 for an SNR of 10 dB but only 0.11 for an SNR of 50 dB. Moreover, the mean correlation improvement is of 0.19 and 0.17 for the iFOD2 and

PTT algorithms, respectively. At last, the average maximum reached is *r* = 0.83 for PTT and *r* = 0.82 for iFOD2. There are no major differences between the three dataset. For the second experiment using high resolution DW-MRI data, we showed an average correlation improvement of 0.15 for the PTT algorithm and 0.17 for the iFOD2 algorithm, in figure 6.b. Incorporating anatomical priors calculating on high resolution data during the tracking process of low resolution data increases the Pearson correlation coefficient with the ground-truth connectivity matrix.

### B. HCP data

Figure 5 displays the segmentation of AF, CG, CST, OR and SLF I obtained with iFOD2 and PTT, with and without prior. The results for the Valid streamlines and the Dice score are presented, respectively, in Figure ?? and Figure ???. For the Valid streamlines score, the addition of priors appears to always increase, on average, the numbers of streamlines that connect both end regions. However, enhanced tractography seems to also increase the variability of the results, as we found a mean variance across all bundles of 7.19% without the prior and 19.84% with. We also noticed that incorporating the priors in the PTT algorithm improve more the Valid streamlines than with the iFOD2 algorithm. Indeed, the average gain for iFOD2 is of 25.08% and 31.37% for PTT. The same observation can be made for the dice measurement, but as opposed to the Valid streamlines measurement, the variability seems to be decreased with the addition of prior. We notice that incorporating the anatomical priors improve the spatial coverage and that a higher fraction of streamline reach the endpoints of the bundle

## V. DISCUSSIONS

In this work, we proposed a novel method for creating and incorporating anatomical priors to any ODF-based tractography algorithms. We showed that incorporating anatomical priors on two state-of-the-art tractography algorithms improve the overall quality of tractograms when priors are extracted from a ground truth, but also from high quality tractograms. Indeed, first, on the DiSCo dataset, the proposed prior-based tractography methods obtained better correlation scores between the reconstructed bundles and the ground truth for every level of noise than the standards probabilistic algorithms, but this improvement decreases with the increasing SNR. Thus, the addition of priors to ODF images improve the quality of tractograms in presence of noise. However, due to the actual form of the extracted TOD it may not be useful on high SNR data. More particularly, in the DISCO experiment, we tried to mimic real clinical context with low spacial resolution data and an SNR of 20. In this case, our results showed a better correlation for two methods incorporating priors between the estimated fiber and the ground-truth, reducing false positive fibers. This could allow us to consider using this method on clinical data of average quality, informed by prior extracted on high quality data. On the HCP data, the same results were obtained with an increase of the fraction of valid streamlines numbers when anatomical priors were incorporating during tracking process. In the same time, due to the increase of
<figure><img src="image_7.png" /><figcaption>Fig. 6. Pearson correlation for the DiSCo experiments: (a) Fiber pathways ground-truth experiment with respect to SNR (b) Multi-Shell ODF experiment for the 3 DiSCo dataset. First column for PTT and second for iFOD2. In all cases, blue points denotes measure without prior and the orange ones, with prior.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. On the left, Dice score result for the HCP experiment. On the right, VS number result for the HCP experiment. From top to bottom: AF, CG, CST, OR, SLF I</figcaption></figure>

the dice score when adding prior, the overall shape of the bundle seems to be also modified. We also reported an increase in variability across subjects of the tractography with priors compared with the ones without priors.

The effect of the priors seems to be more effective in complex fibers configurations, like fanning and crossing fibers or in large curvature radius. Indeed, our method is able to a better delineation and estimation of precise details. For the CST, it is well known that the upper part, that fans into the cortex, is hard to estimate. Without prior, the iFOD2 and PTT algorithms only manage to estimate a portion of the CST. While, with the addition of anatomical prior, they achieve a more complete reconstruction of the fanning portion. Another example is the OR bundle, where the Meyer's loop, with its highly curving fibers, is still a tractography open challenge [51]. When guided by EODF, both tested algorithms show results that, while containing spurious streamlines, allow a better recovery of the anterior extent of the Meyer's loop. This result can be easily explained by the fact that in regions with crossing fibers, more information is taken from the prior. See the right column in Fig. 5 where the brightest regions

<figure><img src="image_9.png" /><figcaption>Fig. 8. Comparison with the BST algorithm</figcaption></figure>
correspond to the most complex fibers micro structure regions and denotes a higher usage of the prior.

It is important to note that the results presented in this paper use iFOD2 and PTT algorithms as proof of concept, but the main idea behind this work is to develop a method to improve the spatial coverage and the numbers of streamline reaching the endpoints of the bundle of any ODF-based tractography algorithm. In fact, future works should study the contribution of prior addition with other algorithm but also generalize this method to other diffusion model.

In [35], authors generated sharp TOD by projecting Dirac in the SH basis. However, our experiences have shown that this also generate a number of negative values which increase with the truncated SH basis order, due to the Gibbs phenomenon. Thus, to overpass this problem and remove the negative values, they use apodized delta function [52] proposed in [38], but at the cost of a loss of angular resolution. Using this PSF, we would have to truncated at a higher order the TOD SH basis to obtain the same angular resolution that we have with the 8th order of the Gaussian PSF [38]. So we decided to keep the Gaussian PDF and the ODF estimated without CSD, at the expense of less sharp ODF peaks, our interest here being to compare the tracking connectivity with and without prior. In the same way, in [35], priors and data are combined by an element-wise multiplication between the two sets of SH coefficients, followed by a normalization. However, our approach (i.e. using a Riemannian framework), allows us a better control over the weighting between the TOD and ODF and minimizes the potential angular smoothing. It also enables, for future works, more complex operation on TOD, such as interpolation between different set of priors or even TOD filtering to remove unnecessary information, in order to improve the prior quality, at the cost of more approximation during the TOD estimation and an increased complexity.

Other works that studies the idea of guiding tractography algorithms, like [32], [33], [34] or [35], use non-linear registration and atlasing techniques during the atlas, or template, creation process, preferring to do an average of the references images. This could smooth the variability across the atlas and could also be the source of errors that would be accumulated in the final atlas and would, in the end, inject a bias inherent to the atlas in the tractography process. In our approach, using a modified Guimond method [42], the references image are, at each step, iteratively registered with a non-rigid transformation onto the current space that become the reference space in the next step and it is only the average of the transformation that is computed at the end of each step. With this approach, all the variability information is accumulated in the atlas and the errors are not included in the final atlas. Thus, no bias is introduced in the tractography.

Also, using pre-segmented fibers to build the atlas may not be appropriate in order to study the global connectivity of the brain, as this approach only allows the enhanced tracking to be performed within a unique bundle. But, although not showed here, our method also permits, by concatenating the individual bundles priors, the construction of full brain atlases, and thus enable guidance of whole brain tracking. Therefore, it could be used for clinical studies where priors would be constructed

from high resolution dMRI data in order to guide tractography algorithm on poor, clinical, resolution data.

Finally, our method has proved to be able to increase the quality of the estimated fibers, in term of spatial extent and number of valid connection when utilized on healthy subjects. These promising results could already enable advances in an academic context for the study of the healthy brain. In a more clinical context, if used on patients suffering from pathologies that only slightly or moderately modify the white matter, such as psychiatric pathologies, this method could also bring improvements in fiber estimation, perhaps through a more in-depth study of α and β parameters. Research along these lines should be pursued in future work. However, when working on patient suffering from severe brain alteration, such as strokes, this method could produce a tractogram solely guided by the prior in the affected region, thus removing the specific microstructural modification and therefore may not be the most appropriate approach. In such scenarios, a solution might be to import other types of information into the priors. Either in the way the data and priors are weighted, using for example Apparent fiber density [(?)] of the data to potentially detect these brain modifications and thus prevent the fibers from being tracked solely on the priors, or directly in the estimation of the priors using other imaging modalities, such as myelin-sensitive relaxometry in order to obtain priors images describing these brain regions in greater detail.

## VI. CONCLUSION

In this paper, we developed a method of anatomical prior creation and combination, which is usable with any ODF-based tractography algorithm. The prior are computed on fiber atlases and expresses in the form of TOD, in order to characterize the brain variability and multiple fiber directions. Then, our prior is incorporated to guide the ODF-based tractography algorithm. Based on our results on DISCO and HCP data, incorporating our anatomical priors improve the fibers reconstruction, in terms of spatial extent and valid streamlines, especially in crossing fiber regions. Furthermore, we have shown that our method outperforms previous methods using anatomical priors in the tracking of the CST and the AF. Moreover, our approach could also greatly enhance the tractography in the context of clinical data, by incorporating priors estimated on high quality data, which could help for the study of neurological diseases.

## REFERENCES

[1] E.O. Stejskal and J.E. Tanner, "Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient," *The journal of chemical physics*, vol. 42, no. 1, pp. 288-292, 1965.

[2] F Dell'Acqua and JD Tournier, "Reconstructing fiber orientations with diffusion mri," *NMR Biomed*, 2017.

[3] D.K. Jones, *Diffusion MRI: Theory, Methods, and Applications*, Oxford University Press, 2010.

[4] S. Mori and P.C.M. 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*, vol. 15, no. 7-8, pp. 468-480, 2002.

[5] B. Jeurissen, M. Descoteaux, S. Mori, and A. Leemans, "Diffusion mri fiber tractography of the brain," *NMR in Biomedicine*, vol. 32, no. 4, pp. e3785, 2019.
[6] M. Mancini, S. B Vos, V. N Vakharia, A. G O'Keeffe, K. Trimmel, F. Barkhof, C. Dorfer, S. Soman, G. P Winston, C. Wu, et al., "Automated fiber tract reconstruction for surgery planning: extensive validation in language-related white matter tracts," *NeuroImage: Clinical*, vol. 23, pp. 101883, 2019.

[7] Kurt G. Schilling et al., "Tractography dissection variability: What happens when 42 groups dissect 14 white matter bundles on the same dataset?," *NeuroImage*, vol. 243, pp. 118502, 2021.

[8] G. Girard, R. Caminiti, A. Battaglia-Mayer, E. St-Onge, K. S. Ambrosen, S. F. Eskildsen, K. Krug, T. B. Dyrby, M. Descoteaux, J.P. Thiran, and G. M. Innocenti, "On the cortical connectivity in the macaque brain: A comparison of diffusion tractography and histological tracing data," *NeuroImage*, vol. 221, pp. 117201, 2020.

[9] P.J. Basser et al., "In vivo fiber tractography using DT-MRI data," *Magn Reson Med.*, vol. 44, pp. 625-32, 2000.

[10] Y. Gao, A. S Choe, I. Stepniewska, X. Li, M. J Avison, and A. W Anderson, "Validation of dti tractography-based measures of primary motor area connectivity in the squirrel monkey brain," *PloS one*, vol. 8, no. 10, pp. e75065, 2013.

[11] M. Descoteaux, R. Deriche, T.R. Knosche, and A. Anwander, "Deterministic and probabilistic tractography based on complex fibre orientation distributions," *IEEE transactions on medical imaging*, vol. 28, no. 2, pp. 269-286, 2008.

[12] J.D. Tournier, F. Calamante, and A. Connelly, "Robust determination of the fibre orientation distribution in diffusion mri: Non-negativity constrained super-resolved spherical deconvolution," *NeuroImage*, vol. 35, no. 4, pp. 1459-1472, 2007.

[13] T.E.J. Behrens, H. Johansen Berg, S. Jbabdi, M.F.S. Rushworth, and M.W. Woolrich, "Probabilistic diffusion tractography with multiple fibre orientations: What can we gain?," *NeuroImage*, vol. 34, no. 1, pp. 144-155, 2007.

[14] J.D. Tournier, R. Smith, D. Raffelt, R. Tabbara, T. Dhollander, M. Pietsch, D. Christiaens, B. Jeurissen, C.H. Yeh, and A. Connelly, "Mrtrix3: A fast, flexible and open software framework for medical image processing and visualisation," *NeuroImage*, vol. 202, pp. 116137, 2019.

[15] Christoff Gössl, Ludwig Fahrmeir, Benno Pütz, Ludwig M Auer, and Dorothee P Auer, "Fiber tracking from dti using linear state space models: detectability of the pyramidal tract," *NeuroImage*, vol. 16, no. 2, pp. 378-388, 2002.

[16] J. Pontabry, F. Rousseau, E. Oubel, C. Studholme, M. Koob, and J.L. Dietemann, "Probabilistic tractography using q-ball imaging and particle filtering: application to adult and in-utero fetal brain studies," *Medical image analysis*, vol. 17, no. 3, pp. 297-310, 2013.

[17] A. Stamm, O. Commowick, C. Barillot, and P. Perez, "Adaptive multi-modal particle filtering for probabilistic white matter tractography," in *International Conference on Information Processing in Medical Imaging*. Springer, 2013, pp. 594-606.

[18] S. N Sotiropoulos, L. Bai, P. S Morgan, C. S Constantinescu, and C. R Tench, "Brain tractography using q-ball imaging and graph theory: Improved connectivities through fibre crossings via a model-based approach," *Neuroimage*, vol. 49, no. 3, pp. 2444-2456, 2010.

[19] S. Jbabdi, M. W Woolrich, J? LR Andersson, and T.E.J. Behrens, "A bayesian framework for global tractography," *Neuroimage*, vol. 37, no. 1, pp. 116-129, 2007.

[20] I. Aganj, C. Lenglet, N. Jahanshad, E. Yacoub, N. Harel, P. M Thompson, and G. Sapiro, "A hough transform global probabilistic approach to multiple-subject diffusion mri tractography," *Medical image analysis*, vol. 15, no. 4, pp. 414-425, 2011.

[21] P. Poulin, D. Jörgens, P.M. Jodoin, and M. Descoteaux, "Tractography and machine learning: Current state and open challenges," *Magnetic resonance imaging*, vol. 64, pp. 37-48, 2019.

[22] S. Jbabdi and H. Johansen-Berg, "Tractography: where do we go from here?," *Brain connectivity*, vol. 1, no. 3, pp. 169-183, 2011.

[23] K. H Maier-Hein, P.F. Neher, J.C. Houde, M.A. Côté, E. Garyfallidis, J. Zhong, M. Chamberland, F.C. Yeh, Y.C. Lin, Q. Ji, et al., "The challenge of mapping the human connectome based on diffusion tractography," *Nature communications*, vol. 8, no. 1, pp. 1-13, 2017.

[24] J.D. Tournier, F. Calamante, A. Connelly, et al., "Improved probabilistic streamlines tractography by 2nd order integration over fibre orientation distributions," in *Proceedings of the international society for magnetic resonance in medicine*. John Wiley & Sons, Inc. New Jersey, USA, 2010, vol. 1670.

[25] D.B. Aydogan and Y. Shi, "Parallel transport tractography," *IEEE Transactions on Medical Imaging*, vol. 40, no. 2, pp. 635-647, 2021.

[26] J. Rafael-Patino, G. Girard, R. Truffet, M. Pizzolato, E. Caruyer, and J.P. Thiran, "The diffusion-simulated connectivity (disco) dataset," *Data in Brief*, vol. 38, pp. 107429, 2021.

[27] F. Rheault, J.C. Houde, and M. Descoteaux, "Visualization, interaction and tractometry: dealing with millions of streamlines from diffusion mri tractography," *Frontiers in neuroinformatics*, vol. 11, pp. 42, 2017.

[28] A. Daducci, A. Dal Palu, M. Descoteaux, and J.P. Thiran, "Microstructure informed tractography: pitfalls and open challenges," *Frontiers in neuroscience*, vol. 10, pp. 247, 2016.

[29] R.E. Smith, J.D. Tournier, F. Calamante, and A. Connelly, "Sift2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography," *NeuroImage*, vol. 119, pp. 338-351, 2015.

[30] A. Daducci, E.J. Canales-Rodriguez, H. Zhang, T.B. Dyrby, D.C. Alexander, and J.P. Thiran, "Accelerated microstructure imaging via convex optimization (amico) from diffusion mri data," *Neuroimage*, vol. 105, pp. 32-44, 2015.

[31] P.F. Neher, M. Gotz, T. Norajitra, C. Weber, and K.H. Maier-Hein, "A machine learning based approach to fiber tractography using classifier voting," in *International Conference on Medical Image Computing and Computer-Assisted Intervention*. Springer, 2015, pp. 45-52.

[32] A. Yendiki, P. Panneck, P. Srinivasan, A. Stevens, L. Zollei, J. Augustinack, R. Wang, D. Salat, S. Ehrlich, T. Behrens, et al., "Automated probabilistic reconstruction of white-matter pathways in health and disease using an atlas of the underlying anatomy," *Frontiers in neuroinformatics*, vol. 5, pp. 23, 2011.

[33] P.A. Cook, H. Zhang, S.P. Awate, and J?C. Gee, "Atlas-guided probabilistic diffusion-tensor fiber tractography," in *2008 5th IEEE International Symposium on Biomedical Imaging: From Nano to Macro*. IEEE, 2008, pp. 951-954.

[34] P.T. Yap, J. H Gilmore, W. Lin, and D. Shen, "Poptract: population-based tractography," *IEEE transactions on medical imaging*, vol. 30, no. 10, pp. 1829-1840, 2011.

[35] F. Rheault, E. St-Onge, J. Sidhu, K. Maier-Hein, N. Tzourio-Mazoyer, L. Petit, and M. Descoteaux, "Bundle-specific tractography with incorporated anatomical and orientational priors," *NeuroImage*, vol. 186, pp. 382-398, 2019.

[36] J. Wasserthal, P. Neher, and K.H. Maier-Hein, "Tractseg-fast and accurate white matter tract segmentation," *NeuroImage*, vol. 183, pp. 239-253, 2018.

[37] M. Descoteaux, E. Angelino, S. Fitzgibbons, and R. Deriche, "Regularized, fast, and robust analytical q-ball imaging," *Magnetic Resonance in Medicine*, vol. 58, no. 3, pp. 497-510, 2007.

[38] T. Dhollander, L. Emsell, W. Van Hecke, F. Maes, S. Sunaert, and P. Suetens, "Track orientation density imaging (todi) and track orientation distribution (tod) based tractography," *NeuroImage*, vol. 94, pp. 312-336, 2014.

[39] A. Goh, C. Lenglet, P.M. Thompson, and R. Vidal, "A nonparametric riemannian framework for processing high angular resolution diffusion images and its applications to odf-based morphometry," *NeuroImage*, vol. 56, no. 3, pp. 1181-1201, 2011.

[40] J. Rafael-Patino, G. Girard, R. Truffet, M. Pizzolato, E. Caruyer, and J.P. Thiran, "The diffusion-simulated connectivity (disco) dataset," *Data in Brief*, vol. 38, pp. 107429, 2021.

[41] A. L Alexander, J.E. Lee, M. Lazar, and A.S. Field, "Diffusion tensor imaging of the brain," *Neurotherapeutics*, vol. 4, no. 3, pp. 316-329, 2007.

[42] A. Guimond, J. Meunier, and J.P. Thirion, "Average brain models: A convergence study," *Computer vision and image understanding*, vol. 77, no. 2, pp. 192-210, 2000.

[43] V. Arsigny, O. Commowick, X. Pennec, and N. Ayache, "A log-euclidean framework for statistics on diffeomorphisms," in *International Conference on Medical Image Computing and Computer-Assisted Intervention*. Springer, 2006, pp. 924-931.

[44] R.O. Suarez, O. Commowick, S.P. Prabhu, and S.K. Warfield, "Automated delineation of white matter fiber tracts with a multiple region-of-interest approach," *Neuroimage*, vol. 59, no. 4, pp. 3690-3700, 2012.

[45] Maier-Hein Klaus Wasserthal Jakob, Neher Peter, "High quality white matter reference tracts (1.2.0) [data set]," *Zenodo*, 2018.

[46] X. Geng, T.J. Ross, W. Zhan, H. Gu, Y.P. Chao, C.P. Lin, G.E. Christensen, N. Schuff, and Y. Yang, "Diffusion mri registration using orientation distribution functions.," in *IPMI*, 2009, vol. 9, pp. 626-637.

[47] H. Akaike, "Information theory and an extension of the maximum likelihood principle," in *Selected papers of hirotugu akaike*, pp. 199-213. Springer, 1998.

[48] I. Aganj, C. Lenglet, G. Sapiro, E. Yacoub, K. Ugurbil, and N. Harel, "Reconstruction of the orientation distribution function in single- and
multiple-shell q-ball imaging within constant solid angle,” *Magnetic Resonance in Medicine*, vol. 64, no. 2, pp. 554–566, 2010.

[49] B. Jeurissen, J.D. Tournier, T. Dhollander, A. Connelly, and J. Sijbers, “Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion mri data,” *NeuroImage*, vol. 103, pp. 411–426, 2014.

[50] W.R. Crum, O. Camara, and D.L.G. Hill, “Generalized overlap measures for evaluation and validation in medical image analysis,” *IEEE transactions on medical imaging*, vol. 25, no. 11, pp. 1451–1461, 2006.

[51] G.A. Bertani, L. Bertulli, E. Scola, A. Di Cristofori, M. Zavanone, F. Triulzi, P.M. Rampini, and G.G. Carrabba, “Optic radiation diffusion tensor imaging tractography: An alternative and simple technique for the accurate detection of meyer’s loop,” *World Neurosurgery*, vol. 117, pp. e42–e56, 2018.

[52] D. Raffelt, J.D. Tournier, S. Crozier, A. Connelly, and O. Salvado, “Reorientation of fiber orientation distributions using apodized point spread functions,” *Magnetic Resonance in Medicine*, vol. 67, no. 3, pp. 844–855, 2012.