ConvTract: A Convolutional-based Framework for Tractography ConvTract: A Convolutional-based Framework for Tractography

—In recent years, a few data-driven tractography algorithms have been proposed to address the limitations of traditional methods, however, these approaches are time-consuming and lack generalizability. To tackle these shortcomings, we introduce ConvTract, a convolutional-based framework for automatic fiber tractography. In this approach, a Convolutional Neural Network (CNN) model is employed to estimate the fiber Orientation Distribution Function (fODF) from the resampled diffusion-weighted magnetic resonance imaging (DW-MRI) data within a local block around each voxel. The proposed model consists of a 3D CNN feature extractor followed by a 2D convolutional regressor head. Our pipeline represents the extracted fODFs by spherical harmonics coefficients. The predicted fiber ODFs can be used for both deterministic and probabilistic tractography. We showed that the deterministic tracking algorithm SD-Stream performs the most precise whole-brain tractography from the obtained fODFs. ConvTract was tested using the Tractometer tool and outperformed the classical and the data-driven tractography methods, by achieving the best scores in five of the eight metrics, including the percentage of valid connections and number of valid bundles. Moreover, the precision and robustness of the proposed framework in estimating fODFs and performing tractography was evaluated using both simulated and real diffusion data. The qualitative and quantitative assessments illustrate the superior performance of ConvTract over other available approaches.


I. INTRODUCTION
RACTOGRAPHY is a modeling technique to reconstruct the structures of the white matter (WM) streamlines using diffusion magnetic resonance imaging (dMRI) data [1]. Fiber tractography based on dMRI images is an important and practical tool for studying the white matter of the brain, which allows finding and analyzing fiber tracts and studying brain connectivity [2]. Another significant application of tractography is investigating neurological disorders and brain diseases that have changed the structure of nerve fibers from their normal state [3]. Standard methods for tracking fibers are often iterative approaches [4]. In these algorithms, the orientation of fiber streamline is initially estimated in each region of the WM using DWI data. Subsequently, the streamline is followed for a certain step size in the acquired orientation. This process is repeated until some stopping criteria are met, including exceeding a maximum valid streamline length, changing orientation sharper than an angular threshold (e.g., = 45°), and leaving white matter mask.
Estimating the local tissue orientations is a challenging task because a single voxel can have numerous different directions. Tractography algorithms can be performed globally [5][6][7] that reconstructs all brain fibers simultaneously or locally that is subdivided into deterministic [8], [9] and probabilistic [10], [11]. Deterministic approaches follow the strongest orientation, while probabilistic methods sample a direction aligned with the dominant orientation to continue the streamline tracing. Traditional methods use a mathematical model that fits to the diffusion data and estimates the fiber orientations. To name a few of these models, we can mention diffusion tensor models [8], [12], multi-tensor models [13], ball-and-sticks [14], FORECAST [15], CHARMED [16], and DIAMOND [17]. Other methods intend to reconstruct the fODFs that determine the probability of the presence of fibers orienting a specific direction on the sphere. Q-ball imaging [18], [19], diffusion spectrum imaging (DSI) [20], constraint spherical deconvolution (CSD) [21], [22], multi-shell multi-tissue (MSMT) CSD [23], constant solid angle Q-ball model (CSA) [24], orientation probability density function (OPDF) [25], and sparse fascicle model (FSM) [26] are in this class of algorithms. The work in [27] introduced a clustering method to estimate the principal diffusion directions (PDD) from the fODF data.
Choosing the apt model to obtain the fiber direction is indistinct [28], [29]. Other main disadvantages of such models include incorporating specific assumptions about the properties of white matter and DWI signals that can be inherently different across individual subjects and imaging devices [30]. Besides, the tractograms produced by traditional methods usually generate numerous invalid nerve fibers (high false positive) [28]. These approaches require rigorous high-level rules to monitor the generated streamlines and filter the invalid ones to increase the quality of the tractography [31]. A few works [32][33][34] have tried to address these limitations. Anatomically-Constrained Tractography (ACT) introduced in [32] utilizes anatomical information to provide priors for streamline generation and increase the biological plausibility of the tractography. The work in [33] proposed a filter framework for tractography that employs an unscented Kalman filter to fit a mixture of gaussian tensor models and propagate the streamline in a direction most consistent with the previous ones. In [34] a bundle-specific tractography method was introduced that estimates a template of streamlines and uses anatomical and orientation priors to track more accurate streamlines and better spatial coverage. However, fiber tracking from diffusion data is basically ambiguous, since it tries to capture global connectivity from local signals [28]. When several fiber tracts intersect in a bottleneck area and diverge afterward, the diffusion signal using lower order models in that region may represent a single orientation. Therefore, the accurate exit direction for each of the bundles cannot be distinguished correctly [35]. In fact, using the local signal, various directions are plausible in the tracking process, and that is the reason for producing a large amount of false positive streamlines [28].
In recent years, Machine Learning (ML) and Deep Learning (DL) methods have shown a great potential to perform fiber tractography and tackle the limitations of traditional approaches, including generating many false positives, dependency on the data properties, and not covering the whole brain properly. Although data-driven tractography pipelines have demonstrated the superiority of their performance to traditional algorithms, they are mostly time-consuming and fail to generate high-quality tractograms for various DWI data [4]. To address these limitations, we propose ConvTract 1 , a DLbased approach to map diffusion data to the fiber orientation distribution functions, which are used subsequently to perform whole-brain tractography. The proposed algorithm uses a novel CNN model to precisely acquire the spherical harmonics (SH) coefficients from resampled diffusion data. Since this pipeline does not incorporate the CNN model in the tracking process, it significantly reduces the runtime of the tractography, which makes our approach more practical compared to other datadriven methods. Besides, producing detailed fODFs by the model allows for a top-quality fiber tracking. The qualitative and quantitative evaluations indicate that our method outperforms current data-driven tractography pipelines.

