Infrared Image Deblurring via High-Order Total Variation and Lp-Pseudonorm Shrinkage

The quality of infrared images is affected by various degradation factors, such as image blurring and noise pollution. Anisotropic total variation (ATV) has been shown to be a good regularization approach for image deblurring. However, there are two main drawbacks in ATV. First, the conventional ATV regularization just considers the sparsity of the first-order image gradients, thus leading to staircase artifacts. Second, it employs the L1-norm to describe the sparsity of image gradients, while the L1-norm has a limited capacity of depicting the sparsity of sparse variables. To address these limitations of ATV, a high-order total variation is introduced in the ATV deblurring model and the Lp-pseudonorm is adopted to depict the sparsity of lowand high-order total variation. In this way, the recovered image can fit the image priors with clear edges and eliminate the staircase artifacts of the ATV model. The alternating direction method of multipliers is used to solve the proposed model. The experimental results demonstrate that the proposed method does not only remove blurs effectively but is also highly competitive against the state-of-the-art methods, both qualitatively and quantitatively.


Introduction
Thermal infrared imagers can create images of targets when they are completely dark and far away. Moreover, the disguised targets and high-speed moving targets can also be detected in thick smoke screens and clouds. Thus, thermal infrared imagers are widely used in military and civilian applications. Consequently, users prefer higher requirements for the quality of infrared images acquired from thermal infrared imagers. The principle of infrared image generation [1][2][3] is as follows. The object whose temperature is above absolute zero in nature radiates infrared rays outward. Based on the infrared radiation, the infrared system employs a sensor to convert infrared radiation into electricity signals. After signal processing, it presents the corresponding visible light image on the display medium. As the infrared imaging system is more complex than the natural imaging system, the infrared images suffer from relatively more degradation, such as Gaussian blurring, motion blurring, and noise pollution. Therefore, deblurring infrared images plays a significant role in an infrared imaging system. Researchers have proposed some infrared image deblurring methods, for instance, the quaternion and high-order overlapping group sparse total variation model [4], the total variation with overlapping were evaluated, such as the peak signal-to-noise ratio (PSNR) [35], structural similarity (SSIM) [36], and gradient magnitude similarity deviation (GMSD) [37]. The major contributions that have led to significant improvements in the quality of deblurred infrared images are as follows: (1) considering the sparsity of a high-order gradient field based on the ATV model; (2) introducing Lp-pseudonorm shrinkage to express the sparsity of a first-order and high-order gradient field; (3) combining first-order TV and high-order TV with Lp-pseudonorm shrinkage organically; (4) adjusting the parameters in the model separately.
This paper is organized as follows: Section 2 presents the ATV deblurring model. The HTV-Lp model proposed in this paper is described in Section 3. The algorithm for solving the HTV-Lp model is presented in Section 4. Next, the publicly available datasets, evaluation metrics, results of the extensive experiments conducted for evaluation of the six models, and experimental analysis are presented in Section 5. Finally, Section 6 concludes the paper.

Traditional Anisotropic Total Variation (ATV) Model
The ATV deblurring model [38] is described as follows: The Gaussian blurring kernel H ∈ R N×N is represented by a point spread function [5], F ∈ R N×N represents the original image, and G ∈ R N×N represents the observed image with noise. 1 2 H * F − G 2 2 is the fidelity term, R ATV (F) is the prior term, and the balance factor u is used to balance the prior and fidelity terms.
where K h = [−1, 1], K v = −1 1 represents the two-dimensional convolution kernels, K h * F describes information in the horizontal direction of the image, while K v * F describes the information in the vertical direction of the image.

