Multichannel Color Image Denoising via Weighted Schatten p-norm Minimization

1 Introduction

Noise corruption is inevitable during the image acquisition process and may heavily degrade the visual quality of an acquired image. Image denoising is thus an essential preprocessing step in various image processing and computer vision tasks [Chatterjee and Milanfar, 2010; Ye et al., 2018; Wang et al., 2018b]; moreover, it is also an ideal test platform for evaluating image prior models and optimization methods [Roth and Black, 2005]. As a result, image denoising remains a challenging yet fundamental problem. Early denoising algorithms were mainly devised on the basis of filter and transformation [Dabov et al., 2007b], such as wavelet transform and curvelet transform [Starck et al., 2002]. State-ofthe-art denoising methods are mainly based on sparse representation [Dabov et al., 2007b], low-rank approximation [Gu et al., 2017; Xie et al., 2016], dictionary learning [Zhang and Aeron, 2016; Marsousi et al., 2014; Mairal et al., 2012], nonlocal self-similarity [Buades et al., 2005; Dong et al., 2013a; Hu et al., 2019] and neural networks [Liu et al., 2018; Zhang et al., 2017a; Zhang et al., 2017b; Zhang et al., 2018].

∗Co-corresponding Author.

Current methods for RGB color image denoising can be categorized into three classes. The first kind involves applying grayscale image denoising algorithms to each channel in a channel-wise manner. However, these methods ignore the correlations between R, G and B channels, meaning that unsatisfactory results may be obtained. The second method is to transform the RGB color image into other color spaces [Dabov et al., 2007a]. This transform, however, may change the noise distribution of the original observation data and introduce artifacts. The third type of method involves making full use of the correlation information across each channel and conduct the denoising task on R, G and B channels simultaneously [Zhang et al., 2017a].

Similar to the case of grayscale images, noise from each channel can be generally, regarded as additive white Gaussian noise. However, the noise levels of each channel are diverse due to camera sensor characteristics and the imagery environment, such as fog, haze, illumination intensity, etc. Moreover, the on-board processing steps in digital camera pipelines may also introduce different noise strength assigned to different channels [Nam et al., 2016]. This reveals the difficulties and challenges faced by RGB color image denoising. Intuitively, if we apply the grayscale image denoising algorithm to color images in a band-wise manner without considering the mutual information and noise difference between each channel, artifacts or false colors could be generated [Mairal et al., 2008]. An example of this is presented in Fig.1. Here, we apply a representative low-rank-based method, WSNM [Xie et al., 2016], and our proposed method, MCWSNM (see Section 2), to the Kodak PhotoCD Dataset 1 . It can be seen from the figure that WSNM retains a large amount of noise and, to some extent, introduces some artifacts. Thus, the issue of how noise differences in each channel should be modeled is key to designing a good RGB color denoising algorithm.

One well-known color image denoising method is color block-matching 3D (CBM3D) [Dabov et al., 2007a]. This is a transformation field method that may destroy the noise distribution across each channel and introduce artifacts. The method in [Zhu et al., 2016] concatenates similar patches of RGB channels into a long vector. However, this concatenation operation regards each channel equally and ignores the noise level differences across each channel. For a long time, low-rank theory has been widely used in many fields, such as recommendation system [Wang et al., 2018a], community search [Fang et al., 2020b; Fang et al., 2020a], multi-out task [Liu et al., 2019] and multi-label[Liu et al., 2017], etc.

1 http://r0k.us/graphics/kodak/

Recently, based on the non-local self-similarity of an image and SVD, several low-rank-based denoising algorithms have been proposed, such as weighted nuclear norm minimization (WNNM) [Gu et al., 2017] and WSNM [Xie et al., 2016]. The nuclear norm and Schatten p-norm [Xie et al., 2016] are two surrogates of the rank function [Xu et al., 2017a]. In WNNM and WSNM, singular values are assigned different weights using a weight vector. It has been proven that WNNM is a special case of WSNM and that WSNM outperforms WNNM. Multi-channel WNNM (MCWNNM) [Xu et al., 2017b] is an extension of WNNM from grayscale image to color image that operates by introducing a weight matrix to adjust each channel's contribution to the final results based on noise level.

