High energy ﬂash X-ray image restoration using region extrema and kernel optimization

The quality of high energy ﬂash X-ray images is crucial to the high-precision diagnosis of object density. High energy ﬂash X-ray radiography is susceptible to the system blur, which usually causes the poor quality of static images. In response to this, a novel restoration algorithm using region extrema and kernel optimization (REKO) is presented. Based on the observation that the region extrema distribution of blurred high energy ﬂash X-ray images deviates from opposite ends of image grey domain, the sparseness-inducing prior for regularizing image region extrema is applied to construct the restoration model. Considering the sparse characteristics of blur kernels, the sparseness-inducing regularization is incorporated to constrain blur kernels in the restoration model. The non-convex and non-linear objective function is gradually minimized through energy alternating minimization and dually linear approximation. Furthermore, a continuity enforced kernel optimization algorithm is proposed to estimate more accurate blur kernels. The discontinuous kernel elements are suppressed by extracting the main structure of blur kernels and constructing kernel continuity function in cross windows. Experimental results demonstrate that our algorithm can more accurately estimate blur kernels and achieve restoration results with sharper edges on high energy ﬂash X-ray images.


INTRODUCTION
In high energy flash X-ray radiography, the system blur which consists of light source fuzzy and detector blurring has a negative effect on image quality. High energy flash X-ray image restoration aims to estimate the system blur kernel and to recover the sharp image from the blurred image. When the blur kernel is uniform and spatially invariant, a general formulation of a blurred image is expressed as where x is the sharp image, k is the blur kernel, y denotes the blurred and noisy observation, n denotes the additive noise, and * is the convolution operator. Blind image restoration in high energy flash X-ray radiography is a highly ill-posed problem as infinite pairs of x and k can give rise to the same as y. In terms of the involved technologies, blind image restoration algorithms can be divided into three categories, namely, strong edge prediction-based, neural network-based and Bayesian framework-based algorithms.
Generally speaking, strong edge prediction-based models can extract strong edges and be instructed for image restoration via the obvious step change of the intensity at image edges. To ensure the preservation of sharp edges, heuristic combination steps are utilized in the restoration procedure for eliminating the noises generated by predicting edges, such as adaptive thresholds [1,2], bilateral filters [2], shock filters [3] and guided filters [4]. Strong edge prediction-based algorithms can achieve accurate restoration results and reduce computational complexity. These algorithms may not always effectively predict edges, especially when the observed image is highly blurred or contains few textures. Instead of predicting strong edges, neural network-based image restoration algorithms [5][6][7] usually need network training and can accurately recover the blur kernel and the real image. The plug-and-play framework has been recently proposed in the neural network-based algorithm [8] or other categories of image restoration algorithms [9,10]. It can implicitly define the prior via off-the-shelf denoiser and solve such ill-posed problem. These algorithms perform well on blurred images including blurred low-textured images [11].
Formulated within the Bayesian framework, maximum a posteriori (MAP) based blind image restoration algorithms have the advantages of good flexibility and low computational complexity because the image dataset training or image pre-processing step is not required. Prior knowledge can be divided into image prior and blur kernel prior, which have been introduced in these algorithms to favour sharp images rather than blurred ones. For image prior, most studies generally utilize the heavy-tailed distribution of image gradient [12,13] or other types of sparse statistics, such as sparse channel prior [14,15], as strong sparsity prior to regularize the restoration model. These algorithms are widely used in visible image restoration to obtain sharp restoration results. As one type of greyscale images, high energy flash X-ray images usually have the property of single edges, low texture, and concentrated grey distribution. To some extent, the characteristics of high energy flash X-ray images are similar to those of text images. Pan et al. [16] studied the sparse distribution of intensity and gradient and utilized the sparseness-inducing norm to constrain image for text image restoration. Unfortunately, the global intensity-based sparseness-inducing prior is not suitable for describing the concentrated grey distribution of high energy flash X-ray images. It is likely to cause the elements dispersion of the estimated blur kernel and the regional grey distortion of the restored image.
With regard to blur kernel prior, the quadratic prior based on the continuity of blur kernels have been widely introduced for visible image restoration [14][15][16][17][18]. These methods have some defects proved by Fang et al. [19] that the quadratic constraint neglects the strong sparsity of blur kernels and complicates the restoration processing without significant performance improvement. The value of blur kernel sparsenessinducing norm is numerically larger than that of the kernel quadratic constraint. This suggests that the sparseness-inducing norm of blur kernels has a higher degree of numerical discrimination. Therefore, the kernel sparseness-inducing prior (KSP) can effectively constrain the blur kernel in the restoration model. Only using the sparseness-inducing constraint of blur kernels may generate discontinuous kernel elements, which would have a negative effect on restoration results. The mixture of the quadratic and sparseness-inducing norm is inefficient because the weight parameters need to be given manually. Therefore, an efficient kernel optimization algorithm is expected.
In this paper, we propose an effective restoration algorithm for high energy flash X-ray images. Within the MAP framework, we build a novel objective function for image restoration. For the choice of image prior, we use the sparsity of image region extrema whose distribution deviates from opposite ends of image grey domain under the influence of blurring to match the properties of high energy flash X-ray images. The blur kernel is regularized by the sparseness-inducing norm with a high degree of numerical discrimination. That can enforce the strong sparsity of blur kernel and improve the restoration quality of images. Then we introduce energy alternating minimization and dually linear approximation to handle this non-convex and nonlinear objective function. Subsequently, a kernel optimization algorithm enforced by continuity is developed. The discontinuous kernel elements are suppressed by extracting the main structure of the blur kernel and defining the kernel continuity function in cross windows, which can estimate more accurate blur kernels. The proposed algorithm can be applied to high energy flash X-ray image restoration, which performs better than the state-of-the-art restoration algorithms.
The rest of this paper is organized as follows. Section 2 introduces a broad-scale overview of the related work. The restoration model and the blur kernel optimization algorithm are presented in Sections 3 and 4, respectively. Section 5 demonstrates the performance of the proposed algorithm. The conclusion is given in Section 6.

