To denoise diffusion MRI (dMRI) data in a model-independent way. The thermal noise in dMRI signals propagates through the estimation process to the diffusion model parameters of interest. Therefore, reducing the noise variance in the dMRI signals increases the precision of the diffusion parameter estimators. We here adopt the idea of noise removal by means of transforming redundant data into the Principal Component Analysis (PCA) domain and preserving only the components that contribute to the signal¹. Unfortunately, an objective criterion to discriminate between the “signal-carrying” and “noise-only” components has been lacking. Indeed, the number of signal-carrying principal components is unknown and is expected to depend on factors such as resolution and SNR. Commonly used criteria include thresholding of the eigenvalues associated to the principal components by an empirically set threshold value². Here, we objectify the threshold for PCA denoising by exploiting the fact that the noise-only eigenvalues are expected to obey the universal Marchenko-Pastur (MP) distribution³.

Oversampled dMRI data is represented by a number of parameters (P) that is significantly smaller than the number of diffusion encodings (N). Therefore, a M×N data matrix covering $$$M\gg 1$$$ voxels within a sliding window (e.g. $$$[5\times 5\times 5]$$$) is highly redundant. The PCA transformation of redundant data compresses the noise-free signal variance into P components, whereas the noise is spread over all components. The bulk of the PCA eigenspectrum of the data matrix is expected to show the spectral behavior of noise. According to Random Matrix Theory, the MP distribution describes the universal noise eigenspectrum³ (Fig.1):

$$ p(\lambda)=\frac{\sqrt{(\lambda_+ -\lambda)(\lambda -\lambda_-)}}{2\pi \gamma \sigma^2} \quad\mathrm{if}\quad\lambda_-\leq \lambda\leq\lambda_+\quad (p(\lambda)=0\:\mathrm{ otherwise)}$$

with $$$\lambda_\pm=\sigma^2(1\pm\sqrt{\gamma})^2$$$, $$$\sigma$$$ the noise level and $$$\gamma=M/(N-P)$$$. Distribution fitting, i.e. minimizing the error between the MP distribution $$$p(\lambda)$$$ and the histogram $$$n(\lambda)$$$ of the lowest eigenvalues, yields an accurate estimate of $$$\sigma$$$, P, and an objective threshold to idenitify significant signal components, i.e. $$$\lambda_+$$$ . We use the P significant eigenvalues $$$(\lambda >\lambda_+)$$$ to estimate the expected value of the signal in each voxel and diffusion gradient direction. Since the noise level and signal are estimated simultaneously, both parameters can be corrected for the nc-$$$\chi$$$ noise bias afterwards⁴. We compare our technique (MPPCA) with adaptive non-local mean⁵ (ANLM) and total generalized variation⁶ (TGV).

Simulation: A M×N matrix was generated by projecting M, ranging from 100 to 2500, axially symmetric diffusion tensors onto N=30, 60, and 90 diffusion gradient directions with $$$b=1\,\mathrm{ms/\mu m^2}$$$. Monte Carlo simulations showed that the accuracy as function of M and N the method was >99.3%, whereas the noise standard deviation (SD) was reduced to $$$σ\sqrt{P/N+P/M}$$$.

Real data ($$$5\times b=0$$$, and 3 repetitions of $$$ 90\times b=\{1,\, 2.5\}\,\mathrm{ms/\mu m^2}$$$): We applied MPPCA, ANLM, and TGV to a subset of the $$$b=2.5\,\mathrm{ms/\mu m^2}$$$-shells (1 repetition, N=60). Fig. 2 shows the denoised images for qualitative comparison. The blurring for ANLM and TV is quantitatively confirmed by comparing the k-energy density as function of the distance to the k-space center (Fig. 3). Indeed, the suppression of high frequencies explains the low-pass filtering effect. The suppression of high frequencies not strictly reflecting thermal noise, but rather anatomical details, also results in residuals, normalized by σ, that follow a distribution with SD exceeding one (meaning that they are contaminated by signal), whereas MPPCA shows a distribution that matches the simulations, i.e. zero-centered normal distribution with SD=$$$\sigma\sqrt{1-P/N-P/M}$$$(Fig. 4). Hence, MPPCA preserves the actual signal better than competing methods. The same trends are observed in the fractional anisotropy (FA) maps derived from the $$$ b=\{0,\,1\} \,ms/\mu m^2$$$ subsets (N_b>0=30, 60, and 90). We observe a significant increase in precision in the estimation of FA, determined using bootstrapping with replacement based on the 3 repetitions (Fig. 5). Although ANLM might outperform MPPCA in terms of precision, ANLM is shown to be less accurate when comparing the denoised results to a ground truth, derived from $$$3\times (5\times b=0+ 90 \times b=\{1,\, 2.5\} \,\mathrm{ms/\mu m^2})$$$.

Unlike competing methods, MPCCA is a selective denoising technique that reduces signal fluctuations solely rooting in thermal noise. The selective nature of the proposed technique is created by exploiting data redundancy in the PCA domain using universal properties of the eigenspectrum of random covariance matrices. Indeed, Typical dMRI data has shown to exhibit sufficient redundancy ($$$P\ll M,N$$$), which allows to distinguish between the noise and signal components in a model-independent way, and without the need to empirically set the PCA threshold value. The result show improved precision in the estimation of diffusion model parameters whereas the accuracy is preserved. Similar results have been obtained for other all diffusion parameters.