Samuel St-Jean

Pierrick Coupé

Maxime Descoteaux

Non Local Spatial and Angular Matching : Enabling higher spatial resolution diffusion MRI datasets through adaptive denoising

Keywords: Diffusion MRI, Denoising, Block Matching, Noise bias, Dictionary learning

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Diffusion magnetic resonance imaging (MRI) is an imaging technique that allows probing microstructural features of the white matter architecture of the brain. Due to the imaging sequence used, the acquired images have an inherently low signal-to-noise ratio (SNR), especially at high b-values. Acquiring data at high bvalues contains relevant information and is now of great interest for connectomics [START_REF] Van Essen | The WU-Minn Human Connectome Project: An Overview[END_REF] and microstructure [START_REF] Alexander | Orientationally invariant indices of axon diameter and density from diffusion MRI[END_REF] studies. High noise levels bias the measurements because of the non-Gaussian nature of the noise, which in turn prohibits high resolution acquisition if no further processing is done. This can also lead to a false and biased estimation of the diffusion parameters, which impacts on the scalar metrics (e.g. fractional anisotropy (FA)), or in the fitting of various diffusion models (e.g. diffusion tensor imaging (DTI) and high angular resolution diffusion imaging (HARDI) models). This can further impact subsequent tractography and connectivity analysis if the spatially variable noise bias is not taken into account. Therefore, high SNR diffusion weighted images (DWIs) are crucial in order to draw meaningful conclusions in subsequent data or group analyses [START_REF] Jones | White Matter Integrity, Fiber Count, and Other Fallacies: The Do's and Don'ts of Diffusion MRI[END_REF]. This paper focuses on denoising techniques since they can be used on all existing data, without adding to the scanning time. They also can be readily applied to any already acquired dataset just like motion and eddy current corrections that are commonly applied on acquired datasets. One possible way to acquire higher quality data is to use better hardware, but this is costly and not realistic in a clinical setting. One can also use a bigger voxel size in order to keep the relative SNR at the same level, but at the expense of a lower spatial resolution or acquiring fewer directions to keep an acceptable acquisition time [START_REF] Descoteaux | Diffusion-Weighted Imaging[END_REF]. Averaging multiple acquisitions also increases the SNR, but this should be done either using Gaussian distributed noisy data [START_REF] Eichner | Real diffusion-weighted MRI enabling true signal averaging and increased diffusion contrast[END_REF] or in the complex domain to avoid the increased noise bias [START_REF] Jones | White Matter Integrity, Fiber Count, and Other Fallacies: The Do's and Don'ts of Diffusion MRI[END_REF].

With the advance of parallel imaging and acceleration techniques such as the generalized autocalibrating partially parallel acquisitions (GRAPPA) or the sensitivity encoding for fast MRI (SENSE), taking into account the modified noise distribution is the next step [START_REF] Aja-Fernández | Noise estimation in parallel MRI: GRAPPA and SENSE[END_REF][START_REF] Dietrich | Influence of multichannel combination, parallel imaging and other reconstruction techniques on MRI noise characteristics[END_REF].

The noise is usually modeled with a Rician distribution when SENSE is used and a non central Chi (nc-) distribution with 2N degrees of freedom (with N the number of receiver coils) when a Sum of Squares (SoS) reconstruction is used. If GRAPPA acceleration is also used with a SoS reconstruction, the degrees of freedom of the nc-distribution will vary between 1 and 2N [START_REF] Aja-Fernández | Noise estimation in parallel MRI: GRAPPA and SENSE[END_REF]. Some techniques have been specifically adapted by the medical imaging community to take into account the Rician nature of the noise such as non local means algorithms [START_REF] Coupé | An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images[END_REF][START_REF] Manjón | Adaptive non-local means denoising of MR images with spatially varying noise levels[END_REF][START_REF] Tristán-Vega | DWI filtering using joint information for DTI and HARDI[END_REF], Linear Minimum Mean Square Error estimator [START_REF] Descoteaux | Noise and signal estimation in magnitude MRI and Rician distributed images: a LMMSE approach[END_REF] Email address: samuel.st-jean@usherbrooke.ca (Samuel St-Jean) 1 Sherbrooke Connectivity Imaging Laboratory, Computer Science Department, Université de Sherbrooke 2 Unité Mixte de Recherche CNRS (UMR 5800), Laboratoire Bordelais de Recherche en Informatique, Bordeaux, France [START_REF] Brion | Noise correction for HARDI and HYDI data obtained with multi-channel coils and sum of squares reconstruction: an anisotropic extension of the LMMSE[END_REF], generalized total variation [START_REF] Liu | Generalized total variation-based MRI Rician denoising model with spatially adaptive regularization parameters[END_REF], a majorize-minimize framework with total variation denoising [START_REF] Varadarajan | A Majorize-Minimize Framework for Rician and Non-Central Chi MR Images[END_REF], maximum likelihood [START_REF] Rajan | Nonlocal maximum likelihood estimation method for denoising multiple-coil magnetic resonance images[END_REF] or block matching [START_REF] Maggioni | Nonlocal transform-domain filter for volumetric data denoising and reconstruction[END_REF]. Some methods [START_REF] Bao | Structure-adaptive sparse denoising for diffusion-tensor MRI[END_REF][START_REF] Becker | Adaptive smoothing of multi-shell diffusion weighted magnetic resonance data by msPOAS[END_REF][START_REF] Brion | Noise correction for HARDI and HYDI data obtained with multi-channel coils and sum of squares reconstruction: an anisotropic extension of the LMMSE[END_REF][START_REF] Gramfort | Denoising and fast diffusion imaging with physically constrained sparse dictionary learning[END_REF][START_REF] Lam | Denoising diffusion-weighted magnitude MR images using rank and edge constraints[END_REF][START_REF] Manjón | Diffusion Weighted Image Denoising Using Overcomplete Local PCA[END_REF][START_REF] St-Jean | Non Local Spatial and Angular Matching : a new denoising technique for diffusion MRI[END_REF][START_REF] Tristán-Vega | DWI filtering using joint information for DTI and HARDI[END_REF] have also been specifically designed to take advantage of the properties of the diffusion MRI signal such as symmetry, positivity or angular redundancy. Since the data acquired in diffusion MRI depicts the same structural information, but under different sensitizing gradients and noise realization, these ideas take advantage of the information contained in the multiple acquired diffusion MRI datasets.

We thus propose to exploit the structural redundancy across DWIs through a common sparse representation using dictionary learning and sparse coding to reduce the noise level and achieve a higher SNR.

Our method can be thought of a Non Local Spatial and Angular Matching (NLSAM) with dictionary learning. To the best of our knowledge, most recent state-of-the-art denoising algorithms either concentrate on modeling the nc-noise bias or the spatially varying nature of the noise in a Rician setting. Our method thus fills the gap by being robust to both of these aspects at the same time, as seen in Table 1. We will compare our method against one structural MRI method and two other publicly available algorithms : the Adaptive Optimized Non Local Means (AONLM) [START_REF] Manjón | Adaptive non-local means denoising of MR images with spatially varying noise levels[END_REF], which is designed for 3D structural MRI, the Local Principal Component Aanalysis (LPCA) [START_REF] Manjón | Diffusion Weighted Image Denoising Using Overcomplete Local PCA[END_REF] and the multi-shell Position-Orientation Adaptive Smoothing (msPOAS) algorithm [START_REF] Becker | Adaptive smoothing of multi-shell diffusion weighted magnetic resonance data by msPOAS[END_REF], both designed for processing diffusion MRI datasets. More information on each method features and parameters will be detailed later.

The contributions of our work are : i) Developing a novel denoising technique specifically tailored for diffusion MRI, which takes into account spatially varying Rician and nc-noise.

ii) Quantitatively comparing all methods on common diffusion MRI metrics.

iii) Quantifying the impact of denoising on local reconstruction models.

iv) Analyzing the impact of denoising on tractography with a synthetic phantom and a high spatial resolution dataset.

Noise type

AONLM LPCA msPOAS NLSAM Stationary Rician X X X X nc- X X X X Variable Rician X X X X nc- X X X X Use 4D angular information X X X X
Table 1: Features of the compared denoising algorithm, see Section 3.3 for an in-depth review of each method. The NLSAM algorithm is the only technique robust to both the spatially varying nature of the noise and the nc-bias at the same time.

