An Efficient Image Reconstruction Framework Using Total Variation Regularization with Lp-Quasinorm and Group Gradient Sparsity

The total variation (TV) regularization-based methods are proven to be effective in removing random noise. However, these solutions usually have staircase effects. This paper proposes a new image reconstruction method based on TV regularization with Lp-quasinorm and group gradient sparsity. In this method, the regularization term of the group gradient sparsity can retrieve the neighborhood information of an image gradient, and the Lp-quasinorm constraint can characterize the sparsity of the image gradient. The method can effectively deblur images and remove impulse noise to well preserve image edge information and reduce the staircase effect. To improve the image recovery efficiency, a Fast Fourier Transform (FFT) is introduced to effectively avoid large matrix multiplication operations. Moreover, by introducing accelerated alternating direction method of multipliers (ADMM) in the method to allow for a fast restart of the optimization process, this method can run faster. In numerical experiments on standard test images sourced form Emory University and CVG-UGR (Computer Vision Group, University of Granada) image database, the advantage of the new method is verified by comparing it with existing advanced TV-based methods in terms of peak signal-to-noise ratio (PSNR), structural similarity (SSIM), and operational time.


Introduction
Owing to the influence of sensing equipment, imaging systems, environmental conditions, and human factors, digital images are usually subject to a certain degree of degradation such as image blur, image noise, and partial image information loss.It is widely assumed that the convolution of a clear original image with a blur kernel gives a blurred image, which will further turn into a degraded image when interfered with by noise.Common types of blur kernels are Gaussian blur kernel, motion blur kernel, and average blur kernel [1].According to the distribution condition of probability density function (PDF) of their amplitude, noises are classified into Gaussian noise, impulse noise, Rayleigh noise, exponential noise, gamma noise, etc. [1].Designing image optimization models and proposing efficient algorithms by using degraded images and some prior information properly to reconstruct clear images is of great importance to the advanced postprocessing of images.
Reconstructing a degraded image to the original image is an inverse problem.Seeking solutions to an inverse problem is ill-posed in the sense that there is a lack of any solution or a lack of a unique, stable solution [2,3].To deal with such a type of ill-posed problem, image reconstruction based on regularization is an effective method.In a regularization model, the image reconstruction function consists of two parts: a fidelity term and a regularization term.The fidelity term can be obtained from a priori information of the image.Given that an imaging process may be affected by different types of degradation factors, the fidelity terms may be modeled differently.However, the use of a different regularization term in a regularization function has a large effect on the quality of reconstructed images.The relative significance of fidelity and regularization terms can be adjusted by regularization parameters.
The selection of a regularization function is of great research interest in the field of image processing, and many effective models have been proposed by scholars for different applications [4][5][6][7].Phillips [8] and Tikhonov [9] proposed smooth regularization terms, which were derived from the L2-norm and referred to as the Tikhonov regularization method.As a classic regularization method, the Tikhonov regularization method has the advantage of calculation simplicity but usually leads to oversmoothing of the edges of images.In 1992, Rudin, Osher, and Fatemi [10] pioneered the total variation (TV) regularization method, also known as the ROF model, which can well preserve image edge features and has attracted widespread attention in the image denoising field.
The TV models are divided into the anisotropic total variation (ATV) model and the isotropic total variation (ITV) model [11][12][13].In first-order TV models, the image is piecewise smooth, which makes the models obviously advantageous in preserving image edges, but the models are likely to produce a staircase effect.Meanwhile, given the nondifferentiability of total variation functionals, the solution is difficult to find.Since then, some extended models and their algorithms have been proposed on the basis of the TV models and have been widely used in image reconstruction [14][15][16], photoelectric detection [17], geological exploration [18], remote sensing [19], and medical image analysis [20], among others [21].For example, the total generalized variation (TGV) model can effectively approximate a polynomial function of arbitrary order [22][23][24].
The TGV model has proven, theoretically and practically, to be a more effective method to remove the staircase effect [23], and has attracted wide attention from scholars.For example, Knoll and Bredies et al. [25] applied TGV regularization terms to magnetic resonance imaging, and Kong and Peng et al. [26] applied TGV regularization terms to the denoising of seismic signals.Nonlocal TV (NLTV) models accomplish the removal of the staircase effect by introducing local information of the image [27,28].Although NLTV models are superior to TV models in removing the staircase effect, their computational complexity is higher than that of traditional TV models.Fractional-order TV (FrTV) models [29] extend integer-order TV models to the fractional order, taking into account the nonlocal characteristics of the image while considering the local characteristics of the image, so that the staircase effect is effectively inhibited.In 2011, Wang et al. [30] proposed Fast Total Variation deconvolution (FTVd) to solve the TV_L1 problem for image deblurring.
In 2011, Sakurai et al. [31] proposed the four-directional total variation (4-TV) model to extend gradient information into four directions in contrast to the two vertical directions or horizontal directions considered by traditional variation methods, to improve the denoising performance.In recent years, Selesnick and Chen [32,33] proposed an overlapping group sparsity TV (OGSTV) model in which a regularization term based on overlapping group sparsity is introduced, and applied it with the L2-norm fidelity term to the denoising of one-dimensional Gaussian noise.The OGSTV regularization term considers the neighborhood information of an image gradient and retrieves the neighborhood structural similarity of the image gradient.Liu et al. [34] generalized this to a two-dimensional OGS_L1 regularization term and applied it with an L1-norm-constrained fidelity term to the removal of pulse noise from images, and Liu et al. [35] applied it to speckle noise removal.Their experiments showed that the regularization term can effectively suppress the staircase effect.However, the L1-norm, which is the convex relaxation of the L0-norm [36], has a limited ability to induce sparsity.Meanwhile, given the increased computational complexity owing to the neighborhood information of an image gradient, a longer computational time is required for images with complicated edge information.This study proposes a new regularization model that incorporates the Lp-quasinorm x p p is named as Lp-quasinorm for simplicity) in the fidelity term and the group gradient sparsity in the regularization term to denoise and deblur images containing impulse noise.Recently, Woodworth and Chartrande [37] proved that the Lp-quasinorm shrinkage is superior to the soft-thresholding shrinkage when recovering sparse signals.In this study, the regularization term of the group gradient sparsity is tentatively combined with an Lp-quasinorm constraint to overcome the drawback of traditional variation methods-namely, focusing on the gradient at a given pixel in the image without considering the neighborhood information.By retrieving the gradient information of the surrounding pixels as a reference to form group gradients.This decreases the group gradient between a single noisy pixel and its neighborhood while still allowing large group gradients between edge pixels and their neighborhood.
In this study, by setting a reasonable threshold, the edge information of the image is well preserved while the noise is efficiently removed.To find a solution, the alternating direction method of multipliers (ADMM) [38,39] and the majorization-minimization (MM) [40] algorithm are used in the proposed model.To improve the image reconstruction efficiency, a Fast Fourier Transform (FFT) is introduced to transform the image difference operation in the time domain into the frequency domain to effectively avoid large matrix multiplication operations.In particular, accelerated ADMM that allows a faster restart of the optimization process is further introduced in the proposed model, thus greatly improving the operational speed.In numerical experiments, the FTVd, FrTV, OGS_L1, and TGV models are compared in terms of objective performance indicators [such as the peak signal-to-noise ratio (PSNR) [41], structural similarity (SSIM) [42] and operational time] with the proposed model in this study.The results show that the proposed model achieves higher PSNR, SSIM and better visual effects for reconstructed images.The proposed model uses a generalized regularization term applicable to other image reconstruction problems, such as magnetic resonance imaging (MRI) process, optical coherence tomography (OCT) images noise reducing [43], etc.
The rest of this paper is organized as follows.In the next section, prerequisite knowledge related to the algorithms is presented, including the definition of the OGS_L1 regularization term and the MM algorithm.The third section presents the model proposed in this paper, and then details the problem-solving process in the ADMM framework, after which the section elaborates on an accelerated ADMM for faster restart of the optimization process and faster calculation speed.Numerical experimental results are given in the fourth section.Finally, the fifth section draws the conclusion.

