Image Deblurring under Impulse Noise via Total Generalized Variation and Non-Convex Shrinkage

Image deblurring under the background of impulse noise is a typically ill-posed inverse problem which attracted great attention in the fields of image processing and computer vision. The fast total variation deconvolution (FTVd) algorithm proved to be an effective way to solve this problem. However, it only considers sparsity of the first-order total variation, resulting in staircase artefacts. The L1 norm is adopted in the FTVd model to depict the sparsity of the impulse noise, while the L1 norm has limited capacity of depicting it. To overcome this limitation, we present a new algorithm based on the Lp-pseudo-norm and total generalized variation (TGV) regularization. The TGV regularization puts sparse constraints on both the first-order and second-order gradients of the image, effectively preserving the image edge while relieving undesirable artefacts. The Lp-pseudo-norm constraint is employed to replace the L1 norm constraint to depict the sparsity of the impulse noise more precisely. The alternating direction method of multipliers is adopted to solve the proposed model. In the numerical experiments, the proposed algorithm is compared with some state-of-the-art algorithms in terms of peak signal-to-noise ratio (PSNR), structural similarity (SSIM), signal-to-noise ratio (SNR), operation time, and visual effects to verify its superiority.


Introduction
Image deblurring under impulse noise is important in digital image processing research. It is widely applied in many aspects of science and engineering such as space detection, photodetection, remote sensing technology, medical image analysis, and geological exploration. In the process of acquiring, storing, and transmitting images, owing to limitations of imaging equipment and conditions, images are prone to quality degradation, such as image blurring, noise pollution, contrast decline, and loss of details. Image restoration of degraded images is an inverse problem whose solution is often ill-posed [1]. For such types of problems, regularization proved to be an effective method. This involves the proper use of some prior information of the degraded image as a fidelity term and an image optimization model that contains an efficient solution algorithm as a regularization term, because recovering a clear image is of great importance for advanced image processing in the later stages.
The total variation (TV) regularization proposed by Rudin, Osher, and Fatemi [2], also known as the ROF model, proved to be one of the most advanced regularizations in the field of image reconstruction. Based on the TV regularization, some extended models and their algorithms were widely applied to image restoration. For example, Wang et al. [3] proposed a fast total variation deconvolution (FTVd) model for image deblurring. The model uses the total variation (TV) regularization based on the assumption that the image is piecewise smooth; the method can remove image noise sufficiently, but it is prone to producing a "staircase effect". Furthermore, the FTVd model adopted the L1 norm to express the sparsity of impulse noise. However, it cannot depict the statistic characteristics of impulse noise.
Recently, some TV extensions emerged. Sakurai et al. [4] proposed a four-directional total variation model to extend the vertical and horizontal gradient information in the TV regularization into four directions. Ren et al. [5] proposed a fractional-order TV (FTV) model and tried to extend the integer-order TV model to the fractional order. Selesnick and Chen [6] proposed an overlapping group sparsity TV (OGSTV) model, in which the neighborhood information of the image gradient was taken into account and applied in conjunction with an L2-norm fidelity term to the Gaussian-noise denoising of one-dimensional signals. Liu et al. [7] extended the model to a two-dimensional model based on an L1-norm fidelity term and applied the extended model to image deconvolution and impulse noise removal, achieving good results [8]. Bredies et al. [9] proposed a total generalized variation (TGV) model based on total variation, which simultaneously constrains the first-order and second-order gradients of the image and has many excellent mathematical properties such as convexity, lower semi-continuity, and rotation invariance; it is also capable of approximating any number of items. The TGV model proved-in theory and practice-to be able to remove the staircase effect effectively while preserving the edges and details of the image. TGV regularization terms attracted wide attention in the academic community and were applied to the reconstruction of nuclear magnetic resonance images [10] and seismic signal denoising [11,12]. Inspired by the advantages of TGV regularization, we incorporate it here in the image deblurring process.
The aforementioned image deblurring methods mostly employ the L1 norm to constrain the fidelity term. However, the L1 norm is a convex relaxation of the L0 norm, with limited ability to induce sparsity [13]. Recently, Woodworth and Chartrand [14] proved that Lp-pseudo-norm shrinkage is a method with better sparsity compared with soft-threshold shrinkage, which attracted extensive attention in academia. Zheng et al. [15] theoretically compared the advantages of Lp-minimization and L1-minimization, and tested them in a Gaussian noisy setting. Zhang et al. [16] applied Lp shrinkage [17] to the reconstruction of tomographic images to obtain an effect superior to the L1-norm constraint. Chen et al. [18] applied an Lp-pseudo-norm constraint to the sparse time-frequency representation of the seismic signal spectrum to obtain a time-frequency spectrum with higher time-frequency resolution. Wong et al. [19] proposed a model combining an Lp-pseudo-norm fidelity term and OGS regularization, which uses neighborhood structure similarity of image gradients for image denoising.
Considering the superiority of Lp-pseudo-norm and TGV, here, we employ Lp-pseudo-norm to describe the sparsity of impulse noise and adopt the TGV regularization term to fit the image sparsity. On the one hand, the Lp-pseudo-norm constraint increases the degrees of freedom of the fidelity term and better represents the sparsity of the impulse noise. In this way, the deblurring process is more robust to noise pollution. On the other hand, the TGV regularization term simultaneously puts constraints on the information of the first-order and second-order gradients of the image, thereby effectively suppressing the staircase effect.
To solve the problem, the alternating direction method of multipliers (ADMM) [20,21] is employed to decompose the complex problem into a number of sub-problems that are relatively easy to solve. At the same time, to improve the efficiency of the algorithm, it is assumed here that the image satisfies a periodic boundary condition. Following a fast deconvolution method [3,22,23], fast Fourier transform (FFT) is employed based on the convolution law to transform image differences in the time domain to differences in the frequency domain to avoid large matrix multiplication operations. In numerical experiments, the proposed new method is compared with FTVd, TGV, and OGSTV in terms of the following objective indicators: peak signal-to-noise ratio (PSNR), structural similarity (SSIM), signal-to-noise ratio (SNR), and computation time.
The rest of this paper is organized as follows: Section 2 presents prerequisite knowledge related to the algorithms used in this study, including the modeling of the TGV regularization term and the Lp Algorithms 2019, 12, 221 3 of 24 norm. Section 3 firstly presents the model proposed in this study and then elaborates on an efficient method for solving the model in the ADMM framework. Results of the numerical experiment are presented in Section 4, and conclusions are given in Section 5.

