Non-Blind Image Deblurring Method Using Shear High Order Total Variation Norm

: In this paper, we propose a shear high-order gradient (SHOG) operator by combining the shear operator and high-order gradient (HOG) operator. Compared with the HOG operator, the proposed SHOG operator can incorporate more directionality and detect more abundant edge information. Based on the SHOG operator, we extend the total variation (TV) norm to shear high-order total variation (SHOTV), and then propose a SHOTV deblurring model. We also study some properties of the SHOG operator, and show that the SHOG matrices are Block Circulant with Circulant Blocks (BCCB) when the shear angle is π 4 . The proposed model is solved efficiently by the alternating direction method of multipliers (ADMM). Experimental results demonstrate that the proposed method outperforms some state-of-the-art non-blind deblur-ring methods in both objective and perceptual quality.


Introduction
Image deblurring is a classical problem in low level vision.A blurred image can be considered as a convolution between the latent image and the point spread function (PSF) of the image capture device.The blurring process can be modeled as   f Ku n (1)  where is an original image without any form of degradation, n n   R K is a degradation matrix, f is an observed image, and n is the additive noise.
The goal of image deblurring is to estimate the original image u from an observed image f degraded by blur kernel and noise [1] .In this paper, we mainly perform image restoration processing on blur kernel and Gaussian noise.The function of Gaussian noise can be expressed as where  and  represent the standard derivative and average value of the noise distribution n, respectively.Image deblurring is an ill-posed inverse problem.In other words, its solution neither exists nor is unique, and small disturbances in the data may cause large errors in reconstruction.Hence, in order to overcome this ill-posedness, regularization techniques need to be considered to produce reasonable approximate solutions from noisy data.Rudin-Osher-Fatemi (ROF) model [2] is one of the most successful regularization methods due to its ability to preserve sharp edges [3] .The corresponding minimization task is as follows:  denotes the Euclidean norm, TV  ‖ ‖ is the discrete TV regularization term.The first term of (3) is called the fidelity term, and the second term is called the regularization term (or penalty) which describes prior knowledge of the image.0  ∨ is the regularization parameter which controls the balance between fidelity term and regularization term.During the last few decades, a series of methods based on total variation and its variants have been developed.
The well-known first-order total variation methods such as FTVd (fast total variation deconvolution) [4] and ADMM (alternating direction method of multipliers) TV [5,6] have also achieved satisfactory results in the image restoration.FTVd is a total variation model based on splitting techniques.In Ref. [7], Chan et al proposed a frame-based deblurring problem using the ADMM algorithm.In Ref. [5], Jiao et al used the ADMM iterative method to solve the TV deblurring problem (ADMM TV).Compared with FTVd, the ADMM TV method can achieve comparable experimental performance without the need for parameter selection.
However, it is well known that TV-based restoration methods suffer from the staircase artifact in flat regions.There are many image restoration methods that have been developed to overcome these shortcomings.A well-known method is to replace the first-order total variation in the regularization term with a second-order (or higher) total variation.
In 2000, Chan et al [8] proposed the high-order TV (HOTV) model to improve the restoration quality.Contrary to TV regularization that penalizes gradient magnitude, regularizations using second-order derivatives do not penalize ramp regions whose intensity varies linearly, and therefore do not force the intensity of the region to remain constant.In other words, these regularizations have ability to potentially recover a broader class of images that consist of more than just piecewise-constant regions [9] .In Ref. [10], a second-order Laplace model for image deblurring was proposed by You and Kaveh.In 2003, a second-order derivative regularization [11] (Hessian regularization) model was proposed, called the Lysaker, Lundervold and Tai (LLT) model, which can effectively suppress the staircase artifacts and maintain the smoothness of flat areas of the image.In 2012, Hu et al [12] improved the model for image recovery and designed the maximum-minimum (MM) algorithm to solve the model.Their numerical experiments show that second-order TV has an advantage over classical TV regularization in avoiding staircase artifacts.In Ref. [13], the authors pointed out that high-order TV's ability to preserve edges became weaker compared with first-order TV regularization and tended to smoothen edges and other small details.Thus, hybrid models that combine higher-order TV with other regularization were proposed in many literatures [13][14][15] .
Among the hybrid models, by replacing the first-order TV with a second-order TV, Zhang et al [16] obtained better results than the model with strongly convex terms presented in Ref. [17].Furthermore, in Refs.[18]  and [19], a combination of first-order and second-order TV priors is used to recover blurred images and the minimization model is solved by the ADMM algorithm.To further improve the restoration performance, Liu et al [20] investigated a spatially adaptive regularization parameter updating scheme.As an adaptive balancing scheme between first-order and second-order TV, Wang et al [21] developed a Poisson denoising framework based on an iterative weighted total generalized variational (TGV) model.In Ref. [22], Liu et al pointed out that TGV can maintain sharp edges and perform well in areas of gradual intensity slopes.
TV or HOTV norm only contains finite difference in horizontal and vertical directions, and this leads to losing some detail information in other directions.In order to incorporate more direction information in gradient domain, in this paper we propose a novel shear high-order TV (SHOTV) norm which contains finite difference in different directions.Compared with higher order TV regularization, the proposed SHOTV regularization method provides a more meaningful description of the gradient information as it contains information about the different directions of the image gradients.
where the ˆ( , ) i j -th entry of shr X is given, for ˆ1, , i m   and ˆ1, , j n   , by shr ˆˆ, , :  1.2 Shear Gradient and Shear HOTV Norm TV regularization was introduced by Rudin et al [2] , for image denoising and further extended to other image restoration tasks.For the gray-scale image , if , if TV norm only contains finite difference in horizontal and vertical directions, and this leads to losing some detailed information in other directions.
TV regularization preserves fine features and sharp edges, but it usually produces staircase artifact that is not presented in real images [24,25] .In many applications [26,27] , some authors have proposed a higher-order TV to reduce staircase artifact.HOTV of u is defined as ) j n i   -th entry of the vector u .By combining the higher-order finite difference operator and the shear operator, we propose the shear high-order gradient (SHOG) operator by the following way: Then for the gray-scale image u , its discrete shear higher-order total variation (SHOTV) norm is The shear gradient has following properties: .In fact, for , by the definition of gradient operator and shear operator, we have , , where j is defined as in Eq. ( 7).