RELATED WORK
Blind image restoration is a highly ill-posed inverse problem, the aim of which is to estimate the blur kernel and to recover the sharp image from the given blurred images. Most of the existing blind image restoration algorithms are based on the Bayesian framework with flexibility and interpretability. Choosing appropriate priors for images and blur kernels is crucial to the restoration performance of such algorithms. Therefore, in this section, we focus on the related work of Bayesian-based restoration algorithms. As for image prior, Fergus et al. [20] utilized mixed Gaussian function to learn the image gradient prior within variational Bayesian (VB) framework, while the work in [21] shows that the sparse gradient prior-based algorithm still makes image blurred rather than sharp. To overcome this limitation, Thakur et al. [22] introduced a 1 −norm data fidelity term combining with the total variation filtering under the VB framework, which can restore more image details.
Compared with the VB-based algorithms, MAP-based algorithms show greater flexibility and faster computation efficiency. According to the observation that gradients of visible images obey the heavy-tailed distribution, various algorithms [12,13] are widely proposed with the interactions model between pixel pairs. The non-local self-similarity prior is introduced in blind image deblurring domain, which can reduce the inaccuracy in modelling adjacent image pixels and has been widely used in image denoising problems [23,24]. In the work by Ren et al. [17], the low-rank prior is combined with the gradient prior to estimate an accurate blur kernel by using a weighted nuclear norm minimization, which assigns different singular values with different weights.
Apart from the low-rank prior and the gradient prior, sparse representation is developed and incorporated into the regularization for blind image restoration [14][15][16]. Pan et al. [14] proposed the dark channel prior (DCP) by applying 0 −norm regularization at the dark channel, and the method generates latent images with more details and clear textures. The DCP algorithm was first applied in single-image dehazing by He et al. [25]. It is based on the observation that in most visible images, at least one of the colour channels possesses some pixels whose value tends to zero intensity. In the work by Pan et al. [14], the DCP model is modified that the dark channel of visible images is sparse instead of zero. By enforcing the sparsity of the dark channel, this work has excellent restoration results on benchmark datasets of visible images. Furthermore, Yan et al. [15] introduced the bright channel prior (BCP) to eliminate the limitation, which the DCP may not effectively restore the intermediate latent image if no dark pixels at the image. These channel priors-based algorithms are robust on visible images. They fully consider that each local region of natural images exists some pixels with extreme intensities.
High energy flash X-ray image usually has the properties of single edges, low texture, and concentrated grey distribution, which is similar to text image. Cao et al. [26] adopted text-specific multiscale dictionaries to introduce image prior on the text fields, which can effectively deal with the different scales of strings at a blurred image. Liu et al. [27] presented a surface-aware strategy based on the intrinsic geometry for image restoration, which restores saturated text images well. Pan et al.
[16] further formulated the sparseness-inducing constraint for both intensity and gradient distribution of images. However, the global intensity-based sparseness-inducing prior is simple in structure and is not suitable for describing the greyscale distribution characteristic of high energy flash X-ray image, which would lead to kernel dispersion and regional grey distortion in the restoration results.
Blur kernel prior also plays an important role in blind image restoration. Most existing algorithms research the properties of blur kernels and constrain the blur kernel via the quadratic prior [13][14][15]. The analysis in [19] shows the value of the blur kernel quadratic constraint is much smaller than that of the sparsenessinducing constraint. Therefore, the contribution of the kernel quadratic regularization for the restoration model convergence and the restored image quality can be negligible. Accordingly, Fang et al. [19] introduced the sparseness-inducing prior of blur kernels to improve the quality of the restoration results. Even so, only using the sparseness-inducing prior for blur kernels still generates discontinuous kernel elements, which would generate bad points in the restoration results.
To achieve high blur kernel estimation accuracy for high energy flash X-ray image, the contributions of this paper are as follows: (1) According to the sparsity of blur kernels and region extrema of high energy flash X-ray image, a novel restoration model based on the sparseness-inducing prior is introduced to improve the accuracy of restoration results. (2) Energy alternating minimization and dually linear approximation are adopted to effectively deal with the non-convex and non-linear restoration model. (3) An efficient blur kernel optimization algorithm enforced by continuity is proposed, which is applicable to deal with Gaussian and motion blur kernels. By extracting the main structure of blur kernels and constructing the kernel continuity function in cross windows, the discontinuous kernel elements are suppressed and the blur kernel estimation accuracy is better guaranteed.