L1-and L2-Norm Modeling in TV Regularization
The goal of image reconstruction is to estimate the original image by using the observed degraded image.An image degradation model is as follows [1]: G ∈ R N×N represents an observed blurred image with noise, F ∈ R N×N represents the image reconstructed by the model, H ∈ R N×N represents a blur kernel function, and the symbol * represents the convolution operator.N ∈ R N×N represents additive noise.When the noise is Gaussian, the TV regularization (ROF model) is expressed as [10] In Formula (2), arg min is the fidelity term using the L2-norm, and µR(F) is the sparsity-constrained regularization term, with µ as the regularization parameter measuring the weight of the fidelity term vs. the regularization term.
In the ATV model, R(F) is defined as follows [11]: In the ITV model, R(F) is defined as follows [2]: where T represent the horizontal and vertical differential convolution operators, respectively; • 1 represents the Euclidean L1-norm, and • 2 represents the Euclidean L2-norm.
The noise distribution might conform to a distribution other than the Gaussian distribution, e.g., the Laplace distribution, and in such case the fidelity term of the L1-norm will be used to replace that of the L2-norm.This paper is primarily intended to explore image reconstruction in the presence of impulse noise.The impulse noise is additive noise, which is mainly accounted for by the bright and dark noise generated by sensors, transmission channels, and decoding processing.In 2004, Nikolova proposed the use of a L1 data-fidelity term for impulse noise related problems [44].When the noise is impulse noise, consider the following L1-norm model:

Lp-Quasinorm
Noise that conforms to a Laplace distribution is usually modeled with the L1-norm.The L1-norm is the convex relaxation of the L0-norm [36].In recent years, Woodworth and Chartrand demonstrated that the Lp shrinkage is superior to the soft-thresholding shrinkage when recovering sparse signals [37], which has attracted widespread attention in academia [4,45].The Lp-norm is defined as , and the Lp-quasinorm is defined as The L1-norm and L2-norm are special cases of the more general Lp-norm at p = 1 and p = 2, respectively.Figure 1a shows the Lp-quasinorm contour, with p values of 0.25, 0.5, 1, 2, and 4 in the order of the inner to outer contour lines.Figure 1b shows a schematic diagram, assuming that the image is contaminated by impulse noise with standard deviation σ.The solid line in Figure 1b indicates the fidelity term for image reconstruction.Owing to noise interference, the fidelity term will fluctuate between the dashed lines.As seen in the figure, the fidelity term intersects more with the contours of 0 < p < 1 at the axis of higher sparsity.The Lp-quasinorm is more likely than the L1-norm to induce sparse solutions to the impulse noise removal problem.
of the inner to outer contour lines.Figure 1b shows a schematic diagram, assuming that the image is contaminated by impulse noise with standard deviation σ .The solid line in Figure 1b indicates the fidelity term for image reconstruction.Owing to noise interference, the fidelity term will fluctuate between the dashed lines.As seen in the figure, the fidelity term intersects more with the contours of 0 1 p < < at the axis of higher sparsity.The Lp-quasinorm is more likely than the L1-norm to induce sparse solutions to the impulse noise removal problem.

