# Image Denoising Using Nonlocal Means with Shape-Adaptive Patches and New Weights

Authors: Chenglin Zuo; Jun Ma; Hao Xiong; Lin Ran
Journal: Shock and Vibration
 (2021)

---

## Abstract

Digital images captured from CMOS/CCD image sensors are prone to noise due to inherent electronic fluctuations and low photon count. To efficiently reduce the noise in the image, a novel image denoising strategy is proposed, which exploits both nonlocal self-similarity and local shape adaptation. With wavelet thresholding, the residual image in method noise, derived from the initial estimate using nonlocal means (NLM), is exploited further. By incorporating the role of both the initial estimate and the residual image, spatially adaptive patch shapes are defined, and new weights are calculated, which thus results in better denoising performance for NLM. Experimental results demonstrate that our proposed method significantly outperforms original NLM and achieves competitive denoising performance compared with state-of-the-art denoising methods.


---

## Body

## 1. Introduction

Digital imaging devices such as digital cameras and camera phones are ubiquitous in our daily life, which use complementary metal oxide semiconductors (CMOS) or charged coupled devices (CCD) image sensors to acquire images. However, since the CMOS and CCD image sensors are subject to noise from two notable sources, i.e., electronic instruments and the photo-sensing devices [1, 2], the quality of captured images is usually not satisfactory, especially when images are taken in low light condition, which leads to degraded imaging results. Hence, denoising has become a fundamental image restoration problem in image signal processor (ISP).During the past decades, image denoising has been widely studied. However, until now, how to remove the noise efficiently while preserving significant image details has remained a challenge. Early smoothing methods, such as Gaussian filter [3], anisotropic filter [4], total variation [5], and bilateral filter [6], perform noise removal solely based on the information provided in a local neighborhood, thereby resulting in disturbing artifacts around edges. Later, transform domain-based denoising methods have been developed and extensively studied as well [7–19]. Wavelet transform (WT) [7] decomposes the image into multiple frequency components, where the noise is removed with thresholding [8, 9] or statistical modeling [10–12]. By transforming back the processed wavelet coefficients into spatial domain, denoising is accomplished. Late development of WT denoising includes ridgelet [13] and curvelet [14] methods. In [15], adaptive principal components-based denoising method was proposed. Compared with WT that uses a fixed wavelet basis, it computes the locally fitted basis to decompose the image. In [16, 17], the highly overcomplete dictionary was trained by using K-SVD algorithm for sparse and redundant image representation. In [18, 19], discrete cosine transform (DCT) was applied to the local neighborhood, which achieves very sparse representation of the image and hence leads to effective denoising performance.Recently, nonlocal methods, that exploit the image nonlocal self-similarity, have achieved outstanding denoising performance. In [20, 21], Buades et al. first proposed the nonlocal principle-based denoising method, called nonlocal means (NLM). In this method, noise-free pixel is estimated as a weighted average of all pixels in the image, where the weights are determined based on the similarity between the patch centered at the pixel being estimated and the patches centered at other pixels. Since NLM exploits the fact that similar patches appear abundantly in the image and can contribute for denoising, it obtains high quality denoising performance. Subsequently, numerous extensions of NLM have been developed. In [22], the optimal neighborhood for each pixel was chosen during the iteration procedure to balance the accuracy of approximation and the stochastic error. In [23], the noisy image was classified into several region types, according to which the patch size was adaptively adjusted to match the local property. In [24], adaptive patch size and bandwidth were selected pixel-wise, relying on the feature metric that can provide a quantitative measure of local geometric structures. In [25], the smoothing parameter was chosen automatically based on noise estimation, and then the two-stage NLM with adaptive smoothing parameter was performed. In [26], quadtree-based NLM was proposed, which employs quadtree decomposition on each image patch to obtain more subpatches of various sizes. In addition to considering the patch sizes, some methods try to handle variable patch shapes. In [27], the adaptive binary shape for each patch was estimated by thresholding the difference between the central pixel and other patch pixels. In [28], shape-adaptive patches were constructed to match more homogeneous pixels successfully, especially in textured areas. In [29], several types of patches with various shapes were predefined and applied in NLM, respectively. Then, local estimates associated with these shapes were combined using Stein’s unbiased risk estimate (SURE). To improve the accuracy of similarity measure, nonlocal similarity of residual image structures in method noise was further exploited in [30, 31]. Besides, rotation invariant patch comparison, that can handle rotational similarity existing in the image, was also studied in [32–36]. Analogously in [37], affine invariant similarity measure was applied to find more similar patches. Although NLM and its extensions have achieved significant denoising results, only exploiting the spatially nonlocal redundancy still limits their performance. Therefore, some methods combine the nonlocal principle with other techniques [38, 39], resulting in state-of-the-art denoising performance [40–43].In this paper, we address these issues and propose an efficient denoising method, as shown in Figure1. First, the original NLM is employed to obtain an initial estimate of the noisy image. However, due to inaccurate weight computation with noise interference, the initial estimate does not contain complete image details, which means the method noise still contains residual image information. For well preserving the residual image in method noise, the wavelet thresholding is used to smooth noise as much as possible. Then, the preserved residual image is combined with initial estimate to obtain a basic denoising result, based on which spatially adaptive patch shapes are defined using LPA-ICI and new weights are calculated. Finally, NLM denoising is implemented again but with the shape-adaptive patches and new weights.Figure 1
Flowchart of the proposed image denoising method.The remainder of this paper is structured as follows: in Section2, original NLM is briefly reviewed. In Section 3, our proposed denoising method is described in detail. In Section 4, we present and analyse the comparative experimental results. Finally, Section 5 concludes this paper.