FTVd Model
We assumed that the linear mathematical model of an observed degraded image is expressed as follows [24]: For the sake of simplicity, it was assumed, in this study, that the processed image is an N × N square matrix, which can be easily generalized to an image of size M × N. G ∈ R N×N is a degraded image matrix, F ∈ R N×N represents an original unknown clear image matrix, H is a blur kernel matrix, * represents the two-dimensional convolution operator, and N ∈ R N×N represents additive noise, such as Gaussian noise and impulse noise. For impulse noise, a discrete ROF model is constructed using the L1 norm, where the FTVd model is expressed as where · 1 represents the Euclidean L1 norm. In Equation (2), H * F − G 1 is the fidelity term using the L1-norm, and R(F) is the regularization. In FTVd mode, 1] T are the horizontal and vertical convolution operators, respectively, and µ is a regularization parameter higher than zero, which serves to balance the fidelity term and the regularization term. As is shown in Equation (1), FTVd employs the ATV as the regularization and the L1 norm to depict the sparsity of impulse noise. However, ATV assumes the processed image to be piecewise constant. As a result, staircase artefacts occur. In addition, the L1 norm is the convex relaxation of the L0 norm. It cannot express the sparsity of impulse noise precisely. Considering the disadvantages of FTVd, the TGV regularization is employed to replace the ATV regularization, and the Lp-pseudo-norm is adopted to depict the statistic characteristics of impulse noise.
where α 0 and α 1 serve as the balance coefficients between the first-order total variation and V h , V v ∈ R N×N represent the horizontal and vertical gradients of the processed image.
Since the TGV regularization puts sparse constraints on not only the first-order gradients but also the second-order gradients, the staircase artefacts are relieved. Thus, we adopt the TGV regularization instead of ATV regularization in Equation (2) to further improve the quality of the restored image.