OGS_L1 Model
Based on the traditional ATV model, when the noise is impulse noise, the OGS_L1 model is expressed as follows [46,47]: ϕ(K h * F) and ϕ(K v * F) represent the horizontal and vertical group gradient, respectively.The regularization term considering the group gradient sparsity is expressed as If the information of K × K pixels of the image neighborhood is considered, where K is the size of the group, one can define where V ∈ R N×N represents the image gradient.The variable V i,j,K,K composed of K × K pixel matrix of the 2D image is defined as follows: where x is the floor function of x, which takes the value of the largest integer equal to or less than x; the element V i,j is at the center of the matrix V i,j,K,K , and the center is enclosed by the selected group.When K = 1, the OGS_L1 model degenerates into an ATV model.
As shown by Formula (8), the group gradient has fully considered the gradient information of the pixel neighborhood, and the regrouping of the gradient information through the two norms increases the difference between the smooth regions and the image edge regions.
To illustrate this problem in a more straightforward manner, a schematic diagram is given as shown in Figure 2, where Figure 2a depicts the vertical gradients in the smooth image regions with the two pixels 1 and 2 heavily contaminated by noise, while Figure 2b shows the vertical gradients in the image edge regions without noise contamination, with the pixel gradients shown for the two pixels 3 and 4 .The hollow circles in the figure denote pixels of low grayscale values, while the solid circles denote pixels of high grayscale values.For convenience of discussion, it is assumed that the grayscale values of the hollow circles are 0, and the grayscale values of the solid circles are 1.
Information 2019, 10, x 6 of 22  Given that the probability that a number of consecutive pixels in the smooth regions are all contaminated by noise is extremely small, this structural similarity can be used for denoising.
According to (8), when 3 K = , the group gradient of pixels ○ 1 and ○ 2 is 2 , while the group gradient of pixels ○ 3 and ○ 4 is 5 .As long as the threshold is set to 2 , the noise of pixels ○ 1 and ○ 2 can be well removed, while the group gradient of pixels ○ 3 and ○ 4 shrink to 5-2 .Thus, the image edge is well preserved.
In summary, by using group gradients to calculate gradients at pixels, it is possible to highlight the difference between the highly noisy pixels in the smooth regions and the pixels in the edge regions to achieve image reconstruction in a more robust manner.

MM Algorithm
To seek the optimal solution containing the group gradient sparsity ( )  ϕ V , the MM algorithm is adopted here [40].This is a progressive method focused on finding a multivariate auxiliary function to construct an iterative series to approximate the solution for a given problem.When solving the problem of ( 9), the iterative cycle is expressed as (10) [34] ( ) N N × ∈ I denotes the identity matrix, k denotes the number of iteration steps.For any is a diagonal matrix with the elements on the diagonal defined as [ ] In (11), the diagonal elements of D can be calculated using MATLAB's built-in function conv2.

Propose Model and Solution
In this section, a new image reconstruction model is proposed based on the Lp-quasinorm and the regularization term of group gradient sparsity as follows: As mentioned above, the regularization term of the group gradient sparsity is beneficial to If the conventional TV regularization methods are applied to Figure 2a, pixels 1 and 2 are very likely to be mistaken for and preserved as image edge pixels since they have grayscale values similar to the gradients at pixels 3 and 4 , so that the noise cannot be effectively removed.In conventional TV models, the four pixels have an identical horizontal difference and vertical difference, so that no matter how much the threshold is, the gradients of the four pixels show consistent variation trends.
Given that the probability that a number of consecutive pixels in the smooth regions are all contaminated by noise is extremely small, this structural similarity can be used for denoising.According to (8), when K = 3, the group gradient of pixels 1 and 2 is √ 2, while the group gradient of pixels 3 and 4 is √ 5.As long as the threshold is set to √ 2, the noise of pixels 1 and 2 can be well removed, while the group gradient of pixels 3 and 4 shrink to √ 5 − √ 2. Thus, the image edge is well preserved.
In summary, by using group gradients to calculate gradients at pixels, it is possible to highlight the difference between the highly noisy pixels in the smooth regions and the pixels in the edge regions to achieve image reconstruction in a more robust manner.

MM Algorithm
To seek the optimal solution containing the group gradient sparsity ϕ(V), the MM algorithm is adopted here [40].This is a progressive method focused on finding a multivariate auxiliary function to construct an iterative series to approximate the solution for a given problem.When solving the problem of ( 9), the iterative cycle is expressed as (10) [34] 2 denotes the identity matrix, k denotes the number of iteration steps.For any X ∈ R N×N , D ∈ R N 2 ×N 2 is a diagonal matrix with the elements on the diagonal defined as Information 2019, 10, 115 7 of 21 In (11), the diagonal elements of D can be calculated using MATLAB's built-in function conv2.

Propose Model and Solution
In this section, a new image reconstruction model is proposed based on the Lp-quasinorm and the regularization term of group gradient sparsity as follows: As mentioned above, the regularization term of the group gradient sparsity is beneficial to improving the difference between the smooth regions and the edge regions in the image after noise contamination.Meanwhile, with the ability of Lp-quasinorm in characterizing the image gradient sparsity, the proposed model is able to effectively preserve the image edge information and remove the staircase effect while blurring and denoising the image.
Next, a method for seeking a solution to the proposed model in the ADMM framework is introduced.For convenience, the following intermediate variables are introduced: . According to the ADMM framework, an augmented Lagrangian multiplier (ALM) method [38] is adopted to solve the constrained optimization problem of (12).By doing so, the problem is transformed into an unconstrained problem whose augmented Lagrangian objective function can be expressed as In ( 13), the dual variable ) is introduced as a Lagrangian multiplier, β i > 0(i = 1, 2, 3) is a penalty parameter which controls the linear constraint, and X, Y represents the inner product of two matrices X and Y.
It is difficult to solve the maximum value problem for the joint variable X 1 -X 3 .To solve this problem, the original optimization problem is decomposed into several smaller local subproblems in the framework of ADMM [38,48], that is, X 1 -X 3 are separately used in the solution-seeking process, followed by alternating iteration until the method converges to the optimal solution of the original model.The subproblems after the decomposition of ( 13) are as follows, where k denotes the number of iteration steps: Information 2019, 10, 115 8 of 21

