Non-Convex High Order Total Variation With Overlapping Group Sparsity Denoising Model Under Cauchy Noise

It is widely known that the total variation regularization model preserves the edges well in the restored images but has some staircase effects. We consider using non-convex high-order total variation and overlapping group sparsity as a hybrid regularization to present a new denoising model. The proposed model can well preserve edges and reduce the staircase effect in the smooth region of the restored images. In order to solve the proposed hybrid model, we develop an efficient alternating minimization method. Compared with other models for removing Cauchy noise, numerical experimental results demonstrate that the superiority of the proposed model and algorithm, both in terms of visual and quantitative measures.


I. INTRODUCTION
Image restoration mainly includes image denoising and deblurring, which is an important problem in image processing. We all known that image restoration is a typically linear inverse problem [1]. The objective of image restoration is to establish an appropriate image degradation model, and then make a reasonable estimate of the original image by solving the inverse problem. In mathematics, the image degradation model can be described as follows: where u ∈ R n×n is the original image, f is the observed image, τ denotes the noise. The objective of noise reduction is to estimate the true clean image u from its noisy observation f . To deal with this problem, adding some regularization terms to energy is usually considered. Tikhonov regularization [2] is one of the commendable regularization method, where Tikhonov et al. use a quadric functional of the L 2 norm of the gradient magnitude ∇u i.e. ∇u 2 2 as the regularization term. However, this model produces over-smoothness The associate editor coordinating the review of this manuscript and approving it for publication was Yilun Shang . while removing noise. In order to overcome this drawback, Rudin et al. [3] proposed the following total variation (TV) regularization model(also known as ROF model) arg min u 1 2 u−f 2 2 + λ ∇u 1 , where λ > 0 is a regularization parameter. Although the classical TV regularization model has the ability of preserving the edges, it also suffers from the well-known staircase artifacts. In order to reduce the staircase artifacts, one choice is using high-order total variation (HTV) [4], [5] as a regularization. Although HTV regularization can alleviate the staircase effect, it may reduce the ability of edges-preserving. More recently, Liu et al. proposed a new model using the overlapping group sparsity (OGS) regularization [6] to reduce staircase artifacts. By combining TV and HTV regularization, some hybrid regularization models are proposed [7], [8] to overcome the staircase effects and preserve edges simultaneously.
The above models have a good effect in processing Gaussian noise, but have poor processing ability for non-Gaussian noise, such as Rician noise, Gamma noise and Cauchy noise. In this paper, we focus on restoring images with Cauchy noise, which usually appeared in radar and sonar, biomedical images, atmospheric and underwater acoustic [9]- [12]. In recent years, many methods have been proposed to remove Cauchy noise. In [13], Chang et al. considered an image restoration method with Cauchy noise based on the recursive recovery algorithm of Markov random field model. Achim et al. [14] removed Cauchy noise using a bivariate maximization posterior estimation in the complex wavelet domain. In [15], Wan et al. proposed a novel segmentation method for color images to remove Cauchy noise. Based on TV regularization and non-convex data-fidelity term, Mei et al. [16] proposed a non-convex model for Cauchy noise removal. Due to the nonconvexity and nonsmoothness of the proposed model, they developed a specific alternating direction method of multiplier to solve it. Although, they proved that their method converged to a stationary point, the restored results by their method strongly depends on the initializations. To overcome this disadvantage, Sciacchitano et al. [17] proposed the following convex variational model by adding a quadratic penalty function to the nonconvex data-fidelity term arg min where u 0 denotes the image obtained by applying the median filter to the noisy image f , ·, · X denotes the standard inner product, 1 is an n × n matrix of ones, µ > 0, α > 0 are the parameters. In order to reduce staircase artifacts and preserve edges, Yang et al. [18] replaced TV regularization in the model (2) by TV and HTV regularization and proposed a hybrid convex variational regularization model (''HTVAM'' for short). Moreover, they also introduced an efficient alternating minimization algorithm which can adaptively select regularization parameters. By combining total variation with overlapping group sparsity, Ding et al. [19] proposed the following convex variational model (''OGSTVL1'' for short) for eliminating Cauchy noise: where, φ(∇u) is the OGS-TV regular function. They utilized the alternating direction method of multiplies(ADMM) and majorization minimization(MM) method to propose an efficient algorithm to solve the model (3). The numerical experiments showed that their method maintains the edge-preserving property of the TV method and overcomes staircase effects effectively. On the other hands, the l 1 norm regularization terms in TV-based image restoration gives rise to over-smoothing problem. Compared with l 1 norm, l p (0 < p < 1) norm can measure sparsity better. In order to overcome this disadvantage, many non-convex l p norm regularization are introduced. In [20], Chartrand first proposed a non-convex optimization problem using l p norm as the objective function and he gave the theoretical l p reconfigurable conditions for arbitrary sparse signals. In addition, Wen et al. [21] further theoretically demonstrated the superiority of the l p norm-based method. In [22], Zhang et al. proved that non-convex l p norm regularization in TV-based image restoration lead to even better edge preservation when compared to the convex l 1 norm regularization in their numerical experiments.
Motivated by the above studies, we propose the following hybrid model combining non-convex HTV and OGS-TV regularization for removing Cauchy noise: where α > 0, µ > 0 and ω > 0 are regularization parameters, p ∈ (0, 1). To the best of our knowledge, this is the first time that non-convex HTV and OGS-TV regularization as a hybrid regularization term is considered for Cauchy noise. The main contributions of this paper are as follows:(1) We first combine non-convex HTV and OGS-TV regularization together in Cauchy noise removal model. This hybrid regularization model is beneficial to improve the OGS-TV sparsity and effectively remove the Cauchy noise in the degraded image.
(2)Using the technique of variable substitution, we develop a new algorithm under the framework of alternating direction method of multipliers (ADMM) to solve the proposed model. (3)We conduct extensive experiments to demonstrate that our proposed method has superior performance over the state-of-the-art methods in Cauchy denoising. By adding OGS-TV term to non-convex l p regularization, the edges and texture of the restored image are enhanced significantly. At the same time, this method greatly reduces the staircase artifacts caused by total variation. The structure of the paper is as follows. In Section II, we give the definitions of second-order TV function and OGS-TV function. In Section III, we propose an efficient alternating minimization algorithm for solving the proposed model. In Section IV, we give some numerical results to demonstrate the effectiveness of the proposed model and algorithm. Finally, we conclude this paper in Section V.