Theory

We now define two important terms used throughout the present work. Firstly, a patch is defined as a 3D region of neighboring spatial voxels, i.e. a small local region of a single 3D DWI. Secondly, a block is defined as a collection of patches taken at the same spatial position, but in different DWIs, i.e. a block is a 4D stack of patches which are similar in the angular domain. The reader is referred to Fig. 1 for a visual representation of the process.

(a) A block is made of the b0 and some angular neighbors (b) Angular neighbors position on the sphere. The Block Matching Algorithm. Reusing the key ideas from the non local means, the block matching algorithm [START_REF] Dabov | Image denoising by sparse 3-D transform-domain collaborative filtering[END_REF] further exploits image self-similarity. Similar 2D patches found inside a local neighborhood are stacked into a 3D transform domain and jointly filtered via wavelet hard-thresholding and Wiener filtering. Combining these filtered estimates using a weighted average based on their sparsity leads to superior denoising performance than the non local means filter. The idea has been extended in 3D for MRI image denoising in [START_REF] Maggioni | Nonlocal transform-domain filter for volumetric data denoising and reconstruction[END_REF] and an adaptive patch size version for cardiac diffusion MRI image denoising was successfully employed by [START_REF] Bao | Structure-adaptive sparse denoising for diffusion-tensor MRI[END_REF].

The Dictionary Learning Algorithm. Dictionary learning has been used in the machine learning community to find data driven sparse representations [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF][START_REF] Mairal | Online dictionary learning for sparse coding[END_REF]. Typically, a set of atoms (called the dictionary) is learned over the data, providing a way to represent it with a basis tailored to the signal at hand [START_REF] Olshausen | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF]. This is analogous to using an off-the-shelf basis like the discrete cosine transform or wavelets, but in a data-driven manner which gives better results than using a fixed, general purpose basis. The main difference is that this is not necessarily a basis in the sense it can also be overcomplete, i.e. it can have more atoms than coefficients. Given a set of input data X = [x 1 , . . . , x n ] 2 R m⇥n organized as column vectors, the process is expressed as

min D,↵ 1 n n X i=1 ✓ 1 2 kx i D↵ i k 2 2 + k↵ i k 1 ◆ s.t. kDk 2 2 = 1 , (1)
where D 2 R m⇥p is the learned dictionary, is a trade-off parameter between the data fidelity term and the penalization on the coefficients ↵ = [↵ 1 , . . . , ↵ n ] 2 R p⇥n . A higher value of promotes sparsity at the expense of the similarity with the original data. The columns of D are also constrained to be of unit `2 norm in order to avoid degenerated solutions [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF][START_REF] Gramfort | Denoising and fast diffusion imaging with physically constrained sparse dictionary learning[END_REF][START_REF] Mairal | Online dictionary learning for sparse coding[END_REF]. The key is to devise a sparse representation to reconstruct structural information and discard noise, since the latter does not typically allow a sparse representation in any basis. Using a penalization on the `1 norm of the coefficients promotes sparsity, hence providing denoising through the regularized reconstruction.

This idea has led to inpainting and denoising applications from the machine learning community [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF][START_REF] Mairal | Online dictionary learning for sparse coding[END_REF] or even to accelerated acquisition process in the diffusion MRI community for diffusion spectrum imaging (DSI) [START_REF] Gramfort | Denoising and fast diffusion imaging with physically constrained sparse dictionary learning[END_REF].

Adjusting for various noise types. Although the original formulations of Eq. ( 1) and Eq. ( 4) assume stationary, white additive Gaussian noise, this is usually not true in diffusion MRI data, especially at high b-values and low SNR. The noise is usually modeled as following a Rician distribution or a nc-distribution when used with parallel imaging depending on the reconstruction algorithm and the number of coils N used during the acquisition [START_REF] Aja-Fernández | Noise estimation in parallel MRI: GRAPPA and SENSE[END_REF][START_REF] Dietrich | Influence of multichannel combination, parallel imaging and other reconstruction techniques on MRI noise characteristics[END_REF]. This introduces a bias which depends on the intensity of the signal that must be taken into account to recover the expected value of the original signal as shown in Fig. 2. Note, though, that other common preprocessing corrections, such as motion correction or eddy current correction, require interpolation and could thus change the theoretical noise distribution [START_REF] Veraart | Comprehensive framework for accurate diffusion MRI parameter estimation[END_REF]. The key idea lies in the fact that the nc-distribution is actually made from a sum of Gaussians, from which the Rician distribution is a special case with N = 1. By making the hypothesis that each of the 2N

Gaussian distributions share the same standard deviation G (Koay et al., 2009a), one can map a value m from a nc-distribution to an equivalent value m from a Gaussian distribution. We first compute estimates for G and ⌘ (which is an estimate of the signal value in a Gaussian setting). If ⌘ is below the noise floor due to a low local SNR, that is when ⌘ < G p ⇡/2, we set ⌘ = 0 instead of being negative as suggested by [START_REF] Bai | A framework for accurate determination of the T2 distribution from multiple echo magnitude MRI images[END_REF]. The next step uses the cumulative distribution function (cdf) of a nc-distribution and the inverse cumulative distribution function (icdf) of a Gaussian distribution to find the equivalent value m between the two distributions. This effectively maps a noisy nc-distributed signal m to an equivalent noisy Gaussian distributed signal m. See Fig. 3 for a synthetic example with a visual depiction of the process for mapping nc-signals to Gaussian distributed signals and (Koay et al., 2009a) for the original in-depth details.

Using a variance stabilization means considering the noise as additive white Gaussian noise, which allows any already designed technique for Gaussian noise to be used without any modification. The author of [START_REF] Foi | Noise estimation and removal in MR imaging: The variance-stabilization approach[END_REF] has shown that techniques with a Rician noise adaptation performed equally well as their Gaussian noise version through the use of a noise stabilization approach. The same idea has been directly applied with block matching [START_REF] Dabov | Image denoising by sparse 3-D transform-domain collaborative filtering[END_REF] for structural MRI in [START_REF] Maggioni | Nonlocal transform-domain filter for volumetric data denoising and reconstruction[END_REF]. The classical solution to remove the noise bias is to include the noise model into the denoising algorithm itself, as for example done in [START_REF] Descoteaux | Noise and signal estimation in magnitude MRI and Rician distributed images: a LMMSE approach[END_REF][START_REF] Becker | Adaptive smoothing of multi-shell diffusion weighted magnetic resonance data by msPOAS[END_REF][START_REF] Lam | Denoising diffusion-weighted magnitude MR images using rank and edge constraints[END_REF][START_REF] Manjón | Adaptive non-local means denoising of MR images with spatially varying noise levels[END_REF]. The drawback with this solution is that each method has to be rethought to account for any other noise type not considered in its original formulation.

Method

Adjusting for various noise types. In this paper, we will deal with both the Rician and nc-noise model on a voxelwise basis through the noise stabilization technique of (Koay et al., 2009a). This indeed makes our algorithm easily adaptable for any noise type by simply changing the pre-applied transformation as needed. We will use the Probabilistic Identification and Estimation of Noise (PIESNO) (Koay et al., 2009b) to estimate the stationary noise standard deviation. PIESNO works on a slice by slice basis, assuming the background noise as stationary along the selected slice, and is designed to find the underlying standard deviation of the Gaussian noise given its Rician or nc-nature. Voxels that are considered as pure background noise are found automatically by the method, using the fact that the squared mean of those voxels follows a Gamma distribution. Once automatically identified, the standard deviation G of those voxels can be computed and a new estimation of the Gamma distribution is made with the updated G until convergence.

In the case of spatially varying noise, we will use a method similar to [START_REF] Manjón | Adaptive non-local means denoising of MR images with spatially varying noise levels[END_REF], where the noise is estimated locally as

2 i = min ||u i u j || 2 2 , 8i 6 = j, (2)
with u i a noisy patch computed by subtracting a patch to a low-pass filtered version of itself and applying the local Rician correction factor of [START_REF] Koay | Analytically exact correction scheme for signal extraction from noisy magnitude MR signals[END_REF]. If the background was masked automatically by the scanner or is unreliable due to the scanner preprocessing for statistical estimation, we use a similar idea by computing the local standard deviation of the noise field as