Subproblem Solving
(1) Subproblem F When solving subproblem F of ( 14), to effectively avoid the computational complexity caused by large matrix multiplication operations, a fast two-dimensional Fourier transform is introduced to transform the operation of time-domain to that of frequency-domain.The frequency domain of the subproblem F is expressed as In ( 15), x denotes the frequency spectrum of x, and the symbol • denotes element-wise matrix multiplication.Given that the right side of the formula is a smooth convex function, its first-order difference result is set to zero to find the optimal solution [5], which is obtained by using an inverse two-dimensional Fourier transform, as shown below: where FFT −1 2D represents the two-dimensional inverse Fourier transform, and the division symbol represents element-wise matrix division.
(3) Subproblems X 2 and X 3 Subproblems X 2 and X 3 in ( 14) can be solved by using the following iterative loop based on the MM algorithm and (10) In (18), I ∈ R N 2 ×N 2 denotes the identity matrix, and ) denotes the n-th iteration of the MM algorithm in the k + 1-th outer loop.
Information 2019, 10, 115 9 of 21 Finally, according to the ADMM solution rule, the objective function of the dual variable This can be updated by the ascend gradient method as ( 20) where γ is a step length parameter, also known as a relax parameter.The convergence of the ADMM algorithm has been proven when γ ∈ (0, √ 2 ) [48].The entire algorithm proposed in this study is summarized as Algorithm 1 and is named GGS_Lp.

Fast ADMM Solver with Restart
The convergence rate of the traditional ADMM is O(1/k), while GGS_Lp involves calculation of the neighborhood gradient, which increases the computational burden.By introducing accelerated ADMM, the convergence rate can be increased to O(1/k 2 ) [49] to improve the computational efficiency.In this section, an accelerated ADMM framework is used to solve the proposed model at a faster convergence speed.To do this, it is necessary to introduce auxiliary variables X i (i = 1, 2, 3) and Λ i (i = 1, 2, 3), the acceleration step ε i (i = 1, 2, 3), and the sum of primal-dual residuals c i (i = 1, 2, 3).In the accelerated ADMM, the sum of the primal-dual residuals is monitored.With an increase in the sum, the iterative result is considered unable to converge, and thus the process is restarted to ensure high computational efficiency.
In the accelerated ADMM framework, X i (i = 1, 2, 3) is updated as where X 2 2 , and X 3 The sum of primal-dual residuals is calculated as When inequality (24) is not satisfied, the optimization process is restarted.η is a number close to 1, and is set to 0.96 to prevent the optimization process from being frequently restarted.
i , the auxiliary variable and the acceleration step are updated as When a restart occurs, the variables are updated according to the following formula: The algorithm incorporating the accelerated ADMM is summarized as Algorithm 2 and is named GGS_LP_Fast in this study.

Numerical Experiments
In this section, standard test images of various styles (Figure 3) are selected as experimental objects to evaluate the model proposed in this study.In particular, the image "Satellite" is downloaded from Emory University image database [50], while other images are downloaded from CVG-UGR (Computer Vision Group, University of Granada) image database [51], with each image consisting of 512 × 512 pixels.The downloaded images are converted to an image size of 256 × 256 pixels before the image reconstruction test.To verify the rationality and effectiveness of the proposed models, Gaussian blur and impulse noise are added to the experimental objects in the simulation experiment, whereby a large number of numerical experimental results are obtained to demonstrate the advantages of the proposed method.

Evaluation Criteria and Running Environment
A comparison is performed in terms of two evaluation indicators commonly used in image processing: the peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) [52].The definitions of PSNR, and SSIM are as follows: In (28), X denotes the original image and Y denotes the reconstructed image; µ X and µ Y denote the mean of X and Y, respectively; µ 2 X and µ 2 Y denote the variance of X and Y, respectively; µ XY is the covariance of X and Y; and constants k 1 and k 2 ensure a nonzero denominator in (28), with k 1 = 0.01 and k 2 = 0.03 in the experiment.Generally, the larger the PSNR of the reconstructed image, the better the reconstructed image quality.The SSIM lies in the range (0, 1), and the closer the value is to 1, the higher the similarity, the better the image reconstruction effect.
In this section, standard test images of various styles (Figure 3) are selected as experimental objects to evaluate the model proposed in this study.In particular, the image "Satellite" is downloaded from Emory University image database [50], while other images are downloaded from CVG-UGR (Computer Vision Group, University of Granada) image database [51], with each image consisting of 512 512 × pixels.The downloaded images are converted to an image size of 256 256 × pixels before the image reconstruction test.To verify the rationality and effectiveness of the proposed models, Gaussian blur and impulse noise are added to the experimental objects in the simulation experiment, whereby a large number of numerical experimental results are obtained to demonstrate the advantages of the proposed method.