BLIND IMAGE RESTORATION MODEL
In this section, the proposed restoration model is described, which includes three parts: (1) the restoration model formulation, (2) blur kernel estimation, and (3) intermediate image estimation. We introduce the quadratic norm to constrain the data fidelity term and the sparseness-inducing norm to constrain image region extrema, image gradients and blur kernels. Then energy alternating minimization and dually linear approximation are utilized to handle this restoration model.

Restoration model formulation
From the probabilistic point of view, simultaneously estimating the intermediate image and the blur kernel by Equation (1) is equivalent to solving the standard MAP estimation problem [28]. We proposed a novel model based on the MAP framework for solving x and k: where P (k) and P (x) denote kernel prior and image prior probability; p(k) = − log P (k) and p(x) = − log P (x) denote the prior distribution of k and x, respectively; , are weight parameters; and the term ‖x * k − y‖ 2 2 is the data fidelity term, which can ensure that the convolution result of the restored image and blur kernel should be similar to the observed blurred image.
The image prior and the blur kernel prior are crucial to the regularization which is an ill-posed inverse problem, also crucial to blur kernel estimation and latent image recovery. Different types of prior constraint represented by different matrix norms have been proposed [29]. Although the features of high energy flash X-ray images are similar to those of text images, the global intensity-based sparseness-inducing constraint has poor accuracy in estimating the blur kernel of high energy flash X-ray images. In order to deal with this problem, we introduce image region theory, analyse the feasibility of the greyscale form of extreme channel prior [15] at high energy flash X-ray images, and define the region extrema prior in greyscale.
The region extrema of an image x, including minimum and maximum, are formulated respectively as where p denotes pixel location and Ω(p) denotes an image region patch. For most sharp images with continuous greyscale, region extrema tend to be the maximum and minimum of image greyscale domain. The region maximum becomes smaller, while the region minimum becomes bigger after blurring.
Referring to Equation (1), a blurred image y is obtained by the convolution of the sharp image x and the blur kernel k. For discrete signals (images), the convolution is defined as the sum of the product of one and another signal after inversion and shifting: where z is the pixel location of the blur kernel k, Ω k denotes the kernel domain, s is the size of the kernel, and ⌊⋅⌋denotes the floor rounding operator. In addition, blur kernels should meet the constraints of k(z ) ≥ 0 and ∑ z∈Ω k k(z ) = 1, which are called blur kernel normalization. Therefore, the blurred image y can be regarded as the weighted sum of local pixels at the sharp image x. The following formula can be obtained: We conclude that the region minimum of blurred images is higher than that of sharp ones, and the region maximum of blurred images is lower than that of sharp ones: For the convenience of calculation, T is defined as the maximum value at image greyscale domain and the region maximum in pixel p can be accordingly calculated byT − R max (x, p). The non-zero intensity pixels in image region extrema are counted by 0 −norm to show the strong sparsity of region extrema in sharp images: The following conclusion can be concluded that the value of the sparseness-inducing constraint of region minimum at a sharp image is lower than that at a blurred one. The corresponding region maximum is shown as follows: We further validate the above inference by computing the region extrema distribution of French test object (FTO) optical path image and high energy flash X-ray static images under different blur kernels. As shown in Figure 1, region extrema of sharp images are concentrated at the opposite ends of image grey domain, while those of blurred images deviate. We apply the region extrema sparsity as priors by using ‖R min (x)‖ 0 and ‖T − R max (x)‖ 0 in the restoration model.
As discussed above, the sparsity of image region extrema can be regularized by 0 −norm. The sparseness-inducing regularization for image gradients is applied to retain large gradients and remove tiny ones: where , are weight parameters; In terms of blur kernel prior, most existing algorithms often use the quadratic prior, but such a regularizer has less assistance. To ensure the accuracy of kernel estimation, the mixture of the quadratic and sparseness-inducing norm is introduced but the weight parameters which are given manually may lead to an additional optimization budget and make the deblurring iteration difficult to converge. The kernel continuity using the quadratic norm has been already partially included in the data fidelity term. Compared with the value of the kernel quadratic norm, the value of the kernel sparseness-inducing norm is bigger. Therefore, the kernel sparseness-inducing prior is used in our restoration model to take the kernel strong sparsity into consideration.
Taking these priors together, Equation (2) is rewritten as It is difficult to solve Equation (15) directly because unknown variables are more than the known ones. We use an efficient alternating minimization algorithm based on half-quadratic splitting. Intermediate latent image and blur kernel are estimated alternatively by assuming one of them is known. Then the optimization problem becomes two sub-problems: The sub-problem Equation (16) is composed of terms with variable k, which assumes variable x is known. By fixing the value of variable k, Equation (17) denotes the sub-problem which comprises all x-included terms.