II. RELATED WORK
The first ML-based tractography method was introduced in 1 Code available at https://github.com/smh-hosseiny/ConvTract 2015. A random forest (RF) classifier was used to obtain direction proposals for streamline tracking from raw diffusion data [36]. Several initial direction candidates are acquired in a local neighborhood of each voxel to find the best directions for a deterministic tractography. The model was evaluated using Tractometer tool [30] and showed promising results comparing to other submissions. The authors of [37] proposed a GRUbased Recurrent Neural Network (RNN) to map diffusion data to streamline orientation. The model exploits diffusion measurements in each voxel as well as the previous voxels along the streamline. The training was performed as a regression problem. The model was tested on the 2015 ISMRM challenge dataset [28] and scored using the tractometer tool. The scores were competitive with the best submissions. In [38], a RNN model was used to solve a multi-class classification problem instead of the conventional regression approach to find the fiber direction in each voxel. In this method, the diffusion data is initially resampled into 100 gradient directions evenly distributed on the unit sphere. Then, the three-dimensional ground-truth vectors are resampled on 724 directions in order to allow solving the direction estimation as a classification problem. The classification framework permits the usage of the predicted fODFs for both deterministic and probabilistic tractography. This model was also tested using tractometer tool and showed encouraging performance. The method introduced in [39] used a probabilistic MLP model to find the diffusion direction in each voxel based on the Fischer von Mises distribution in order to perform a probabilistic tractography. The model receives the diffusion data in a 3×3×3 cube around each voxel as well as the streamline direction in the previous voxel. The model was trained on a subject of the TractSeg dataset [40] and tested using the tractometer tool.
Unlike other data-driven approaches, Wasserthal et al. proposed a bundle-specific tracking algorithm using an encoder-decoder CNN [41]. The model generates bundlespecific fODF peaks called Tract Orientation Maps (TOM) and employs a deterministic tractography method to create the tractogram. The MSMT-CSD algorithm was utilized to obtain three direction peaks in each voxel and fed them into a U-net CNN to predict the bundle-specific TOMs for 20 bundles. The TOMs can be manipulated in different ways to perform tractography. The authors mentioned that deterministic tractography methods yield the best performance results. The model was trained on 105 HCP datasets [40]. The study in [42] introduced HAMLET (hierarchical harmonic filters for learning tracts from diffusion MRI). The model takes SH coefficients as input to predict bundle-specific spherical tensors that can be served as input to classical tractography methods. The authors didn't perform any evaluations and comparisons to other methods. Recently, Théberge et al. introduced the first successful deep reinforcement learning framework for tractography [43]. The agent, which is a neural network, finds the streamline directions at each step and performs fiber tracking based on the environment state, including the SH coefficients of the order 6 for each voxel and its six immediate neighbors (up, down, left, right, front, back), the tracking masks of these neighbors, and the last four streamline directions. A reward is calculated based on the alignment of the agent's action, which is the estimated direction of the streamline, and the ground truth peak direction of the fiber ODF obtained from diffusion data. This novel method was evaluated using the tractometer tool and demonstrated competitive performance to other approaches. Besides, the algorithm showed a decent generalization to unseen datasets.
Other works focus more on fODF or diffusion tensor estimation rather than the tractography process. Li et al.
proposed SuperDTI, a framework that uses a deep encoderdecoder network to map the DWI signal to diffusion tensor parameter maps [44]. Using as few as 6 DWI volumes for this purpose, significantly reduces the scan time and offers a feasible application in clinical routines. The work in [45] introduced a CNN model to estimate fODFs. It was illustrated that the proposed pipeline outperforms the ball-and-stick model in terms of angular error. The study in [46] used an MLP model with residual blocks to derive fODF from diffusion data, which was compared with some classical approaches using angular correlation coefficient (ACC). The study in [47] implemented an autoencoder model to generate priors for regularizing CSD algorithm to estimate fODFs. Lin et al. proposed a CNN model to generate fODF from down-sampled DWI data in a feasible acquisition runtime [48]. The model was tested using simulated and in vivo datasets and demonstrated superior performance over the MSMT-CSD algorithm. More recently, a CNN model was introduced to predict precise fODF that can be utilized for brain tractography as well as connectivity analysis [49]. Extensive experiments were carried out to assess the framework performance in comparison to classical and data-driven approaches using weighted average angular error (WAAE) and Jensen-Shannon divergence (JSD) metrics as well as evaluating the accuracy of fiber tract reconstruction by expert ratings.
For the sake of validation of fiber tractography studies and assessing their accuracy, the tractometer tool has been a prevalent choice for benchmarking the data-driven tractography pipelines. According to the reported scores from the tractometer tool [4], the aforementioned learning-based algorithms have some edges and downsides compared to one another, and no method outperforms others in all aspects. Hence there is no state-of-the-art automatic pipeline for fiber tractography in general.