In this paper, based on the low-rank property of the nonlocal self-similar patches and the different noise strength across each channel, we propose a multi-channel weighted Schatten p-norm minimization (MCWSNM) model for RGB color image denoising. More specifically, considering a small local RGB patch in a noisy image, we first find its nonlocal similar cubic patches in a search window with an appropriate size. These similar cubic patches are then vectorized and concatenated to construct a noisy low-rank matrix; then Schatten p-norm minimization method is performed to recover the noise-free low-rank matrix. Moreover, in order to make full use of the redundant information across each channel, a weight matrix is introduced to balance each channel's contribution to the final denoising results. This model can be solved via alternating direction method of multipliers (ADMM) framework, and each variable can be updated with a closed-form solution. In addition, the convergence property of our proposed algorithm is also provided.

2 Our ProposedMethod

2.1 Method Formulation

The image denoising task involves recovering a clean image xc from its noisy observation data yc, where c = {r, g, b} is the index of the R, G and B channels. Here, we assume that the noisy image is corrupted by additive white Gaussian noise nc with σc = (σr, σg, σb), where σr, σg, σb denote the noise standard deviation in the R, G, B channels, respectively. Mathematically, yc = xc + nc. Low-rank-based denoising methods utilize the low-rank property of the grouped nonlocal similarity patches; this group procedure is illustrated in Fig.2. Given a noisy RGB color image yc, each local patch of size s × s × 3 is extracted and stretched into a vector y = (yT r , yT g , yT b )T ∈ R3s2 , where yr, yg, yb ∈ Rs2 are the corresponding patches in the R, G and B channels. For each vectorized local patch y, we find its M most similar patches (including y itself) by Euclidean distance in a proper size search window around it. We then stack the M similar patches vector column by column to obtain the noisy patch matrix Y = X + N ∈ R3s2×M , where X and N are the corresponding clean and noise patch matrix, respectively. It is clear that Y is a low-rank matrix and X can be recovered using a low-rank approximation algorithm.

As a nonconvex surrogate of the rank function, the weighted Schatten p-norm of a matrix Z ∈ Rm×n is defined as ‖Z‖w,Sp = (∑min{m,n} i=1 wiσp i ) 1p , 0 < p < 1, where w = (w1,w2, . . . ,wmin{m,n})T is a non-negative weight vector and σi is the i-th singular value of Z. The weighted nuclear norm [Gu et al., 2017] is a special case of the weighted Schatten p-norm when the power p = 1. It has been proven that the Schatten p-norm with 0 < p < 1 is a better approximation of rank function than the nuclear norm [Zhang et al., 2013].

Intuitively, when denoising a noisy color image, if a channel is heavily polluted, the contribution of this channel to the final denoised results should be less, and vice versa. Thus, a weight matrix W assigned to the noise strength of each channel is introduced. Based on the low-rank property of the non-local self-similarity, we propose the following multichannel WSNM (MCWSNM) model to recover the local lowrank patch X:

where λ > 0 is a tradeoff parameter used to balance the

Algorithm 1 minX ‖Y -X‖2F + λ‖X‖pw,Sp via GST

Input: Y , weight {wi}ri=1 , λ, p, K

• 1: Singular value decomposition Y = UΣV T , Σ =

diag{σ1, σ2, . . . , σr}; 2: for i = 1 : r do 3: τ GST p (wi, λ) = (2λwi(1-p)) 1 2-p +λwip(2λwi(1p)) p-1 2-p ; 4: if |σi| ≤ τ GST p (wi) then 5: δi = 0; 6: else 7: δ(k) i = |σi|; 8: for k = 0, 1, . . . , K do 9: δ(k+1) i = |σi| -wip(δ(k) i )p-1 ; 10: k = k + 1; 11: end for 12: δi = sgn(σi)δ(k) i ; 13: end if 14: end for 15: ∆ = diag(δ1, δ2, . . . , δr); Output: X∗ = U∆V T

data fidelity term and the regularization term, while W = diag(τrI, τgI, τbI), I ∈ Rs2×s2 is the identity matrix, (τr, τg, τb) = min{σr, σg, σb}/(σr, σg, σb).

2.2 Model Optimization

We utilize the variable splitting method to the MCWSNM model (1). By introducing an augmented variable Z, model (1) can be reformulated as a linear constrained optimization problem. Thus, model (1) can be expressed as follows:

Because the two variables X and Z in problem (2) are separable, this programming problem can be solved using the ADMM framework. The augmented Lagrangian function of (2) is

where L is the augmented Lagrangian multiplier and ρ > 0 is the penalty parameter. By taking derivatives of L with respect to X and Z and setting the derivative function to be zero, the variables can be alternatively updated as follows: (1) Xk+1

(2) Zk+1