Blur kernel estimation
In this subsection, we estimate the blur kernel k according to Equation (16). We use the image gradient instead of the image intensity with the 2 −norm constraining data fidelity term to employ more accurate kernel estimation. Therefore, we can transform Equation (16) intô We solve this sub-problem by alternating minimization which is described in Algorithm 1. By introducing auxiliary variable h, the alternative objective function can be rewritten as Given the auxiliary variable h, we use the alternating minimization algorithm to solve this problem. Equation (19) can be rewritten aŝ Equation (20) can be regarded as the least square problem to be solved. Therefore, the blur kernel k can be calculated by fast Fourier transform (FFT) at each iteration: where  (⋅) and  −1 (⋅) denote the FFT and inverse fast Fourier transform (IFFT), and• is the complex conjugate operator. Calculating  (∇x) (∇y) and  (∇x) (∇x) by FFT in advance, we can efficiently solve the blur kernel k. By fixing the variable k, Equation (19) can be rewritten aŝ The above function can be handled by pixel-wise minimization algorithm as

Intermediate image estimation
In this subsection, we fix the blur kernel k and optimize the latent sharp image x. By introducing new auxiliary variables e, l and g = (g h , g v ) T , which are corresponding toR min (x), T − R max (x) and ∇x, respectively, we can rewrite the sub-problem Equation (17) as the following energy minimization function: where , and are positive penalty parameters. When , and tend to be infinity, solving Equation (24) is equivalent to solving Equation (17). By using a half-quadratic splitting strategy to calculate the variable x, e, l and g, we can obtain the solution of Equation (24).
Given the intermediate latent image x, the sub-problems with the variables g, e and l are solved respectively. The corresponding sub-optimization function is All the three functions mentioned above belong to the pixelwise minimization problem. We compute the optimal solutions as As R min (x) and T − R max (x) are non-linear, we introduce M min and M max as equivalent linear operators to calculate these difficult problems. M min and M max are fundamentally mapping matrices, which map the pixels to the region minimum and maximum, respectively: For a sharp image, M min x = R min (x) and M max x = R max (x). In the process of computing the sub-problem Equation (24) 1 Initialize blur kernel k 0 ; 2 Compute the max pyramid level L max according to [2]; 3 for ilevel = 1 ∶ L max 4 Down-Sampling y to the current image pyramid to gety iiter ; 5 for iiter = 1 ∶ N max 6 Solve for x by Algorithm 2; 7 Solve for k by Algorithm 1; 8 end for; 9 Up-Sampling k to fit the next level; 10 end for; 11 Extracting the main-structure S (k) of blur kernel k; 12 Solve for kernel continuity function C (h, v) using (37); 13 Suppress discontinuous pixels in sliding cross windows; 14 Update and normalize k.
15 Restore x by a non-deblurring algorithm with y and k.
known, we modified the sub-problem in Equation (24) as follows to calculate image x: (31) Similar to solving the kernel k, we solve the intermediate image x by FFT at each iteration: where ; ∇ h and ∇ v denote horizontal and vertical differential operations, respectively. The main steps of our intermediate image estimation are shown in Algorithm 2.

BLUR KERNEL OPTIMIZATION
In Section 3, the blur kernel is constrained by the sparsenessinducing prior instead of the common quadratic prior, which can ensure that the blur kernel tends to be sparse and stable at the iterative process. The constraint is used because the value of the kernel quadratic norm is much smaller than that of the kernel sparseness-inducing norm with a higher degree of numerical discrimination. However, due to the blur kernel continuity without fully considered, the solution process of the restoration model may give rise to some discontinuous kernel elements. Therefore, we extract the main structure of blur kernels with considering continuity constraint and construct a continuity function in sliding cross windows to suppress discontinuous elements, which can obtain a more accurate blur kernel.

