Total Variation with Overlapping Group Sparsity for Image Deblurring under Impulse Noise

The total variation (TV) regularization method is an effective method for image deblurring in preserving edges. However, the TV based solutions usually have some staircase effects. In order to alleviate the staircase effects, we propose a new model for restoring blurred images under impulse noise. The model consists of an ℓ1-fidelity term and a TV with overlapping group sparsity (OGS) regularization term. Moreover, we impose a box constraint to the proposed model for getting more accurate solutions. The solving algorithm for our model is under the framework of the alternating direction method of multipliers (ADMM). We use an inner loop which is nested inside the majorization minimization (MM) iteration for the subproblem of the proposed method. Compared with other TV-based methods, numerical results illustrate that the proposed method can significantly improve the restoration quality, both in terms of peak signal-to-noise ratio (PSNR) and relative error (ReE).


Introduction
Image deblurring and denoising problems have been widely studied in the past decades. In the literature, it is widely assumed that observed images are the convolution of standard linear and space invariant blurring functions with true images plus some noise. Let g denote the blurred and noisy image, h the blur kernel, f the original image and η the noise. The image f is assumed to be a real function defined on a bounded and piecewise smooth open subset Γ of R 2 . In general, the image formation process can be modeled as: g = h ? f + η, where "?" denotes the two-dimensional convolution operation. Image deblurring is to estimate the true image f from the blurred and noisy image g. As is well known, image deblurring and denoising is a typically illposed problem [1,2]. To handle this problem, regularization technique is usually considered to obtain a stable and accurate solution. In this way, we need to solve the following problem where the first term is called the regularization term, the second term is called the fidelity term (' 2 -fidelity), μ > 0 is the regularization parameter, and ψ is the regularization functional. Without loss of generality, a discretized image may have n 1 × n 2 pixels. Simply, in our work, we assume that n 1 = n 2 = n, then f, g and η are vectors of length n 2 . It is easy to extend to any n 1 × n 2 image. Let H be the corresponding blurring matrix of n 2 × n 2 from h [3]. Then the discretized form of the minimization problem (1) is equivalent to the following matrix-vector form where k Á k 2 denotes the Euclidean ' 2 norm. Notice that H is a matrix of block circulant with circulant blocks (BCCB) structure when periodic boundary conditions are applied [3]. How to choose a good regularization functional is an active area of research in the imaging science. In the early 1960s, D. L. Phillips [4] and A. N. Tikhonov [5] proposed the definition of ψ as an ' 2 -type norm (academically called Tikhonov regularization), that is, c ¼ kLf k 2 2 with L an identity operator or difference operator. Although the functional ψ of this type has the advantage of facilitating the calculations, it is rarely used in current practice because it has the drawback of penalizing discontinuities in resulting solutions, for instance, over-smoothing edges. Therefore, this is not a good choice since natural images have many edges.
To overcome this drawback, many different types of regularization functionals have been proposed. One popular model was introduced by Rudin, Osher and Fatemi (ROF) in [6]. They proposed a total variation (TV) regularization with an ' 2 -fidelity term (' 2 -TV) for image restoration. Its corresponding minimization task is: where BV(Γ) denotes the space of functions of bounded variation. kfk TV is defined by kf k TV ≔ X 1⩽i;j⩽n kðrf Þ i;j k 2 ¼ X 1⩽i;j⩽n ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi j ðr x f Þ i;j j 2 þ j ðr y f Þ i;j j 2 q which is called as isotropic TV, or kf k TV ≔ X 1⩽i;j⩽n kðrf Þ i;j k 1 ¼ X 1⩽i;j⩽n j ðr x f Þ i;j j þ j ðr y f Þ i;j j which is named as anisotropic TV, where k Á k 1 denotes the Euclidean ' 1 norm. Operator r: R n 2 ! R 2×n 2 denotes the discrete gradient operator (under periodic boundary conditions) which is defined by (rf) i,j = ((r x f) i,j , (r y f) i,j ), with ðr x f Þ i;j ¼ ( f iþ1;j À f i;j if i < n; f 1;j À f n;j if i ¼ n; ðr y f Þ i;j ¼ ( f i;jþ1 À f i;j if j < n; f i;1 À f i;n if j ¼ n; for i, j = 1, 2, Á Á Á, n, where f i,j refers to the ((j − 1)n + i)th entry of the vector f (it is the (i, j)th pixel location of the n × n image, and this notation remains valid throughout the paper unless otherwise specified). More details for the definition of kfk TV can be referred to [7,8].
Many methods have been proposed to solve the restoration model (3) such as the fast TV deconvolution (FTVd) [9,10], the augmented Lagrangian method (ALM) [11][12][13][14], the dual methods [13,15], and the split Bregman method [11,16]. These studies provided many efficient iterative algorithms to make the computation on TV related models be more convenient for Gaussian noise removal. Furthermore, in many cases, the noise does not satisfy the Gaussian assumption, for instance, the noise may follow a Laplace distribution [17]. There has been a growing interest in using an ' 1 -fidelity term instead of the ' 2 -fidelity term for image restoration in many literature such as [8,[18][19][20][21][22][23][24][25][26][27] for considering another non-Gaussian noise-impulse noise. We consider the corresponding regularization model with an ' 1 -fidelity term as: A classic approach is to use the TV regularizer by ψ(f) = kfk TV , which we call ' 1 -TV. Recently, Yang et al. [20] used the FTVd method to solve the ' 1 -TV model fast. Guo et al. [8] proposed a fast algorithm for image restoration in the ' 1 -TV model. Their method was to add a penalty term by using the variable substitution method, which belongs to penalty methods in optimization. They employed an alternating minimization method to solve it. They also proved the convergence of their method. Numerical tests showed that their method got better restorations and faster than FTVd. Wu et al. [25] used ALM to solve the ' 1 -TV model. More recently, Chan et al. [26] proposed a constrained total variation (TV) regularization method for image restoration for ' 1 -TV. Their method used a box constrained projection to ensure the restored images stay in a given dynamic range. They used the alternating direction method of multipliers (ADMM) [11][12][13][14][28][29][30][31][32][33] to solve the model based on augmented Lagrangian method. They got better results than other methods such as FTVd and ALM. Their numerical results showed that for some images where there are many pixels with values lying on the boundary of the dynamic range, the gain could get very high numerical superiority in the peak signal-tonoise ratio. As far as our knowledge goes, the box constrained projection in image restoration is significant both in practice and theory.
Although the TV regularization using in the restoration problems can recover sharp edges of a degraded image, it also gives rise to some undesired effects and transforms smooth signal into piecewise constants, the so-called staircase effects [34,35]. To overcome this deficiency, one usual method is to replace the original TV norm by a high-order TV norm. The highorder TV regularization schemes have been studied so far mainly for overcoming the staircase effects while preserving the edges in the restored image. More details can be referred to [35][36][37][38]. The high-order TV based methods may have some other behaviors. For example, it may transform the smooth signal to over-smoothing and take more time to compute.
More recently, Bredies et al. [39] proposed total generalized variation for image restoration with overcoming the staircase effects while preserving the edges. Pryre and Fadili [40] considered group sparsity and overlapping group sparsity (OGS) for image processing, such as denoisng, compressed sensing. Their numerical experiments showed that the PSNR (defined in Section 5) gain by OGS is a consistent improvement on comparing with group sparsity across a wide range of natural images. Selesnick and Chen [41] proposed an OGS TV regularizer to one-dimensional signal denoising. They applied the majorization minimization (MM) method to solve their model. Their numerical experiments showed that their method can overcome staircase effects effectively. Their method has the disadvantages of the low speed of computation and the difficulty to directly extend to the two-dimensional case. However, all these methods only considered the ' 2 -fidelity term under Gaussian noise. To our knowledge, the methods for image deblurring and denoising under impulse noise with these regularization terms are still missing in the literature. Particularly, Liu-Huang-Liu in [27] proposed a hybrid model for image deblurring and denoising under impulse noise, but their model is not convex and it may take more time since its regularization term is the combination of TV and high-order TV with adapted parameter selection.
In this paper, inspired by the works from [40] and [41], we propose a new model for images deblurring under impulse noise by setting ψ in (4) to be the OGS-TV functional. In this model we first extend the OGS-TV functional in [41] to general two-dimensional case as a new regularization term, then consider an ' 1 -fidelity term for application in images deblurring under impulse noise. Moreover, we impose a box constraint to the proposed model to obtain more accurate solutions similarly as [26]. Using the technique of variable substitution, we provide an efficient algorithm to solve the model under the framework of ADMM. We use an inner loop which is nested inside the majorization minimization (MM) iteration for the subproblem of the proposed method. Our main contributions are firstly combining the former three parts together in one model in application and making the 2D TV with OGS case be more easily solved than the 1D case in [41]. According to a brief explanation of our model and our numerical results, we can observe that our new model using the OGS-TV regularizer coincides to maintain the edge-preserving property of TV methods and overcome the staircase effects, which is similar as total generalized variation [39]. In addition, the numerical results also show that our method is very effective and competitive with other TV-based methods, such as Chan et al.'s method [26] and Guo et al.'s method [8], especially on getting higher PSNR and better visual quality.
The outline of the rest of this paper is as follows. In the next section, we will briefly introduce the definition of the OGS regularization functional. We will also review the MM method and ADMM, which are used in our proposed method. In Section 3, we propose the OGS-TV based model for recovering images under blur and impulse noise. We also provide the efficient solving algorithm in this section. The numerical results are given in Section 4. Finally, we conclude this paper in Section 5.

OGS-TV
In [41], the authors denoted a K-point group (K denotes the group size) of the vector t 2 R n by Note that t i,K can be seen as a block of K contiguous samplings of t staring at index i. With the notation (5), a group sparsity regularizer for one-dimensional case is defined in [41] as Similarly, we can define a K-square-point group of a two-dimensional signal such as images considered in this work v 2 R n 2 , where vector v is stacked in column-wisely, in other words, the (i, j)th entry of a matrix is assigned to be the ((j − 1)n + i)th entry of the vector v. Clearly, . . . v iþK r ;jÀK l v iþK r ;jÀK l þ1 Á Á Á v iþK r ;jþK r where K l ¼ ½ KÀ1 2 , K r ¼ ½ K 2 and [x] denotes the largest integer less than or equal to x. The group size is denoted by K 2 . Note thatṽ i;j;K;K can be seen as a square block of K × K contiguous samplings of v with the center at index (i, j). Here we choose a group entries around the objective point rather than a group following the objective point like one-dimensional in [41] because of the faster and easier computation in the experiments. Moreover, to the best of our knowledge, the former is much better than the latter in image restoration because the pixels in image are related to or depended on all the ambient pixels rather than partial surrounding pixels. Let v i,j, K,K be a K 2 -vector obtained by arranging the K × K elements ofṽ i;j;K;K in lexicographic order. This notation also remains valid throughout the paper unless otherwise specified. Then the overlapping group sparsity functional of the two-dimensional array can be defined by From the definition above, we can easily get that this function is convex. Consequently, we define the regularization functional ψ in (4) to be the form We call the regularizer ψ in (9) as the OGS anisotropic TV functional because we handle the r x f and r y f separately (while the isotropic TV is defined differently by , and call the corresponding convex minimization model (4) L1-OGS-ATV.

The MM method
The MM method is an asymptotical method in solving optimization problems, for instance, a minimization problem of the form as follows, where α is a positive parameter and the functional φ is defined in (8). The point of the MM method is that, instead of directly solving the difficult minimization problem P(v), the MM approach solves a sequence of easier optimization problems Q(v, v k ) (k = 0, 1, 2, . . .) firstly and then manages to get the minimizer of P(v). Generally, an MM iterative algorithm for minimiz- When P(v) is convex, then under the former conditions, the sequence v k produced by (11) converges to the minimizer of P(v) [42,43]. We aim to solve the special problem (10), and it is obvious that P(v) in (10) is convex. Therefore, the MM approach is available for solving the problem (10). While we were concluding this manuscript, we became aware of very recent related work in [44]. The authors in [44] have studied the problem (10) elaborately, which is a subproblem of our method. Moreover, for the sake of completeness, we briefly introduce the solving method here by our way independently.
First of all, to derive an efficient algorithm with the MM scheme for solving the problem (10), we want to find a majorizor of P(v). Here, we only need to find a majorizor of φ(v) because of the simple enough quadratic term of the first term in (10). Note that for all v and u 6 ¼ 0 (u, v 2 R n 2 ) with equality when u = v. Substituting each group of φ(v) into (12) and summing them, we get a majorizor of φ(v) with Sðv; uÞ ! φðvÞ; Sðu; uÞ ¼ φðuÞ; provided kv i,j,K,K k 6 ¼ 0 for all i, j. After simple calculation, S(v, u) can be rewritten as where C(u) is independent of v, and Λ(u) is a diagonal matrix with each diagonal component with m = 1, 2, Á Á Á, n 2 . The entries of Λ can be easily computed by using Matlab built-in function conv2. Then a majorizor of P(v) can be easily given by with the solutionv where I is an identity matrix with the same size of Λ(v k ). We can easily get that Λ 2 (v k ) is also a diagonal matrix with each diagonal component [Λ 2 (v k )] m,m equaling to the form of removing the out root of the right term of (16). Moreover, the inversion of the matrix I þ 1 a L 2 ðv k Þ can be computed very efficiently since it only requires simple componentwise calculation. Therefore, we obtain the Algorithm 1 for solving the problem (10).

Algorithm 1
The MM method for solving (10)

Variable splitting and ADMM
Consider an unconstrained optimization problem in which the objective function is the sum of two functions as where ϕ i : R n i ! R are closed proper convex functions, χ i R n i are closed convex sets, A i 2 R l×n i , and b 2 R l is a given vector. The augmented Lagrangian function ( [28]) of (20) where λ 2 R l is the Lagrange multiplier, β is a penalty parameter which controls the linear constraint, and C does not depend on x 1 , x 2 . The idea of ADMM is to find a saddle point of L. Usually, ADMM consists in alternated minimizing L on x 1 , x 2 , λ, for instance, minimizing L with respect to x 1 by fixing x 2 and λ. That delivers to the following simple but powerful algorithm classic ADMM.
Algorithm 2 Classic ADMM for the minimization problem (20) initialization: ; until a stopping criterion is satisfied.
According to [12], we can see classic ADMM is convergent because of the nonexpansive and absolute summable properties of the x 1 and x 2 subproblems. More details can be found in [12,[28][29][30][31]. However, the convergence speed is not too fast. In order to speed the convergence, we can introduce a step length parameter γ for updating the multiplier [29,32,33]. The algorithm framework is outlined as follows called general ADMM.

Proposed Method
With the definition of (9), we will consider a minimization problem of the form (L1-OGS-ATV) Similarly as [26], we will solve the problem We refer to this model as CL1-OGS-ATV. Obviously, this model is also convex. Particularly, each term of the model (24) has the properties of additivity and separability. Therefore, we can rewrite it as follows (under the periodic boundary conditions), From the Eq (25), we can observe that our model (24) can be seen as a combination of n × n coupled subproblems, which are all (i, j) terms of the last line in (25). Each subproblem is approximate to the original TV regularization model (4) with ψ(f) = kfk TV in a local region, which has the edge-preserving property. More details please refer to [45]. Therefore, our model (24) can coincide to 1) maintain the edge-preserving property of TV methods, and 2) have the property of smoothing the local regions, which can be seen as overcoming the staircase effects. However, we did not solve our model decoupled as (25) in this work, so we choose the following algorithm to solve (24) based on ADMM.
For the model (24), by introducing new auxiliary variables v x , v y , z, w, we transform the minimization problem (24) to the equivalent constrained minimization problem Note that the constraint is now imposed on w instead of f. The augmented Lagrangian function of (26) is where β 1 , β 2 , β 3 > 0 are penalty parameters and λ 1 , λ 2 , λ 3 , λ 4 2 R n 2 are the Lagrange multipliers.
According to the scheme of general ADMM mentioned above (Algorithm 3), for a given (v k x ; v k y ; z k ; w k ; f k ; l k 1 ; l k 2 ; l k 3 ; l k 4 ), the next iteration ðv kþ1 x ; v kþ1 y ; z kþ1 ; w kþ1 ; f kþ1 ; l kþ1 1 ; l kþ1 2 ; l kþ1 3 ; l kþ1 4 Þ is generated as follows: (27) with respect to v x and v y . The minimizers are obtained by v kþ1 It is obvious that problems (28) and (29) match the framework of the problem (10), thus the solutions of (28) and (29) can be obtained by using Algorithm 1, respectively. 2. Compute z k+1 .
The minimization with respect to z can be given by the well-known Shrinkage [20] explicitly: where j Á j, sgn and "" represent the componentwise absolute value, signum function, and componentwise product, respectively.
The minimizer is given explicitly by 4. Compute f k+1 by solving the normal equation where " Ã " denotes the conjugate transpose, see [13] for more details. Since all the parameters are positive, the coefficient matrix in (32) are always invertible and symmetric positive definite. In addition, note that H, r x , r y and their conjugate transpose have BCCB structure under periodic boundary conditions. We know that the computations with BCCB matrix can be very efficient by using fast Fourier transforms (FFTs).

Update the multipliers via
Based on the discussions above, we present the algorithm on ADMM using inner MM iteration for solving the convex CL1-OGS-ATV model (24) shown as Algorithm 4.

Algorithm 4 CL1-OGS-ATV-ADM4 for the minimization problem (24)
initialization: CL1-OGS-ATV-ADM4 is a special form of general ADMM for the case with two blocks of variables (v y , v x , w, z) and f. If the Step (1) of Algorithm 4 can be solved exactly, the convergence for CL1-OGS-ATV-ADM4 can be guaranteed. In this case, if the relax parameter g 2 ð0; ffi ffi 5 p þ1 2 Þ, Algorithm 4 is convergent, more details can be referred to [29,30,32,33]. Besides, although step (1) of Algorithm 4 can not be solved exactly, we can find a convergent series to ensure the convergence as [12]. Particularly, our numerical experiments verify the convergence of Algorithm 4.

Numerical Results
In this section, we present several numerical results to illustrate the performance of the proposed method. We compare our method CL1-OGS-ATV-ADM4 ("Ours" for short) with other state-of-the-art methods, Chan et al.'s ADM2CTVL1 proposed in [26] ("CTY" for short, Algorithm 2 in [26] for the constrained TV-L1 model), Guo et al.'s fast ' 1 -TV proposed in [8] ("GLN" for short) and a high-order method derived by ourselves.
All experiments are carried out on Windows 7 32-bit and Matlab 2010a running on a desktop equipped with an Intel Core i3-2130 CPU with 3.4 GHz and 3.4 GB of RAM.
The quality of the restoration results is measured quantitatively by using the peak signal-tonoise ratio (PSNR) in decibel (dB) and the relative error (ReE): where f and f denote the original and restored images respectively, and Max I represents the maximum possible pixel value of the image. In our experiments, Max I = 1. The stopping criterion used in our work is set to be where F k is the objective function value of the respective model in the kth iteration, which is We only consider the salt-and-pepper noise in our experiments, since the variation method is easy to extend to the random value noise case.

Study on the rest parameters
Firstly, we set the group size parameter K = 3 to find a good maximum inner iterations NIt.
Our experiments are on the image "Cameraman" blurred by Gaussian blur kernel with 7 × 7 and standard deviation 5 and corrupted by 40% salt-and-pepper noise. The results are shown in Table 1. From Table 1, we can choose maximum inner iterations NIt = 5 for the best. Here we find that the larger the inner iteration NIt (> 5) is, the lower PSNR values are, which is a drawback for the selection of NIt. This is because we do not tune the other parameters to be best for different NIt. After we test more similar tests, and due to the CPU time and little error between different selections of NIt, we choose NIt = 5 for balance with the following experiments. Then we fix NIt = 5 and repeat more experiments for choosing a good group size  parameter K. We operate the three 256-by-256 images (a) "Cameraman", (b) "Satellite", and (c) "House" for this best option of parameter K. The results are shown in Fig 2. From the figure, we can see that K = 3 is better for all the tests both on CPU time and PSNR. Therefore, we fix that NIt = 5 and K = 3 for balance in our work. Then, we test how to select a good regularization parameter μ for different images. We will point out several important advantages of our method in the following experiments. For the sake of simplicity, we focus on the above three test 256-by-256 images (a) (b) and (c). Under the Gaussian blur with 7 × 7 window size and standard deviation 5, the images corrupted by added salt-and-pepper noise from 30% to 60% are tested. In Fig 3, we plot PSNR, ReE, and Time for our algorithm against different values of the regularization parameter μ. Each row in Fig 3 corresponds to the four salt-and-pepper noise levels. In fact, for all μ, our method always gives high PSNR values. Moreover, the PSNR curves of our method are very flat, which shows that our method is stable for a wide range of μ, which is wider than that in [26]. That is to say, our method is more robust than the method in [26]. This is the first advantage of our method.
In addition, as far as we know, a good image restoration algorithm should satisfy the following two properties. a) It is fast and can reach good results in term of both numerical values and high visual quality. b) It is not sensitive to parameters. Our method meets the requirement of these two properties, as it obtains good restoration results with the same parameters for different images under the same blur and noise. This is the second advantage of our method. Here and in the following experiments, under the same blur and noise, we choose the same parameters for all the test images. Particularly, for the images under the Gaussian blur with 7 × 7 window size and standard deviation 5 and corrupted by salt-and-pepper noise from 30% to 60%, we set μ = 100, 80, 60, 40 respectively. After similar tests as Fig 3, we list all the selection rule of μ: for the images under the Gaussian blur with 15-by-15 window size and standard deviation 5 and corrupted by salt-and-pepper noise from 30% to 60%, μ = 120, 110, 100, 90 respectively; for the images under the average blur with 7 × 7 window size and corrupted by salt-and-pepper noise from 30% to 60%, μ = 100, 80, 60, 40 respectively.