Evaluation Criteria and Running Environment
A comparison is performed in terms of two evaluation indicators commonly used in image processing: the peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) [52].The definitions of PSNR, and SSIM are as follows: In the experiment, image reconstruction effectiveness is compared between the FTVd model (The matlab could be download from the website https://www.caam.rice.edu/~{}optimization/L1/ftvd/ Version 4.1.),FrTV model, TGV model, OGS_L1 model, and the two models proposed in this study (GGS_Lp and GGS_LP_Fast), both of which are frequency-domain models.To ensure the objectivity and fairness of the evaluation, the iterative conditions of the above algorithms are terminated when satisfying (29), where F (k) denotes the objective function of a model after the k-th iteration.In the experiment, by setting the relax parameter γ = 1.618 and adjusting the regularization parameters of each algorithm, each algorithm is allowed to operate at its best performance to ensure the objectiveness of the test results.
Hardware environment: The processor is an Intel®Core™ i7-6700 CPU @ 3.4 GHz with 16.0 GB of memory.Simulation platform: MATLAB R2018a.
Each blur kernel used in the test is generated by the built-in MATLAB function fspecial ('gaussian', N1, N2), which generates a Gaussian blur kernel with a standard deviation of N2 and window size of N1 × N1, with the kernel hereinafter referred to as a GN1 × N1, σ = N2 kernel for simplicity.Impulse noise is added at a noise level of 30%, 40%, 50%, and 60%.

Number of Inner Iteration Steps
When solving the X 2 , X 3 subproblems, it is necessary to consider the number of inner iteration steps of the MM algorithm.In general, the greater the number of iteration steps, the more accurate the solution of a subproblem, but the more time-consuming the computational process.To find a proper tradeoff, the values taken by the parameter n are first explored while other parameters are fixed.Each image is subject to multiple tests at different noise levels, as exemplified below.
The images of (a) Lena, (d) plane, and (f) fish are convoluted with a G7 × 7, σ = 5 kernel at an impulse noise level of 30%.The experimental results are presented in Figure 4.As indicated by Figure 4a,b, both PSNR and SSIM are improved with an increasing number n of inner iterations, but when n exceeds 5, the promotional effect of the iteration step number on PSNR and SSIM becomes less significant.With an increase in the iteration step number, the computational time increases, Figure 4c.As shown in Figure 4a,b, some of the values have little dropped as the graph is progressed, which is because the other parameters do not be tuned to the best for different n.Based on a comprehensive consideration of the computational time and optimal PSNR and SSIM values, the number of inner iteration steps is set to 5 in the subsequent test.

Group Gradient Parameter K
The group gradient parameter K of the proposed algorithm is tested and compared to evaluate its effect on the overall performance of the algorithms.Given a different image and a different noise level, K is continuously changed from 1 to 10, and at each value, multiple tests are performed.Take the image of (e) seabird as an example, which is convoluted with a 7 7  5 blur kernel at noise levels of 30%, 40%, 50%, and 60%.The PSNR and SSIM results are recorded under these different scenarios, as shown in Figure 5.When 1 K = , the algorithm degenerates into a traditional ATV model.
As shown in Figure 5, an increase in K has different effects on PSNR and SSIM at different noise levels.For example, at a 30% noise level with 5 K = , both PSNR and SSIM reach their maximum values, and this trend is invariant.It is evident that given a low noise level, the neighborhood information of the image has positive effects on the algorithm performance.However, when the noise level is 60%, the structural characteristics of the neighborhood gradient are seriously damaged, and in this circumstance, a greater K may lead to decreased algorithm performance.
In addition, the greater the parameter K , the more information is used from the neighborhood during the computation and the higher the computational complexity, thus increasing the computational time.With K set to a proper value, the proposed algorithm has a greater capability to preserve image edges and to resist high-level noise.Based on a comprehensive consideration of both computational time and the optimal PSNR and SSIM values, K is set to 3 in the subsequent test.

Group Gradient Parameter K
The group gradient parameter K of the proposed algorithm is tested and compared to evaluate its effect on the overall performance of the algorithms.Given a different image and a different noise level, K is continuously changed from 1 to 10, and at each value, multiple tests are performed.Take the image of (e) seabird as an example, which is convoluted with a G7 × 7, σ = 5 blur kernel at noise levels of 30%, 40%, 50%, and 60%.The PSNR and SSIM results are recorded under these different scenarios, as shown in Figure 5.When K = 1, the algorithm degenerates into a traditional ATV model.
As shown in Figure 5, an increase in K has different effects on PSNR and SSIM at different noise levels.For example, at a 30% noise level with K = 5, both PSNR and SSIM reach their maximum values, and this trend is invariant.It is evident that given a low noise level, the neighborhood information of the image has positive effects on the algorithm performance.However, when the noise level is 60%, the structural characteristics of the neighborhood gradient are seriously damaged, and in this circumstance, a greater K may lead to decreased algorithm performance.

Group Gradient Parameter K
The group gradient parameter K of the proposed algorithm is tested and compared to evaluate its effect on the overall performance of the algorithms.Given a different image and a different noise level, K is continuously changed from 1 to 10, and at each value, multiple tests are performed.Take the image of (e) seabird as an example, which is convoluted with a 7 7  5 blur kernel at noise levels of 30%, 40%, 50%, and 60%.The PSNR and SSIM results are recorded under these different scenarios, as shown in Figure 5.When 1 K = , the algorithm degenerates into a traditional ATV model.
As shown in Figure 5, an increase in K has different effects on PSNR and SSIM at different noise levels.For example, at a 30% noise level with 5 K = , both PSNR and SSIM reach their maximum values, and this trend is invariant.It is evident that given a low noise level, the neighborhood information of the image has positive effects on the algorithm performance.However, when the noise level is 60%, the structural characteristics of the neighborhood gradient are seriously damaged, and in this circumstance, a greater K may lead to decreased algorithm performance.
In addition, the greater the parameter K , the more information is used from the neighborhood during the computation and the higher the computational complexity, thus increasing the computational time.With K set to a proper value, the proposed algorithm has a greater capability to preserve image edges and to resist high-level noise.Based on a comprehensive consideration of both computational time and the optimal PSNR and SSIM values, K is set to 3 in the subsequent test.In addition, the greater the parameter K, the more information is used from the neighborhood during the computation and the higher the computational complexity, thus increasing the computational time.With K set to a proper value, the proposed algorithm has a greater capability to preserve image edges and to resist high-level noise.Based on a comprehensive consideration of both computational time and the optimal PSNR and SSIM values, K is set to 3 in the subsequent test.

Regularization Parameter µ and Quasinorm p
Next, the effects of regularization parameter µ on the test results are evaluated.Given a different test image at various impulse noise levels, µ is gradually increased from a small value to a high value while the values of PSNR and SSIM of image reconstruction are recorded.Here, take the test images of (c) boat and (g) house as examples, where the images are convoluted with a G7 × 7, σ = 5 blur kernel at noise levels of 30%, 40%, 50%, and 60%, respectively, while the values of PSNR at various values of µ are recorded, as shown in Figure 6.The experimental results show that given a test image, PSNR is relatively stable when µ is greater than a certain value, suggesting that the proposed algorithm is not sensitive to the regularization parameter µ.That is, the optimal values of µ are relatively stable, and the proposed algorithm is robust.In the subsequent experiment, µ is set to 100, 90, 80, and 60 at noise levels of 30%, 40%, 50%, and 60%, respectively.
Information 2019, 10, x 14 of 22 PSNR at various values of μ are recorded, as shown in Figure 6.The experimental results show that given a test image, PSNR is relatively stable when μ is greater than a certain value, suggesting that the proposed algorithm is not sensitive to the regularization parameter μ .That is, the optimal values of μ are relatively stable, and the proposed algorithm is robust.In the subsequent experiment, μ is set to 100, 90, 80, and 60 at noise levels of 30%, 40%, 50%, and 60%, respectively.
( With the above-determined parameters, p is allowed to vary from 0 to 1 in steps of 0.05 in light of the characteristics of the impulse noise, and a loop statement is used for a thorough search for the p value that optimizes PSNR, as shown in Table 1.Test images are convoluted with a 7 7 5 G σ × = ， blur kernel at noise levels of 30%, 40%, 50%, and 60%, respectively.The results show that to a certain test image, the p value that optimizes PSNR is relatively stable.To further test the effectiveness of the proposed algorithm, tests are subject to in-depth comparisons with other state-of-the-art algorithms.The group of images in Figure 3 is used as the test objects, which are convoluted with a 7 7  5 G σ × = ， blur kernel at impulse noise levels of 30%, 40%, 50%, and 60%, respectively.The test results for the different images are listed in Table 2, with the optimal indicators highlighted in boldface and black font.With the above-determined parameters, p is allowed to vary from 0 to 1 in steps of 0.05 in light of the characteristics of the impulse noise, and a loop statement is used for a thorough search for the p value that optimizes PSNR, as shown in Table 1.Test images are convoluted with a G7 × 7, σ = 5 blur kernel at noise levels of 30%, 40%, 50%, and 60%, respectively.The results show that to a certain test image, the p value that optimizes PSNR is relatively stable.Based on the data in the above table, the following conclusions can be drawn: 1.
With the addition of blurring and different levels of noise in various images, the proposed models in this study are superior to other excellent models for image reconstruction in terms of PSNR and SSIM.This indicates that the proposed models have good deblurring and denoising performance, while allowing the reconstructed images to be more similar to the original images.

2.
When the proposed models are employed to reconstruct the eight images, the PSNR values are higher than those of the TGV model by 0.204-2.177dB, and the superiority of the proposed models is more evident at higher levels of noise.For example, PSNR of the proposed models (28.081 dB) is 2.177 dB higher than that of the TGV model (25.904 dB) for the image of Lena at a noise level of 60%, and the PSNR difference becomes 2.174 dB for the image of the boat at the same noise level (26.797 dB vs. 24.623dB).

3.
Compared with OGS_L1, the proposed models increase the Lp-quasinorm shrinkage and improve the ability to characterize image gradient sparsity.When recovering the eight images, the proposed models have PSNR values that are 0.328-1.964dB higher than those of the OGS_L1 model.For the butterfly image at the noise level of 60%, PSNR of the proposed models (26.077 dB) is 1.964 dB higher than that of the OGS_L1 model (24.113 dB).
To better understand the test results, some reconstructed images are illustrated below: Figure 7 depicts the image reconstruction results of several models for the butterfly, boat, house, and plane images blurred by convolution with a G7 × 7, σ = 5 blur kernel at impulse noise levels of 30-60%.The proposed model of this study shows good image reconstruction performance at each noise level, with PSNR values greater than those of the other advanced models.It is noteworthy that even in the presence of high-level noise contamination, the proposed algorithm is still able to provide the reconstructed images with good visual effect.8 depicts zoomed-in seabird images recovered by several models after image blurring owing to a G5 × 5, σ = 5 blur kernel at a 40% level of impulse noise.As shown by the visual effects of the image recovery results, the FTVd method leads to a staircase effect in the recovered image, as shown in Figure 8d.The protectory effect of the FrTV method on the image edges needs to be improved, as shown in Figure 8e.The TGV and OGS_L1 methods have satisfactory inhibitory performance on the staircase effect, but local areas of the TGV-reconstructed image still contain some highly noisy pixels, as shown in Figure 8f,g.The OGS_L1 method reconstruction leads to relatively smooth image reconstruction results for smooth image regions, but some image edges are oversmoothed.The proposed model is able to reconstruct the image pixels of similar grayscale values in the smooth image regions well, and avoids the staircase effect while preserving the image edges, thereby providing a satisfactory overall reconstruction performance.blur kernel at a 40% level of impulse noise.As shown by the visual effects of the image recovery results, the FTVd method leads to a staircase effect in the recovered image, as shown in Figure 8d.The protectory effect of the FrTV method on the image edges needs to be improved, as shown in Figure 8e.The TGV and OGS_L1 methods have satisfactory inhibitory performance on the staircase effect, but local areas of the TGV-reconstructed image still contain some highly noisy pixels, as shown in Figure 8f,g.The OGS_L1 method reconstruction leads to relatively smooth image reconstruction results for smooth image regions, but some image edges are oversmoothed.The proposed model is able to reconstruct the image pixels of similar grayscale values in the smooth image regions well, and avoids the staircase effect while preserving the image edges, thereby providing a satisfactory overall reconstruction performance.blur kernel at an impulse noise level of 40%.To further compare the results of various models, the 100th-channel signals are extracted from the image.It is evident that the TGV-reconstructed image still contains some noisecontaminated pixels in some local areas, as shown in Figure 9d.In Figure 9e, the OGS_L1 algorithm results in thorough noise removal, but the image edges are oversmoothed to some extent.The proposed algorithm can remove the noise entirely while preserving the image edges well, as shown in Figure 9f.Finally, to further test the timeliness of the proposed models, the mean time spent for image reconstruction is compared between the FTVd, TGV, and OGS_L1 models and the proposed models under the premise that the abovementioned test procedures are free of bias.Here, the proposed models in this study are named the GGS_Lp model and the GGS_LP_Fast model, with the latter being a model for accelerated restart of the former.In the test, images are blurred with a G7 × 7, σ = 5 blur kernel at impulse noise levels of 30-60%.The test results are partially presented in Table 3, with the optimal results highlighted in boldface and black font.

Comparison of Image Reconstruction Time for Several Algorithms
Finally, to further test the timeliness of the proposed models, the mean time spent for image reconstruction is compared between the FTVd, TGV, and OGS_L1 models and the proposed models under the premise that the abovementioned test procedures are free of bias.Here, the proposed models in this study are named the GGS_Lp model and the GGS_LP_Fast model, with the latter being a model for accelerated restart of the former.In the test, images are blurred with a 7 7 5 G σ × = ， blur kernel at impulse noise levels of 30-60%.The test results are partially presented in Table 3, with the optimal results highlighted in boldface and black font.

Amplitude
Amplitude Amplitude Amplitude Amplitude Amplitude Analysis of the data in Table 3 leads to the following conclusions: 1.
The GGS_Lp model and its accelerated restart version GGS_LP_Fast, both proposed in this study, lead to similar PSNR values when the two models are used to reconstruct images, but their PSNR values are all greater than those of the FTVd, TGV, and OGS_L1 models.

2.
As a classic algorithm, FTVd is still advantageous in terms of computational speed.TGV and OGS_L1 give better image reconstruction results than FTVd, but they are more time-consuming.In this study, neighborhood gradient information is thoroughly considered, and the Lp-quasinorm is introduced as a regularization constraint, which leads to a greater computational complexity and thereby greater time consumption for image reconstruction compared with the FTVd model.By contrast, the GGS_LP_Fast model increases computational efficiency to a great extent, and its performance is especially prominent for images with a high level of noise.

Figure 2 .
Figure 2. Two-dimensional diagram of sparse groups: (a) smooth region contaminated by noise and (b) edge region not contaminated by noise.

Figure 2 .
Figure 2. Two-dimensional diagram of sparse groups: (a) smooth region contaminated by noise and (b) edge region not contaminated by noise.

Figure 4 .
Figure 4. Effects of iteration steps number on the results: (a) PSNR, (b) SSIM, and (c) processing time.

Figure 5 .
Figure 5. Effects of group gradient parameter K on experimental results: (a) PSNR, (b) SSIM, and (c) processing time.4.2.3.Regularization Parameter μ and Quasinorm p Next, the effects of regularization parameter μ on the test results are evaluated.Given a

Figure 4 .
Figure 4. Effects of iteration steps number on the results: (a) PSNR, (b) SSIM, and (c) processing time.

Figure 4 .
Figure 4. Effects of iteration steps number on the results: (a) PSNR, (b) SSIM, and (c) processing time.

Figure 5 .
Figure 5. Effects of group gradient parameter K on experimental results: (a) PSNR, (b) SSIM, and (c) processing time.4.2.3.Regularization Parameter μ and Quasinorm p Next, the effects of regularization parameter μ on the test results are evaluated.Given a different test image at various impulse noise levels, μ is gradually increased from a small value to a high value while the values of PSNR and SSIM of image reconstruction are recorded.Here, take the test images of (c) boat and (g) house as examples, where the images are convoluted with a 7 7 5 × = , G σ

Figure 5 .
Figure 5. Effects of group gradient parameter K on experimental results: (a) PSNR, (b) SSIM, and (c) processing time.

Figure 7 .
Figure 7. Random examples of images reconstructed by several advanced models and proposed model.Top row depicts the blurred and noisy images of the butterfly, boat, house, and plane owing to convolution with a 7 7 5 G σ × = ， blur kernel and impulse noise levels of 30%, 40%, 50%, and 60%.

Figure 7 .
Figure 7. Random examples of images reconstructed by several advanced models and proposed model.Top row depicts the blurred and noisy images of the butterfly, boat, house, and plane owing to convolution with a G7 × 7, σ = 5 blur kernel and impulse noise levels of 30%, 40%, 50%, and 60%.

Figure
Figure8depicts zoomed-in seabird images recovered by several models after image blurring owing to a G5 × 5, σ = 5 blur kernel at a 40% level of impulse noise.As shown by the visual effects of the image recovery results, the FTVd method leads to a staircase effect in the recovered image, as shown in Figure8d.The protectory effect of the FrTV method on the image edges needs to be improved, as shown in Figure8e.The TGV and OGS_L1 methods have satisfactory inhibitory performance on the staircase effect, but local areas of the TGV-reconstructed image still contain some highly noisy pixels, as shown in Figure8f,g.The OGS_L1 method reconstruction leads to relatively smooth image reconstruction results for smooth image regions, but some image edges are oversmoothed.The proposed model is able to reconstruct the image pixels of similar grayscale values in the smooth image regions well, and avoids the staircase effect while preserving the image edges, thereby providing a satisfactory overall reconstruction performance.

Figure 8
Figure 8 depicts zoomed-in seabird images recovered by several models after image blurring owing to a 5 5 5 G σ × = ，blur kernel at a 40% level of impulse noise.As shown by the visual effects of the image recovery results, the FTVd method leads to a staircase effect in the recovered image, as shown in Figure8d.The protectory effect of the FrTV method on the image edges needs to be improved, as shown in Figure8e.The TGV and OGS_L1 methods have satisfactory inhibitory performance on the staircase effect, but local areas of the TGV-reconstructed image still contain some highly noisy pixels, as shown in Figure8f,g.The OGS_L1 method reconstruction leads to relatively smooth image reconstruction results for smooth image regions, but some image edges are oversmoothed.The proposed model is able to reconstruct the image pixels of similar grayscale values in the smooth image regions well, and avoids the staircase effect while preserving the image edges, thereby providing a satisfactory overall reconstruction performance.

Figure 8 .
Figure 8.Comparison of the zoomed-in views of reconstructed images by several advanced models and proposed model: (a) selected area, (b) zoomed-in part of selected area, (c) blurred and noisy zoomed-in part owing to convolution with a 5 5 5 G σ × = ， blur kernel at a impulse noise level of 40%, and (d-h) zoomed-in parts of reconstructed images by FTVd, FrTV, TGV, OGS_L1, and proposed model, respectively.

Figure 9
Figure 9 depicts a comparison of grayscale values in single channels of the satellite image as recovered by several algorithms after image blurring with a 9 9 5 G σ × = ，blur kernel at an impulse noise level of 40%.To further compare the results of various models, the 100th-channel signals are extracted from the image.It is evident that the TGV-reconstructed image still contains some noisecontaminated pixels in some local areas, as shown in Figure9d.In Figure9e, the OGS_L1 algorithm results in thorough noise removal, but the image edges are oversmoothed to some extent.The proposed algorithm can remove the noise entirely while preserving the image edges well, as shown in Figure9f.

Figure 8 .
Figure 8.Comparison of the zoomed-in views of reconstructed images by several advanced models and proposed model: (a) selected area, (b) zoomed-in part of selected area, (c) blurred and noisy zoomed-in part owing to convolution with a G5 × 5, σ = 5 blur kernel at a impulse noise level of 40%, and (d-h) zoomed-in parts of reconstructed images by FTVd, FrTV, TGV, OGS_L1, and proposed model, respectively.

Figure 9
Figure9depicts a comparison of grayscale values in single channels of the satellite image as recovered by several algorithms after image blurring with a G9 × 9, σ = 5 blur kernel at an impulse noise level of 40%.To further compare the results of various models, the 100th-channel signals are extracted from the image.It is evident that the TGV-reconstructed image still contains some noise-contaminated pixels in some local areas, as shown in Figure9d.In Figure9e, the OGS_L1 algorithm results in thorough noise removal, but the image edges are oversmoothed to some extent.The proposed algorithm can remove the noise entirely while preserving the image edges well, as shown in Figure9f.
Figure9depicts a comparison of grayscale values in single channels of the satellite image as recovered by several algorithms after image blurring with a G9 × 9, σ = 5 blur kernel at an impulse noise level of 40%.To further compare the results of various models, the 100th-channel signals are extracted from the image.It is evident that the TGV-reconstructed image still contains some noise-contaminated pixels in some local areas, as shown in Figure9d.In Figure9e, the OGS_L1 algorithm results in thorough noise removal, but the image edges are oversmoothed to some extent.The proposed algorithm can remove the noise entirely while preserving the image edges well, as shown in Figure9f.

4. 2 . 5 .
Comparison of Image Reconstruction Time for Several Algorithms

Figure 9 .
Figure 9.Comparison of the proposed model and several advanced models in terms of the 100thchannel data of the reconstructed satellite images.(a) refers to the grayscale values of the 100th pixel channel in the original satellite image, with the yellow vertical line on the figure indicating the 100thchannel position.(b) refers to a blurred image owing to convolution with a 9 9 5 G σ × = ， blur kernel at impulse noise level of 40%.(c-f) depict images with grayscale values of the 100th channel recovered by the FTVd, TGV, OGS_L1, and proposed algorithms, respectively.

Figure 9 .
Figure 9.Comparison of the proposed model and several advanced models in terms of the 100th-channel data of the reconstructed satellite images.(a) refers to the grayscale values of the 100th pixel channel in the original satellite image, with the yellow vertical line on the figure indicating the 100th-channel position.(b) refers to a blurred image owing to convolution with a G9 × 9, σ = 5 blur kernel at impulse noise level of 40%.(c-f) depict images with grayscale values of the 100th channel recovered by the FTVd, TGV, OGS_L1, and proposed algorithms, respectively.

Table 1 .
The p value that optimizes PSNR for different images.

Table 2 .
PSNR (dB) and SSIM values generated by various algorithms in deblurring and denoising different images.

Table 1 .
The p value that optimizes PSNR for different images.Comparative Testing of Image Reconstruction Effectiveness between Several Algorithms To further test the effectiveness of the proposed algorithm, tests are subject to in-depth comparisons with other state-of-the-art algorithms.The group of images in Figure3is used as the test objects, which are convoluted with a G7 × 7, σ = 5 blur kernel at impulse noise levels of 30%, Information 2019, 10, x 16 of 22 even in the presence of high-level noise contamination, the proposed algorithm is still able to provide the reconstructed images with good visual effect.