Proposed Model
Because the pixel information is not fully considered in the traditional ATV model, only the information of the first-order gradient is considered. To enhance the deblurring effect, the traditional ATV model is extended to the high-order TV model, which considers the first-order gradient information, and adds the high-order gradient information to the prior term. Figure 1 depicts the contour maps that highlight the advantage of Lp-pseudonorms. The contours of the L1-norm, L2-norm and Lp-pseudonorm are depicted in Figure 1a-c. The L1-norm and L2-norm can be expressed as N j=1 X 2 ij , respectively; whereas, the Lp-pseudonorm is defined in the following paper. In Figure 1, the dotted lines represent the fidelity term, while the solid blue lines represent the contours of the prior item. The red dots represent the intersection of the dotted line and the solid line, while the intersections with the axes represent the better sparseness of the image gradients. It is also observed that the Lp-pseudonorm provides a greater degree of freedom compared with L1-norm and L2-norm, which can better depict the sparsity of the image gradients. Therefore, the Lp-pseudonorm was added to the prior terms to improve the robustness of the image gradients. Appl. Sci. 2019, 9, x FOR PEER REVIEW 4 of 20 The high-order gradient information and Lp-pseudonorm were introduced to increase the accuracy of the prior knowledge, thereby preserving the edge of the image while suppressing the influence of the small edge on the estimation of the blurring core. Therefore, the proposed high-order total variation with Lp-pseudonorm shrinkage (HTV-Lp) model is defined as follows: , the Lp-norm is defined as while the Lp-pseudonorm is defined as To resolve the HTV-Lp model defined by Equation (3), the ADMM framework [33,39] is used. The complex problem was transformed into several simple decoupled sub-problems using variable substitution.
The split variables can be defined as According to the ADMM principle, dual variables ( 1,2, ,5) (4) can be converted into an unconstrained extended Lagrange function: The high-order gradient information and Lp-pseudonorm were introduced to increase the accuracy of the prior knowledge, thereby preserving the edge of the image while suppressing the influence of the small edge on the estimation of the blurring core. Therefore, the proposed high-order total variation with Lp-pseudonorm shrinkage (HTV-Lp) model is defined as follows: where , and u i (i = 1, 2, · · · , 5) are the balance parameters.

The Proposed Model Solution Based on the Alternating Direction Method of Multipliers
To resolve the HTV-Lp model defined by Equation (3), the ADMM framework [33,39] is used. The complex problem was transformed into several simple decoupled sub-problems using variable substitution. The split variables can be defined as The original problem is transformed into the following equations: where According to the ADMM principle, dual variables Z i (i = 1, 2, · · · , 5) are introduced. The problem shown in Equation (4) can be converted into an unconstrained extended Lagrange function: where X, Y represents the inner product of X and Y, and β i is the penalty factor of the quadratic penalty function.

F Sub-Problem Solving
The sub-function of the sub-problem of F is written as follows, A supplement to Equation (6) can be formulated as follows, The sub-function then becomes: The Fourier transform is performed on the above expression to obtain the frequency domain representation and transform Equation (8) into the following equation: where the symbol • represents the point multiplication operation and X represents the spectrum of X.
The derivative of Equation (9) to F is as follows, Let Lhs and Rhs be defined as Equations (11) and (12), respectively, Then, according to Equation (10), F (k+1) is defined as follows, where i f f t2 denotes fast inverse Fourier transform.
The objective function of the Z 1 subproblem is: We adopted Lp-pseudonorm shrinkage to solve the Equation (14). In contrast to Wang et al., [21] where the value of p was the same, we separated the p values of different gradients and directions of Appl. Sci. 2020, 10, 2533 6 of 19 the image. The advantage of this operation lies in the fact that it can improve the performance of the deblurring image as described in Section 5.3. Lp-pseudonorm shrinkage is defined as shrink p i (ξ, τ) = max |ξ| − τ 2−p i |ξ| p i −1 , 0 (i = 1, 2, · · · , 5). According to the Lp-pseudonorm shrinkage rule, Z 1 can be updated as follows, Similarly, we have: The objective function of the sub-problem of Z i can be written as follows, According to the mountain climbing method [40], where γ is the learning rate. Table 1 presents the pseudocode of HTV-Lp for infrared image deblurring.

HTV-Lp Pseudo-Code
Input: observed image G and fuzzy kernel H; where tol represents the threshold.