S D S S D S S D S D D D
By Theorem 1, the shear gradient operator is essentially a generalization of gradient operator.
The comparison between HOG operator and SHOG operator applied on Barbara image are shown in Fig. 2. It can be observed that more direction information is captured and more detailed edge information can be retained by the SHOG operator.From the local comparison graph, it can be seen that the SHOG operator preserves more information in the texture part of the image.Comparison between 2nd-gradient operator and shear 2nd-gradient operator applied on Barbara image u

BCCB Matrices
Under the periodic boundary conditions, both the blurring matrix and gradient matrix are Block-Circulant-with-Circulant-Blocks (BCCB), and thus are diagonalizable by the 2D discrete Fourier transform.In this subsection, we will show that the SHOG matrices are BCCB when the shear angle is π 4 .

D
It can be noticed that the above matrix is not a BCCB matrix.

Proposed Model
Considering the advantages of high-order TV and shear operator, we propose the following optimization model where the term 1  ‖ ‖ represents convex 1  norm regularization.The variables , 0   ∨ are regularization parameters that control the data fidelity term and the shear high-order regularization.
( ) I  u is an indication function to impose hard constraints on the objective function defined as 0 , if ( ) , if . To minimize problem ( 16) robustly, we hereby adopt framework of alternating direction method of multipliers (ADMM) [28] .Equation ( 16) can be converted into a constrained optimization problem: The corresponding augmented Lagrangian function of (17) can be written as ) 2 where 1 2 ,   are penalty parameters.Variables 1 μ and 2 μ are the Lagrange multipliers associated with the constraints S  w D u and  z u , respectively.ADMM minimizes the Lagrangian function on two sets of variables u , w and z .Since w and z are decoupled from each other, they can be solved separately.
ADMM is an algorithm using proximal splitting techniques for solving the following convex optimization problem: ) where :  are nonempty closed convex sets, and c is a given vector.
Part of the appeal of ADMM is the fact that the algorithm lends themselves to parallel implementations.The algorithm is given in Algorithm 1. 2 Iteration: 1) Until a stopping criterion is satisfied.
In this section, we discuss the optimization strategies for solving each sub-problem of (18) one by one.