## 2. Nonlocal Means

Given a noise-free imageu defined on a discrete grid I, the noisy observation of u at pixel i∈I is defined as(1)vi=ui+ni,where n(i) is the zero-mean white Gaussian noise perturbation at pixel i. Let Ni denote the patch centered at pixel i, and its noisy observation is defined as(2)vNi=vj|j∈Ni.In classical NLM [20], for a pixel i, the estimated value of its noise-free version, u^i, is calculated as the weighted average of all noisy observations of the pixels in the image:(3)u^i=∑j∈Iwi,jvj∑j∈Iwi,j.The weightwi,j depends on the similarity of patch Ni and patch Nj, which is defined as(4)wi,j=e−vNj−vNj22/λ2,where ⋅2 denotes the Euclidean norm to measure the similarity between the patches, and the parameter λ acts as a degree of filtering, therefore controlling the decay of the weight. Since the classical NLM only estimates a single pixel, it is referred as the pixel-wise NLM.Later, the patch-wise NLM does also exist [21]. Here, the whole patch is estimated as follows:(5)u^Ni=∑j∈Iwi,jvNi∑j∈Iwi,j.As the patches are overlapped, each pixel in a patch achieves multiple estimates. Therefore, for a pixeli, its different estimates are aggregated to obtain the final result:(6)u^i=1Ai∑j∈Aiu^Nj,i,where Ai=j|i∈Nj and u^Nj,i denote the estimated value of noise-free patch Nj at pixel i. Thanks to the multiple estimations for each pixel, the patch-wise NLM achieves better denoising performance than the pixel-wise one. In this paper, when NLM is mentioned, it refers to the patch-wise NLM.Based on the denoising result, method noise is defined as follows:(7)n^=v−u^.It can be inferred that, if NLM performs well, the method noise must look like a noise and should contain as little structure as possible. In Figures2(a) and 2(b), an example image and its noisy version (σ=20) are shown, respectively. Figure 2(c) shows the NLM denoising result, while Figure 2(d) shows the corresponding method noise. As can be seen, obvious image structure appears in the method noise, which means that some image details are removed from the denoised image. Therefore, we make use of the residual image in method noise to exploit nonlocal self-similarity further.Figure 2
NLM denoising performance. (a) Original image; (b) noisy image (σ = 20); (c) initial estimate; (d) method noise; (e) method noise thresholding with wavelet; (f) basic estimate combining initial estimate and filtered method noise.
(a)(b)(c)(d)(e)(f)


## 3. Proposed Denoising Method

### 3.1. Method Noise Thresholding