III. MATERIALS AND METHODS
The proposed method has a similar pipeline to [41], where a DL-based model processes the input dMRI data and generates an output map containing the diffusion orientation information, which is eventually utilized to run a classical tractography algorithm. The pipeline of [41] takes diffusion direction peaks in each voxel as input and predicts bundle-specific fODF peaks. In contrast, we resample and normalize the input DWI data, then feed them into a CNN model to obtain the fODF in spherical harmonic bases. The predicted fiber ODFs can be used for either deterministic or probabilistic tractography. Fig.  1 shows the pipeline of end-to-end whole-brain tractography. In the following, we describe (1) our dataset, (2) introduced model, (3) implementation details, and (4) fiber tracking process.

A. Model
ConvTract contains a deep CNN model that maps the resampled diffusion data to the fiber ODF represented in SH of  Fig. 2 and consists of two components: a 3D CNN feature extractor and a proceeding 2D convolutional network that serves as the regressor head. The DWI data are initially resampled into 100 directions evenly distributed on a unit sphere. The resampled values from a 3×3×3 cube around each voxel are considered as the input of the model, so the input data is four-dimensional (3×3×3×100). The 3D CNN feature extractor bifurcates into two branches. The first branch acquires the spatial diffusion information in the entire 3×3×3 block and determines the principal diffusion directions in the local area. The second branch processes the diffusion data in a single voxel at the center of the block for which the SH coefficients are going to be predicted. Residual blocks [50] are utilized in both branches to enhance the performance of the network. The first component of the network creates a 40×40×3 feature map by concatenating the produced feature maps of both branches. The obtained representation comprises the diffusion information for the voxel of interest as well as for its surrounding voxels. Table  I describes the details of the 3D feature extractor architecture. The regressor head, which is shown in Fig. 3, is basically an image classifier with a regression output layer containing 45 neurons. This component of the network is composed of a few 2D convolutional layers to estimate the SH coefficients using the input feature map.