Lp-Pseudo-Norm
This study focused on exploring image restoration in the presence of impulse noise pollution. An impulse noise is often the noise of alternating bright (white) and dark (black) dots resulting from a number of factors, such as a faulty pixel in the imaging sensor, a wrong memory location in the hardware, a transmission channel noise, and a decoding process; it is considered an additive noise. In recent years, with the introduction of optimization methods, sparsity was widely applied to image restoration. However, the non-convexity of the L0 norm makes solution optimization (NP-hard Algorithms 2019, 12, 221 4 of 24 problem) challenging. In 2004, Nikolova [26] discovered that, when the noise is an impulse noise, the use of an L1 norm for constructing a fidelity term model allows better removal of outliers. Later, many studies used an L1-norm fidelity term for image denoising related to impulse noise or outliers. Image deconvolution in the presence of impulse noise is often treated as solving a minimization problem involving an L1-norm fidelity term and a regularization function. In fact, the sparsity of the L1 norm, which is the convex relaxation of the L0 norm, is not as good as that of the L0 norm. In 2016, Woodworth and Chartrand [14] pointed out that the Lp norm is a better convex approximation of the L0 norm than the L1 norm is, and proved that Lp shrinkage has better sparsity than soft-threshold shrinkage. The Lp norm is defined as To better elaborate on this issue, the contours of the Lp-pseudo-norm are shown in Figure 1. Figure 2 presents a group of diagrams, where the contours represent the fidelity terms for image restoration, while the dashed lines represent the regularization terms. As shown in the figure, it is more likely for the regularization term to intersect with the fidelity term of p < 1 than that of p = 1 on the coordinate axis, thereby making the likelihood of obtaining a sparse solution greater. Therefore, by using the Lp-pseudo-norm as a constraint, it is possible to better extract the sparse characteristics of signals in the transform domain to more precisely restore the image in the optimization problem. restoration. However, the non-convexity of the L0 norm makes solution optimization (NP-hard problem) challenging. In 2004, Nikolova [26] discovered that, when the noise is an impulse noise, the use of an L1 norm for constructing a fidelity term model allows better removal of outliers. Later, many studies used an L1-norm fidelity term for image denoising related to impulse noise or outliers. Image deconvolution in the presence of impulse noise is often treated as solving a minimization problem involving an L1-norm fidelity term and a regularization function. In fact, the sparsity of the L1 norm, which is the convex relaxation of the L0 norm, is not as good as that of the L0 norm. In 2016, Woodworth and Chartrand [14] pointed out that the Lp norm is a better convex approximation of the L0 norm than the L1 norm is, and proved that Lp shrinkage has better sparsity than soft-threshold shrinkage. The Lp norm is defined as To better elaborate on this issue, the contours of the Lp-pseudo-norm are shown in Figure 1. Figure 2 presents a group of diagrams, where the contours represent the fidelity terms for image restoration, while the dashed lines represent the regularization terms. As shown in the figure, it is more likely for the regularization term to intersect with the fidelity term of 1 p  than that of 1 p  on the coordinate axis, thereby making the likelihood of obtaining a sparse solution greater. Therefore, by using the Lp-pseudo-norm as a constraint, it is possible to better extract the sparse characteristics of signals in the transform domain to more precisely restore the image in the optimization problem.

Proposed Model
To use the FFT to achieve a faster algorithm computation speed, it is assumed that the image satisfies a periodic boundary condition. Moreover, the model is expressed as a matrix, and the  restoration. However, the non-convexity of the L0 norm makes solution optimization (NP-hard problem) challenging. In 2004, Nikolova [26] discovered that, when the noise is an impulse noise, the use of an L1 norm for constructing a fidelity term model allows better removal of outliers. Later, many studies used an L1-norm fidelity term for image denoising related to impulse noise or outliers. Image deconvolution in the presence of impulse noise is often treated as solving a minimization problem involving an L1-norm fidelity term and a regularization function. In fact, the sparsity of the L1 norm, which is the convex relaxation of the L0 norm, is not as good as that of the L0 norm. In 2016, Woodworth and Chartrand [14] pointed out that the Lp norm is a better convex approximation of the L0 norm than the L1 norm is, and proved that Lp shrinkage has better sparsity than soft-threshold shrinkage. The Lp norm is defined as To better elaborate on this issue, the contours of the Lp-pseudo-norm are shown in Figure 1. Figure 2 presents a group of diagrams, where the contours represent the fidelity terms for image restoration, while the dashed lines represent the regularization terms. As shown in the figure, it is more likely for the regularization term to intersect with the fidelity term of 1 p  than that of 1 p  on the coordinate axis, thereby making the likelihood of obtaining a sparse solution greater. Therefore, by using the Lp-pseudo-norm as a constraint, it is possible to better extract the sparse characteristics of signals in the transform domain to more precisely restore the image in the optimization problem.

Proposed Model
To use the FFT to achieve a faster algorithm computation speed, it is assumed that the image satisfies a periodic boundary condition. Moreover, the model is expressed as a matrix, and the

Proposed Model
To use the FFT to achieve a faster algorithm computation speed, it is assumed that the image satisfies a periodic boundary condition. Moreover, the model is expressed as a matrix, and the gradient difference is expressed in the form of a convolution to avoid large-scale matrix multiplication operations, thereby facilitating the subsequent discussion. Combining the Lp-pseudo-norm with the TGV regularization, this study proposes a new regularization model for image restoration to more thoroughly extract the structural characteristics of the first-order and second-order gradient matrices of the image, and obtain a sparser solution. The model is named the TGV-Lp model, described below.
where µ is the regularization parameter, · p p represents the Lp-pseudo-norm, α 0 and α 1 serve as the balance coefficients, and V h , V v represent the horizontal and vertical gradients of the processed image.