Since the method noise contains obvious image structure, the residual image in it is estimated firstly. Here, we use the BayesShrink wavelet thresholding method [44] to suppress the noise as much as possible.BayesShrink is an adaptive, data-driven thresholding strategy via soft-thresholding which derives the threshold in a Bayesian framework, assuming a generalized Gaussian distribution for the wavelet coefficients. This method is adaptive to each sub-band because it depends on data-driven estimates of the parameters. The threshold for a given sub-band is derived by minimizing Bayesian risk as follows:(8)T=σn2σw,where σn2 is the noise variance estimated from sub-band HH1 by a robust median estimator, given by(9)σ^n=MedianYi,j0.6745,Yi,j∈HH1,and σw2 is the variance of wavelet coefficients in that sub-band, whose estimate is computed using(10)σ^w2=maxσ^y2−σ^n2,0,where σ^y2=1/MN∑i,j=1M,NYi,j2.Figure2(e) shows the filtered method noise thresholding with wavelet. As can be seen, the residual image is preserved well while the noise is smoothed efficiently. Then, we combine the initial estimate, denoted as u^initial, and the filtered method noise, denoted as n^f, together to obtain a basic estimate:(11)u^basic=u^initial+n^f.In Figure2(f), we show the combined result. It can be seen that the basic estimate preserves more image details than the initial estimate, such as the regions marked by green boxes, which means that it will be more accurate to estimate the spatially adaptive patch shapes and to calculate the weights based on the basic estimate.


### 3.2. Spatially Adaptive Patch Shape Estimating

The anisotropic local polynomial approximation- (LPA-) intersection of confidence intervals (ICI) technique is used to estimate the spatially adaptive shape for each patch in the image.Figure3 shows the implementation of the LPA-ICI-based patch shape estimating. For a pixel i, eight directions are first predefined. For every specified direction θk=k−1/4π, k=1,…,8, a varying-scale family of narrow “linewise” directional LPA convolution kernels gh,θkh∈H is used to obtain a corresponding set of directional varying-scale estimates u^h,θkh∈H, u^h,θk=v⊗gh,θk, h∈H, where H∈R+ is the set of scales. Then, for each estimate, a confidence interval is built as follows:(12)Dh,θk=u^h,θk−Γσ,u^h,θk+Γσ,where Γ>0 is a tuning parameter and σ is the noise standard deviation. Based on ICI rule, an adaptive scale h+i,θk∈H is defined for every direction θk. Finally, the shape-adaptive patch Ni+ is constructed as the polygonal hull of suppgh+i,θk,θkk=18.Figure 3
Implementation of the LPA-ICI-based patch shape estimating. “Linewise” one-dimensional directional LPA kernels are used for 8 directions. The shape-adaptive patchNi+ is constructed as the polygonal hull of the adaptive-scale kernel supports.In Figure4, we show some examples of the estimated shape-adaptive patches in the noise-free, noisy, initially estimated, and basically estimated images, respectively. It can be seen that, due to the influence of strong noise, the estimated patch shapes in the noisy image are incorrect. The same goes for those in the initially estimated image but because of the loss of image details during NLM denoising. By contrast, patch shapes in the basically estimated image are more accurate.Figure 4
Some examples of the estimated shape-adaptive patches in (a) noise-free image, (b) noisy image (σ = 50), (c) initially estimated image by NLM, and (d) basically estimated image, respectively.
(a)(b)(c)(d)


### 3.3. New Weight Calculating

Based on the basic estimate, weights between the patches are calculated again. Since we estimate the spatially adaptive shape for each patch, the new weight is calculated as follows:(13)wi,j+=e−u^basicNi+−u^basicNji,+22/γ2,where Ni+ denotes the estimated shape-adaptive patch, Nji,+denotes the patch using the same shape with Ni+, and γ is the filtering factor.


### 3.4. NLM with Shape-Adaptive Patches and New Weights