i = std(u i low pass(u i )) (3)
If a noise map was acquired during the scanning session, it can be sampled directly to estimate the parameters of the noise distribution. In the event that such a map is unavailable, a synthetic one can be constructed by subtracting the image from its low-pass filtered counterpart (see Eq. 3). Since the noise is assumed as independent and identically distributed across DWIs, we apply a median filter on the 4D dataset to get a 3D noise field. Finally, a Gaussian filter with a full-width at half maximum of 10 mm is applied to regularize the noise field, which is then corrected for the more general nc-bias with the correction factor of [START_REF] Koay | Analytically exact correction scheme for signal extraction from noisy magnitude MR signals[END_REF]. A similar approach based on extracting the noise field with a principal component analysis was used by [START_REF] Manjón | Diffusion Weighted Image Denoising Using Overcomplete Local PCA[END_REF].

Locally Adapting the Dictionary Learning. In order to locally adapt the method to spatially varying noise, we add some more constraints to the classical formulation of Eq. ( 1). Firstly, since the measured signal in diffusion MRI is always positive, we use this assumption to constrain the positivity of the global dictionary D and the coefficients ↵, i.e. D 0, ↵ 0 as done in [START_REF] Gramfort | Denoising and fast diffusion imaging with physically constrained sparse dictionary learning[END_REF]. We fixed the regularization parameter for Eq. ( 1) in the same fashion as [START_REF] Mairal | Online dictionary learning for sparse coding[END_REF], that is = 1.2/ p m, with m = ps 3 ⇥ an, ps is the patch size and an the number of angular neighbors. Secondly, once D is known, we use Eq. (4) (see the next paragraph) iteratively until convergence with the constraint ↵ 0 and i = 2 i (m + 3 p 2m), where 2 i is the local noise variance found either with PIESNO or Eq. (3). In accordance with [START_REF] Candès | Enhancing sparsity by reweighted l1 minimization[END_REF], i is an upper bound on the `2 norm of the noise. We set the convergence as reached for ↵ i at iteration j when max |↵ i,j ↵ i,j 1 | < 10 5 or until a maximum of 40 iterations is realized.

Adaptive and Iterative `1 Minimization. While Eq. (1) will both construct the dictionary D and find the coefficients ↵, there are specialized iterative algorithms for solving `1 problems in order to yield sparser solutions [START_REF] Candès | Enhancing sparsity by reweighted l1 minimization[END_REF][START_REF] Daubechies | Iteratively reweighted least squares minimization for sparse recovery[END_REF]). An equivalent constrained formulation for solving

each column i of ↵ is min ↵i kw j ↵ i k 1 s.t. 1 2 kx i D↵ i k 2 2  i , (4)
where w j is a weighting vector penalizing the coefficients of ↵ i at iteration j. Eq. ( 4) can thus be iterated to further identify non zero coefficients in ↵ i by setting w j+1 = 1 |↵i|+✏ for the next iteration. The algorithm is then started with w 0 = 1 and ✏ = max|D T ⇠|. As suggested by [START_REF] Candès | Enhancing sparsity by reweighted l1 minimization[END_REF], ⇠ ⇠ N (0, 2 ) is set as a random Gaussian vector, which gives a baseline where significant signal components might be recovered.

While similar in spirit to Eq. (1), Eq. ( 4) provides a way to find the sparser representation for ↵ i while bounding the `2 reconstruction error.

To the best of our knowledge, our paper is also the first to use the noise variance as an explicit bound on the `2 reconstruction error. This yields a sparse representation while controlling at the same time the fidelity with respect to the original data, while the classical way is to use the variance as a regularized penalization factor.

The proposed algorithm

Our new NLSAM algorithm combines ideas from block matching and sparse coding. We will use the same kind of framework, but by replacing the thresholding part in the block matching with a step of dictionary denoising instead, allowing the penalization on the sparsity of the signal to regularize the noisy blocks. We also take explicit advantage of the fact that diffusion MRI data is composed of multiples volumes of the same structure, albeit with different noise realizations and contrasts across DWIs. This allows sparser estimates to be found, further enhancing the separation of the data from the noise [START_REF] Olshausen | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF]. Our method is thus composed of three steps :

1. Correct the noise bias if needed.

2. Find angular neighbors on the sphere for each DWI.

3. Apply iterative local dictionary denoising on each subset of neighbors.

Step 1. In case the noise is not Gaussian distributed, we first correct for the noise bias by finding the Gaussian noise standard deviation with PIESNO (Koay et al., 2009b). If the background is masked, we instead use Eq. (3). We then transform the DWIs into Gaussian distributed, noisy signals using the correction scheme of (Koay et al., 2009a).

Step 2. We find the angular neighbors for each of the DWIs. In this step, the local angular information is encoded in a 4D block of similar angular data, as seen in Fig. 1. The gradients are symmetrized to account for opposite polarity DWIs, which share similar structure to their symmetrized counterpart. The search is also made along all the shell at the same time, since structural information (such as sharp edges) is encoded along the axial part of the data. This encodes the similar angular structure of the data along the 4th dimension in a single vector.

Step 3. The dictionary D is constructed with Eq. ( 1) and the blocks are then denoised with Eq. ( 4). This step can be thought of finding a linear combination with the smallest number of atoms to represent a block. In order to adapt to spatially varying noise, each block is penalized differently based on the local variance of the noise. This enables the regularization to adapt to the amount of noise in the block, which is usually stronger as the acquired signal is farther from the receiver coils. Since each overlapping block is extracted, each voxel is represented many times and they are recombined using a weighted average based on their sparsity as in [START_REF] Maggioni | Nonlocal transform-domain filter for volumetric data denoising and reconstruction[END_REF][START_REF] Manjón | Diffusion Weighted Image Denoising Using Overcomplete Local PCA[END_REF]. For each voxel i with intensity v i contained in multiple overlapping blocks V j in neighborhood V , we set the final value of v i as

v i = P j2V v j (1 + ||V j || 0 ) P j2V 1 + ||V j || 0 , ( 5
)
where V is the same spatial position for voxel i across multiple blocks V j . This assumes that more coefficients in block V i also means more noise in the reconstruction. The `0 norm thus penalizes reconstructions with more coefficients and assigns a lower weight in that case for the overlapping weighted average. This third step is then repeated for all the DWIs. Since each DWI will be processed more than once with a different set of neighbors each time (see Fig. 1 for the block formation process), we obtain multiple denoised volumes of exactly the same data, but denoised in a different angular context. Once all the DWIs have been processed, we average the multiple denoised versions obtained previously in order to further reduce any residual noise. See Appendix A for an outline of the NLSAM algorithm as pseudocode. The result will be a denoised version of the input, through both dictionary learning and spatial and angular block matching.

Datasets and acquisition parameters

Synthetic phantom datasets. The synthetic data simulations are based on the ISBI 2013 HARDI challenge phantom3 and were made with phantomas4 . We used the given 64 gradients set from the challenge at b-values of 1000 and 3000 s/mm 2 . For simplicity, we will now refer to these datasets as the b1000 and the b3000 datasets. The datasets were generated with Rician and nc-noise profile, both stationary and spatially varying, at two different signal-to-noise ratios (SNR) for each case. In total, we thus have 8 different noise profiles for each b-value. The stationary noise was generated with SNR 10 and 20 and the spatially varying noise was generated with SNR varying linearly from 5 to 15 and from 7 to 20. The noise distributions were generated for each SNR by setting N = 1 for the Rician noise and N = 12 for the ncnoise. The noisy data was generated according to

Î = v u u t N X i=0,j=0 ✓ I p N + ✏ i ◆ 2 + ✏ 2 j , where ✏ i , ✏ j ⇠ N (0, 2 ), ( 6
)
where Î is the resulting noisy volume, N (0, 2 ) is a Gaussian distribution of mean 0 and variance 2 with = mean (b0) / SNR and mean(b0) is the mean signal value of the b = 0 s/mm 2 image. is a mask to create the noise distribution set to 1 in the constant noise case and as a sphere for the spatially varying noise case. For the spatially varying noise experiments, has a value of 1 on the borders up to a value of 3 at the middle of the mask, thus generating a stronger noise profile near the middle of the phantom than for the stationary (constant) noise case. As shown on Fig. 4, this results in a variable SNR ranging from approximately SNR 5 and SNR 7 in the middle of the phantom up to SNR 15 and SNR 20 for the spatially varying noise case. This noise mimics a homogeneous noise reconstruction as implemented by some scanners while still having a spatially varying noise map. The acquisition was made on the same scanner, but during another scanning session. No further processing nor denoising was done on this dataset for the reasons mentioned above. This can be thought of having a higher angular resolution at the cost of a lower spatial resolution for a comparable acquisition time.