Main structure extraction of blur kernel
As a simple greyscale image, the blur kernel can be represented by a two-dimensional matrix. Assuming that the real blur kernel is known, we can establish the following relation of the estimated blur kernel by location information: wherek(h, v) denotes the estimated blur kernel, k(h, v)denotes the real blur kernel, and d (h, v) denotes the estimation error. The more similar the estimated blur kernel is to the real one, the smaller the gap of position integral between them: The background of blur kernels is single and effective pixels are mainly concentrated in the foreground, which is different from complex images. Therefore, the difference in the position integral between the estimated blur kernel and the real one can be expressed by the difference of each pixel intensity and the blur kernel shape in the foreground. We can optimize the estimated blur kernel by minimizing the gap of the main structure S (k) between the two blur kernels: where S (k) and S (k) contain the shape and strength of the real blur kernel and the estimated one, respectively: where In order to minimize Equation (35), we need to obtain the main structure which is most similar to the real blur kernel. We rely on the skeleton extraction algorithm [30] to effectively obtain the trajectory of the estimated blur kernel. For Gaussian blur kernels and square blur kernels, the use of the skeleton extraction algorithm alone cannot completely extract the main structure as shown in Figure 2(b). To this end, we propose a simple algorithm to obtain the boundary of the blur kernel. The non-zero pixels in the blur kernel are be traversed. We judge whether some pixels with intensities of zero (i.e., the background of blur kernels) exist in the neighbourhood of the traversed pixel. If there exist, define it as blur kernel boundary, otherwise continue to traverse. The algorithm can compensate for the deficiency that the skeleton extraction algorithm cannot effectively extract the main structure of Gaussian and square kernel. The skeleton and boundary of the kernel together constitute the main structure S (k) of the kernel. It is beneficial for extracting the main structure and suppressing the discontinu- (a) Intermediate estimated blur kernel, and discontinuous kernel pixels are in the red patch, (b) the process of discontinuous pixels suppression ity of blur kernel elements to integrating the two algorithms above. Figure 3 shows the extraction of the main structure for the motion and Gaussian blur kernel. For the motion kernel, both the proposed algorithm and KSP by Fang et al. [19] extract the main structure of the blur kernel. For Gaussian kernels, the extraction effect of KSP is not as good as the proposed algorithm. Fang et al. [19] minimize the gap of the blur kernel structure in the camera shake time only by skeleton extraction algorithm, and extract the skeleton similar to the main structure of motion kernel. Therefore, it is only effective for motion blur kernel. We minimize the difference of blur kernel structure in different pixel positions, and add a boundary acquisition step based on the skeleton extraction. That is suitable for different types of blur kernels, including Gaussian kernels, motion kernels and so on.

Discontinuous kernel elements suppression
In this subsection, the pixels of the main structure in blur kernel are traversed in sliding cross windows, and then the kernel continuity function is constructed. By taking the nonzero pixels number of blur kernel in each window as a continuity measure, we delete pixels whose continuous function is smaller than the width w of the sliding cross window until all pixels are processed to update for the main structure of the blur kernel. The blur kernel continuity function is defined as where C (h, v) denotes the number of non-zero pixels in the crossing window; h and v denote the horizontal and vertical directions, respectively. Assuming the width of the cross window w = 2, the pixel with C (h, v) = 6 is retained, which is larger than w , while the pixel with C (h, v)= 1 is removed as shown in Figure  3(b). After discontinuous suppression of shape(k), strength(k) is updated and normalized. The whole suppression process of kernel optimization is shown in Figure 4. In Figure 4(b), the green pixels denote the main structure, the orange patch denotes the cross window on the current red pixel, and the yellow arrows denote cross traversal. In Figure 4(c), the kernel before optimization consists of pixels in orange and blue. The kernel after optimization consists of orange pixels. Algorithm 3 summarizes the specific steps of our algorithm based on region extrema and kernel optimization. The overall framework of the proposed algorithm is depicted in Figure 5. Our algorithm performs the following steps: (i) half-quadratic splitting of the proposed restoration model, (ii) multi-scale iterative process, (iii) kernel optimization, and (iv) non-blind deblurring.

EXPERIMENTAL RESULTS
Because of the confidentiality of most high energy flash X-ray images, our experiments are mainly carried out on FTO image and non-confidential high energy flash X-ray images to validate the effectiveness of our algorithm (abbreviated as REKO).

FIGURE 5
The procedure of the proposed algorithm

FIGURE 6 Experimental blur kernels. (a)-(d) Motion A-D, (e)-(h) Gaussian A-D, (i)-(j) square A, B
Since the proposed constraint of image region extrema and blur kernel sparsity are also suitable for describing the blurring process of visible images. We further evaluate the restoration quality of our algorithm in comparison with other existing algorithms on blurred visible images. Experiments are carried out on both synthetic datasets and real blurred images. The images of the synthetic datasets are generated by ten blur kernels illustrated in Figure 6 and added by Gaussian noise with variance 0.01 2 ‖x true * k‖ 2 2 , where x true is the real sharp image.
In all experiments, some parameters are empirically set as follows: = = = 0 .004, = 0 .002, max = max = max = 1, max = 8. The maximum number of iterations in our restoration model is set as: N max = 5. The sensitivity analysis with respect to these parameters will be discussed in Section 5.5.

Effectiveness of region extrema constraint
As mentioned in Section 3.1, we consider the features of high energy flash X-ray images and utilize image region extrema prior to obtain more accurate restoration results. In this case, the effectiveness of the region extrema constraint is mainly evaluated by comparing our algorithm with KSP. For the experimental images, the text image dataset in [19] is applied to the KSP and our algorithm to carry out a fair comparison.
We analyse the convergence rate and the estimation accuracy by using the kernel similarity proposed by Hu et al. [31] as the metric. The kernel similarity utilizes the maximum response of where k andk denote the real blur kernel and the estimated one, respectively; is the possible shift between the two blur kernels; and represents element coordinates. Experiments are carried out within a multi-scale framework with five iterations set for each scale. The blur kernel gradually converges during each iteration and serves as the initial state value of the next scale. The curve of kernel similarity is shown in Figure 7(a), which can relatively reflect the ability of restoration algorithms in describing the convergence rate and the blur kernel estimation accuracy. Figure 7(b) calculates the average similarity of each blur kernel. Our algorithm can estimate more accurate blur kernels than KSP. Average peak signal to noise ratio (PSNR) and structural similarity index measure (SSIM) of the restoration results can be seen in Figure7(c) and (d). These evaluation metrics exhibit the fact that our algorithm acquires a better quality in most images except for some text images containing sparse toning. The reason is that the binarization characteristics of these text images are more obvious and the use of global intensity prior has completely satisfied the restoration needs for such images. We construct the restoration model with image region constraint and extremely emphasize the detail of images, resulting in a poor restoration of these images. The reason why our algorithm can be applied to the text image restoration and achieve great performance is that the proposed image prior is the same as that in KSP when image region length r is set to 1. When r ≥ 1, our restoration model more fully considers the regional properties of images. It can be understood that the region extrema prior can reasonably constrain the restoration model and obtain a fast convergence rate and high accuracy of restoration.