With shape-adaptive patches and new weights, we implement the NLM denoising again to remove the noise in the noisy image. For a shape-adaptive patchNi+, the estimated value of its noise-free version, u^Ni+, is calculated as follows:(14)u^Nj+=∑j∈Iwi,j+vNji,+∑j∈Iwi,j+.Similarly, for a pixeli, its different estimates are aggregated to obtain the final result:(15)u^finali=1Ai+∑j∈Ai+u^Nj+,i,where Ai+=j|i∈Nj+ and u^Nj+,i denote the estimated value of noise-free patch Nj+ at pixel i.


### 3.5. Computational Complexity

For the image of sizeN×N, the computational complexity of the original NLM for initial estimate is ΟnpN2, where np denotes the patch size used by NLM. In our practical implementation, a limited window of size w1×w1 is used to restrict the search of similar patches, which reduces the complexity to Οnpw1N. Besides, by using the moving average filter together with weight symmetry, the complexity can be further brought down to Οw1N. Then, the computational complexity of method noise thresholding with BayesShrink wavelet is usually ΟN. For the spatially adaptive patch shape estimation with LPA-ICI technique, since it is based on convolutions against one-dimensional kernels for a very limited number of directions, its computational overhead for the whole noise removal processing is negligible. In the final estimating procedure, for the search window of size w2×w2, the computational complexity is Οw2N.


## 3.1. Method Noise Thresholding

Since the method noise contains obvious image structure, the residual image in it is estimated firstly. Here, we use the BayesShrink wavelet thresholding method [44] to suppress the noise as much as possible.BayesShrink is an adaptive, data-driven thresholding strategy via soft-thresholding which derives the threshold in a Bayesian framework, assuming a generalized Gaussian distribution for the wavelet coefficients. This method is adaptive to each sub-band because it depends on data-driven estimates of the parameters. The threshold for a given sub-band is derived by minimizing Bayesian risk as follows:(8)T=σn2σw,where σn2 is the noise variance estimated from sub-band HH1 by a robust median estimator, given by(9)σ^n=MedianYi,j0.6745,Yi,j∈HH1,and σw2 is the variance of wavelet coefficients in that sub-band, whose estimate is computed using(10)σ^w2=maxσ^y2−σ^n2,0,where σ^y2=1/MN∑i,j=1M,NYi,j2.Figure2(e) shows the filtered method noise thresholding with wavelet. As can be seen, the residual image is preserved well while the noise is smoothed efficiently. Then, we combine the initial estimate, denoted as u^initial, and the filtered method noise, denoted as n^f, together to obtain a basic estimate:(11)u^basic=u^initial+n^f.In Figure2(f), we show the combined result. It can be seen that the basic estimate preserves more image details than the initial estimate, such as the regions marked by green boxes, which means that it will be more accurate to estimate the spatially adaptive patch shapes and to calculate the weights based on the basic estimate.


## 3.2. Spatially Adaptive Patch Shape Estimating

The anisotropic local polynomial approximation- (LPA-) intersection of confidence intervals (ICI) technique is used to estimate the spatially adaptive shape for each patch in the image.Figure3 shows the implementation of the LPA-ICI-based patch shape estimating. For a pixel i, eight directions are first predefined. For every specified direction θk=k−1/4π, k=1,…,8, a varying-scale family of narrow “linewise” directional LPA convolution kernels gh,θkh∈H is used to obtain a corresponding set of directional varying-scale estimates u^h,θkh∈H, u^h,θk=v⊗gh,θk, h∈H, where H∈R+ is the set of scales. Then, for each estimate, a confidence interval is built as follows:(12)Dh,θk=u^h,θk−Γσ,u^h,θk+Γσ,where Γ>0 is a tuning parameter and σ is the noise standard deviation. Based on ICI rule, an adaptive scale h+i,θk∈H is defined for every direction θk. Finally, the shape-adaptive patch Ni+ is constructed as the polygonal hull of suppgh+i,θk,θkk=18.Figure 3
Implementation of the LPA-ICI-based patch shape estimating. “Linewise” one-dimensional directional LPA kernels are used for 8 directions. The shape-adaptive patchNi+ is constructed as the polygonal hull of the adaptive-scale kernel supports.In Figure4, we show some examples of the estimated shape-adaptive patches in the noise-free, noisy, initially estimated, and basically estimated images, respectively. It can be seen that, due to the influence of strong noise, the estimated patch shapes in the noisy image are incorrect. The same goes for those in the initially estimated image but because of the loss of image details during NLM denoising. By contrast, patch shapes in the basically estimated image are more accurate.Figure 4
Some examples of the estimated shape-adaptive patches in (a) noise-free image, (b) noisy image (σ = 50), (c) initially estimated image by NLM, and (d) basically estimated image, respectively.
(a)(b)(c)(d)