Other denoising algorithms for comparison

We now present the various features and cases covered by the denoising algorithms studied in this paper. The Adaptive Optimized Nonlocal Means (AONLM) method [START_REF] Manjón | Adaptive non-local means denoising of MR images with spatially varying noise levels[END_REF] is designed for Rician noise removal in a 3D fashion and works separately on each DWIs volume. It also includes a Rician bias removal step and is able to spatially adapt to a varying noise profile automatically. We used the recommended default parameter of a 3D patch size of 3 x 3 x 3 voxels with the Rician bias correction in all cases. The Local Principal Component Analysis (LPCA) method [START_REF] Manjón | Diffusion Weighted Image Denoising Using Overcomplete Local PCA[END_REF] is also made to take into account the Rician noise bias and is spatially adaptive, but also uses the information from all the DWIs in the denoising process. We used the automatic threshold set by the method with the Rician noise correction for all experiments. Both AONLM and LPCA can be downloaded from the author's website5 .

The multi-shell Position-Orientation Adaptive Smoothing (msPOAS) algorithm [START_REF] Becker | Adaptive smoothing of multi-shell diffusion weighted magnetic resonance data by msPOAS[END_REF] was designed for both Rician and nc-noise, while also taking into account the angular structure of the data for adaptive smoothing. We discussed with the authors of msPOAS6 for their recommendations and using their suggestion, we set k ? = 12 and = 18. We also supplied the correct value for N and used the implemented automatic detection of the noise standard deviation from msPOAS. For the NLSAM algorithm, we used a patchsize of 3 x 3 x 3 voxels with 4 angular neighbors, which corresponds to the number of angular neighbors at the same distance on the sphere for each selected DWI. The value of N was given to the algorithm and the number of atoms was set to two times the number of voxels in a block for the dictionary learning part, which was repeated for 150 iterations. The other parameters were set as described in Section 3. As shown in Table 1, our method is designed to work on both stationary and spatially variable Rician and nc-noise.

The NLSAM algorithm is implemented in python and is also freely available. 7Finally, we quantitatively assess the performance of each method by comparing them against the noiseless synthetic data using i) The peak signal-to-noise ratio (PSNR) in dB and the structural similarity index (SSIM) on the raw data intensities [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF].

ii) The dispersion of the FA error, computed from a weighted least-square diffusion tensor model.

iii) The mean angular error (AE) in degrees and the discrete number of compartments (DNC) error for a region of crossings [START_REF] Daducci | Quantitative comparison of reconstruction methods for intra-voxel fiber recovery from diffusion MRI[END_REF][START_REF] Paquette | Comparison of sampling strategies and sparsifying transforms to improve compressed sensing diffusion spectrum imaging[END_REF].

iv) The Tractometer [START_REF] Côté | Tractometer: Towards validation of tractography pipelines[END_REF] ranking platform on deterministic tractography algorithms for the synthetic datasets. This platform computes global connectivity metrics, giving an insight on the global coherence of the denoised datasets in a tractography setting.

v) Tracking some known bundles on the high resolution in-vivo dataset and qualitatively comparing them to their lower spatial resolution counterpart.

Local models reconstruction and fiber tractography

The weighted least-square diffusion tensors were reconstructed using the default parameters of Dipy [START_REF] Garyfallidis | Dipy, a library for the analysis of diffusion MRI data[END_REF] to compute the FA values. We used the Constrained Spherical Deconvolution (CSD) [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF] with a spherical harmonics of order 8 to reconstruct the fODFs and extract the peaks subsequently used for the deterministic tracking. To compute the fiber response function (frf), we used all the voxels in the white matter that had an FA superior to 0.7. If less than 300 voxels meeting this criterion were found, the FA threshold was lowered by 0.05 until the criterion was met. See Sections 4.2 and 5.3 for more information about the bias introduced in the FA. For the synthetic datasets, the tracking was done inside the white matter mask and the seeding was done from the bundles extremities to mimic seeding from the white-gray matter interface [START_REF] Girard | Towards quantitative connectivity analysis: reducing tractography biases[END_REF]. We used 100 seeds per voxels to allow sufficient bundle coverage, a stepsize of 0.2 mm and a maximum angle deviation of 60 degrees. The other parameters used were the defaults supplied by the tractometer pipeline [START_REF] Côté | Tractometer: Towards validation of tractography pipelines[END_REF]. The in-vivo datasets deterministic tracking was made with the technique of [START_REF] Girard | Towards quantitative connectivity analysis: reducing tractography biases[END_REF] by seeding from the white matter and gray matter interface with the particle filtering and generating approximately 1 million of streamlines. White matter masks were created by segmenting a T1 image with FSL FAST 8 from the same subject and then registered with ANTs 9 to each in-vivo dataset. The bundles were finally automatically segmented using the White Matter Query Language (WMQL) [START_REF] Wassermann | On describing human white matter anatomy: The white matter query language[END_REF] Tract Querier tool with regions obtained from a T1 white matter and gray matter parcellation. This atlas-based automatic dissection method extracts fiber bundles automatically using anatomical definitions in a reproducible manner for all methods, as opposed to the traditional way of manually defining including and excluding ROIs to define bundles.

Visualization of fODFs, peaks and tractography was made using the fibernavigator All denoised methods were applied on the high spatial resolution 1.2 mm dataset. We also show an acquisition of the same subject at 1.8 mm for visual comparison. Our NLSAM algorithm is able to recover more crossings from the 3 way junction of the SLF, the CST and the CC as shown in the yellow and white boxes. While the 1.8 mm dataset is less noisy, its lower spatial resolution also means that each voxel contains more heterogeneous tissues and mixed diffusion orientations. The 1.2 mm denoised dataset shows more crossings without the averaging effect of the larger voxel size. For a comparable acquisition time, the denoised high resolution dataset has more information than its lower resolution counterpart without processing.

Fig. 7 shows the PSNR and SSIM for the SNR 10 (stationary noise) and SNR 15 (spatially varying noise) synthetic datasets. The LPCA algorithm performs best in term of PSNR on the Rician noise case, but attains a lower score for nc-noisy datasets. The same trend is seen for AONLM and msPOAS algorithms, where the SNR 15 nc-case is the hardest test case. In contrast, our NLSAM technique is above 30 dB for the PSNR and 0.9 for the SSIM in most cases, with a relatively stable performance amongst the majority of tested cases. We also note that even though msPOAS is made to adjust itself to nc-noise, the fact that the algorithm does not account for the intensity bias makes the perceptual metrics drop for the nc-noise cases.

(a) PSNR b1000 (b) PSNR b3000 (c) SSIM b1000 (d) SSIM b3000

Bias introduced in the FA

As shown by the FA difference map on Fig. 8, our NLSAM method commits a small FA error locally with a smaller maximum error than the compared methods. Voxelwise underestimation is denoted in blue and overestimation in red, where white means the computed value is close to the expected value. The noisy data largely overestimates the FA values for the synthetic datasets, while other denoising methods underestimate the real FA value most of the time. On the b1000 datasets, NLSAM has the smallest spread of FA error. The effect of stabilizing the data prior to denoising can also be seen by the stable FA median error committed by NLSAM across all noise types. For the b3000 datasets, the need to correct the intensity bias caused by the noise becomes more important, as seen by the increased error in underestimating the correct FA value for most methods. For the spatially varying Rician noise case, our method commits the lowest overestimation, as opposed to AONLM and LPCA, which are developed for this particular noise case. It is also important to note that in contrast to the other methods, msPOAS does not explicitly correct for the intensity bias by design, but rather leaves this correction to subsequent processing steps. The SNR 15 nc-noise case is where all the methods make the biggest error, as they reduce the variance but still suffer from a large bias in FA. Overall, our method restores the value of the FA for large bundles more accurately. We also see that most methods make their largest error near the partial volume ball mimicking cerebrospinal fluid (CSF).

Impact on angular and discrete number of compartments (DNC) error