Analysis of blur kernel optimization algorithm
To analyse the effectiveness of blur kernel optimization, we compare the restoration results involving our kernel optimization algorithm. The experiment is carried on seven images from the dataset in [32]. Figure 8 ?? depicts the PSNR and SSIM of restoration results involving optimization of different blur kernels. Through three types of blur kernels, the restoration results exhibit the fact that our algorithm with kernel optimization performs better than without kernel optimization, because of the main structure extraction and discontinuous kernel elements suppression in cross windows. It can be seen that our blur kernel optimization algorithm has a satisfactory improvement on restoration results especially under the influence of Gaussian blur.  For further qualitative analysis, Figure 9 shows the restoration results of four images involving blur kernel optimization. Our algorithm before blur kernel optimization basically removes blur and restores image contour except for image details. Our blur kernel optimization algorithm assists in restoring sharp edges and abundant details. It is clear that our algorithm provides greater improvement of restoration results than that without blur kernel optimization.

Synthetic datasets
To have an insight into the effectiveness of our algorithm, the experiment is carried on two different visible image datasets: the dataset of Sun et al. [33] and the dataset of Lai et al. [34]. Each image is blurred by Motion A, Gaussian C, Square B in Figure  6 with additive Gaussian noise of variance 0.01 2 ‖x true * k‖ 2 2 . These images are utilized to compare the restoration quality of our algorithm with those of Levin et al. [28] (EMLO), Hu et al. [31] (GR), Sroubek et al. [35] (MC), Ren et al. [17] (ELR), Pan et al. [14] (DCP), Yan et al. [15] (ECP), Liu et al. [36] (SGF), Zhang et al. [37] (TSI) and our algorithm before kernel optimization.
We use the error ratio proposed by Levin et al. [38] to evaluate restored images. The error ratio computes the difference between the restored image x and the ground truth image x gt over the difference between the non-deblurred image x k with the ground truth kernel k: where represent element coordinates. Figure 10 plots the cumulative distribution of error ratio (e.g., the success rate at horizontal axis 3 denotes the percentage of test examples achieving error ratio below 3). It is clear that the restoration success rate of our algorithm without kernel optimization is always higher than the other eight algorithms at the error ratio below 4.5. After kernel optimization, the success rate with our algo-

FIGURE 11
Average kernel similarity of restoration results on each category of images rithm increases by about 3%. Obviously, images restored by our algorithm where the error ratio is less than 1.5 account. Furthermore, to observe the accuracy of blur kernel estimation, we compare the average kernel similarity on the visible image datasets. The similarity of estimated blur kernels for different categories of images and kernels is shown in Figure 11 and Table 1, respectively. The metrics reveal the positive effects of our algorithm in estimating the blur kernel. Figure 12 specifically shows the images and blur kernels restored by the mentioned nine algorithms. The typical difference can be seen from the results of four visible images degraded by Motion A, Gaussian A and Square B. The typical difference can be seen from the restoration results on four visible images degraded by Motion A, Gaussian A and Square B. GR selected regions from the blurred image for blind restoration via fast motion deblurring [1]. Due to the uncertainty of good region selection, the restoration results by GR are barely satisfactory in Figure 12. EMLO generates the relatively poor quality of restored images. The main reason is that though the shape of the estimated blur kernel is accurate, most pixels of the estimated blur kernel deviate