## 3.3. New Weight Calculating

Based on the basic estimate, weights between the patches are calculated again. Since we estimate the spatially adaptive shape for each patch, the new weight is calculated as follows:(13)wi,j+=e−u^basicNi+−u^basicNji,+22/γ2,where Ni+ denotes the estimated shape-adaptive patch, Nji,+denotes the patch using the same shape with Ni+, and γ is the filtering factor.


## 3.4. NLM with Shape-Adaptive Patches and New Weights

With shape-adaptive patches and new weights, we implement the NLM denoising again to remove the noise in the noisy image. For a shape-adaptive patchNi+, the estimated value of its noise-free version, u^Ni+, is calculated as follows:(14)u^Nj+=∑j∈Iwi,j+vNji,+∑j∈Iwi,j+.Similarly, for a pixeli, its different estimates are aggregated to obtain the final result:(15)u^finali=1Ai+∑j∈Ai+u^Nj+,i,where Ai+=j|i∈Nj+ and u^Nj+,i denote the estimated value of noise-free patch Nj+ at pixel i.


## 3.5. Computational Complexity

For the image of sizeN×N, the computational complexity of the original NLM for initial estimate is ΟnpN2, where np denotes the patch size used by NLM. In our practical implementation, a limited window of size w1×w1 is used to restrict the search of similar patches, which reduces the complexity to Οnpw1N. Besides, by using the moving average filter together with weight symmetry, the complexity can be further brought down to Οw1N. Then, the computational complexity of method noise thresholding with BayesShrink wavelet is usually ΟN. For the spatially adaptive patch shape estimation with LPA-ICI technique, since it is based on convolutions against one-dimensional kernels for a very limited number of directions, its computational overhead for the whole noise removal processing is negligible. In the final estimating procedure, for the search window of size w2×w2, the computational complexity is Οw2N.


## 4. Results and Discussion