Solution Method
To solve the model described by Equation (4), the complex optimization problem is decomposed into several sub-problems that are solved using the ADMM framework. Then, alternating iterations of the solutions of the sub-problems are performed until the solutions converge to the optimal solution of the original optimization problem. To this end, a series of intermediate variables X i are introduced, with X i ∈ R N×N (i = 0 ∼ 5), and they are defined as The original problem is now transformed into the following constraint problem: The augmented Lagrangian multiplier (ALM) method [20] is used to transform the constraint problem described by Equation (6) into an unconstrained augmented Lagrangian function as follows: Algorithms 2019, 12, 221 6 of 24 where µ 0 = µα 0 and µ 1 = µα 1 ; the dual variable Λ i ∈ R N×N (i = 0 ∼ 5) is introduced as a Lagrangian multiplier, β j > 0 ( j = 0 ∼ 2) is the penalty parameter, and X,Y represents the inner product of two matrices X and Y. Upon adding the term to Equation (7) to complete the formula, Equation (8) is obtained after rearrangement.

Solving Sub-Problems F, V h , and V v
According to the ADMM principle, F, V h , and V v are decoupled with the intermediate variables and the dual variables, while there is coupling between F, V h , and V v . It is necessary to establish Equation (9) to obtain a simultaneous solution.
It is assumed here that the image to process satisfies a periodic boundary condition. According to the convolution law, the convolution result of two matrix spaces corresponds to the point multiplication of two matrix spectra in the frequency domain. To effectively avoid the computational complexity introduced by large-scale matrix operations, FFT is used for convolutional computation to solve the Algorithms 2019, 12, 221 7 of 24 problem described by Equation (9). To this end, F, V h , and V v in Equation (9) are re-expressed in the frequency domain as follows: where X denotes the frequency spectrum of X, and • denotes the element-wise matrix multiplication operator. The right side of the equation is a smooth convex function. To seek the optimal solution of Equation (10), the partial differentiation of J F , J V h , J V v with respect to F, V h , V v is set to zero as follows: For the sake of reading convenience, Equation (12) is re-formulated into simultaneous equations of the three variables F, V h , and V v .
By re-organizing Equation (11), we obtain Algorithms 2019, 12, 221 where 1 denotes a matrix of ones. Now, F, V h , and V v can be solved using FFT based on Cramer's rule.
In Equation (16), F −1 2D denotes a two-dimensional inverse Fourier transform, and the division symbol / denotes element-wise matrix division. According to Cramer's rule, || * in Equation (16) refers to the following:

Solving Sub-Problems Containing Dual Variables
For the sub-problems containing dual variables Λ i (i = 0, 1, 2, · · · , 5), the objective functions are expressed as In the ADMM framework, dual variables can be updated using the ascending gradient method as follows: where γ denotes the learning rate or the slack variable. The convergence of the algorithm was proven in the case of γ ∈ (0, √ 5+1 2 ) [21]. Hitherto, all the sub-problems of the model proposed in this study were solved, and the algorithm is herein referred to as TGV_Lp, which is summarized in Algorithm 1.

Numerical Experiments
To verify the appropriateness and validity of the proposed model, standard images of different styles were selected as experimental objects. The images were degraded with Gaussian blur, average blur, motion blur, and different degrees of impulse noise so that a large number of numerical simulation experiments could be performed to make indicator-based comparisons and evaluations in an objective manner. In particular, the image "Satellite" was downloaded from the Emory University image database [27], and the other images were downloaded from CVG-UGR (Computer Vision Group, University of Granada) image database [28]. These test images were widely used in other relevant publications, thereby allowing for objective and fair comparison of the test results. All calculations in this study were performed using MATLAB R2018a on a Windows 10 system with an Intel ®Core ™ i7-7700 central processing unit (CPU) (4.2 GHz) with 16.0 GB of memory.