Experimental Environment
To demonstrate the superiority of the proposed HTV-Lp model, Figure 2 shows eight different test images downloaded from the publicly available datasets found in http://adas.cvc.uab.es/elektra/datasets/ far-infra-red/, and http://www.dgp.toronto.edu/~{}nmorris/data/IRData/. Infrared images of http: //adas.cvc.uab.es/elektra/datasets/far-infra-red/were obtained using an infrared camera (PathFindIR Flir camera) with 19-mm focal length lens. The size of the "Store.BMP" is 506 × 408 and that of the rest images are 384 × 288. The Gaussian blur and motion blur kernels are generated by MATLAB functions. For example, MATLAB function "fspecial" ('gaussian', B × B × σ) generates a B × B Gaussian blur kernel with a standard deviation of σ. For convenience, the kernel is referred to as (G, B, B, σ).
"fspecial" ('motion', L, θ) generates a motion blur kernel with a motion displacement of L and motion angle of θ, which is referred to as (M, L, θ) for convenience. Furthermore, Gaussian noise standard deviations of 1, 3, and 5 are used.
http://www.dgp.toronto.edu/~nmorris/data/IRData/. Infrared images of http://adas.cvc.uab.es/elektra/datasets/far-infra-red/were obtained using an infrared camera (PathFindIR Flir camera) with 19-mm focal length lens. The size of the "Store.BMP" is 506 × 408 and that of the rest images are 384 × 288. The Gaussian blur and motion blur kernels are generated by MATLAB functions. For example, MATLAB function "fspecial" ('gaussian', B × B × σ) generates a B × B Gaussian blur kernel with a standard deviation of σ. For convenience, the kernel is referred to as (G, B, B, σ). "fspecial" ('motion', L, θ) generates a motion blur kernel with a motion displacement of L and motion angle of θ, which is referred to as (M, L, θ) for convenience. Furthermore, Gaussian noise standard deviations of 1, 3, and 5 are used.

Evaluation Metrics
The main comparison evaluation metrics considered in this study are the PSNR, SSIM, and GMSD, which are expressed in Equations (19) (20) where X and Y refer to the original image and the restored image, respectively. u X and u Y are the average of the sum of the images, and 2 σ X , 2 σ Y represent the variance of X and Y , respectively. Here, k1 and k2 are used to ensure that the SSIM expression is non-zero; in this experiment, k1 is set as 0.01 and k2 as 0.03. where,

Evaluation Metrics
The main comparison evaluation metrics considered in this study are the PSNR, SSIM, and GMSD, which are expressed in Equations (19)-(21), respectively: where X and Y refer to the original image and the restored image, respectively. u X and u Y are the average of the sum of the images, and σ 2 X , σ 2 Y represent the variance of X and Y, respectively. Here, k 1 and k 2 are used to ensure that the SSIM expression is non-zero; in this experiment, k 1 is set as 0.01 and k 2 as 0.03. where, where m r and m d refer to the gradient amplitude of the image in the horizontal and vertical directions, respectively, and c is a constant of small value to guarantee that the denominator is a non-zero number. Larger PSNR values and smaller GMSD values correspond to better image recovery quality. The range of SSIM is (0, 1), where higher values indicate better deblurring performance. The iterative stop condition in the algorithm can be expressed as follows, where tol is set to be 10 −4 .

The Sensitivity of the Parameters
The parameters u i , β i , p i , γ(i = 1, 2, · · · , 5) should be properly selected. We try each parameter tentatively until the iterated image is properly recovered. The ranges of each parameter are determined empirically as u i , β i = 0-20, p i , γ = 0-1(i = 1, 2, · · · , 5). The parameters are adjusted by traversing these ranges with a step size of 0.01. When PSNR achieves a maximum value, the corresponding parameter values are selected as the optimal ones. Figure 3 depicts the effect of u i on PSNR during the adjustment. Notably, when the PSNR reaches its maximum, the optimal selections of u i (i = 1, 2, · · · , 5) vary, hence, it is effective to adjust these parameters separately. Figure 4 depicts the effect of variations in p i on PSNR during the adjustment. When the p i is too large, the image gradient sparsity is not sufficiently strong. However, if p i is too small, the image gradients become too sparse. Thus, the value of p i should be adjusted until the image recovery is satisfactory. As shown in Figure 4, the values of p i are different when PSNR reaches the maximum value under the same experimental conditions. Therefore it is effective to separate the p i of different gradients and directions of image.