In this section, we compare our proposed method with original NLM [21] and other two state-of-the-art denoising methods: shape-adaptive patches-based NLM (NLM-SAP) [29] and shape-adaptive PCA-based BM3D (BM3D-SAPCA) [41]. In the experiments, a set of 8 natural images commonly used in the literature of image denoising are used for the comparison, and their noisy versions are simulated by adding independent white Gaussian noise with varying noise levels. The results of our proposed method are generated using the scales H=1,2,3,4,5,6 and a search window of 21×21. The parameter γ is selected experimentally as γ=0.4σ. The results of other three methods are obtained by using the codes available online with recommended parameters. To evaluate the quality of denoised images, the popular peak signal to noise ratio (PSNR) and the structural similarity index (SSIM) are calculated.The results of the experiments are shown in Tables1 and 2, where the best results among the four methods are highlighted. It can be seen that our proposed method significantly outperforms NLM and NLM-SAP and achieves competitive denoising performance compared with BM3D-SAPCA. Particularly, in some cases of high noise levels, our method performs even slightly better than BM3D-SAPCA. In terms of SSIM results, our proposed method is quite close to BM3D-SAPCA, and its superiority exists in all cases with respect to NLM.Table 1
The PSNR results by different denoising methods. In each cell, the results of the four denoising methods are presented in the following order: top left, NLM [21]; top right, NLM-SAP [29]; bottom left, BM3D-SAPCA [41]; bottom right, our proposed denoising method.
σ203550100C. man29.8529.7527.0926.9825.3724.8221.4220.8630.9129.9728.1727.5426.5926.0222.8722.50House32.4032.6229.7929.0027.7426.3723.1622.7933.8933.2831.3730.9029.5228.9725.0824.77Peppers30.1730.6427.0827.5025.1625.2621.0221.0231.5730.8228.7528.2226.9826.5123.2422.90Lena31.5831.9728.9529.0027.3927.1223.9323.9933.2032.6230.7230.2829.0628.6725.3625.37Barbara30.4730.5027.6026.8325.7224.6922.1422.0232.3731.3629.6128.3427.6826.4023.2222.65Boats29.8029.6426.9926.8025.2825.0422.1922.3131.0230.3428.5127.9626.8826.4723.6823.42Man29.7529.5827.0326.9225.4125.3922.5422.9530.8330.3228.3828.0426.9326.5923.9623.98Hill29.7729.4527.1426.9225.4925.5522.8423.2830.8530.3228.6128.2127.1926.8624.2624.29Average30.4730.5127.7027.4925.9425.5322.4022.4031.8331.1229.2628.6827.6027.0623.9523.73The bold values represent the best results among the four methods, which has been explained in the second paragraph of Section4.Table 2
The SSIM results by different denoising methods. In each cell, the results of the four denoising methods are presented in the following order: top left, NLM [21]; top right, NLM-SAP [29]; bottom left, BM3D-SAPCA [41]; bottom right, our proposed denoising method.
σ203550100C. man0.8400.8450.7630.7730.7140.7060.5410.4930.8860.8540.8270.8010.7870.7580.6430.587House0.8310.8490.7750.7870.7330.7150.5530.5310.8760.8700.8380.8360.8070.7880.6760.639Peppers0.8400.8630.7590.7940.7030.7300.5180.5600.8860.8770.8330.8260.7920.7720.6660.631Lena0.8300.9180.7640.8580.7240.8010.5740.6440.8800.8790.8370.8190.8010.7720.6740.652Barbara0.8610.9250.7790.8430.7110.7680.5100.6080.9120.8940.8630.8150.8110.7370.6000.531Boats0.7860.8760.6890.7890.6250.7200.4710.5610.8280.8100.7610.7370.7080.6790.5780.531Man0.7930.8710.6920.7850.6250.7210.4730.5820.8400.8250.7630.7450.7100.6840.5790.555Hill0.7650.8510.6590.7570.5890.6940.4480.5720.8090.7900.7290.7090.6750.6480.5490.535Average0.8180.8740.7350.7980.6780.7310.5110.5680.8640.8490.8060.7860.7610.7290.6200.582The bold values represent the best results among the four methods, which has been explained in the second paragraph of Section4.Let us then focus on the visual quality of the denoised images by the four methods. In Figures5 and 6, we show the denoising results on two typical images with moderate and strong noise corruption, respectively. It can be seen that our proposed method is very effective in reconstructing both the smooth and the texture/edge regions. When the noise level is not very high, as shown in Figure 5, our proposed method performs better than NLM and NLM-SAP on edge preservation, and BM3D-SAPCA achieves the best visual output. When the noise level is high, as shown in Figure 6, however, details in the denoised images by NLM and NLM-SAP become blurred, and BM3D-SAPCA tends to generate many visual artifacts. By contrast, our proposed method performs much better, which preserves the image details well and generates much less artifacts.Figure 5
Denoising performance on the C. man image with moderate noise corruption. (a) Original image; (b) noisy image (σ = 20); denoised image (c) by NLM [21]; (d) NLM-SAP [29]; (e) BM3D-SAPCA [41]; (f) our proposed method.
(a)(b)(c)(d)(e)(f)Figure 6
Denoising performance on the Lena image with strong noise corruption. (a) Original image; (b) noisy image (σ = 100); denoised image (c) by NLM [21]; (d) NLM-SAP [29]; (e) BM3D-SAPCA [41]; (f) our proposed method.
(a)(b)(c)(d)(e)(f)


## 5. Conclusions

In this work, we have presented an efficient image denoising method. By exploiting the residual image in the method noise, spatially adaptive patch shapes are defined, and new weights are calculated to improve the denoising performance of NLM further. Experimental results demonstrate that our proposed method is effective in noise removal and texture/edge preservation and can achieve competitive denoising performance compared with state-of-the-art denoising methods.



---

*Source: 9532702-2021-07-27.xml*