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.
Image denoising must be performed as the first step of the image-processing pipeline. Interpolation or smoothing in other processing steps, such as motion and distortion correction, may alter the noise characteristics and thus violate the assumptions upon which MP-PCA is based.
Typical use will be:
dwidenoise dwi.mif out.mif -noise noise.mif
dwi.mif contains the raw input DWI image,
out.mif is the denoised
DWI output, and
noise.mif is the estimated spatially-varying noise level.
We always recommend eyeballing the residuals, i.e. out - in, as part of the quality control. The lack of anatomy in the residual maps is a marker of accuracy and signal-preservation during denoising. The residuals can be easily obtained with
mrcalc dwi.mif out.mif -subtract res.mif mrview res.mif
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
-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
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:
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.
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.