Constrained Minimization Problem for Image Restoration Based on Non-Convex Hybrid Regularization

It is widely known that the classic total variation(TV) model has been proven to be very effective in preserving sharp edges. However, the TV model suffers from the staircase effects which produce blocking artifacts in the restored images. In this paper, we propose a new hybrid regularization model by combining non-convex second order total variation with wavelet transform to restrain staircase effects and protect some details of the images. To compute the new model effectively, we propose an alternating minimization method for recovering images from the blurry and noisy observations. The new model is first transformed into several sub-problems, and the generalized iterated shrinkage algorithm, the Fourier transform method and projection method are used to solve these sub-problems, respectively. Numerical experiments show that the proposed model can restrain blocking artifacts while projecting sharp edges, and the restoration quality outperforms several state-of-the-art methods.


I. INTRODUCTION
With the development of computer technology, image restoration which is a fundamental problem in the field of image processing has attracted more and more attention. Image restoration refers to restoring clean images from degraded images which contaminated by noise and other factors, such as camera shake, illumination and so on. Mathematically, the image degradation can be described by the following mathematical model: where u ∈ R s is the original image and f ∈ R s is an observed blurred and noisy image, η ∈ R s denotes the additive Gaussian white noise, and matrix K ∈ R s×s is a linear operator representing the convolution, here s ∈ M × N . It is well known that image restoration problem is a typical ill-conditioned problem, and it is quite difficult to solve it directly. Hence, to solve the ill-conditioned problem, we first The associate editor coordinating the review of this manuscript and approving it for publication was Wei Zhang. need to transform it into a well-posed problem. However, for now, regularization is the best way which transforms ill-conditioned problem into well-posed problem. The linear system (1) can be replaced by the following optimization problem min u 1 2 where · 2 is the Euclidean norm, the the positive parameter µ is the regularization parameter and φ(u) refers to some regularization functions, such as, Tikhonov regularization [1] and total variation regularization [2]. It is generally known that Tikhonov regularization has a defect that tends to make images overly smoothed and fails to protect some details such as sharp edges. A good regularization function can improve the effect of image restoration. In order to effectively preserve edges and some details in the restored image, Rudin, Osher and Fatemi [2] introduced TV regularization that can protect the sharp edges of images. The image restoration model with total VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ variation regularization can be described as where ∇u 1 is a discretization of the total variation about u, ∇ ∈ R s×s denotes the local finite difference operator. Due to the high nondifferentiability and nonlinearity of the objective function, the total variation model (3) is difficult to solve. In order to address this problem, many efficient methods have been proposed to solve the total variation model [3]- [7]. With the deepening of research, scholars have found that total variation regularization also has a disadvantage, which causes staircase effects in piecewise affine signal regions in the presence of noise [8]. To overcome this shortcoming, several methods have been proposed to restrain the staircase effects. One common method is to use high-order regularization schemes which can deal with the staircase effects while preserving the edge information in the restored image. Lysaker, Lundervold, Tai in [10] proposed second-order TV regularization model for additive noise removal in medical magnetic resonance images. Another more effective method to deal with staircase effects is the total generalized variation (TGV) [14], [15] that is a generalization of TV. To simultaneously reduce staircase artifacts and preserve edges, some hybrid regularization models are proposed [12], [13].
All the models mentioned above are convex, which are easy to solve and can guarantee the convergence [16]. In addition, they can be efficiently solved by many convex optimization algorithms. On the other hand, many studies [17]- [25] have shown that non-convex regularizations are superior to convex ones, due to their edge preserving property. There are two classes of non-convex regularization methods for image restoration problems. The first class is using a non-convex function to construct a non-convex restoration model. For example, Zhang et al. [18] defined two non-convex function: the Exponential-Type (ET) function [19] and the Geman function [20], and proposed a non-convex model for image restoration, which is composed of a non-convex data fitting term and a TV regularization term for blur and noise removal. In [21], Nikolova et al. proposed a amily of cost functions for signal and image recovery. Their idea is to use the l 1 data fitting term combined with non-convex regularization. By combining nonconvex first-order TV and nonconvex second-order Laplace, Tang et al. [22] proposed a generalized hybrid nonconvex variational regularization model. The second class is to replace the l 1 norm regulation with the l p (0 < p < 1) norm regularization. In [23], the authors replaced the l 1 norm of gradient by the nonconvex l p (0 < p < 1) norm to improve the quality of TV denoising results. Using a convex combination of non-convex TV( u p 2 (0 < p < 1)) and nonconvex higher-order TV( 2 u p 2 (0 < p < 1)), Oh et al. [24] propose a novel hybrid variational model for the image denoising problem. Later, Liu et al. [25] extend the nonconvex hybrid regularizer model to the restoration of noisy blurred images. At present, there are three main methods for solving l p norm-based problems: the iterative re-weighted l 1 (IRL1) [26] minimization algorithm, the iterative re-weighted least square (IRLS) [27] algorithm and generalized iterated shrinkage algorithm (GISA) [28].
In general, most natural images have complex texture structures. The traditional restoration model can not meet the requirements of some details recovery. Recently, Image restoration models based on total variation and wavelet transform [29] have been proposed and widely used in the domain of compressed sensing. These models can make full use of the sparsity prior information and image restoration is more effective in details. In addition, wavelet transform can suppress the staircase effects created by TV-based models. The basic idea of wavelet transform based approaches is that the images can be sparsely approximated by properly designed wavelet transform. In most cases, the regularization of model based wavelet transform is l 1 norm.
Inspired by the above mentioned advantages of non-convex regularization, second-order total variation regularization and wavelet transform, we propose a hybrid non-convex regularization image restoration model based on second-order total variation and wavelet transform. The model for deblurring is to solve the following optimization problem: In the practical applications, due to the noises in the observed image f , it is not exactly to make Ku = f . In this case, instead of (4), we consider the following constrained optimization problem: where the parameter δ > 0 is an estimate of the noise level, W represents separately the wavelet transform which defined as Haar wavelet in this model, ∇ 2 u = (∇ xx u, ∇ xy u, ∇ yx u, ∇ yy u) denotes the second-order discrete gradient of u. For more details about the second-order discrete gradient refers to [6]. In order to solve (5), we design an alternate minimization algorithm based on the classical alternating minimization method and projection method. The fast iterative technique is applied to the projection method to accelerate the efficiency of the algorithm. To effectively handle the l p norm, we adopt a generalized iterated shrinkage algorithm [28]. The effectiveness of the performance of the proposed method is compared with other methods.
The paper is organized as follows. In Section 2, we use the proposed alternate minimization algorithm to solve the proposed model (5). In Section 3, we present numerical results and performance comparisons. Finally, we conclude the paper in Section 4.

