A MULTI-STAGE ALGORITHM FOR IMAGE DENOISING BASED ON PCA AND ADAPTIVE TV-REGULARIZATION

In this paper, we present an image denoising algorithm comprising three stages. In the ﬁrst stage, Principal Component Analysis (PCA) is used to suppress the noise. PCA is applied to image blocks to characterize localized features and rare image patches. In the second stage, we use the Gaussian curvature to develop an adaptive total-variation-based (TV) denoising model to effectively remove visual artifacts and noise residual gener-ated by the ﬁrst stage. Finally, the denoised image is sharpened in order to enhance the contrast of the denoising result. Experimental results on natural images and computed tomography (CT) images demonstrated that the proposed algorithm yields denoising results better than competing algorithms in terms of both qualitative and quantitative aspects.


Introduction
Imaging devices always produce some noise during the image acquisition process, and denoising is an essential step in most of image processing pipelines. Let u : Ω ⊂ R 2 → R be a clean image, and let f be a noisy image (i.e., the degraded version of u). The commonly used model of image degradation assumes that the observed image f is perturbed by a Gaussian noise of the form where η(0, σ 2 ) denotes an additive white Gaussian noise with zero mean and a variance σ 2 . Reconstruction of the clean image u is an inverse problem, which is generally ill-posed. Adding a regulariza-tion term, which provides a priori information, enforces a unique approximation of the solution. Rudin-Osher-Fatemi (ROF) model is one of the most famous works in this direction that uses the TV functional as a regularization energy [Rudin et al., 1992]. The main drawback of the ROF model is to produce the staircase effect, which appears at the extrema and at the boundary of an image. To improve the denoising performance of the ROF model, various types of weighted functions are incorporated to the TV regularization in order to control the amount of smoothing [Zhou and Li, 2013;Liu et al., 2016;Wang and Qi, 2018;Bazhanov et al., 2018;Lai et al., 2019;Thang et al., 2019;Pham et al., 2020].
Among filter-based algorithms, non-local means filtering (NLMF) was proven to be asymptotically optimal under a generic statistical image model [Buades et al., 2005b]. NLMF takes advantage of the redundancy of natural images. By this, some image blocks at different positions in a natural image have self-similar features, which means they have a certain correlation. NLMF estimates the denoised image by computing the average of all pixels in the image, weighted by the similarity between these pixels and the target pixel. The similarity is measured as a decreasing function of the weighted Euclidean distance [Buades et al., 2005a].
The image denoising problem can also be solved by exploiting the observation that natural images can be approximated by images with sparsity [Deledalle et al., 2011;Zhang et al., 2012;Ji et al., 2010;Gu et al., 2014;Gu et al., 2017]. In the work of Zhang et al. [Zhang et al., 2010], the authors presented an image denoising scheme comprising two stages. In the first stage, the noise is reduced in the PCA domain. In the second stage, the denoising procedure is iterated one more time to further enhance the denoising performance of the first stage. Experimental results showed that the proposed scheme achieves very competitive denoising performance compared with state-of-the-art algorithms. However, this approach also produces visual artifacts due to thresholding in the transformed domain.
Recently, deep learning networks have achieved great success on tasks of image restoration such as Lefkimmiatis, 2018;Zhang et al., 2018]. These denoising models yield significant gain of denoising capability, especially for noisy images which have low SNR. Yet, it requires a large amount of resources to train and deploy models. This is a crucial limitation of these learning models, especially when they are integrated in the pre-processing block of a computer vision system.
In this paper, we investigate a denoising algorithm for images corrupted by Gaussian noise. Our contribution is threefold. First, we propose a multi-stage denoising algorithm which comprises a PCA-based denoiser, an adaptive TV-based denoiser and a sharpening procedure. The proposed algorithm can achieve the denoising performance better than that of widely used algorithms. Second, we propose an adaptive TV-based model using Gaussian curvature for the second stage that allows the model to adjust the edge-preservation and the noise reduction abilities automatically. Third, we modify the thresholding technique of the PCA-based denoiser to reduce resulting artifacts.
The rest of the paper is organized as follows. The proposed algorithm is presented in Sections 2. Section 3 discusses experimental results.

Proposed algorithm
The proposed algorithm comprises three stages as shown in Fig.1. The first stage yields the initial estimation of the denoised image by using a modified patchbased PCA method with the Stein thresholding. The choice of PCA-based denoising is based on the fact that it can produce denoising results competive with the stateof-the-art denoising algorithms [Deledalle et al., 2011;Zhang et al., 2010]. Most of the noise is suppressed in this stage. Yet, visual artifacts are produced due to thresholding in the PCA domain and the influence of the noise. In the second stage, the output of the first stage will be further refined by using an adaptive TV-based denoiser to reduce visual artifacts while preserving edges. Since the TV-based denoiser reduces the contrast of the denoised image, a sharpening procedure is applied in the third stage to enhance the contrast of the image. Details of the proposed algorithm are described as follows.