Comparison with CTY and GLN for the test image "Cameraman"
In this subsection, we mainly compare our proposed method to CTY and GLN for deblurring problems under salt-and-pepper noise. We use image "Cameraman" for the experiments in this subsection. Our purposes are (1) to show the improvement of PSNR and to demonstrate the efficiency of our proposed method mainly via a comparison with CTY and GLN, and (2) to illustrate that our proposed method can overcome the staircase effects effectively and get better visual quality than CTY and GLN, which is the third advantage of our method. Firstly, we generate the blurred images by two Gaussian blurs (i) and (ii) with periodic boundary conditions as mentioned above, and then corrupt the blurred images by salt-andpepper noise from 30% to 60%. For CTY and GLN, we have tuned the parameters manually to give the best PSNR improvement. The numerical results by the three methods are shown in Table 2. From the table, we see that both our proposed method and CTY are much faster and can get higher PSNR than GLN. Our proposed method needs the fewest iterations than the other two methods, and the time is always close to CTY. Particularly, the iterations of the GLN method always reach the maximum number of iterations which we set to be 200.
We also show the images restored by the three methods. We display the degraded images and the restored images by three methods under two Gaussian blurs (i) and (ii) and 50% level The tests are on the images "Cameraman", "Satellite", "House" that are blurred by 7 × 7 Gaussian blur kernel with standard deviation 5 and corrupted by salt-and-pepper noise from 30% to 60%. From left to right, results on PSNR, ReE, CPU time respectively. From top to bottom, results on noise level 30%, 40%, 50%, 60% respectively. We can easily see the third advantage of our proposed method that our method can overcome the staircase effects effectively and get better visual quality than others. Moreover, we also plot the evolution of the PSNR over time and iterations for the three different methods in Fig 5 for the image blurred by 7 × 7 Gaussian blur and corrupted by 40% level salt-and-pepper noise.
From the expriments and the description in [8], the GLN method has three sensitive parameters that depend on blur, noise level and test images rather than only one sensitive regular parameter μ in CTY and our proposed method. Besides, the results of GLN is nearly same as FTVd, while CTY is much better than FTVd in [26]. From the above tests, we also observe that both our proposed method and CTY can get better results than GLN. Moerever, the staircase Table 2. Numerical comparison of the fast ℓ 1 -TV method (GLN) [8], the ADM2CTVL1 method (CTY) [26], and our proposed method (Ours) under two Gaussian blurs (Bls) (i) fspecial('gaussian', 7, 5) and (ii) fspecial('gaussian', 15,5) and corrupted by salt-and-pepper noise from 30% to 60%.  The PSNR results can be found clearly from Table 2.
doi:10.1371/journal.pone.0122562.g004 effects by GLN method is also existent. Therefore, we omit the following comparison with GLN and only list the comparison with CTY. Remark 2. Here and in the following tests for CTY, we tune the regularization parameter μ to be optimal by checking the highest PSNR and the lowest ReE under a "for loop" of Matlab code from 1 to 70 by step length 1 for different images under different blurs and noise levels. Besides, in the experimentsunder the 15 × 15 Gaussion blur with standard deviation 5, we find that if we increase the inner penalty parameter β 2 we will get higher PSNR. We do not change the parameter in our proposed method because of the good properties of a good restoration algorithm introduced above. However, in Table 2 for the CTY method under 15 × 15 Gaussion blur, we set the inner penalty parameter β 2 = 50 instead of β 2 = 20 in [26] for higher PSNR. For other tests following we also set this inner penalty parameter of CTY to be as default (β 2 = 20) in [26].