We now study the angular error and the mean relative error in the discrete number of compartments (DNC) [START_REF] Daducci | Quantitative comparison of reconstruction methods for intra-voxel fiber recovery from diffusion MRI[END_REF][START_REF] Paquette | Comparison of sampling strategies and sparsifying transforms to improve compressed sensing diffusion spectrum imaging[END_REF]. The mean relative discrete number of compartments error is defined as DN C i = 100 ⇥ |P itrue P iest |/P itrue for voxel i, P itrue and P iest is the number of crossings respectively found on the noiseless dataset and on the compared dataset. All metrics were computed on the voxels containing at least two crossings fibers on the noiseless dataset shown previously in Fig. 5.

Fig. 9 shows the distribution of the angular error and of the DNC error found in the region studied in addition to the mean angular error. All of the denoising algorithms have a lower median and mean angular error than the noisy datasets. In addition, the NLSAM denoised datasets have an almost equal or lower angular error than the other denoising methods, but with a lower maximum error most of the time as

shown by the smaller whiskers. For the b1000 dataset DNC error, all three of AONLM, LPCA and NLSAM improve on the noisy dataset for the Rician noise case as they are devised for this kind of data. LPCA also has a better performance than the other two for the spatially varying Rician noise case, while NLSAM has a lower mean DNC error for both of the nc-noise case. The effect of the intensity bias is also seen on msPOAS, where the DNC error is always lower than the noisy dataset, but also higher than all the other methods which take into account the intensity bias. The b3000 dataset is much harder, where no method seems to have a clear advantage in all cases over the others. One interesting thing to note is that the noisy dataset has a low DNC error for both of the Rician noise case, but the confidence interval indicates it is in the same range as the denoised datasets.

Impact on tractography

We now show how denoising techniques impact tractography by evaluating the number of valid bundles (VB), invalid bundles (IB) [START_REF] Côté | Tractometer: Towards validation of tractography pipelines[END_REF] and the valid connection to connection ratio (VCCR) [START_REF] Girard | Towards quantitative connectivity analysis: reducing tractography biases[END_REF] found by the tracking algorithm. A valid bundle is defined as connecting two ROIs in the ground truth data while an invalid bundle is a connection made between two ROIs which is not supported by the ground truth data. The valid connection to connection ratio is the total of valid connections (VC) over the sum of valid and invalid connections (IC), i.e. V CCR = V C V C+IC . A good denoising algorithm should find a high number of valid bundles, a low number of invalid bundles and a high percentage of valid connection to connection ratio.

Deterministic tractography on the synthetic phantom. Table 2 shows the results of deterministic tractography on the SNR 10, 15 and 20 synthetic datasets for both b1000 and b3000. The noiseless b1000 dataset had 25/27 valid bundles, 55 invalid bundles and a valid connection to connection ratio of 65% and the noiseless b3000 dataset had 27/27 valid bundles, 40 invalid bundles and a valid connection to connection ratio of 68%. One of the first thing to note is that even though the noisy dataset always has a high number of valid bundles, it is at the price of a huge number of invalid bundles. Moreover, the valid connection to connection ratio is systematically lower for the SNR 10 datasets than any of the denoising methods. This indicates that only looking at the number of valid and invalid bundles does not show how many streamlines reached each region since only at least one streamline is required to make a connection, thus counting as a valid bundle. Another observation is that denoising helps controlling the number of invalid bundles and gives a better valid connection to connection ratio in most cases over the noisy data. For the SNR 15 cases, NLSAM has the highest number of valid bundles in almost all cases, but at the price of a larger number of invalid bundles at lower SNR. Another interesting trend is the tradeoff between valid bundles and invalid bundles :

AONLM and LPCA both manage to get a lower number of invalid bundles, but also tend to have a lower number of valid bundles than msPOAS or NLSAM overall.

For the SNR 20 stationary noise cases, all methods are close in valid bundles with some difference in the number of invalid bundles. This shows that tractography could benefit from variable tracking parameters instead of fixed values depending on the preferred trade-off for the task at hand [START_REF] Chamberland | Real-time multi-peak tractography for instantaneous connectivity display[END_REF]. Tracking the real data. We now look at tractography on the in-vivo high spatial resolution dataset and its clinical spatial resolution counterpart of the same subject previously shown on Fig. 6. The high spatial resolution dataset at 1.2 mm isotropic has 40 unique gradient directions while the lower spatial resolution dataset at 1.8 mm isotropic has 64 unique gradient directions for a comparable acquisition time. The background is masked by the scanner and has a spatially varying Rician noise profile due to the SENSE reconstruction, which is the specific noise case covered by the AONLM and LPCA denoising algorithm. We use the deterministic tractography algorithm from [START_REF] Girard | Towards quantitative connectivity analysis: reducing tractography biases[END_REF], which considers anatomical constraints for more anatomically plausible tractography. Fig. 10 (which is displayed in landscape to present all the bundles at once) shows from top to bottom the left arcuate fasciculus (AF), the inferior fronto-occipital fasciculus (IFOF) and the corticospinal tract (CST) as dissected automatically by the Tract querier [START_REF] Wassermann | On describing human white matter anatomy: The white matter query language[END_REF]. The noisy 1.2 mm AF stops prematurely in the frontal part of the bundle, while the 1.8 mm noisy AF misses the temporal lobe. In contrast, the streamlines from the NLSAM denoised bundle go further into the temporal lobe. Also note how the right IFOF has a better coverage for all the 1.2 mm datasets and more fanning near the front of the brain than the noisy 1.8 mm dataset. We also see that the left IFOF is thinner than its right counterpart, but most of the bundles tracked from the denoised datasets produces less spurious tracks while keeping the anatomical details. The LPCA denoised IFOF stops prematurely for the left posterior part of the bundle, possibly because of a lost crossing region along the fibers during the denoising process. The CST does show some commissural fibers through the pons in the noisy 1.8 mm dataset, while they are present but look like spurious fibers on the noisy 1.2 mm dataset. AONLM can recover some of those commissural fibers, while NLSAM is the only algorithm which recovers them in addition to richer fanning near both sides of the motor cortex. We also show a T1 weighted image aligned in the diffusion space for anatomical reference. Top : The left arcuate fasciculus. Note how the denoised NLSAM arcuate fasciculus goes further into the frontal and temporal region than both of its noisy 1.2 mm and 1.8 mm counterparts. Middle : The inferior fronto-occipital fasciculus. The AONLM denoised bundle has a denser part for the right IFOF while the LPCA bundle stops prematurely for the left IFOF, possibly due to a missing crossing along the bundle. Bottom : The corticospinal tract. We see that NLSAM recovers the commissural fibers in the pons from the noisy 1.2 mm dataset, which are not even present on the noisy 1.8 mm dataset nor on the other denoising algorithm's bundles. NLSAM also recovers more fanning to both sides of the brain than all the compared methods.

Discussion

Enhancing the raw data

We quantitatively showed in Fig. 7 that denoising restores perceptual information when compared to the unprocessed noisy data. Taking the spatially varying aspect and the particular nature of the noise into account is also important since modern scanners implement parallel imaging which changes the nature of the noise [START_REF] Dietrich | Influence of multichannel combination, parallel imaging and other reconstruction techniques on MRI noise characteristics[END_REF], leading to a lower performance for denoising methods not fully taking into account the introduced bias. Fig. 6 shows that this is also qualitatively true for in-vivo data, where denoising visually restores information in regions heavily corrupted by noise. While perceptual metrics might indicate the performance of an algorithm, one must remember that the relative signal difference is of interest in diffusion MRI, which is not fully captured by perceptual metrics like the PSNR or the SSIM. One is also usually interested in diffusion MRI metrics as opposed to perceptual information brought by the raw diffusion MRI datasets. AONLM is able to remove most of the noise, but still shows some residuals near the inferior part of the brain, possibly due to only considering the 3D volumes separately, which means that the algorithm can not benefit from the additional angular information brought by diffusion MRI. LPCA can restore visual information and sharp edges from the noisy dataset, but the region in the pons, where the noise level is higher and crossing fibers are more complex, also seems to be piecewise constant. This might arise from the fact that the algorithm uses all of the DWIs at once for its PCA decomposition step and treats all intensities at the same in the noise removal step. As for NLSAM, the algorithm only works in the local angular domain, thus exploiting similar contrast and redundant edge structure under different noise realization. msPOAS also uses a similar idea, where the angular similarity is weighted according to the Kullback-Leibler divergence to control the importance of dissimilar intensities in the denoising process.