Comparison of the Deblurring Performance
The parameters of the six algorithms are adjusted by the traversal to achieve the best deblurring indicators. The test results for the different images are presented in Tables 2 and 3, where the values in bold refer to the optimal indicator values. Table 2 presents the deblurred results obtained by the six different models for degraded images with Gaussian blur (G, 5, 5, 7) and Gaussian noise standard deviations of 1, 3, and 5, respectively. Table 3 presents the results obtained by the six different models for degraded images with motion blur (M, 10, 10) and Gaussian standard deviations of 1, 3, and 5, respectively. Figure 5 depicts the changes in the methods in terms of the evaluation metrics with respect to the increase in the iteration number.
Appl. Sci. 2020, 10, 2533 9 of 19 is too large, the image gradient sparsity is not sufficiently strong. However, if i is too small, the image gradients become too sparse. Thus, the value of i p should be adjusted until the image recovery is satisfactory. As shown in Figure 4, the values of i p are different when PSNR reaches the maximum value under the same experimental conditions. Therefore it is effective to separate the i p of different gradients and directions of image.  (e) Figure 3. Effects of variations on the peak signal-to-noise ratio (PSNR) for the "Car" images degraded by the motion blur kernel (M, 10, 10) with Gaussian noise standard deviations 1, 3, and 5 respectively.

Comparison of the Deblurring Performance
The parameters of the six algorithms are adjusted by the traversal to achieve the best deblurring indicators. The test results for the different images are presented in Tables 2 and 3, where the values in bold refer to the optimal indicator values. Table 2 presents the deblurred results obtained by the six different models for degraded images with Gaussian blur (G, 5, 5, 7) and Gaussian noise standard deviations of 1, 3, and 5, respectively. Table 3 presents the results obtained by the six different models for degraded images with motion blur (M, 10, 10) and Gaussian standard deviations of 1, 3, and 5, respectively. Figure 5 depicts the changes in the methods in terms of the evaluation metrics with respect to the increase in the iteration number.   The performance of the proposed method, illustrated intuitively in Figure 5, and Tables 2 and 3, can be summarized as follows: (1) The performance indicators obtained by the HTV-Lp model are higher than all the other models, thereby demonstrating that the proposed method has better The performance of the proposed method, illustrated intuitively in Figure 5, and Tables 2 and 3, can be summarized as follows: (1) The performance indicators obtained by the HTV-Lp model are higher than all the other models, thereby demonstrating that the proposed method has better deblurring and denoising effects. (2) Observed from Table 2 when restoring the eight degraded images with the motion blur kernel (M, 10, 10) and Gaussian noise standard deviation of 1-5, the HTV-Lp model shows average PSNR values 0.584 dB higher than those of the ATV method, 0.536 dB higher than those of the ITV method, 0.389 dB higher than those of the ATpV method, 0.575 dB higher than those of the FAV4Lp method, and 0.100 dB higher than those of the Hybrid TV method. (3) Observed from Table 3 when restoring the eight degraded images with Gaussian blur kernel (G, 5, 5, 7) and Gaussian noise standard deviation of 1-5, the HTV-Lp model shows average PSNR values 0.531 dB higher than those of the ATV method, 0.391 dB higher than those of the ITV method, 0.300 dB higher than those of the ATpV method, 0.589 dB higher than those of the FAV4Lp method, and 0.175 dB higher than those of the Hybrid TV method. Hence, we conclude that the high-order gradients sparsity of the image is helpful for image recovery and Lp-pseudonorm is better than L1-norm.