This is a weighted Schatten p-norm minimization problem. It was noted in [Xie et al., 2016] that if the weights follow the non-descending permutation, (5) can be equivalently transformed into independent nonconvex ℓp-norm subproblems,

Algorithm 2 Solve MCWSNM via ADMM

Input: data Y , weight W, p, K1, µ > 1, tol > 0;

• 1: Initialization: X0 = Z0 = 0, L0 = 0, ρ0 > 0, flag = False, k = 0;
• 2: while flag == False do
• 3: Update X by using (4);
• 4: Update Z by solving (5);
• 5: Update L: Lk+1 = Lk +ρk(Xk+1 -Zk+1);
• 6: Update ρ: ρk+1 = µρk;
• 7: k = k +1;
• 8: if (Convergence condition is satisfied) or (k ≥ K1) then
• 9: flag = True;
• 10: end if
• 11: end while

Output: X∗ .

the global optima of which can be efficiently solved by the generalized soft-thresholding (GST) algorithm. This is summarized in Algorithm 1.

1. ρk+1 = µρk, (µ > 1).

The above alternative updating steps are repeated until either the convergence condition is satisfied or the number of iterations exceeds a preset K1. The ADMM algorithm converges when 1) ‖Xk+1 - Zk+1‖F ≤ tol, 2) ‖Xk+1 - Xk‖F ≤ tol and 3) ‖Zk+1-Zk‖F ≤ tol are simultaneously satisfied; here, tol is a small tolerance value. The complete updating procedures are summarized in Algorithm 2.

2.3 Convergence and Complexity

Although the proposed MCWSNM (1) is nonconvex, the global optimum also can be obtained for the limitations of weights permutation in GST. From Fig.3, it can be seen that ‖Xk+1 -Xk‖2, ‖Zk+1 -Zk‖2 and ‖Xk+1 -Zk+1‖2 simultaneously approach 0 during the iteration process. This curve is based on a synthetic experiment on Kodak PhotoCD Dataset and the test image is "kodim01" (see Experiments Section in detail).

We next provide a theorem that theoretically guarantees the convergence property of Algorithm 2.

Theorem 1. Assume that the weights w are sorted in a non-descending order and the parameter ρk is unbounded; then the sequences {Xk}, {Zk} and {Lk} generated in Algorithm 2 satisfy: (a) limk→+∞‖Xk+1 - Zk+1‖F = 0; (b) limk→+∞ ‖Xk+1 - Xk‖F = 0; (c) limk→+∞ ‖Zk+1 - Zk‖F = 0.

Proof. We first prove that the sequence {Lk} generated by Algorithm 2 is upper bounded. UkΛkV T k denotes the SVD of matrix 1 ρk Lk + Xk+1 in the k + 1 iteration, where Λk is the diagonal singular value matrix. By using the GST algorithm for WSNM, we have Zk+1 = Uk∆kV T k , where ∆k = {diag(δ1 k , δ2 k , . . . , δn k )} is the diagonal singular value matrix after generalized soft-thresholding operation. Then:

Hence, the sequence {Lk} is upper bounded.