Nevertheless, these perceptual metrics show that denoising improves upon the noisy data, but one should also look at metrics derived from the studied object of interest i.e. tensor or fODF derived metrics, since high perceptual metrics might also reflect blurring of diffusion features, which is the main interest in this type of acquisition rather than the perceived quality.

Impact of the stabilization algorithm on the compared denoising methods

Fig. 11 shows the FA map when the compared denoising algorithm are applied on the stabilized data with the algorithm of (Koay et al., 2009a). For this experiment, we consider a voxel as being degenerated if its FA is exactly 0. The first thing to note is that the algorithm only reprojects the noisy data on plausible Gaussian distributed values and does not do any denoising. While we used here the algorithm of (Koay et al., 2009a) to correct the noise bias, another interesting approach consist of producing realvalued datasets as shown in [START_REF] Eichner | Real diffusion-weighted MRI enabling true signal averaging and increased diffusion contrast[END_REF]. This approach does not require estimation of 2 or N , but instead use information contained in the phase of complex-valued acquisitions. Secondly, all of the other compared denoising algorithms produce some invalid voxels on the raw dataset, while having less degenerated voxels on the stabilized dataset as shown in Table 3. Nevertheless, only our NLSAM algorithm does not produce any degenerated FA voxel on the in-vivo dataset. As tractography might rely on a thresholded FA mask [START_REF] Chamberland | Real-time multi-peak tractography for instantaneous connectivity display[END_REF], any missing white matter voxel will end the tractography early and produce anatomically invalid tractography. In the same way, computing FA based statistics in search of group differences inside a white matter mask might lead to erroneous conclusions when degenerated voxels are present. This undesirable side effect should be avoided when possible by choosing a method producing a low number of invalid voxels, such as NLSAM. The percentage of degenerated voxels is indicated in parenthesis for each mask, where a voxel is considered degenerated if its FA value is exactly 0.

Reducing the diffusion metrics bias

Fig. 8 shows that knowing where errors are committed gives a better view of how denoising improves upon the noisy data. We see that our NLSAM algorithm actually has a smaller maximum error in underestimating the FA most of the time while other methods both over and underestimate the real FA value and make larger errors near CSF or at borders with the background. This could indicate that they are subject to problems with partial volume effect, which seems less important for NLSAM.

While stabilizing the data alleviates the FA underestimation problem in most cases, it also helps to reduce the number of degenerated voxels in the in-vivo data as shown in Fig. 11. Both AONLM and msPOAS produce less degenerated FA voxels on the stabilized dataset as shown in Table 3, while NLSAM does not produce any degenerated voxel. In contrast, the noisy data and LPCA have an increased number of degenerated FA voxels, which might be caused by the diffusion signal being near the noise floor, thus producing a flat profile which is not properly recovered in this case. Reducing the FA bias and avoiding degenerated voxels by including denoising in the processing pipeline could improve the statistical analysis in along-tract metrics [START_REF] Colby | Along-tract statistics allow for enhanced tractography analysis[END_REF] when looking for group differences.

Restoring the coherence of local models

The CSD algorithm relies on the estimation of the fiber response function (frf), which in turn relies on the diffusion tensor. To estimate the frf, one must select voxels containing only a single fiber population. One way to do this is to estimate it from voxels with a high FA, usually with FA > 0.7 [START_REF] Descoteaux | Deterministic and probabilistic tractography based on complex fibre orientation distributions[END_REF][START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF]. We observed that for the SNR 10 dataset with nc-noise, the noisy dataset, AONLM and LPCA could not find as much single fiber voxels based on the FA threshold method as msPOAS or NLSAM since their reconstructed tensors have an inherently lower FA. This in turns impacts deconvolution since the estimates used for the deconvolution kernel are less stable, and lowering the FA threshold too much might violates the single fiber assumption, which is crucial for the CSD algorithm. One way to circumvent this could be by using the method of [START_REF] Tax | Recursive calibration of the fiber response function for spherical deconvolution of diffusion MRI data[END_REF], which is based on a peak amplitude criterion instead of an FA threshold to identify single fiber voxels.

Fig. 9 shows that msPOAS and NLSAM have larger angular error than AONLM or LPCA, but this does not seem to impact much the number of valid bundles found by deterministic tractography. Indeed, the noisy data has the largest angular error in all cases, but still has a high number of valid bundles in most cases. This also suggests that a large overestimation or underestimation of fiber crossings (as reflected by the DNC error) has a higher impact on tractography. Both LPCA and msPOAS have a lower number of valid bundles than AONLM or NLSAM, which both have a rather symmetric under and over estimation of the number of peaks. This means that an overall estimation or underestimation of the number of crossings bias tractography, as it tends to follow false crossings or stops prematurely due to a lack of crossings, while a distributed error is not skewed toward these effects.

In Fig. 6, we see that denoising restores coherence in regions of crossing fibers that were lost on the noisy dataset or not even present in the lower spatial resolution 1.8 mm dataset due to a smaller voxel size. We also see that NLSAM restores more coherent crossings than the other denoising methods in the junction of the SLF and the CST, whereas the noisy dataset only shows incoherent crossings. This actually enables tractography algorithms to reconstructs tracts which are in agreement with the expected anatomy. In the same amount of acquisition time, one can thus acquire higher spatial resolution DWIs and get better angular information by post-processing the acquired data with denoising.

Limiting spurious fibers from tractography

We studied the impact of denoising techniques on deterministic tractography on a synthetic dataset in Section 4.4. One often has to choose between having a high number of valid bundles and invalid bundles, or having less valid bundles and at the same time reducing the amount of invalid bundles. The noisy dataset always reaches a high number of valid bundles, but also at the price of having the highest number of invalid bundles most of the time. Our NLSAM algorithm shows a good balance between the number of valid and invalid bundles at low SNR, especially for the spatially varying noise case. This is always a tradeoff that one has to choose as seen in the ISMRM 2015 tractography challenge11 .

For example, the LPCA algorithm has always a low number of invalid bundles, but also the lowest number of valid bundles for the spatially varying noise case. In opposition, NLSAM has a high number of valid bundles, but also a high number of invalid bundles most of the time.

Regarding the deterministic tractography, changing the stepsize or the maximum curving angle would give different results in terms of connectivity metrics, indicating that the tractography algorithm and chosen tractography parameters have a non negligible influence on the end results [START_REF] Chamberland | Real-time multi-peak tractography for instantaneous connectivity display[END_REF][START_REF] Girard | Towards quantitative connectivity analysis: reducing tractography biases[END_REF]. We also used a seeding strategy of 100 seeds per voxel from the ROIs at each bundles endpoints to ensure a maximal number of valid bundles, which promotes a high number of valid bundles for each dataset. This shows that the missed bundles are hard to recover or not supported by the data itself, as opposed to being missed because of inadequate seeding [START_REF] Côté | Tractometer: Towards validation of tractography pipelines[END_REF]. On the other hand, this can artificially increase the number of invalid bundles, which could be reduced by reducing the number of seeds per voxel. Since automatic tractography pruning techniques such as [START_REF] Côté | Cleaning up the mess: tractography outlier removal using hierarchical QuickBundles clustering[END_REF] might help reduce the number of spurious tracks, this would indicate that having a higher number of valid bundles would be preferable since invalid bundles could be potentially removed afterward. In contrast, a low number of valid bundles cannot be circumvented with further postprocessing. Nevertheless, denoising increases the valid connection to connection ratio and reduces the number of invalid bundles, thus bringing confidence in the validity of the tractography results when compared to the noisy datasets.