Comparison with CTY for other test images
In this subsection, we focus on comparisons between our proposed method and CTY for deblurring problems under salt-and-pepper noise. Fifty-six degraded test images are generated in the way similar to that in the last subsection. That is, we first generated the blurred images operating on images (b)-(h) with the periodic boundary condition by two blurs Gaussian blur (i) (also as G) and average blur (iii) (also as A), then corrupted the blurred images by salt-andpepper noise from 30% to 60%. The parameters of our proposed method are set as above description and the parameters of CTY as Remark 2.
Conclusions similar to those in the last subsection can be made based on the results in Table 3 and Table 4. For example, our proposed method is always more accurate, with a possible improvement of more than 2.90 dB in PSNR (see image (c) with Gaussian blur and a 30% level of noise). For all images, the lower the noise level is, the better the improvement of PSNR will be. Even in high noise level, our method is more accurate than CTY by sacrificing partial time. For a further step, the iterations of our method are fewer than CTY for almost all test.
Finally, in Fig 6 and Fig 7, we display the zoom parts of the degraded images examples and the zoom parts of the restored test images by two methods respectively for Gaussian blur and average blur with noise level from 30% to 60%. We can easily see the visual improvement in the images by using our method. More specifically, from the third column in Fig 7, although the numerical result does not improve too much, the image visual quality of our proposed method is much better than CTY. This superiority is obvious for almost all the test images in our work.