We then prove that the sequence of the Lagrange function {L(Xk+1,Zk+1,Lk+1, ρk+1)} is also upper bounded. The inequality L(Xk+1, Zk+1, Lk, ρk) ≤ L(Xk, Zk, Lk, ρk) always holds since we have the globally optimal solution of X and Z in their corresponding subproblems (step 3 and step 4 in Algorithm 2. Based on the updating rule of L, it yields

Since {Lk} is upper bounded, the sequence {Lk+1 -Lk} is also upper bounded. Denote by a the upper bound of {Lk+1Lk} for all k ≥ 0, i.e. {Lk+1 - Lk} ≤ a, ∀k ≥ 0. We therefore conclude that

Thus, {L(Xk+1, Zk+1, Lk+1, ρk+1)} is upper bounded. We next prove the sequences of {Xk} and {Zk} are upper bounded. Since {L(Xk,Zk,Lk, ρk)} and {Lk} are upper bounded and

and the sequence {W(Y -Xk)} and {Zk} are upper bounded. Since Lk+1 = Lk + ρk(Xk+1 -Zk+1), {Xk} is also upper bounded. Thus, there exists at least one accumulation point for {Xk, Zk}. More specifically, we obtain that

and that the accumulation point is a feasible solution to the objective function. Thus, equation (a) is proved.

Finally, we prove that the change of sequence {Xk} and {Zk} in adjacent iterations tends to be 0. For Xk+1 and Xk = 1 ρk-1 (Lk - Lk-1) + Zk, we have

This completes the proof.

The main computational cost in a single iteration of Algorithm 2 consists of two parts. The first part involves updating X, in which the time complexity is O(max{s4M,M 3 }). The second part involves updating Z. The predominant cost of this updating is the SVD and the calculation of X∗ . Its complexity is O(min{s4M, s2M 2 } + s2 r4M). The costs for updating L and ρ can be ignored. Therefore, the total time complexity of our MCWSNM for solving problem (1) is about O(K1M 3 ).

3 Experiments

To illustrate the performance of the proposed MCWSNM, we implement MCWSNM on synthetic and real color noisy images. We also compare the proposed method with several recent proposed methods, including NC [Lebrun et al., 2015], NCSR [Dong et al., 2013b], PGPD [Xu et al., 2015], MCWNNM [Xu et al., 2017b], DnCNN [Zhang et al., 2017a], FFDNet [Zhang et al., 2018] and IRCNN [Zhang et al., 2017b]. In more detail, NC is an online blind image denoising platform, NCSR and PGPD are designed for grayscale images and MCWNNM is a color image denoising method while DnCNN, FFDNet, and IRCNN are three methods that operate on the basis of convolutional neural networks. All the parameters of the comparison algorithms are either optimally assigned, or chosen as described in the reference papers.

Noise level of MCWSNM and most comparative methods should be provided as parameters. In the synthetic experimental case, the noise (σr, σg, σb) in the R, G and B channels are assumed to be known. In the case involving real noise, we utilize the noise estimation method outlined in [Chen et al., 2015] to estimate the noise level of the noisy image for each channel. We implement the grayscale denoising methods, NCSR and PGPD, on the color noisy images in a bandwise manner with the corresponding channel noise level.

3.1 Synthetic Noisy Color Image Experiments

The Kodak PhotoCD Dataset is first utilized in our synthetic experiments. It includes 24 color images, each of which is either 768 × 512 or 512 × 768 in size. The noisy image is generated by adding zero mean Gaussian noise with σr = 40, σg = 20 and σb = 30. In MCWSNM, we set the local search Figure 5: Denoised images of "kodim1" in Kodak PhotoCD Dataset window size of each patch as 20, the similar patch number M = 70, each patch size s = 6, K1 = 10, K = 4, c = 2√2, λ = 0.6, p = 0.999 and ρ = 3.

Figure 6: Denoised images of cropped "kodim3" in the Kodak PhotoCD Dataset with σr = 40, σg = 20 and σb = 30.

The line graph of the PSNR results for MCWSNM and the comparison methods are presented in Fig.4. It can be seen from this graph that our proposed method outperforms the other competing methods in most cases. Moreover, Fig.5 and Fig.6 show the denoised results of "kodim1" and "kodim3"

in Kodak PhotoCD Dataset, respectively. Visual inspection reveals that MCWSNM can obtain the best visual result and clear texture information. Color artifacts are generated by NCSR and PGPD, while NC over-smooths the image. While MCWNNM, DnCNN, FFDNet and IRCNN can get a clean image, some of the detailed information is lost compared with the original image.

3.2 Real Noisy Color Image Experiments

We implement the WCWSNM on a real noisy color image dataset, CC [Nam et al., 2016], to evaluate the method's performance. This dataset includes 11 static scenes. The noisy images were captured in an indoor environment with different cameras and camera settings. For each scene, 500 images were taken using the same camera and camera settings. The mean image of each scene was then calculated to generate the noisy-free images, which were roughly regarded as the ground truth (GT). As the size of the original image is very large, 60 cropped smaller images with size 512 × 512 were provided in [Nam et al., 2016]. In this paper, we select 10 cropped images (see Fig.7) to conduct the experiment.

For MCWSNM, we set the local search window size for each patch as 20, the similar patch number M = 60, each patch size s = 5, K1 = 10, K = 4, c = 2√2, λ = 0.6, p = 0.999 and ρ = 3. Because the GT are provided, a quantitative assessment of each method is accessible. PSNR results for different cameras and camera settings are provided in Table

1. It can be seen from the table that the proposed MCWSNM obtains the highest PSNR value in most noisy images. Fig. 8 and Fig. 9 show the visual results of #7 and #10 images in CC dataset, respectively. Compared with other methods, MCWSNM gets the best visual effect. It not only remove the noise completely, but also preserve the detailed information effectively.

damental Research Funds for the Central Universities under Grant 413000092 and 413000082.