II. PRELIMINARY
In this section, we introduce some preliminaries that will be used in this paper.

A. TOTAL VARIATION OF THE MATRIX
u ∈ R n×n is a gray image. We introduce the following first-order forward and backward difference operators: According to the above operators, we can defined the second-order difference operators as where u i,j represents (i, j)th pixel in the image, i, j = 1, · · · , n. Then, we can obtain the discrete TV and HTV of u as:

B. OGS-TV
For any vector t ∈ R n , Selesenick and Chen [23] define a K-point group as And then, they use t i,K to define a group sparsity regularizer for one-dimensional case as Similarly, for the two-dimensional image matrix u ∈ R n×n , Liu [6] define a K -square-point group as x denotes the largest integer not more than x. Arranging the K × K of u i, j K in lexicographic order, we obtain a vector u(i, j) K . Then, we can defined the overlapping group sparsity functional of the two-dimensional array as Consequently, the regularization functional φ(∇u) in (3) can be given as:

III. PROPOSED METHOD
In this section, we develop an efficient alternating minimization method to solve our proposed nonconvex model (4). By introducing three auxiliary variables t, q, w, the minimization problem (4) can be changed into the following equivalent form arg min Firstly, we define the following augmented Lagrangian function of (6): where β 1 , β 2 , β 3 > 0 and λ i , i = 1, 2, 3 are penalty parameter and the lagrange multipliers, respectively. Based on the classical ADMM, starting at , the iterative scheme is implemented via the following subproblems: In the following, we solve the subproblems in (8) in the detail.
(a) u-subproblem By using the direct differential method, the Euler-Lagrange equation of (9)is given as follows:

VOLUME 9, 2021
Under the periodic boundary conditions for u, the solution of equation (10) can be obtained by using fast Fourier transform operation F and inverse fast Fourier operation F −1 where It is easy to see that we can use majorization minimization (MM) method [6] to solve the problem (12).
In order to minimize problem (15), we use the iteratively re-weighted l 1 algorithm (IRL1) [24]. First, we use the similar strategies in [24] to approximate the l p minimization problem (13) to the following weighted l 1 problem where the weight η i is given by the following updating scheme with ε being a small positive number to avoid division by zero. Let x k+1 = ∇ 2 u k+1 + λ k 2 β 2 , then minimization problem (14) can be given by one-dimensional shrinkage If we denote The gradient and Hessian matrices of G(v) are given by the following equations, The solution of problem (16) can be obtained by Newton's method where k + 1 and l + 1 are the numbers of outer iteration and Newton iteration, respectively. (e) Updating multipliers via λ 1 , λ 2 , λ 3 Finally, three Lagrange multipliers in (8) are updated as follows: The whole algorithm for removing Cauchy noise is given in Algorithm 1.

IV. EXPERIMENTS AND DISCUSSION
In this section, we present several numerical results to demonstrate the effectiveness of the proposed model and method. All test images are shown in FIGURE 1, which are gray images with the size of 256 × 256. All the numerical experiments are implemented under Windows 10 and Matlab R2018a running on a laptop equipped with 2.60 GHz Inter(R) Core(TM) i5-3230M and 4GB memory.

A. EXPERIMENT SETTING
In the following experiments, we use the following degradation equation to generate the noisy image f where ξ > 0 represents the noise level, τ 1 and τ 2 are independent random variables and follow from the Gaussian  distribution with mean 0 and variance 1. It is the same as the parameter setting in [19], we set γ = √ ξ for all the experiments in this paper.
The peak signal-to-noise (PSNR) and the structural similarity index measurement (SSIM) are used to measure the quality of the restored images. The PSNR is given by the following formulation: where u and u are the ideal image and the recovered image, respectively. Generally, the larger PSNR values imply that the restored images are better. The SSIM [25] is defined as follows: , where µ u , µ u , σ u and σ u represent the mean values and the standard variance of images u and u, σ u u is the covariance of u and u, C 1 and C 2 are two positive constants. Note that  The stopping criterion for all of the tested algorithms was set to where u k is the restored image in the k-th iteration. Throughout all the experiments, we set the penalty parameters β 1 = 30, β 2 = 2.5, β 3 = 10 for the noise level ξ = 0.02; and β 1 = 50, β 2 = 9.5, β 3 = 10 for the noise level ξ = 0.04.