Evaluation Indicators and Stopping Criteria
For objective evaluation of image quality, PSNR, SSIM, and SNR were used. Their definitions are given below [29].
where X represents an original image, and Y represents a restored image; L represents the maximum gray level of the image, which was set as 255 in this study. Generally, a larger PSNR value of a restored image denotes less image distortion and greater similarity between the restored image and the original image.
where µ X and µ Y represent the means of images X and Y, respectively, σ 2 X and σ 2 Y represent the variances of X and Y, respectively, σ XY represents the covariance of X and Y, and the constants k 1 and k 2 are very small positive numbers introduced to prevent the denominator of Equation (22) from being zero; they were set as 0.01 and 0.03 in this study, respectively. The maximum value of SSIM is 1, which indicates that the original image and the restored image are 100% identical; as the value of SSIM approaches 1, the similarity between the original and restored images increases.
Images with low noise have higher SNR values, making SNR an important metric for characterizing the performance of the restoration algorithm.
In the experiments, the iteration stopping criterion for comparing various algorithms was Equation (24), where F (k) and F (k+1) denote the objective functions after the k-th iteration and (k + 1)-th iteration, respectively.
The blur kernels used in the test were Gaussian blur, average blur, and motion blur, all of which were generated by built-in functions of MATLAB. The function fspecial ('gaussian', S,σ) generates an S × S Gaussian blur kernel standard deviation of σ, and the kernel is hereinafter referred to as (G, S, σ) for convenience; fspecial ('average', S) generates an S × S average blur kernel, which is referred to as (A, S); fspecial ('motion', L,θ) generates a motion blur kernel with a length of L and angle of θ, which is referred to as (M, L, θ). Noise effects were generated by built-in function imnoise (Ib, 'salt & pepper', level), where "Ib" is the blurred image and "level" is the noise density, i.e., the percentage of noise in the total number of pixels.

Parameter Selection and Sensitivity Analysis
Parameters involved in the algorithm of this study are µ 0 , µ 1 , β i (i = 0, 1, 2), p, and γ. In particular, µ 0 = µα 0 and µ 1 = µα 1 are regularization parameters with values greater than zero, which serve to balance the fidelity term and regularization term. Based on our previous experience with testing, we assumed here that α 0 = 2α 1 while gradually increasing the value of µ between 0.1 and 10 to find the optimal effect of image restoration. Given the characteristics of impulse noise, the value of the Lp-pseudo-norm was gradually increased from 0.1 to 1 by a step of 0.05, and a loop statement was adopted to traverse all paths to find the p value for the optimal PSNR. In all tests, the penalty parameter β i (i = 0, 1, 2) of the ALM method was set to have ratios of β 0 : β 1 : β 2 = 50 : 1 : 5, and the slack variable was set to 1, i.e., γ = 1.