Comparison with other TV-based methods
For image restoration under Gaussian noise, there has been many literature for high order TV or other TV-based methods, such as Lysaker-Lundervold-Tai (LLT) model [36], Wu-Tai's work [13], total generalized variation [39]. However, to our knowledge, the methods for image deblurring and denoising under impulse noise with these TV-based methods are still missing in the literature. More recently, Liu-Huang-Liu in [27] proposed a hybrid model for image deblurring and denoising under impulse noise, but their method might spend more time due to the combination.
In order to get a fair comparison of our method and other TV-based methods, we take the high-order TV model by changing the ' 2 -fidelity term in LLT to ' 1 ("HOTV" for short) and the total generalized variation method by changing the ' 2 -fidelity term in [39] to ' 1 ("TGV" for short) as examples. Other TV-based methods can be similarly treated. Moreover, we also impose the box constraint in these examples and solve them by ADMM to compare them with our method under Gaussian blur (i) and 30% to 40% impulse noise for diffident images. We only list the results for three images, "Cameraman", "House" and "Weatherstation" for examples. The codes for HOTV and TGV are projected by ourselves based on ADMM. Furthermore, we think that the superiority of our method can be seen clearly from these tests. The numerical results are summarised in Table 5. The numerical improvement by our method is obvious from the table. Moreover, TGV usually performs better than HOTV but spends more time. Particularly, the time of each iteration by HOTV is less than other two methods ours and TGV. In Fig 8, we display the zoom parts of the degraded images examples and the zoom parts of the restored test images respectively for Gaussian blur with noise level 40% by three methods. We can see the visual improvement in the images by using our method. In addition, after test more experiments, we find that TGV preforms better than our method on the smooth region of the image, and our method performs better on the edges. Remark 3. In this work, we do not compare our method with other two stage methods such as Cai et. al [23] and Yan [24], because our method is a variational method by imposing a new regularization term. Our contribution is that we proposed a new convex model for handling the impulse noise without distinction of salt-and-pepper noise or random-value noise. We could easily apply our new regularization term in the former two stage methods and would consider it in future.