For the in-vivo dataset tracking shown in Fig. 10, we see that tractography benefits from higher spatial resolution acquisitions, but that the produced tracts are slightly more noisy. Combining the high spatial resolution, highly noisy dataset with a denoising algorithm at the beginning of the processing pipeline gives more anatomically plausible tracts in the end. The AF and CST produced by the NLSAM denoised dataset are both more anatomically plausible than their noisy or lower spatial resolution counterpart, which have less fanning fibers in the case of the CST. This shows that high resolution DWIs exhibits more anatomical information thanks to the smaller voxel size which might not be easily discernible at a lower spatial resolution [START_REF] Sotiropoulos | Advances in diffusion MRI acquisition and processing in the Human Connectome Project[END_REF]. Acquiring at higher spatial resolution could also help resolve complicated fiber configurations such as crossings fibers from fanning fibers or disentangle small structures like the optic chiasm [START_REF] Roebroeck | High-resolution diffusion tensor imaging and tractography of the human optic chiasm at 9.4 T[END_REF], which is not possible at lower spatial resolution [START_REF] Calabrese | Investigating the tradeoffs between spatial resolution and diffusion sampling for brain mapping with diffusion tractography: Time well spent? Human Brain Mapping[END_REF][START_REF] Jones | White Matter Integrity, Fiber Count, and Other Fallacies: The Do's and Don'ts of Diffusion MRI[END_REF].

Other methods for high spatial resolution acquisitions

We have shown in Fig. 6 that high spatial resolution acquisitions which are noisy at first can reveal improved anatomical details when they are subsequently denoised. This indeed suggests that high resolution acquisition can now be acquired on clinical scanners. Recently, other algorithms enabling a high spatial resolution at the acquisition level have been suggested [START_REF] Ning | A joint compressed-sensing and super-resolution approach for very high-resolution diffusion imaging[END_REF][START_REF] Scherrer | Accelerated High Spatial Resolution Diffusion-Weighted Imaging[END_REF]. These methods both rely on smartly fusing the (complementary) data of multiple acquisitions acquired at a lower spatial resolution to obtain a single high resolution volume. While the approach we suggest here is to acquire a single volume using a standard sequence, both techniques are fundamentally exploiting different aspects to increase the available spatial resolution. As such, it would be possible to combine our denoising algorithm with the reconstruction algorithms presented in [START_REF] Ning | A joint compressed-sensing and super-resolution approach for very high-resolution diffusion imaging[END_REF][START_REF] Scherrer | Accelerated High Spatial Resolution Diffusion-Weighted Imaging[END_REF].

Current limitations and possible improvements

Although most models assume a Rician or nc-noise distribution, this does not take into account the noise correlation between each coils or the effect of acceleration techniques that subsample the k-space [START_REF] Aja-Fernández | Noise estimation in parallel MRI: GRAPPA and SENSE[END_REF]. The development of correction factors for existing algorithms relies on computing the effective values for the noise standard deviation and the number of degrees of freedom of the ncdistribution, which is expected to be smaller than 2N. These values can be used to take into account the correlation introduced between the coils in parallel imaging acquisitions [START_REF] Brion | Noise correction for HARDI and HYDI data obtained with multi-channel coils and sum of squares reconstruction: an anisotropic extension of the LMMSE[END_REF]. To consider the fact that the noise distribution nature might vary spatially in addition to the noise variance, one can use a priori information obtained from the scanner through the SENSE sensitivity maps or the GRAPPA weights and need to estimate the correlation between each of the receiver coils. We could explicitly add such a correction to our algorithm since we work locally with the stabilization algorithm, Eq. (1) and Eq. (4).

Using multiband acceleration also modifies the noise properties due to the introduced aliasing, which further strengthen the idea that spatially adaptive denoising algorithms should be used on modern scanners and sequences [START_REF] Ugurbil | Pushing spatial and temporal resolution for functional and diffusion MRI in the Human Connectome Project[END_REF]. Nevertheless, obtaining the needed map for a SENSE reconstruction or the required GRAPPA weights might not be easily feasible in a clinical setting. We also intend to revisit the order in which preprocessing algorithms (motion correction, eddy currents correction, distortions correction) should be applied since these steps require interpolation, which could also introduce spatial correlation in the noise profile. This also makes the noise distribution deviate slightly from its theoretical distribution, where parameters vary spatially instead of being fixed constants for the whole volume [START_REF] Aja-Fernández | Noise estimation in parallel MRI: GRAPPA and SENSE[END_REF]. Nevertheless, we have observed experimentally that our NLSAM algorithm is robust to small subject motion thanks to the local neighborhood processing. In case where artifacts (such as epi distortions) might undermine the denoising process, one can first correct for these artifacts using a nearest neighbor interpolation, which should not modify the noise distribution. Subsequent corrections can then be performed after denoising using other kinds of interpolation as needed.

While developing the NLSAM algorithm, we found that using a bigger 3D patchsize did not significantly improve the denoising quality, while augmenting both computing time and memory requirements. Our implementation also allows one create the smallest subset of angular neighbors covering all DWIs through a greedy set cover algorithm. This option (named "NLSAM fast" in Table 4) leads to a speedup of 3 to 4

times, but at the cost of slightly reducing the denoising performance since some DWIs might be denoised only once instead of multiple times. We used the fully covered version for our experiments, which were run on a machine running Ubuntu Linux 12.04 with a quad core Intel i7 930 at 2.8 GHz and 18 GB of RAM.

Conclusion