Regularization Parameter µ
To test the effect of the regularization parameter µ on the image restoration results, different levels of blur and impulse noise were introduced into the test images (a)-(i), and PSNR values were recorded by gradually increasing the value of µ in the range of 0.1-10 while fixing other parameters. The experimental results show that the proposed algorithm is not sensitive to µ; that is, during the recovery process of different images, the value of PSNR does not undergo significant changes with changes in µ, indicating that the proposed algorithm is robust. The test results of the "Boat" image of Figure 3d and "Hill" image of Figure 3f treated with the Gaussian blur kernel at 30%, 40%, 50%, and 60% noise levels are presented in Figure 4. The PSNR value in Figure 4 was not optimal for the proposed algorithm, as other fixed parameters remained to be adjusted optimally. In subsequent experiments, the regularization parameter was set to 1, i.e.,µ = 1.
in an objective manner. In particular, the image "Satellite" was downloaded from the Emory University image database [27], and the other images were downloaded from CVG-UGR (Computer Vision Group, University of Granada) image database [28]. These test images were widely used in other relevant publications, thereby allowing for objective and fair comparison of the test results. All calculations in this study were performed using MATLAB R2018a on a Windows 10 system with an Intel ® Core ™ i7-7700 central processing unit (CPU) (

Evaluation Indicators and Stopping Criteria
For objective evaluation of image quality, PSNR, SSIM, and SNR were used. Their definitions are given below[29].
where X represents an original image, and Y represents a restored image; L represents the maximum gray level of the image, which was set as 255 in this study. Generally, a larger PSNR value of a restored image denotes less image distortion and greater similarity between the restored image and the original image.
where  X and  Y represent the means of images X and Y, respectively, 2  X and 2  Y represent the variances of X and Y, respectively,  XY represents the covariance of X and Y, and the constants 1 k and 2 k are very small positive numbers introduced to prevent the denominator of Equation (22) from being zero; they were set as 0.01 and 0.03 in this study, respectively. The maximum value of SSIM is 1, which indicates that the original image and the restored image are 100% identical; as the value of SSIM approaches 1, the similarity between the original and restored images increases.  imnoise (Ib, 'salt & pepper', level), where "Ib" is the blurred image and "level" is the noise density, i.e., the percentage of noise in the total number of pixels.

Parameter Selection and Sensitivity Analysis
Parameters involved in the algorithm of this study are 0 , 1 , ( = 0,1,2), , and . In particular, 00    and 11    are regularization parameters with values greater than zero, which serve to balance the fidelity term and regularization term. Based on our previous experience with testing, we assumed here that 0 1 =2  while gradually increasing the value of  between 0.1 and 10 to find the optimal effect of image restoration. Given the characteristics of impulse noise, the value of the Lp-pseudo-norm was gradually increased from 0.1 to 1 by a step of 0.05, and a loop statement was adopted to traverse all paths to find the p value for the optimal PSNR. In all tests, the penalty parameter The experimental results show that the proposed algorithm is not sensitive to  ; that is, during the recovery process of different images, the value of PSNR does not undergo significant changes with changes in  , indicating that the proposed algorithm is robust. The test results of the "Boat" image of Figure 3d and "Hill" image of Figure 3f treated with the Gaussian blur kernel at 30%, 40%, 50%, and 60% noise levels are presented in Figure 4. The PSNR value in Figure 4 was not optimal for the proposed algorithm, as other fixed parameters remained to be adjusted optimally. In subsequent experiments, the regularization parameter was set to 1, i.e., =1  .

Value of p in Lp-Pseudo-Norm
When it came to selecting the p value of the Lp-pseudo-norm, given the characteristics of the impulse noise, the p value was gradually increased from 0.1 to 1 in steps of 0.05, and a loop statement was adopted to traverse all paths so as to find the p value that could give the optimal PSNR value of image restoration. For restoration of the nine images treated with ( ,7,5) G at 30-60% noise levels, the best p values were obtained and are listed in Table 1, with the results showing that the p value was relatively stable for the same image. (a) Effect of µ on peak signal-to-noise ratio (PSNR) at noise levels 30%, 40%, 50%, and 60% with image "Boat". (b) Effect of µ on PSNR at noise levels 30%, 40%, 50%, and 60% with image "Hill".

Value of p in Lp-Pseudo-Norm
When it came to selecting the p value of the Lp-pseudo-norm, given the characteristics of the impulse noise, the p value was gradually increased from 0.1 to 1 in steps of 0.05, and a loop statement was adopted to traverse all paths so as to find the p value that could give the optimal PSNR value of image restoration. For restoration of the nine images treated with (G, 7, 5) at 30-60% noise levels, the best p values were obtained and are listed in Table 1, with the results showing that the p value was relatively stable for the same image. In addition, to better show the best selection of the p and µ values, we used a three-dimensional (3D) figure. The test results of the "Hill" image treated with the Gaussian blur kernel at 30% noise level are presented in Figure 5.  In addition, to better show the best selection of the p and  values, we used a three-dimensional (3D) figure. The test results of the "Hill" image treated with the Gaussian blur kernel at 30% noise level are presented in Figure 5.

Numerical Performance Comparison
To test the effectiveness of the proposed algorithm, the following four algorithms were compared: an FTVd model, a TGV model, an OGSTV model, and the proposed frequency-domain TGV-Lp model. The algorithms were adapted from the authors' web pages or shared by the authors. It is noteworthy that all these algorithms are state-of-the-art algorithms in the field of image restoration. In the experiments, regularization parameters were adjusted until each algorithm achieved its best performance for image restoration to ensure an objective and fair comparison.

Numerical Performance Comparison
To test the effectiveness of the proposed algorithm, the following four algorithms were compared: an FTVd model, a TGV model, an OGSTV model, and the proposed frequency-domain TGV-Lp model. The algorithms were adapted from the authors' web pages or shared by the authors. It is noteworthy that all these algorithms are state-of-the-art algorithms in the field of image restoration. In the experiments, regularization parameters were adjusted until each algorithm achieved its best performance for image restoration to ensure an objective and fair comparison.

Numerical Performance Comparison for Image Recovery in the Presence of Gaussian Blur
The results of the image restoration tests performed on nine test images degraded with Gaussian blur kernel (G, 7, 5) and 30-60% impulse noise are listed in Table 2. In the table, the optimal indicators are highlighted in bold. Table 2. Numerical results of PSNR (dB), structural similarity (SSIM), and SNR (dB) for the recovery of nine test images degraded with a Gaussian blur kernel (G, 7,5) and 30-60% impulse noise. FTVd-fast total variation deconvolution; TGV-total generalized variation; OGSTV-overlapping group sparsity total variation. The results of image restoration tests performed on nine test images degraded with a Gaussian blur kernel (G, 15,5) and 30% noise are shown in Table 3. Table 3. Numerical results of PSNR (dB), SSIM, and SNR (dB) for the recovery of nine test images degraded with a Gaussian blur kernel (G, 15,5) and 30% impulse noise. An in-depth analysis of the data in Tables 2 and 3 leads to the following conclusions:
All four methods can effectively recover images that were degraded by Gaussian blur and different degrees of impulse noise. The PSNR, SSIM, and SNR values of the proposed method are generally higher than those of other competitive methods, which indicates that the proposed method has better deblurring and denoising effects. In particular, for images with many lines and complex textures, such as "Fingerprint" in Figure 3a, "Butterfly" in Figure 3b, and "Baboon" in Figure 3c, the proposed method performed well, especially at high noise levels. For example, for the "Fingerprint" image of Figure 3a treated with a Gaussian blur kernel (G, 7, 5) and 50% noise, the PSNR value (35.59 dB) of the proposed method was higher by 7.22 dB than that of the FTVd method (28.37 dB), and was higher by 2.55 dB and 1.6 dB than that of the TGV method (33.04 dB) and the OGSTV method (33.99 dB), respectively.

2.
When restoring the nine degraded images treated with the Gaussian blur kernel (G, 7, 5) and 30-60% impulse noise, the proposed method achieved 0.58-2.85 dB higher PSNR values compared with the TGV method. When restoring the nine degraded images treated with the Gaussian blur kernel (G, 15,5) and impulse noise, the proposed method achieved 0.45-1.73 dB higher PSNR values compared with the TGV method. For example, for the "Bird" image of Figure 3e treated with the Gaussian blur kernel (G, 7, 5) and two different levels (50% vs. 60%) of noise, the PSNR values of the proposed method (32.31 dB and 29.69 dB) were 2.85 dB and 2.7 dB higher than those of the TGV method (29.46 dB and 26.99 dB), respectively.

3.
The OGSTV method is an excellent algorithm for image deblurring and impulse noise removal, and it uses the neighborhood gradient information of the image to form a combined gradient for image restoration, which achieves good image restoration performance for smooth image regions but is less satisfactory in recovering image edge regions containing sharp lines and angles. In the test, it was observed that, for smooth images with non-complex textures, such as the "Plane" image of Figure 3h and "Satellite" image of Figure 3i, the OGSTV method had particularly excellent performance, which was slightly better than that of the proposed algorithm in the case of low noise levels. In the case of high noise levels, however, the proposed algorithm was still superior to the OGSTV method. For example, when restoring the "Satellite" image degraded with (G, 7, 5) and 30% noise, the OGSTV method achieved a slightly higher PSNR value of 36.97 dB compared to the proposed method (36.54 dB), with a difference of only 0.43 dB; however, in the cases of 50% and 60% noise, the PSNR values of the OGSTV method were lower than those of the proposed method.
The above analysis indicates that the use of Lp-pseudo-norm shrinkage in the proposed method increases the degree of freedom by one, making it less difficult to describe image gradient sparsity. The performance of the proposed method in removing Gaussian blur and impulse noise is generally better than that of the three competitive methods, especially for images with complex textures. In addition, the comparison reveals that the image restoration performance of the FTVd algorithm is much worse than that of the other three algorithms; thus, the subsequent experimental results of the FTVd algorithm are not listed in some tables.

Numerical Performance Comparison for Image Recovery in the Presence of Average Blur
The image restoration capability of the proposed method for degraded images treated with average blur and impulse noise was investigated, with Table 4 presenting the results of image restoration tests performed on nine degraded test images treated with an average blur kernel (A, 7) and two levels (30% vs. 60%) of noise. With regard to the image restoration performance of the algorithms for different images treated with average blur and different levels of impulse noise, the data in Table 4 lead to conclusions similar to those obtained in the case of Gaussian blur and impulse noise. When restoring the nine degraded test images treated with an average blur kernel (A, 7) and impulse noise, the proposed method achieved 0.48-4.16 dB higher PSNR values compared with the TGV method. For images with many lines and complex textures, such as the "Fingerprint" image of Figure 3a, "Butterfly" image of Figure 3b, and "Baboon" image of Figure 3c, the proposed method achieved 0.8-2.14 dB higher PSNR values compared with the OGSTV method. However, the OGSTV performed well in restoring smooth images in the presence of low noise levels.

Numerical Performance Comparison for Image Recovery in the Presence of Motion Blur
Furthermore, the image restoration capability of the proposed method in the case of motion blur and impulse noise was investigated. Table 5 presents the PSNR results of image restoration tests performed on the degraded "Butterfly" image of Figure 3b, "Baboon" image of Figure 3c, "Bird" image of Figure 3e, and "Truck" image of Figure 3g treated with different motion blur kernels and 30-60% noise. For different images treated with varying degrees of motion blur and impulse noise, the data in Table 5 reveal that the proposed method outperformed the TGV and OGSTV methods in removing motion blur and impulse noise. When restoring the four images, the proposed method achieved 0.52-5.50 dB and 0.21-2.83 dB higher PSNR values than the TGV method and the OGSTV method, respectively.

Comparison of Visual Effects
To better observe the test results, some restored images are presented below for visual comparison purposes.

Comparison of the Visual Effects of Restored Images in the Case of Gaussian Blur
Figures 6 and 7 present image restorations by the four algorithms of the "Boat" image degraded with a Gaussian blur kernel (G, 7, 5) and 30% impulse noise, and zoomed-in images of the restored images for better comparison, respectively. With respect to the visual effects of the restored images, with the FTVd method (c), a loss of some image details occurred during restoration of the image, and there was an obvious staircase effect. The TGV method (d) and the OGSTV method (e) performed well in suppressing the staircase effect, but slight ringing artefacts were generated when the TGV method was used to restore the image. Furthermore, the OGSTV method performed well in smooth image regions but relatively poorly preserved image edges. The proposed method, however, not only preserved image edges, but also restored smooth image regions well, achieving the best overall visual effects for image restoration.  and 50% impulse noise, as well as the zoomed-in restored images for better comparison. With respect to the visual effects of the restored images, the FTVd method (c) led to a certain degree of staircase effect in the restored image, the TGV method (d) led to blurred image edges in the restored image, and the OGSTV method (e) needed to be improved in preserving image edges. In contrast, the proposed method performed well in restoring the details of the image lines.  and 50% impulse noise, as well as the zoomed-in restored images for better comparison. With respect to the visual effects of the restored images, the FTVd method (c) led to a certain degree of staircase effect in the restored image, the TGV method (d) led to blurred image edges in the restored image, and the OGSTV method (e) needed to be improved in preserving image edges. In contrast, the proposed method performed well in restoring the details of the image lines.  7) and 50% impulse noise, as well as the zoomed-in restored images for better comparison. With respect to the visual effects of the restored images, the FTVd method (c) led to a certain degree of staircase effect in the restored image, the TGV method (d) led to blurred image edges in the restored image, and the OGSTV method (e) needed to be improved in preserving image edges. In contrast, the proposed method performed well in restoring the details of the image lines.