First stage: Patch-based PCA denoising
The first stage of the proposed algorithm relies on the observation that natural images can be approximated by a sparse representation. PCA is applied to decompose an image into uncorrelated variables called principal components. In the PCA domain, the energy of the noiseless image concentrates on a subspace spanned by a small subset of components, while the energy of the noise distributes uniformly over the whole space. Thus, the shrinkage approaches can be applied to achieve a sparse representation in the transformed domain.
In order to perform the PCA transformation, we need a set of training samples of the noisy image f . For this purpose, we use a sliding window with the size of L × L to extract image patches, which are training examples for PCA. To reduce the computational complexity, the sliding window is applied with a stride S. Image patches inside a small image region are grouped into a block to characterize the image locally. Then, PCA is performed for each block. In this way, we can find resulting bases adapted to the local image regions as well as the whole image. A thresholding technique can be applied in the PCA domain to remove the noise. Details of the first stage are described as follows.
We denote by F the sample matrix for a block, which has the following form: where each row comprises n = L × L pixel values of a patch. Each block has m patches.
. . f in be the i-th row of the sample matrix F, which corresponds to a patch of the block. The mean value of each column (j = 1, ..., n) is computed, then placed into the mean vector M as where The sample matrix F is centralized by subtracting the mean vector M from each row of F: where H = 1 1 · · · 1 T is an m × 1 column vector.
The covariance of the centralized sample matrix F is computed as According to the PCA method, the singular value decomposition of the covariance matrix Σ is computed. Let λ 1 ≥ λ 2 ≥ ... ≥ λ n ≥ 0 be the eigenvalues of Σ and X 1 , X 2 , ..., X n ∈ R n be the corresponding eigenvectors. Eigenvectors are also called principal components and form an orthogonal basis. A patch F i of the block can be decomposed as It is clear that for a given patch, only a small subset of components are relevant, but they may vary between different patches. Thus, we denoise patches according to the magnitude of the projection onto the PCA space instead of the variances of principal components. For this purpose, we use the Stein thresholding S λ (·) to reduce resulting artifacts as: where λ is a threshold parameter. The thresholding operator is applied for patches of each block aŝ The overlap of patches creates redundancy of patch estimation. Therefore, for all pixels, the final estimate is the uniform average of candidate estimates.

Second stage: Adaptive TV-based denoising
Most of the noise is suppressed in the first stage, yet the output image contains visually unpleasant artifacts and noise residual. The main reason is due to the fact that the covariance matrix Σ is corrupted by noise, causing the estimation bias of the PCA transformation. In addition, only few components represent an image block so that it is unable to explain rare patches. Therefore, we propose to use a second stage to further refine the denoising performance. Among denoising methods, the ROF model is effective at simultaneously maintaining edges whilst removing noise in flat regions due to the TV-based regularizer. The ROF model is formulated as where ∇ stands for the gradient operator; α is a positive regularization parameter; | · | denotes the L 1 -norm. The first term in (10) measures the fidelity to the noisy image f . The second one, called the TV of the image u, is the regularization term. The TV term allows the existence of edges and is also a smoothing term. The ROF model is not adaptive since the regularization parameter is the same for all pixels. To improve the denoising performance of the ROF model, we modify it by introducing a weighted function into the TV-based regularizer in order to adjust the strength of smoothing for each pixel of the image.
To this end, we propose an adaptive TV-based denoising model as where α(·) is a weighted function, which is expressed by where α 0 -an initial regularization parameter; γ -a contrast parameter; Γ G θ * f is the Gaussian curvature of the smoothed image of f , which is defined as where v x , v y and v xx , v yy , v xy , v yx denote the first-order and the second-order derivatives of v with respect to x, y; denoting the Gaussian filter with a standard deviation θ. The Gaussian filter is used to reduce the noise, which affects the computation of the Gaussian curvature in (13). As can be seen from (11), (12), and (13), the proposed model has the strength of smoothing larger for homogeneous regions compared with that for edges. This property allows the model preserve edges, producing better denoising results compared to the original ROF model.
We use the split Bregman method [Goldstein and Osher, 2009] to solve the problem (11). According to the split Bregman method, we minimize the following discrete unconstrained problem of (11): where d is an auxiliary variable; b is a variable chosen through each iteration; µ is a positive penalty parameter. The problem (14) is solved by iteratively minimizing with respect to u and d separately. We obtain the u and the d subproblems as The u subproblem can be solved via the optimality condition using the Fourier transform as where r ∈ [0, M − 1] and s ∈ [0, N − 1] are the frequencies in the discrete frequency domain; F and F −1 stand for the forward and the inverse Fourier transforms; "-" stands for the point-wise division of matrices.
We solve the d subproblem by using the shrinkage operator [Setzer, 2011] as By iteratively computing the sequence of equations (17), (18) and (19), we can find the solution of the problem (14). Note that instead of computing the weighted function α(·) only once using the observed noisy image f , we iteratively refine it using the denoised image at each iteration. This step ensures that the weighted function is updated with the denoised image, resulting in more pleasant denoising results.