B. Dataset
We used two separate datasets to train the model, (1) five subjects from the Human Connectome Project (HCP) [51] dataset with 288 gradient directions and diffusion weightings of [1000, 2000, 3000] / 2 , and (2) Stanford HARDI dataset [52] with b-value of 2000 and 160 gradient directions. The fiber orientation distribution functions were obtained using the fODF reconstruction in Q-ball imaging with solid angle consideration algorithm implemented in MITK [53].

1) Preprocessing
For preprocessing the diffusion data, similar to [38], we normalized the input DWI volumes by the b0 volume and resampled the diffusion data on 100 directions on the unit sphere. In order to resample the data, the input signal at each of the gradient directions can be written as [54], [18]: where ( , ) are standard polar angles, 0 < < , 0 < < 2 , is the th SH basis, and is the corresponding SH coefficient. The matrix × can be defined as: From (1) and (2), we can see ×1 = × ×1 . Given the input diffusion signal with gradient directions, × is obtained for the 8th SH basis ( = 45). Next, the regularized pseudo-inverse of is calculated to acquire the SH coefficients . A unit sphere with 2 = 100 evenly distributed directions  is considered, then the matrix 2 × is computed for SH of the order = 8. Finally, the diffusion data are resampled on the 100 predefined directions using (3). The preprocessing pipeline for normalizing and resampling DWI data is implemented using DIPY library [55]. Ultimately, the brain and the white matter masks are computed for the tractography phase.

C. Implementation
The last layer of the network implements linear activations to handle the regression task, where rectified linear units (ReLU) [56] are used as nonlinearity for other layers. The model was trained using learning rate 10 −4 , and batch size of 64 on a Tesla K80 GPU. The Nadam optimizer [57] was used to train the weights with the Mean Absolute Error loss function defined in (4).
where is the number of samples and , are respectively the ground-truth and the regressed SH coefficients. The preprocessing stage for both training and tracking is performed for individual data batches separately to avoid memory limit error for high-resolution DWI data. The model is implemented using TensorFlow 2.4 and Python 3.8. Dropout layer is employed in the regressor head to prevent overfitting. Eighty percent of the data was used for training the model and the rest for validation. The validation loss was monitored to reduce the learning rate if the loss did not decrease for two consecutive epochs. The training continued for ten epochs until no further improvement was achieved for the loss function.

D. Fiber tracking
As illustrated in Fig. 1, having obtained the fODFs estimated by the model and the WM mask, both deterministic and probabilistic tractography methods can be performed to generate the tractogram. We tested different classical tractography algorithms for reconstructing fiber streamlines using the produced fODFs. The deterministic algorithm SD-Stream [58] implemented in MRtrix [59] showed the best performance according to the tractometer tool (Table V), so it was selected as the default tractography algorithm. The streamline minimum length was set to 20mm, the number of streamlines to 100K, and other configurations to their default values.