Comparison of Computing Time
Finally, to further test the timeliness of the proposed model, the proposed method was compared with the FTVd, TGV, and OGSTV methods by testing the mean time required to restore

Comparison of Computing Time
Finally, to further test the timeliness of the proposed model, the proposed method was compared with the FTVd, TGV, and OGSTV methods by testing the mean time required to restore the nine images of Figure 3 for different blur kernels at various impulse noise levels. Some test results are shown in Table 6, where the best indicators are highlighted in bold. The calculation times in the table reveal that FTVd, which is a classic algorithm, is still advantageous in terms of computing speed. The proposed method is comparable to the OGSTV method in terms of image restoration time, while they are both faster than the TGV method. The proposed method uses the Lp-pseudo-norm as a regularization constraint, which increases the method's computational burden, thereby leading to a certain degree of increase in the image restoration time. However, the gradient difference is re-expressed in the form of a convolution in this study, where it corresponds to the dot multiplication of two matrix spectra in the frequency domain; moreover, FFT is employed for convolution calculation. These strategies effectively avoid large-scale matrix multiplications, thus speeding up the calculation.

Conclusions
In this study, a regularized method for image restoration based on the Lp-pseudo-norm and TGV was proposed to de-blur images and remove impulse noise. The method constrains the first-order and second-order gradients of the image and takes advantage of the sparse representation ability of the Lp-pseudo-norm for image gradients. In this way, it effectively removes image blur and noise, and reduces undesirable artefacts (such as staircase and ringing artefacts) in image restoration while better preserving the edges and details of the image. A large number of experimental results revealed that the proposed method is consistently superior to several state-of-the-art methods in terms of numerical indicators and visual effects. The regularization parameters µ of the Lp-pseudo-norm used in the experiments were relatively invariant, indicating that the proposed algorithm is robust. In follow-up studies, efforts will be made to further improve the computational efficiency of the proposed method. It should be noted that the proposed algorithm uses a generalized regularization term, which can be easily extended to other algorithms and is suitable for other image restoration problems.