In this paper, we introduced a new denoising method, the Non Local Spatial and Angular Matching (NLSAM), which is specifically designed to take advantage of diffusion MRI data. Our method is based on 1) Correcting the spatially varying Rician and nc-noise bias 2) Finding similar DWIs through angular neighbors to promote sparsity 3) Iteratively denoise similar patches and their neighbors locally with dictionary learning, where the local variance is used as an upper bound on the `2 reconstruction error. We extensively compared quantitatively our new method against three other state-of-the-art denoising methods on a synthetic phantom and qualitatively on an in-vivo high resolution dataset. We also showed that taking into account both the effect of spatially varying noise and non-Gaussian distributed noise is crucial in order to denoise effectively the DWIs. Our NLSAM algorithm is freely available 12 , restores perceptual information, removes the noise bias in common diffusion metrics and produces more anatomically plausible tractography on a high spatial resolution in-vivo dataset when compared to a lower spatial resolution acquisition of the same subject.

Since our NLSAM algorithm can be used on any already acquired dataset and does not add any acquisition time, this shows that denoising the data should be a pre-processing part of every pipeline, just like any other correction method that is commonly applied for artifacts removal. With that in mind, the diffusion MRI community could aim for higher spatial resolution DWIs, without requiring the use of costly new hardware or complicated acquisition schemes. This could in turn reveal new anatomical details, which

are not achievable at the spatial resolution currently used in diffusion MRI.

Figure 1 :

 1 Figure 1: a) A 3D block is made by stacking along the 4th dimension the b0, a DWI and its angular neighbors, which share similar structure, but under a slightly different noise realization. b) Disposition of equidistant angular neighbors on the sphere.

Figure 2 :

 2 Figure 2: Top : a) A noisy acquisition with slowly varying nc-noise and b) the resulting stabilized, Gaussian distributed noisy DWI. c) A noisy acquisition with fast varying Rician noise where the background was masked by the scanner with d) its stabilized counterpart. Bottom : Histogram of the nc-noise distribution in the selected background region of a) before stabilization and b) after stabilization. Note the non-Gaussianity of the noise in a) versus b).

  (a) pdf of nc-distribution. (b) cdf of nc-distribution.(c) icdf of Gaussian distribution.

Figure 3 :

 3 Figure 3: A synthetic example of the stabilization algorithm. a) Given a noisy value m = 678 observed in a nc-distribution with N = 4 and G = 200, the underlying value is estimated as ⌘ = 407. b) The associated probability in the nc-cdf with ⌘ is ↵ = 0.513, c) thus giving from the inverse cdf of a Gaussian distribution with mean µ = 407 and standard deviation G = 200 a new noisy value m = 413.

  Figure 4: Synthetic b1000 datasets with various noise profiles used in the experiments. The top row shows the b0 image, the middle row shows the same DWI across noise types and the bottom row shows the various noise distribution which generated the middle row. From left to right : the noiseless data, SNR 10 with stationary Rician noise, SNR 10 with stationary ncnoise, SNR 15 with spatially variable Rician noise, SNR 15 with spatially varying nc-noise.

Fig. 5 (

 5 Fig.5(displayed in landscape to see all the compared methods on the same figure) shows the b1000 noiseless data, the noisy input data at SNR 10 for nc-stationary noise and the results of the denoising on the synthetic phantom for all compared methods. This is the noise case theoretically covered by msPOAS and our NLSAM algorithm. We also show two zoomed regions of crossings with the reconstructed peaks extracted from fODFs. All perceptual and FA metrics were computed on the slice shown while angular metrics were computed in the zoomed region depicted by the yellow box. Note how the small blue bundle and its crossings are preserved on the NLSAM denoised dataset, while other denoising methods tend to introduce blurring. Fig.6(also in landscape for enhanced viewing) shows the noisy high resolution in-vivo dataset, the denoised version obtained for each algorithm and the low spatial resolution acquisition of the same subject without any denoising. Since our scanner uses a 32 channels head-coil but implements the SENSE reconstruction algorithm, the resulting spatially varying Rician noise distribution is the case covered by AONLM, LPCA and our NLSAM algorithm. We show a coronal slice for the gradient direction closest to (0, 1, 0), the colored FA map and a zoom on two regions of crossings. The yellow region shows the junction of the corticospinal tract (CST) and superior longitudinal fasciculus (SLF) while the white region shows the junction of the corpus callosum (CC) and the CST. While the high resolution dataset is noisier than its lower resolution counterpart, the highlighted crossings regions are well recovered by the denoising algorithms and thus offer an improvement in anatomical details over the lower spatial resolution dataset. We also see in the yellow box that the NLSAM denoised dataset recovers crossings extending from the CC which are almost absent in the compared datasets.

Figure 5 :

 5 Figure5: Phantomas b1000 synthetic dataset at SNR 10 for stationary nc-noise on the y = 24 slice. From top to bottom : Raw diffusion MRI, colored FA map, zoom on extracted peaks from fODF of order 8. Note how NLSAM restores the structure without blurring on the colored FA map and is the only method to restore the peaks from the noisy dataset in the zoomed white box region.

Figure 6 :

 6 Figure6: From top to bottom, the raw high resolution in-vivo data corrupted with spatially varying Rician distributed intensities, the colored FA map and a zoom on two regions of crossings. All denoised methods were applied on the high spatial resolution 1.2 mm dataset. We also show an acquisition of the same subject at 1.8 mm for visual comparison. Our NLSAM algorithm is able to recover more crossings from the 3 way junction of the SLF, the CST and the CC as shown in the yellow and white boxes. While the 1.8 mm dataset is less noisy, its lower spatial resolution also means that each voxel contains more heterogeneous tissues and mixed diffusion orientations. The 1.2 mm denoised dataset shows more crossings without the averaging effect of the larger voxel size. For a comparable acquisition time, the denoised high resolution dataset has more information than its lower resolution counterpart without processing.

Figure 7 :

 7 Figure 7: PSNR and SSIM metrics for the SNR 10 stationary and SNR 15 spatially variable noise cases datasets. All methods can correct the stationary and spatially varying Rician noise case to some extent while only our NLSAM algorithm has the best performance for the nc-noise case, especially for the spatially varying noise case.

Figure 8 :

 8 Figure8: Boxplot of the difference in FA for the synthetic datasets at b1000 (left) and b3000 (right). The whiskers show 1.5 times the interquartile range (1.5 x IQR), where outliers are plotted individually. The bars represent the first quantile, the median and the third quartile. No method performs well on the nc-b3000 spatially varying noise case, which is the hardest test case. NLSAM overall produces less error or is equal to the other methods, but has a lower bias in the FA error along noise type.

Figure 9 :

 9 Figure 9: The mean relative percentage of DNC error for the synthetic datasets. The bar represents the 95% confidence interval on the mean as computed by bootstrapping. The DNC error is the number of peaks found in excess or missing in each voxel with respect to the noiseless dataset.

Figure 10 :

 10 Figure10: Deterministic tractography for selected bundles on the in-vivo dataset. We also show a T1 weighted image aligned in the diffusion space for anatomical reference. Top : The left arcuate fasciculus. Note how the denoised NLSAM arcuate fasciculus goes further into the frontal and temporal region than both of its noisy 1.2 mm and 1.8 mm counterparts. Middle : The inferior fronto-occipital fasciculus. The AONLM denoised bundle has a denser part for the right IFOF while the LPCA bundle stops prematurely for the left IFOF, possibly due to a missing crossing along the bundle. Bottom : The corticospinal tract. We see that NLSAM recovers the commissural fibers in the pons from the noisy 1.2 mm dataset, which are not even present on the noisy 1.8 mm dataset nor on the other denoising algorithm's bundles. NLSAM also recovers more fanning to both sides of the brain than all the compared methods.

Figure 11 :

 11 Figure 11: Effect of the stabilization algorithm on the compared methods. The top row shows an axial slice of the in-vivo FA map computed on the stabilized dataset, where some voxels are degenerated. The bottom row shows their location on a binary brain mask. As shown in Table3, all methods produce degenerated FA voxels on both the regular and stabilized data, with the sole exception of NLSAM.

Table 2 :

 2 Tractometer results for the deterministic tracking.

Table 3 ,

 3 all methods produce degenerated FA voxels on both the regular and stabilized data, with the sole exception of NLSAM.

			AONLM	LPCA	msPOAS NLSAM	Noisy	Mask
	Built-in	Brain mask 83 314 (10.1%) 10 526 (1.3%) 84 319 (10.2%) WM mask 29 664 (5.1%) 1 298 (0.2%) 16 665 (2.9%)	↵ 5 994 (0.7%) 823 068 (100%) ↵ 1 769 (0.3%) 578 418 (100%)
	Stabilization	Brain mask WM mask	10 052 (1.2%) 15 750 (1.9%) 404 (0.1%) 1 468 (0.3%)	29 377 (3.6%) 4 267 (0.7%)	0 (0%) 9 395 (3.6%) 823 068 (100%) 0 (0%) 2 850 (0.5%) 578 418 (100%)

Table 3 :

 3 Number of degenerated FA voxels inside a brain mask and a white matter mask for the in-vivo dataset. All methods were compared with their built-in noise estimation on the stabilized version, but without any additional noise correction factor.

Table 4 :

 4 Table4reports the runtime of the various algorithms in minutes and their RAM usage. While the computing time required by NLSAM is larger than the other methods, our Python implementation is fairly unoptimized and could be sped up to competitive runtimes by various code optimizations or lowering the maximum number of iterations in Eq. (4). Required time and RAM usage for the compared denoising algorithms on the b1000 SNR 10 dataset with stationary Rician noise.

		AONLM LPCA msPOAS NLSAM NLSAM fast
	Time (mins)	22.2	3.7	4.0	37.1	9.8
	RAM usage (MB)	552	640	1543	606	412

http://hardi.epfl.ch/static/events/2013_ISBI/

http://www.emmanuelcaruyer.com/phantomas.php

http://personales.upv.es/jmanjon/denoising/index.htm

http://cran.r-project.org/web/packages/dti/index.html

https://github.com/samuelstjean/nlsam

http://www.tractometer.org/ismrm_2015_challenge/

Acknowledgments

We would like to thank Emmanuel Caruyer for adding our requested features and general help with the phantomas library. We also thank the Tractometer team (Jean-Christophe Houde and Marc-Alexandre Côté, www.tractometer.org) for the help in using the tractography evaluation system. We also extend our thanks to José V. Manjón for making the LPCA publicly available and Karsten Tabelow for the publicly available msPOAS code, helping us in the usage and in the choosing of good parameters for the in-vivo dataset denoising. Another special thanks to Guillaume Gilbert from Philips Healthcare for the help with the in-vivo datasets acquisition. Samuel St-Jean was funded by the NSERC CREATE Program in Medical Image Analysis. Pierrick Coupé was supported by the French State, managed by the French National Research Agency (ANR) in the frame of the Investments for the future Programme IdEx Bordeaux (HL-MRI ANR-10-IDEX-03-02), Cluster of excellence CPU, TRAIL (HR-DTI ANR-10-LABX-57) and the CNRS multidisciplinary project Défi ImagIn.

Appendix A. The NLSAM Algorithm

This appendix outlines the NLSAM algorithm as pseudo code. The original implementation in Python is freely available at https://github.com/samuelstjean/nlsam. foreach DWI in dMRI data do Step 2; Find the closest an angular neighbors; Create 4D block with b0, DWI and its an neighbors; Extract all overlapping patches of size (ps, ps, ps);

Step 3; Apply Eq. ( 1) to find D; Iterate Eq. ( 4) to find ↵ until convergence or max_iter is reached; Average overlapping patches based on sparsity with Eq. ( 5); Return to original shape;