A. fODF estimation
To generate fODFs using the proposed model, the MSMT-CSD, and the SFM algorithm, we used two subjects from the HCP dataset (IDs: 100206,100307), which were not included in   Table II. ConvTract shows a more decent alignment to ground truth directions. the training dataset. MSMT-CSD and SFM are popular fODF reconstruction algorithms employed by many studies to produce ground-truth fODFs. The colored fractional anisotropy (CFA) maps and the fODFs in a 6×6 voxel region are illustrated in Fig. 4. The proposed model has shown competitive performance to MSMT-CSD and SFM algorithms. The HARDI-2013 Phantom dataset [60] was used to quantitatively evaluate our model's precision in fODF estimation. This dataset contains a simulated DWI data with 64 gradient directions, bvalue of 3000, and SNR of 30, along with the ground-truth peak directions. We estimated the SH coefficients and extracted up to five peaks in each voxel to measure the weighted average angular error (WAAE) metric [61] defined in (5). The peak extraction [62] was performed with relative peak threshold of 0.5 and separation angle of 45°. Where is the th groundtruth peak, is its magnitude, and is the predicted peak direction. The scores reported in Table II are obtained from [49], evaluating different fODF reconstruction methods using the WAAE metric. Our approach illustrates one of the best performances in retrieving the peak directions of the Phantom-2013 dataset. The results shown in Table II confirm the accuracy of the introduced model to regress precise fODFs from diffusion data. Fig. 5 represents the ground truth fiber orientations and the direction peaks produced by our model and SFM algorithm with the best performance among classical approaches in Table II. Our model shows a superior acquisition of the ground-truth peaks. Fig. 6. The raw DWI data, the estimated fODFs, and the reconstructed tractograms by ConvTract for (a) an in-house dataset, (b) a simulated data, and (c) a subject from the HCP dataset. DWI data with different configurations and preprocessing pipelines were selected to show the robustness of our pipeline to produce decent tractograms for various DWI datasets.

B. fODF estimation
Three DWI datasets with varied noise, number-of-shells, resolution, and preprocessing pipeline were used to assess the precision and robustness of the implemented framework to estimate the SH coefficients and perform tractography for diverse DWI datasets. Fig. 6 represents the estimated fODFs and the generated tractograms using the proposed pipeline. The test dataset includes (1) a preprocessed high-resolution 3-shell DWI data from the HCP dataset, (2) an in-house data with bvalue of 1000 and 69 gradient directions, and (3) a simulated DWI dataset (obtained from https://nitrc.org) with 108 gradient directions and b-values of [700, 2000] with compensated eddy currents and small head motion at SNR of 20, generated using the framework introduced in [63].

C. Tractometer
Subsequently, we scored our method using the tractometer tool, the high-resolution dataset of ISMRM 2015 tractography challenge was used for tractography. The generated and the gold standard tractograms are shown in Fig. 7. Tractography was performed using the SD-Stream algorithm with the estimated fODFs and the ground truth white matter mask provided along with the dataset. The scores reported from the tractometer tool in Table III are taken from [4] that has provided a comprehensive review of the data-driven tractography methods and their evaluation scores. Tractometer compares any given tractogram to a reference tractogram and outputs the following scores: • The number of valid and invalid reconstructed bundles (VB and IB) to represent the correctness of the generated tractogram. • The percentage of streamlines that belong to valid and invalid bundles (VC and IC), as well as the percentage of streamlines that are either shorter than the minimum valid streamline length or do not connect two regions of the cortex (NC) to evaluate the connectivity of streamlines.
• The percentage of ground truth voxels traversed by the produced streamlines (OL), the ratio of inaccurate voxels covered by reconstructed streamlines over the total number of ground truth voxels (OR), and the harmonic mean of OL and 1-OR (F1-score) to assess the voxel coverage.
According to Table III, ConvTract has shown a high performance, generating a precise 3D reconstruction of the fiber tracts. The produced tractogram recovers 24 of 25 gold standard bundles, 76% of the generated streamlines belonging to valid bundles, and the F1-score reaching 66. All of these results are the best scores achieved by the existing data-driven tractography pipelines. ConvTract reached the highest VB, VC, and F1 in addition to the lowest IC and OR, indicating that the proposed framework produced the most accurate tractogram with the best streamline connectivity and voxel coverage. The other notable advantage of our method is the shorter runtime of the proposed approach compared to other available methods for performing whole-brain tractography from the raw DWI data.

D. Design parameters examination
In this section, we analyze the effect of changing some hyperparameters and structure designs on the performance of our proposed pipeline.

1) Effect of architecture & loss function
We examined the impact of the architecture of the regressor head and the employed loss function on the tractography performance. The proposed 2D CNN with 2.43M parameters and a multi-layer perceptron (MLP) with 2.48M parameters were compared in regressing the SH coefficient from the feature map generated by the 3D network. Fiber tracking was performed using the SD-Stream algorithm with 100k streamlines and a minimum streamline length of 20mm. For each architecture, mean square error (MSE) and mean absolute error (MAE) loss functions were used for training the network on the datasets explained in section III. Models were evaluated using the tractometer tool, and the obtained scores are reported in Table IV