Discussion and Conclusion
In this work, we study a new regularization model by applying TV with OGS in the classic ' 1 -TV model for the image deblurring under impulse noise. We provided the efficient algorithm CL1-OGS-ATV-ADM4 under the framework of general ADMM. In particular, an MM inner iteration is proposed to solve the subproblem instead of Shrinkage [20] in the classic ' 1 -TV model. The numerical results illustrate that our method outperforms CTY [26], GLN [8] and some other TV-based methods both in terms of the PSNR values, ReE, iterations or image visual quality. Our main contributions are firstly combining three existent parts together in one convex model for image deblurring under impulse noise without distinction of salt-and-pepper Several random examples of degraded and restored images of CTY and our method under Gaussian blur. Top row, zoom parts of blurred and noisy images under 7 × 7 Gaussian blur with standard deviation 5 and corrupted by salt-and-pepper noise. Second row, zoom parts of restored images by CTY respectively. Third row, zoom parts of restored images by our proposed method respectively. The PSNR results can be found clearly from Table 3.  Several random examples of degraded and restored images of CTY and our method under average blur. Top row, zoom parts of blurred and noisy images under 7 × 7 average blur and corrupted by salt-and-pepper noise. Second row, zoom parts of restored images by CTY respectively. Third row, zoom parts of restored images by our proposed method respectively. The PSNR results can be found clearly from Table 4.  noise or random-value noise. Furthermore, we could easily apply our new regularization term in the popular two stage methods and would consider it in future.