Comparison of Visual Effects
The deblurring results for the "Store" are depicted in Figure 6. To better exhibit the effects of the six algorithms, a part of "Store" details in the red rectangles was enlarged. It can be seen that the HTV-Lp model minimized the noise while simultaneously reliving the staircase artifacts in slanted and smooth regions of the image during image deblurring. In addition, the Lp-pseudonorm can depict the sparsity of processed variables more precisely compared to the L1-norm. Figure 7 depicts the single columns extracted by original image, degraded image and deblurring images of ATV, ITV, ATpV, FTV4Lp, Hybrid TV, HTV-Lp. It is found that the HTV-Lp achieves the flattest curve among the compared methods.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 20 deblurring and denoising effects. (2) Observed from Table 2 when restoring the eight degraded images with the motion blur kernel (M, 10, 10) and Gaussian noise standard deviation of 1-5, the HTV-Lp model shows average PSNR values 0.584 dB higher than those of the ATV method, 0.536 dB higher than those of the ITV method, 0.389 dB higher than those of the ATpV method, 0.575 dB higher than those of the FAV4Lp method, and 0.100 dB higher than those of the Hybrid TV method. (3) Observed from Table 3 when restoring the eight degraded images with Gaussian blur kernel (G, 5,5,7) and Gaussian noise standard deviation of 1-5, the HTV-Lp model shows average PSNR values 0.531 dB higher than those of the ATV method, 0.391 dB higher than those of the ITV method, 0.300 dB higher than those of the ATpV method, 0.589 dB higher than those of the FAV4Lp method, and 0.175 dB higher than those of the Hybrid TV method. Hence, we conclude that the high-order gradients sparsity of the image is helpful for image recovery and Lp-pseudonorm is better than L1norm.

Comparison of Visual Effects
The deblurring results for the "Store" are depicted in Figure 6. To better exhibit the effects of the six algorithms, a part of "Store" details in the red rectangles was enlarged. It can be seen that the HTV-Lp model minimized the noise while simultaneously reliving the staircase artifacts in slanted and smooth regions of the image during image deblurring. In addition, the Lp-pseudonorm can depict the sparsity of processed variables more precisely compared to the L1-norm. Figure 7 depicts the single columns extracted by original image, degraded image and deblurring images of ATV, ITV, ATpV, FTV4Lp, Hybrid TV, HTV-Lp. It is found that the HTV-Lp achieves the flattest curve among the compared methods.

Discussion and Conclusions
Deblurring infrared images plays a significant role in an infrared imaging system. To improve the performance of deblurring images, the propoed HTV-Lp model combines first-order and highorder gradients of images with Lp-pseudonorm shrinkage. The ADMM algorithm is used to split the proposed model into several decoupled sub-problems. In the process of solving thee problems, the convolution and FFT theorem are applied to image deblurring to avoid large computational complexity. The comparison of HTV-Lp with existing ATV, ITV, ATpV, FTV4Lp, and Hybrid TV models shows that HTV-Lp obtained the average highest PSNR and SSIM and lowest GMSD values of all methods and successfully mitigated staircase artifacts.
A limitation of this study is that non-local patch similarity is not fully considered in the proposed model. Additionally, there is room for further improvement with regard to the speed of the HTV-Lp algorithm within the framework of the accelerated ADMM. Thus, in our future work, we will focus on improving the performance and efficiency of the proposed method.

Discussion and Conclusions
Deblurring infrared images plays a significant role in an infrared imaging system. To improve the performance of deblurring images, the propoed HTV-Lp model combines first-order and high-order gradients of images with Lp-pseudonorm shrinkage. The ADMM algorithm is used to split the proposed model into several decoupled sub-problems. In the process of solving thee problems, the convolution and FFT theorem are applied to image deblurring to avoid large computational complexity. The comparison of HTV-Lp with existing ATV, ITV, ATpV, FTV4Lp, and Hybrid TV models shows that HTV-Lp obtained the average highest PSNR and SSIM and lowest GMSD values of all methods and successfully mitigated staircase artifacts.
A limitation of this study is that non-local patch similarity is not fully considered in the proposed model. Additionally, there is room for further improvement with regard to the speed of the HTV-Lp algorithm within the framework of the accelerated ADMM. Thus, in our future work, we will focus on improving the performance and efficiency of the proposed method.