Hong Hsi Lee^{1,2}, Els Fieremans^{1,2}, and Dmitry S Novikov^{1,2}

Conventional denoising and noise level estimation typically require data redundancy from multiple measurements or prior assumptions, such as smooth image prior or similarity between image patches. Here, we propose a single image denoising algorithm with noise map estimation by identifying the noise-only principle components based on universal properties of random covariance matrices, with the data redundancy created by segmenting data in the Fourier or wavelet domains. The proposed method is applicable to medical and other imaging modalities with spatially-varying noise, and is particularly beneficial to quantitative MRI acquisitions with a limited number of scans.

Denoising is essential for reliable image processing (e.g.,co-registration, super-resolution,etc.)[1]. In particular, estimating the image noise level is necessary for image reconstruction in compressed sensing[2]. To date, denoising and noise level estimation requires either data redundancy from multiple measurements (empirically>30)[3,4] or prior assumptions, such as smooth image prior[1] or similarity between image patches[5,6]. However, neither of these two conditions can be satisfied for some applications.

Here, to achieve single image denosing and noise level/map estimation without prior assumptions, we propose an algorithm to identify the noise-only principal components using universal properties of the random covariance matrix[3,4], for which the required data redundancy is created by stacking re-ordered Fourier transformed (k-space) or wavelet-transformed patches from the single image. This empirical redundancy serves as a ground for the denoising at a single-image level, enabling the separation between image features and the noise. This method is universally applicable to images with spatially-varying noise, not just for medical (e.g.,MRI), but also for other imaging modalities (e.g.,electron microscopy).

$$$\rm\bf{Denoising\,algorithm,\,Figs.1\text{-}2}$$$

Consider a $$$d$$$-dimensional image ($$$n$$$ voxels in total) corrupted by noise (e.g., Gaussian, Poisson, Rician). First, we perform the Fourier transform (FT) of the raw image, and segment the k-space data into $$$m_1$$$,$$$m_2$$$, ...,$$$m_d$$$ equal parts along each dimension respectively. Then we stack the segmented k-space data to create the data redundancy (Fig.2) by a factor of $$$m=\prod_{i=1}^dm_i$$$, and re-shape the stacked data into an $$$\frac{n}{m}$$$×$$$m$$$ redundant matrix $$$\bf{X}$$$ with P (presumably $$$\ll$$$min$$$(\frac{n}{m},m)$$$) linearly independent components. After applying the principal component analysis (PCA) to the matrix $$$\bf{X}$$$, the smallest min$$$(\frac{n}{m},m)$$$-P nonzero eigenvalues are purely contributed by the noise and obey the Marchenko-Pastur (MP) distribution[7]. Finally, we remove the noise components identified by MP distribution in PCA domain, unstack the data, and perform inverse-FT (iFT) to obtain the denoised image. As opposed to non-local means methods[5,6], we do not need ad hoc similarity criterion; rather, MP-PCA sets such a criterion automatically and objectively.

In actual implementations, FT can be replaced by other
orthonormal transformations, such as the wavelet transform (WT). The
FT-approach can estimate only the general noise level for the whole image. In
contrast, the WT-approach, combined with a sliding window (kernel), can
estimate a spatially dependent noise map with lower resolution.

$$$\rm\bf{Residual\,analysis}$$$

To test the performance of our denoising method, we analyzed the normalized residual $$$r$$$ of the raw image $$$I$$$ and the denoised image $$$I^\prime$$$[3-4], defined by

$$r\equiv\frac{I^\prime-I}{\sigma}\quad\quad(1)$$

with $$$\sigma$$$ the estimated noise level/map. The normalized residual should be normal distributed with a zero mean and a variance ≤1 when the denoising algorithm removes only the noise, without corrupting the signal.

$$$\rm\bf{Data\,1:\,T_1\text{-}weighted\,3D\,FGATIR\,MRI\,(single\,image\,in\,3d),\,Fig.3}$$$

We scanned a 20s-year-old female on a 3T-Siemens-Prisma using the T$$$_1$$$-weighted 3D-FGATIR sequence (Fast-Gray-Matter-Acquisition-T$$$_1$$$-Inversion-Recovery, resolution=(0.8mm)$$$^3$$$, FOV=179×230mm$$$^2$$$, TE/TR=4.39/3000ms, TI=409ms, acquisition time=11min), and denoised the brain image using level-2 WT-approach in 3d with a 5×5×5-kernel.

$$$\rm\bf{Data\,2:\,Scanning\,electron\,microscopy\,(SEM)\,(single\,image\,in\,2d),\,Fig.4}$$$

The brain tissue from a female 8-week-old C57BL/6 mouse’s genu of corpus callosum was processed and analyzed with an SEM (Zeiss Gemini-300) in a 6×6nm$$$^2$$$ pixel size[8]. We demonstrated our denoising method on a 400×400 data matrix, using level-3 WT-approach in 2d with a 7×7-kernel.

$$$\rm\bf{Data\,3:\,Diffusion\,MRI\,(4\,images\,in\,3d),\,Fig.5}$$$