FIGURE 12
The restoration results on four challenging images. (a) Blurred image, (b) the magnified region of (a), (c) the ground truth for (b), (d) the results by EMLO [28], (e) the results by GR [31], (f) the results by MC [35], (g) the results by ELR [17], (h) the results by DCP [14], (i) the results by ECP [15], (j) the results by SGF [36], (k) the results by TSI [37], (l) the results by our algorithm without kernel optimization, (m) the results by our algorithm (REKO) from the centre. MC, ELR and TSI produce relatively severe noises and ringing artefacts in restored images. ECP and SGF yield smooth performance but with uneven brightness. In contrast, DCP and our algorithm can well restore these images. However, some elements of blur kernel estimated by DCP are divergent. Our algorithm can obtain more accurate blur kernels and images, especially after kernel optimization. Average PSNR and SSIM for different categories of images and blur kernels can be seen in Figure 13 and Table 2 for further comparison, which reveals the positive effects of our algorithm in restoring images. It can be seen that the restoration results by our algorithm have higher PSNR and comparable SSIM than those by other algorithms.
The gradient change of the FTO transmittance image is not obvious, so we take the negative logarithm of the blurred FTO transmittance image to expand the gradient change and obtain the blurred FTO optical path image for restoration. When applied to the blurred FTO optical path image, TSI cannot effectively converge to the centroid by using the K-means clustering algorithm, hence it cannot determine the middle grey domain of the input image. Thus, TSI is not considered in the  Table 3. With regard to these metrics, our algorithm performs better than other algorithms, effectively alleviating the influence of blur on the optical path image. As shown in Figure 14, the restoration result of the FTO optical path image blurred by Gaussian kernel D is further compared to the other algorithms. In terms of blur the results by EMLO [28], (e) the results by GR [31], (f) the results by MC [35], (g) the results by ELR [17], (h) the results by DCP [14], (i) the results by ECP [15], (j) the results by SGF [36], (k) the results by our algorithm without kernel optimization, (l) the results by our algorithm (REKO) for more than 50%, while those of other algorithms and ours without kernel optimization are less than 50% kernel estimation, the blur kernel estimated by our algorithm is the closest to the real one. From the image restoration results, it can be seen that the results reveal a comparable quality level. The restored image with our algorithm has a slightly sharp boundary between the copper and the background, while the restored image with other algorithms generates a few artefacts and noises. It is clear that our algorithm has a positive effect on restoring blurred FTO optical path images.

High energy flash X-ray static images
High energy flash X-ray image restoration is mainly aimed to estimate the blur kernel of the radiography system and obtain a visually sharp static projection image. To evaluate the effectiveness of our algorithm in practical scenes, high energy flash X-ray static images are used for the restoration experiments in comparison with the eight algorithms mentioned in Section 5.3. Figure 15 shows that the restoration result of our algorithm generates clearer edges and fewer ringing artefacts than that of other algorithms. Due to the unavailability of the estimated blur kernel, ELMO and GR produce severe ringing artefacts in restored images. The blur kernels estimated by ELR and SGF have many interferential elements, which lead to severe artefacts in restored images. The restored image with DCP has relatively sharp edges but still possesses a few artefacts. The reason is that the dark channel prior proposed by DCP cannot effectively constrain the image regions with gentle grey changes as same as constraining the edge with obvious grey changes. ECP and TSI can generate better restoration results through the assistance of extreme channel prior and tri-segment intensity prior, respectively. However, the blur kernel quadratic prior used in ECP and TSI cannot reflect the strong sparsity of blur kernels. As a consequence, some longitudinally interferential elements are generated in the estimated blur kernel, which have a passive influence on image restoration. In contrast, the blur kernel estimated by our algorithm has a symmetric tail, especially after blur kernel optimization. This reflects the properties of system blur kernel in flash X-ray radiography. In this case, the restored images of our algorithm achieve the best quality among those of all mentioned algorithms.
It can be seen in Figure 16 that the blur kernels estimated by these algorithms reveal various estimation quality levels. MC yields a comparable restoration quality but still generates some artefacts of the restored image. Multiple blurred images of the same scenes, which MC requires, cannot be provided in high energy flash X-ray radiography. This may lead to the limitation of MC applied to such image restoration. Yet the blur kernel estimated by SGF has obvious lateral interference elements, which may lead to the worst restoration result on the image. The restoration results of ELR and SGF possess serious artefacts and those of other comparative algorithms have white noises. Compared with these algorithms above, our algorithm (a) Blurred image, (b) the magnified region of (a), (c) the results by EMLO [28], (d) the results by GR [31], (e) the results by MC [35], (f) the results by ELR [17], (g) the results by DCP [14], (h) the results by ECP [15], (i) the results by SGF [36], (j) the results by TSI [37], (k) the results by our algorithm without kernel optimization, (l) the results by our algorithm (REKO)

FIGURE 16
Restoration results on the second group of high energy flash X-ray static images. The yellow dashed lines outline part of the restoration error. (a) Blurred image, (b) the magnified region of (a), (c) the results by EMLO [28], (d) the results by GR [31], (e) the results by MC [35], (f) the results by ELR [17], (g) the results by DCP [14], (h) the results by ECP [15], (i) the results by SGF [36], (j) the results by TSI [37], (k) the results by our algorithm without kernel optimization, (l) the results by our algorithm (REKO) can restore sharper high energy flash X-ray images. As shown in Figure 17, compared with our algorithm without kernel optimization, our algorithm can estimate a more symmetrical tail of the blur kernel and restore fewer artefacts of the image. Benefited from the proposed kernel sparseness-inducing prior and blur kernel optimization algorithm, the boundaries of the restored image are sharper compared with those restored by other algorithms. It is clear that our algorithm assists in providing more reasonable blur kernel and more satisfactory restoration results of high energy flash X-ray image.

Analysis and discussion
In this subsection, we discuss the sensitivity problem of our algorithm from four aspects: convergence property, computational complexity, parameter analysis and the effect of preprocessing.

