DWI denoising

MRtrix includes the command dwidenoise, which implements dMRI noise level estimation and denoising based on random matrix theory. The method exploits data redundancy in the patch-level PCA domain ([Veraart2016a], [Veraart2016b] and [CorderoGrande2019]). The method uses the prior knowledge that the eigenspectrum of random covariance matrices is described by the universal Marchenko-Pastur (MP) distribution.

Advanced options

Patch size

The noise level in MRI is spatially varying, due to the proximity of the coil elements and parallel imaging. Noise level estimation and denoising therefore operates in image patches around each voxel, where the noise can be assumed to be approximately homoscedastic. The patch size can be chosen by the user with the option -extent. For maximal SNR gain (when using Exp2, see below) we suggest to choose \(N \approx M\), where \(M\) is the no. DW volumes and \(N\) is the number of kernel elements. However, larger kernels also extend the required run time, so in large datasets it might be beneficial to select smaller sliding kernels. By default, the command will select the smallest isotropic patch size that exceeds the number of DW images in the input data, e.g., 5x5x5 for data with <= 125 DWI volumes, 7x7x7 for data with <= 343 DWI volumes, etc.

Noise level estimation

The noise level in each patch is experimentally estimated from the eigenvalue spectrum of the local data matrix. Assuming \(M<N\), \(P\) signal-carying components (also estimated), and \(M-P\) noise components, the squared noise level is estimated as:

\[\sigma^2 = \frac{\lambda_{P+1}-\lambda_M}{4\sqrt{\gamma}}\]

where \(\lambda_i\) are the eigenvalues of the covaliance matrix, sorted in decreasing order, and \(\gamma\) is the matrix ratio.

dwidenoise implements two different versions of this estimator, based on a different definition of the matrix ratio \(\gamma\):

  • Exp1 uses the definition used in the original papers [Veraart2016a] and [Veraart2016b], namely \(\gamma = (M-P)/N\).

  • Exp2 uses \(\gamma = (M-P)/(N-P)\) instead, which was shown in [CorderoGrande2019] to improve the estimation. This is now the default.

Complex data

Note that dwidenoise does not correct for non-Gaussian noise biases present in magnitude-reconstructed MRI images. Including MRI phase images can reduce such Rician bias [CorderoGrande2019], and the command now supports complex input data. To this end, users can run:

mrcalc dwi_magnitude.mif dwi_phase_rad.mif -polar dwi_complex.mif
dwidenoise dwi_complex.mif out.mif -noise noise.mif
mrcalc out.mif -abs out_magnitude.mif

Where dwi_magnitude.mif and dwi_phase_rad.mif are the magnitude and phase images respectively. Note that the code above assumes that the phase is scaled in radians.

Whilst using complex data can effectively reduce Rician noise bias, dwidenoise currently does not account for spatial noise correlations introduced due to in-plane accelleration. [CorderoGrande2019] addresses these additional effects at reconstruction, but this is not implemented in MRtrix3.