Third stage: Sharpening
Although the second stage can reduce noise while preserving edges, the TV term causes the loss of contrast in the denoised image [Strong and Chan, 2003;Li et al., 2014]. Since the noise and the artifacts are suppressed significantly after the first two stages, we can utilize the unsharp masking technique to enhance the denoising result in terms of contrast. Without the second stage, the sharpening procedure could not enhance the denoising results due to the noise and artifacts.
Let u(x, y) be the original image. Unsharp masking is created by subtracting the blurred image from the original image as where G ρ is a Gaussian filter with the parameter ρ.
Then a weighted portion of the mask is added to the original image to obtain a sharper image where β is a scaling constant.
We perform experiments on 14 widely used images (Fig. 2). Images are perturbed by additive white Gaussian noise with zero mean and a standard deviation σ. We consider three noise levels, including σ = 10, 20 and 40. We evaluate the performance of algorithms by using three criteria: PSNR, SSIM and visual quality. We, first, demonstrate the denoising performance after each stage of the proposed algorithm. Figure 3 shows the denoised images by each stage for the image Monarch.
The quantitative results are also shown for each denoising result. It can be seen that the first stage produces an image with a large amount of visual artifacts. Although the second stage nearly does not improve the PSNR measure, the SSIM value is gained significantly. One can see from the Fig. 3d that most of visual artifacts are reduced by the second stage while maintaining edges. Yet, its denoised result has lower contrast compared with that of the first stage. The sharpening stage significantly improves the denoising result in terms of both quantitative and qualitative criteria. Table 1 reports the PSNR and SSIM results of different methods. From Table 1 we have the following observations. First, on all noise levels, PCATVS achieves the best results for almost all images. By PSNR, the proposed algorithm surpasses PLPCA, ROF, and NLM by 0.42dB, 1.74db, and 1.3dB, respectively; by SSIM, PCATVS outperforms them by 0.0131, 0.0439, and 0.0448, respectively. Second, PCATVS outperforms competing algorithms by a large margin on images such as Barbara. This is because the image Barbara is dominated by a large amount of repetitive structures, which can be effectively exploited by PCA-based method. Besides, the second stage of the proposed algorithm enhances the denoising performance of the PCA-based stage by reducing visual artifacts. Figure 4 shows denoised images of Barbara by different methods. The original image is corrupted by a Gaussian noise with the standard deviation σ = 10. A portion of images is enlarged for visual evaluation. From Fig. 4, one can see that PCATVS and NLMF yield denoising results, which are more visually pleasing than that of PLPCA and ROF. Compared with PCATVS, lowcontrast textures are smoothed out by NLMF. PLPCA produces visual artifacts, whereas ROF generates blocky image regions and reduces the contrast of the image.
Figures 5 and 6 demonstrate denoised images of Jetplane and Monarch for the noise levels σ = 20 and 40, respectively. It can be seen that features, which are characterized by long and small edges such as characters, are well preserved by PCATVS. NLMF produces noisy pixels along edges. As the noise level increases, PLPCA and NLMF generate more visual artifacts, creating visually unpleasant results. The denoised images by ROF suffers the staircase effect.
We further evaluate methods on a computed tomography (CT) image. CT is one of the major tools in medical diagnostics. The quality of CT images depends on the amount of X-ray. As increasing the X-ray dose, the quality of CT images becomes better. However, highdose of X-ray makes negative impacts on patients. Due to the statistical uncertainty in physical measurements, low-dose CT images are noisy and difficult to observe [Erofeeva et al., 2018]. Thus, denoising low-dose CT images improves the quality of images without increasing the X-ray dose. CT images are often corrupted by speckle noise [Duan et al., 2016;Duan et al., 2015]. Fig.  7 shows denoising results by different methods on the CT image corrupted by speckle noise. It can be seen that artifacts ganerated by PLPCA and ROF are significant, whereas NLMF produce over-smooth image. The proposed algorithm can remove noise while preserving edges and fine details well.

Conclusion
In this paper, we have presented a denoising algorithm for images corrupted by Gaussian noise. The proposed algorithm consists of three stages. First, the noise is suppressed in the PCA domain using the Stein thresholding. The output of the first stage is further refined by a TV-based denoising model, which is adaptive with the reconstructed image. Finally, the denoised image is sharpened using the unsharp masking technique. Experimental results demonstrated that the proposed algorithm outperforms competing algorithm in terms of both qualitative and quantitative aspects, especially for images which have a large amount of repetitive structures. Experimental results on CT images showed the flexibility of the proposed algorithm.