II. ACCELERATED ALTERNATING MINIMIZATION METHOD FOR PROPOSED IMAGE RESTORATION MODEL
In this section, we consider the constrained optimization problem (5), and we develop the accelerated alternating minimization method for solving the proposed image restoration model.
The accelerated alternating minimization method is proposed on the basis of the well-known variable-splitting and penalty technique. Next, we need to introduce three auxiliary variables ω, v, r, and transform the optimization problem (5) into the following constrained problem: Then, the corresponding augmented Lagrangian function is where, λ 1 , λ 2 and λ 3 are the Lagrange multipliers, which control positive penalty parameters β 1 , β 2 and β 3 going to infinity, respectively. Next, each sub-problem is solved alternately. Firstly, it is not difficult to discover that, for fixed u k , v k , λ k 1 , λ k 2 and λ k 3 , the minimization of (7) with respect to ω can be written as The minimizer of ω sub-problem can be solved by the generalized soft-thresholding algorithm (GST) [28], which was proposed by zuo et al. and is very effective for solving the minimization problem with l p norm.
The function T GST p in [28] is defined as where sgn(·) is the signum function, S Secondly, for fixed ω k+1 , v k , λ k 1 , λ k 2 and λ k 3 , the minimization of (7) with u can be written as Then, The optimal value of u k+1 must satisfy the following first-order necessary optimality conditions by using the direct differential method, where (∇ 2 ) T representing the conjugate operator of ∇ 2 . W representing the wavelet transform and W T W = I . Under the periodic boundary condition for u, (∇ 2 ) T ∇ 2 and A T A are all block circulant matrices, more details see [30], and the matrices can be diagonalizable by the 2-dimension fast discrete Fourier transforms [31]. Then, the solution of (11) can be obtained by two discrete Fourier transform DFT (involving one inverse DFT).
Thirdly, for fixed ω k , u k+1 , λ k 1 , λ k 2 and λ k 3 , the minimization of (7) with v can be defined as Due to formula (12) is a l 1 norm optimization issue, therefore, the solution of v can be given by the one dimensional shrinkage operator: Fourthly, for r sub-problem, it can be solved by the following form VOLUME 8, 2020 where = {r ∈ R s | r 2 ≤ δ}, P [·] denotes the projection operator onto the ball under Euclidean norm. To speed up the above iteration algorithm, we add a fast iterative technique which proposed in [32]. Then, the result of r can be rewritten as After using the fast iterative technique, the update iteration process of r is . Finally, we update the Lagrange multipliers λ k 1 , λ k 2 and λ k 3 by λ k+1 where the parameter ξ is a relaxation parameter which control the convergence of the proposed algorithm.
We summarize the solution steps into the following algorithm.

III. EXPERIMENT RESULTS
In this section, we do several sets of numerical experiments to test the effectiveness and efficiency of our method, in comparison with FTVd [3], GISA [28] and ADMM-TV2L1 [33]. The iterations for images restoration by different methods are terminated when the following stop criterion is met, The quality of the recovered images is evaluated by peak-signal-to-noise ratio (PSNR), improvement signal-to-noise ratio (ISNR) and relative error (Rerr), which are defined as: where, MSE is the Mean-Squared-Error per pixel, f , u 0 , u are the observed image, original image, recovered image, respectively.