B. PARAMETER DISCUSSION
In this subsection, we show how to select the best values for parameters α, µ, ω, p, K , Nit-inner for different test images. The images ''Peppers'', ''Boats'' and ''Starfish'' are used to study the choice of these parameters. Cauchy noise with the noise level ξ = 0.02 is added into these three tested images.
In order to find out how the parameters α and µ impact the performance of our proposed method. We first fix other parameters and let α to vary continuously from 3 to 6. Next, we use the same method to tune the parameter µ.
The changes of PSNRs and SSIMs versus different values of parameters α and µ in FIGURE 2. It is easy to see that the performance of the proposed method achieves the best with α in 4.5 nearby. We set α = 4.6 in the following experiments. It can also be observed that the maximum values of PSNR and SSIM are obtained when the parameter µ ∈ [6,12] and µ > 25. From FIGURE 2, we can observe that the values of PSNR and SSIM are stable when µ > 25. However, the calculation time increases with the increase of µ. Then, in the following experiments, we set µ = 6.
Next, we test how to select a good value of the group size K and the number of MM inner iterations Nit-inner. We carry out experiments to compute the values of PSNR, SSIM, CPU Times with respect to K and Nit-inner. From FIGURE 3, we can see that the highest PSNR and SSIM values are obtained at K = 3. Therefore, we choose K = 3 for in the following experiments. From  FIGURE 3, we can also observe that the values of PSNR and SSIM tend to be stable when the iteration number Nit-inner ≥ 10, increasing the MM iteration Nit-inner does not dramatically effect the PSNR, SSIM, but, the algorithm takes more time. Therefore, we set Nit-inner = 10 in our experiments.
Finally, we plot the values of PSNR and SSIM obtained by our algorithm against different values of parameters ω, p in FIGURE 4. It can be observed that that the highest PSNR and SSIM values are obtained when ω = 0.57, p ∈ (0.6, 0.8). Then we empirically choose ω = 0.57, p = 0.7, in the following experiments.

C. DENOISING
To demonstrate the superior performance of the proposed method, we compare it with other two state-of-the-art methods: OGSTVL1 [19] and HTVAM [18]. Experiment 1: We corrupt the images(as shown in FIG-URE 1) adding Cauchy noise with the noise level ξ = 0.02. FIGURE 5 shows the denoising results by three different methods. It is not difficult to see that our proposed method outperforms the others in terms of both noise removal and detail preservation. To intuitively compare the local details of the restored images, we enlarged the red box region of images in FIGURE 5. Comparing to OGSTVL1 and HTVAM, our proposed method can preserve details and edges very well in the restored images while removing most of the noisy pixels. In TABLE 1, we lists the values of PSNR and SSIM for the restored images by different methods. It is easy to see that the values of PSNR and SSIM for the restored images by our method are higher than the other two methods. It is consistent with the visual comparison. Moreover, our proposed method needs less time than the other two methods. Then we conclude that our proposed method behaves much better than the other two methods. Experiment 2: We restore the images ''Babyface'', ''Lena'' by three different methods under the noise level ξ = 0.04. The restored images are shown in FIGURE 7. In order to illustrate the differences between the restoration given by different methods, we shows the zoomed regions of the recovered images in FIGURE 8. From the recovered images, it is easy to see that the HTVAM method and OGSTVL1 method yields staircase artifacts in the denoising results, while our proposed method overcomes this drawback. The restored image obtained by the algorithm in this paper is closer to the original image, and the smoothness of the local area is more prominent. We also report the values of PSNR and SSIM obtained by these three methods in TABLE 2. From the table, we observe that our proposed method can get better restoration results than the other two methods, in terms of  the PSNRs and SSIMs. From FIGURES 7,8, and TABLE 2, we conclude that our proposed method behaves better than the HTVAM method and OGSTVL1 method.

D. THE CONVERGENCE OF THE PROPOSED METHOD
In this subsection, we plot the convergence curve to verify the convergence of our proposed methods experimentally. Here, we tested six images of ''Peppers'', ''Boats'',

V. CONCLUSION
In this paper, we proposed a new Cauchy noise removal model by combining nonconvex l p norm regularization with OGS-TV regularization. The proposed model inherited the advantage of both the HTV and OGS-TV to effectively reduce the stair casing artifacts while preserving sharp edges well. Based on ADMM method, IRL1 method and MM algorithm, we developed an efficient alternating minimization algorithm. Numerical results demonstrated that the proposed model and method performed better both quantitatively and qualitatively.
JIANGUANG ZHU received the Ph.D. degree in applied mathematics from Xidian University, Xián, China, in 2011. He is currently an Associate Professor with the Shandong University of Science and Technology. His current research interests include optimization theory, algorithms and applications in image processing, and sparse representation.