We scanned a 24-year-old female on a 3T-Siemens-Prisma and acquired 1b=0 brain image and 3DWIs with b-value=1000s/mm$$$^2$$$ (resolution=(1mm)$$$^3$$$, FOV=182×184mm$$$^2$$$, TE/TR=68/5200ms, PF=5/8), and denoised the 4 images using level-2 WT-approach in 3d with a 5×5×5-kernel. Further, we acquired 3b=0 images and 30DWIs using the same sequence, and denoised the 33 images using the ordinary MP-PCA (dwidenoise, mrtrix, 5×5×5-kernel)[3,4] as the ground truth. Total acquisition time is 10min. We compared the mean diffusivity (MD) map for raw data, denoised data with proposed method, and the ground truth.

In Fig.3, compared with the raw image acquired by FGATIR, the denoised image shows better contrasts for tracts (dark), e.g., corticospinal tract passing through the pons. The noise map is homogenous over the central region. The normalized residual map shows no structure, and its histogram is close to and below the normal distribution (zero-mean, variance=1), indicating that our method removes most of the noise without corrupting the signal. The SEM image in Fig.4 shows similar results.

In Fig.5, compared with the MD map calculated from the raw data, the MD map calculated from the 4 denoised images using the proposed method is more homogenous and closer to the ground truth. The MD histogram shows that the proposed method increases the precision of the estimated parameter, indicated by the narrower MD distribution than the case of the raw data. The averaged normalized residual is homogenous with no anatomical structure, and its histogram is close to the normal distribution.

1. Liu, Ce, et al. "Noise estimation from a single image." Computer Vision and Pattern Recognition, 2006 IEEE Computer Society Conference on. Vol. 1. IEEE, 2006.

2. Lustig, Michael, et al. "Compressed sensing MRI." IEEE signal processing magazine 25.2 (2008): 72-82.

3. Veraart, Jelle, Els Fieremans, and Dmitry S. Novikov. "Diffusion MRI noise mapping using random matrix theory." Magnetic resonance in medicine 76.5 (2016): 1582-1593.

4. Veraart, Jelle, et al. "Denoising of diffusion MRI using random matrix theory." NeuroImage 142 (2016): 394-406.

5. Buades, Antoni, Bartomeu Coll, and J-M. Morel. "A non-local algorithm for image denoising." Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on. Vol. 2. IEEE, 2005.

6. Manjón, José V., et al. "Adaptive non‐local means denoising of MR images with spatially varying noise levels." Journal of Magnetic Resonance Imaging 31.1 (2010): 192-203.

7. Marčenko, Vladimir A., and Leonid Andreevich Pastur. "Distribution of eigenvalues for some sets of random matrices." Mathematics of the USSR-Sbornik 1.4 (1967): 457.

8. Lee, Hong-Hsi, et al. "Electron microscopy 3-dimensional segmentation and quantification of axonal dispersion and diameter distribution in mouse brain corpus callosum." bioRxiv(2018): 357491.

Figure 1. Single image denosing algorithm. To identify noise components using Marchenko-Pastur (MP) distribution in PCA domain, we segmented and stacked the k-space data to create the required data redundancy (Fig.2). The Fourier transform (FT-approach) and its inverse can be replaced by the wavelet transform (WT-approach) or other orthonormal transformations in actual implementations.

Figure 2. The origin of the created data redundancy. If we applied inverse-FT (iFT) to each k-space data patch, we obtained many small images containing similar structures in the same noise level. By stacking them together, we can apply MP-PCA to denoise the image. The iFT in this figure is unnecessary since iFT is a unitary transformation and has no effect on the PCA afterward.

Figure 3. Data 1: T$$$_1$$$-weighted 3D FGATIR MRI (single image denoising). Upper row shows the raw image, the denoised image (level-2 WT-approach in 3d with a 5×5×5-kernel), the estimated noise map, and the normalized residual map defined in Eq.(1). The denoised image shows better contrasts for the tract (red arrows) passing through the pons. Lower row shows the zoomed in images and the histogram/distribution of the normalized residual (red asterisk), which is close to and below the normal distribution with zero mean and variance=1 (black solid line), indicating that our method removes most of the noise without corrupting the signal.

Figure 4. Data 2: SEM (single image denoising). Upper row shows the raw data (field-of-view=2.4×2.4µm), the denoised data (level-3 WT-approach in 2d with a 7×7-kernel), the noise map, and the normalized residual map defined in Eq.(1). The denoised image shows clearer contours for cell organelles (red arrows in zoomed in image of the lower row). Lower row shows the zoomed in images (field-of-view=0.8×0.8µm) and the histogram/distribution of the normalized residual (red asterisk), which is close to and below the normal distribution with zero mean and variance=1 (black solid line), indicating that our method removes most of the noise, without corrupting the signal.

Figure 5. Data 3: Diffusion MRI (4 image denoising). Upper row shows the MD map (µm$$$^2$$$/ms) of 4 raw images (1b=0+3DWIs), the MD map of 4 denoised images (level-2 WT-approach in 3d with a 5×5×5-kernel), and the MD map of 33 images (3b=0+30DWIs) denoised by using the ordinary MP-PCA (dwidenoise, mrtrix) with a 5×5×5-kernel and serving as the ground truth. Lower row shows the MD histogram before and after denoising, the normalized residual averaged over 4 images, and the histogram of the normalized residual (red asterisk), which is close to the normal distribution with zero mean and variance=1(black solid line).