A. PARAMETERS SETTINGS
In this subsection, we will consider how the parameters α and p in our model can be selected properly. First, to illustrate how to select the suitable parameter p, we test the images ''Barbara'', ''Lena'', ''Cameraman'', and ''Baboon'' which are degraded by Gaussian blur (15,15) and Motion blur (25,35) with Gaussian noise levels δ = 0.001 and δ = 0.01, respectively. We plot the PSNR values for our proposed method against different values of parameter p in Fig. 2. It can be observed that the maximum PSNR values are obtained when the parameter p ∈ (0.5, 0.8). We set p = 0.7 in the following experiments. Next, in order to check how the parameter α impacts the performance of our proposed method, we test the images ''Barbara'', ''Lena'' and ''Cameraman'' which are degraded by Average blur (21,21) and Motion blur(75, 135) with Gaussian noise levels δ = 0.001 and δ = 0.02, respectively. In Fig. 3, we plot the changes of PSNR values for our algorithm against different values of parameter α. It can be seen that the maximum PSNR values are obtained when α ∈ (0.2, 0.8). Then we empirically choose α = 0.5 in the following experiments.

B. THE COMPARISON WITH FTVd, GISA AND ADMM-TV2L1
In the previous statement, our proposed method can suppress the staircase effects and protect details. Hence, in this    subsection, in order to illustrate the superiority of our method, we compare our proposed method with FTVd [3], GISA [28] and ADMM-TV2L1 [33].
The image ''Lena'' which has complex texture structure is used to illustrate the performance of our method. Several kinds of blur kernels are added the image to show that our  proposed method has universal applicability. The blur kernels used are list in Table 1, containing two groups of blur degree: serious blur kernels and slight blur kernels. In this test, the noise level is δ = 0.02.  In Fig. 4, we show the resorted results by four different methods for the ''Lena'' image under slight degree of blur kernels. To make a comparison in detail, we enlarge a region(the red area of Fig.4) of images in Fig.5. From these  figures, we can easy observe that the restored image by our method is much better than other three methods in terms of the staircase effect and edge preservation. In Fig. 6, we plot the convergence results about PSNR and Rerr versus VOLUME 8, 2020 iteration for different methods. We observe that our method is convergent and get higher PSNR values than the other three methods. In order to explain the advantages of the proposed method more comprehensively, we test the three images under slight degree of blur kernels. In Table 2, we list the PSNR and ISNR values for the three test images achieved by FTVd, GISA, ADMM-TV2L1 and the proposed method. The comparison of the PSNR and ISNR values in Table 2 also shows that our method yields a better restoration result. In Fig. 7, we restore the ''Lena'' image degraded by the serious degree of blur kernels. To further illustrate the superiority of the our method, we display the zoom parts of the recovered images in Fig. 8. As can be seen from the figures, our method outperforms the other three methods. One can observe that our method avoids entirely the staircase effect and the details are better recovered. We also report the PSNR and ISNR values for the three test images achieved by four different methods under serious degree of blur kernels. From Table 3, our proposed method can also get better restoration results than the other three methods, in terms of the PSNR and ISNR values.
It is well known that image restoration is sensitive to noise level. Therefore, we consider an experiment on the image ''baboon'' degraded by Average blur with 21 * 21 and the different levels of noise. The restored images produced by the FTVd, GISA, ADMM-TV2L1 and the proposed method are shown in Fig. 9. For a better visual comparison, we have enlarged some details of the restored images in Figs. 10. It is not difficult to see from the visual quality of restored images in Figs. 9 and 10 that our method is better than the other three methods since the image recovered by our method preserves more details and looks more natural. We report the PSNR and ISNR values of the restored images by four methods under the three different levels of noise in Table 3. From the table, we see that our method has higher PSNR and ISNR values than the other three methods. Therefore, we conclude that the proposed method performs better than the FTVd, GISA and ADMM-TV2L1. In Fig. 11, we plot the convergence results about PSNR and Rerr values iteration for different methods under different levels of noise. It also shows that our method is convergent and get higher PSNR values than the other three methods

IV. CONCLUSION
In this paper, we propose a new non-convex hybrid regularization model with constraint based on second-order total variation and wavelet transform. Meanwhile, we propose a new method which combines alternating direction method with the fast iterative technique to solve the proposed model. Our experimental results show that the proposed method for the new model is better than several state-of-the-art methods in terms of the quality of image restoration and the ability to suppress the staircase effects. The PSNR, ISNR and Rerr values obtained by our proposed algorithm are also superior to the recovered results by several state-of-the-art methods. Her current research interests include image processing, pattern recognition, wavelet analysis, and sparse representation. VOLUME 8, 2020