1) u Sub-Problem
The sub-problem of u is a least square problem with the following form From the optimality conditions, we have Considering what was discussed in subsection 1.In this paper, we set the shear angle π 4 d   .T K K is also BCCB matrix under the periodic boundary condition.As a consequence, Eq. ( 21) can be efficiently solved by one FFT operation and one inverse FFT operation as where F represents the two-dimensional discrete Fourier transform and 1 F  is the inverse transform.

2) w Sub-Problem
The w sub-problem is a convex second-order deblurring problem We write for simplicity, and solve problem ( 23) using a two-dimensional shrinkage operator.i.e.,   3) z Sub-Problem The z sub-problem is a projection problem.For an 8-bit image, the range of pixel values can be kept at [0, 255] , and the closed-form solution can be obtained by using the projection operator The Lagrangian multipliers are updated as follows The proposed algorithm is named SHOTV and shown in Algorithm 2.
6) Check the stopping criterion

Numerical Experiments
In this section, we present various numerical results to illustrate the performance of the proposed algorithm for non-blind image deblurring.All experiments are carried out on Windows 10 64-bit and Matlab 2019b running on a desktop equipped with an Intel Core i5-10300 H CPU 2.5 GHz and 16 GB of RAM.The source code for all competing methods was obtained from the original authors and we use the default parameter settings.

Experiment Setting
The test images are shown in Fig. 3, and the pixel values of the images are normalized to [0, 1] for simplicity.Both Peak Signal to Noise Ratio (PSNR) and Structural Similarity metrics (SSIM) [29] are used to evaluate the quality of the reconstructed images.They are defined as follows: where u is the clean image, f is the recovered

Parameter Selection Discussion
In this subsection, we focus on the choice of parameters 1 2  ,   and  .For simplicity, we set the pa-  under three different noise levels.We empirically set the parameters 1 [5,30]     in the subsequent experiments.
Then, we discuss the selection of the parameter  .
Our experiments are tested on Pirate and Barbara images with three different levels of noise and two different blur kernels.For simplicity, we set the parameter 1 2     10. Figure 5 shows the curves of PSNR value with respect to parameter  under three different noise levels.For the parameter  in the subsequent experiments, the empirical value is taken to be in the range [800, 2 000].As a matter of fact, the above empirical parameter settings Fig. 3 The test images Fig. 4 The PSNR results versus the parameter 1  for fixed  Fig. 5 The PSNR results versus the parameter  for fixed
First, we test the 9 9  Gaussian blur with a standard deviation of 9 .The blurred image is generated using the MATLAB function "imfilter" with periodic boundary conditions.Then, the blurred image is corrupted with zero-mean additive white Gaussian noise 2 0.001   .The PSNR and SSIM values for the deblurring experiments are reported in Table 1 and Table 2.The numbers in bold indicate the best performance.One can observe that the proposed method outperforms the other competing methods.We choose Fingerprint, Cameraman, Hill and Shirt images for Gaussian deblurring and denoising comparison, and the recovered images are shown in Fig. 6.From the comparison, we can see that the restored images by our method retain more texture and have fewer artifacts.As can be seen from the Shirt image in Fig. 6, our proposed algorithm retains more texture information through the shear operator, resulting in a significant improvement in image deblurring.
Similarly, we test the motion blur with a len size 35 in the blur kernel and a rotation angle 50º.The blurred image is generated using the MATLAB function "imfilter" with periodic boundary conditions, and the blurred image is corrupted by zero-mean additive white Gaussian noise with 2 0.001   .The PSNR and SSIM results for all competing deblurring methods are reported in Table 2.One can observe that the proposed SHOTV method can obtain better result than other competing mehtods.The blurred images and the deblurred images are shown in Fig. 7. Compared with other methods, our algorithm can preserve more texture structure and the deblurred image contains fewer artifacts compared with first-order TV.