2) Effect of fiber tracking algorithm
We tested different classical fiber tracking algorithms using the fODFs produced by the CNN regressor, which was trained with the MAE loss function. The tractography algorithms were selected from various classes including probabilistic (iFOD1 [58] and iFOD2 [11]), deterministic (SD-Stream [58] and FACT [64]), and global (Gibbs [7]). According to the scores in Table V, the SD-Stream algorithm produced the most precise tractogram with the highest scores in the number of valid reconstructed bundles and valid connections and the lowest IC, NC, and OR. The iFOD1 algorithm also performed high-quality tractography and achieved the highest OL and F1-score. The other probabilistic algorithm, iFOD2, resulted in a tractogram with the least number of invalid bundles. Up to five peaks were extracted [62] in each voxel to serve as input for FACT algorithm, and the obtained tractogram received competitive OR and IC scores. The Gibbs global tractography algorithm generated a tractogram showing the highest number of valid bundles.

V. DISCUSSION
Most of automatic tractography pipelines incorporate the ML/DL model in an iterative fiber tracking process to estimate the streamline directions at each step. Using some data-driven models and a tractography process, we demonstrated that decoupling the fODF estimation from the fiber tracking phase would result in a better performance in regards to (1) the produced tractogram quality and (2) the generalizability of the framework to perform reasonable fODF reconstruction and fiber tractography for various DWI data. The results illustrated in the manuscript validated that compared to the conventional fODF reconstruction methods, the proposed model is well capable of mapping the diffusion measurements to precise fiber orientations. We visually confirmed the robustness of the proposed pipeline to obtain fODFs and generate tractograms for diverse DWI datasets, even of noisy and low-quality. The tractometer tool was employed to test our model and quantitatively compare it to the learning-based and nonlearning tractography approaches. The scores reported show the superior performance of ConvTract in respect of (1) the precision of tractography, (2) streamline connectivity, and (3) the spatial coverage of the whole brain. Finally, we analyzed the effect of different hyperparameters and fiber tracking algorithms on the proposed pipeline performance. Considering the scores obtained in this work, we explored, to some extent, the hyper-parameters of our framework and showed that the introduced design may have the best tractography performance according to the tractometer evaluation tool.

VI. CONCLUSION
In this paper, we presented ConvTract, a DCNN-based framework for automatic whole-brain fiber tractography. The introduced model comprises 3D and 2D convolutional networks to regress fODF coefficients, which are subsequently utilized to perform a deterministic tractography using SD-Stream algorithm. The introduced pipeline was evaluated in terms of fODF estimation accuracy and tractography precision and robustness. Moreover, the framework was scored using the Tractometer tool, and the obtained results showed that our approach outperforms the classical and the current state-of-theart data-driven tractography algorithms. Designing a more competent model to capture more details of the fODF, along with a more intelligent fiber tractography algorithm to exploit the estimated fODF more efficiently, would enhance the quality of the reconstructed tractogram in future.