Convergence property
As the proposed restoration model is non-convex and nonlinear, a natural question is whether the proposed algorithm converges well or not. We quantitatively evaluate the convergence properties of our method on the benchmark dataset by Sun et al. [33]. Figure 18 shows the energy computed from Equation (15) and the average kernel similarity during iterations. It can be seen that our algorithm has converged after 35 iterations. The reason why our algorithm converges well may be that the proposed restoration model Equation (15) is transformed into two convex sub-problems of x and k by half-quadratic splitting. The introduction of quadratic terms and the utilization of energy alternating minimization are helpful to approximate the two subproblems to the global optimal. During the iterative steps, the convexity of the two sub-problems can be further improved by introducing the auxiliary variables h, g, e and l. The auxiliary variables also can quickly jump out from the local (a) Blurred image, (b) the magnified region of (a), (c) the results by EMLO [28], (d) the results by GR [31], (e) the results by MC [35], (f) the results by ELR [17], (g) the results by DCP [14], (h) the results by ECP [15], (i) the results by SGF [36], (j) the results by TSI [37], (k) the results by our algorithm without kernel optimization, (l) the results by our algorithm (REKO) Average Energies

FIGURE 18
The convergence property of our algorithm optimum and converge to the global optimal value during iterations. It can be concluded that our algorithm can obtain high accuracy of blur kernel estimation on the basis of convergence.

Computational complexity
Our algorithm is carried out within a multi-scale framework and iterates from coarse to fine. Intermediate latent images and blur kernels are updated alternatively at each scale. Table 4 shows that the run-time of the evaluated algorithms. The size of blur kernels is initialized as 25×25. All experiments are carried on the machine with Intel Core i7-7700HQ processor and 8GB RAM. The run-time of our algorithm is next only to GR (C++) and MC (MATLAB). However, our algorithm can achieve more satisfactory restoration results than other algorithms. FFT and energy alternating minimization are used to accelerate Equations (20) and (31), which makes our algorithm have super-linear convergence and reduces the total iterations. Overall, our algorithm yields better restoration quality than other algorithms on the basis of short run-time.

Parameter analysis
In the proposed algorithm, some parameters ( , , and ) should be given in advance. By varying one and fixing the others, we evaluate the effects of the parameters on the dataset [33]. The kernel similarity is used to measure the accuracy of estimated kernels. Figure 19 shows that the numerical variations of and have no effect on the kernel similarity. And the kernel similarity is basically constant when the values of and are larger than 0.003. This indicates that our algorithm is insensitive to the value settings of parameters.
The proposed restoration model also relates to the blur kernel size. It is important to explore the effect of the initial kernel size on image restoration. Quantitative metrics of restored images under different sizes of the blur kernel are shown in Figure 20. When the size of the blur kernel is smaller than 10×10, the proposed algorithm cannot converge to the global optimal and cannot obtain satisfactory restoration results. Based on the sparseness-inducing constraint of blur kernels, the discrepancies of the restored image are not conspicuous with the increasing the blur kernel size. Commensurately, the restoration results are qualitatively shown in Figure 21. It can be seen that our

Effect of pre-processing on image restoration
In order to explore the effect of pre-processing on image restoration, we smooth the blurred image by mean filtering or median filtering before the image restoration process. Figure 22 shows that restoration results after pre-processing are not as good as those without pre-processing. It can be concluded that pre-processing would have adverse effects on the restoration accuracy and our restoration algorithm can effectively suppress noises.

CONCLUSION
In this paper, we propose an image restoration algorithm for high energy flash X-ray radiography. The proposed algorithm uses the region extrema prior and blur kernel optimization to accurately restore the blur kernel and sharp image from the blurred image. According to the fact that region extrema distribution of high energy flash X-ray images deviates from opposite ends of image grey domain after being blurred, the sparsity of both region extrema and blur kernel is constrained by the sparseness-inducing norm. we construct the objective function under the MAP framework. The non-convex and non-linear objective function is decomposed into two sub-problems, which are handled by dually linear approximation and energy alternating minimization. In order to further improve the estimation accuracy of the blur kernel, a blur kernel optimization algorithm based on continuity is proposed. The main structure of blur kernels is extracted and the kernel continuity function is constructed in sliding cross windows to suppress discontinuous kernel elements.
Extensive experiments demonstrate that our algorithm can estimate the system blur kernel with a more symmetric trailing on high energy flash X-ray static images, and obtain the restoration results with sharper edges in comparison with stateof-the-art algorithms. The effectiveness of our algorithm is further verified through the analysis of quantitative metrics such as PSNR and SSIM. In addition, the proposed model is suitable for describing the blur characteristics of visible images, excellent restoration results can also be achieved on blurred visible images. Next, we will introduce the clustering algorithm to determine the continuous region where maximum and minimum do not exist in the blurred FTO image. According to the region boundary determined by this algorithm, we will concentrate on the weighted calculation algorithm of the equivalent linear operators in our restoration model, which could improve the discrimination of the interface between tungsten and copper.