Convergence Analysis
In order to verify the convergence of the algorithm numerically, we apply Algorithm 2 to four different images corrupted by different blur kernels and additive white Gaussian noise with deviation 2 0.001

 
, where Gaussian blur is applied to Fingerprint and Shirt images, and Gaussian blur is applied to Barbara and Zebra images.We use the equation to calculate the relative error curve for each iteration of the recovery image of Algorithm 2. As shown in Fig. 8, the relative error decreases as the number of iterations increases, which numerically illustrates the convergence of the algorithm.

S
6) where | .| represents the procedure rounding the real number to the nearest integer in the zero direction.The parameter d indicates the axis of shear (top ' t , bottom 'b , left 'l , right 'r ), and d  determines the degree of shear.For example, d  S X means that the image X is sheared with top axis and d  .In particuis the identity transformation, i.e., no shear.Figure 1 shows the results of shear operator in different directions at an angle of π 4 .
Fig. 1 Effect of shear operator d  S with different angles and directions (a) Original image X; (b) d S X  with u D u D u D u denotes the second-order difference of the (( 1) the 'top'(or 'bottom') and 'left'(or 'right') direction with different angles, respectively.With different types of texture images, we choose different directions based on our experience.

Remark 1 (SS
Matrix form of d  S ) The operator d  is linear, so it can be represented as matrix, which is denoted by and the vertical operator yy D are commutative when 't or b ' d    .i.e.,

S
unfortunately, when the shear angle d has the following form:

D 4 .
D are BCCB matrices only when the shear angle is π While the shear angle is not π 4 , computing the inverse of the left hand of (21) in each iteration is costly.
image, Max u is the largest possible value of the image u .û and f are the mean values of images u and f ,  u and  f are the standard variance of images u and f , respectively, and  uf is the covariance of u and f , 1 2 , 0 c c ∨ are constants.In general, high PSNR and SSIM values imply better image quality.

1 
we investigate the sensitivity of the parameter and  .We test the selection of 1  on both House and Parrot images.They are corrupted by motion blur (fspecial('motion', 35, 50)) and Gaussian blur (fspecial ('gaussian',[9,9], 9)).Figure4shows the change of PSNR values versus the parameter 1 generally valid for other test images.

Fig. 6
Fig. 6 Restoration of four degraded images with Gaussian blur and comparison with other algorithms From left to right: the noisy images, FTVd, ADMM-Frame, ADMM-TV, HOTV, and SHOTV; each value in parentheses represents the corresponding SSIM value of the restored image

Fig. 7
Fig. 7 Restoration of four degraded images with Motion blur and comparison with other algorithmsFrom left to right: the noisy images, FTVd, ADMM-Frame, ADMM-TV, HOTV and SHOTV; each value in parentheses represents the corresponding SSIM value of the restored image

Fig. 8
Fig. 8Relative error values versus iteration number information.We extend the HOTV norm to the SHOTV norm based on the SHOG operator, and then employ SHOTV norm as the regularization term to establish an image deblurring model.We also study some properties of the SHOG operator, and show that the SHOG matrices are Block-Circulant-with-Circulant-Blocks (BCCB) when the shear angle is π 4 .ADMM scheme is exploited to solve the proposed model efficiently.Extensive numerical experiments demonstrate the superiority of the proposed method in terms of visual quality and quantitative metrics.

Table 1 PSNR and SSIM comparison of different methods for Gaussian deblurring
Note: Numbers in bold indicate the best performence