Single-Image Super-Resolution Using Panchromatic Gradient Prior and Variational Model

. Single-image super-resolution (SISR) is a resolution enhancement technique and is known as an ill-posed problem. Motivated by the idea of pan-sharping, we propose a novel variational model for SISR. The structure tensor of the input low-resolution image is exploited to obtain the gradient of an imaginary panchromatic image. Then, by constraining the gradient consistency, the image edges and details can be better recovered during the procedure of restoration of high-resolution images. Besides, we resort to the nonlocal sparse and low-rank regularization of image patches to further improve the super-resolution performance. The proposed variational model is eﬃciently solved by ADMM-based algorithm. We do extensive experiments in natural images and remote sensing images with diﬀerent magnifying factors and compare our method with three classical super-resolution methods. The subjective visual impression and quantitative evaluation indexes both show that our method can obtain higher-quality results.


Introduction
Image super-resolution (SR) is one of the most fundamental problems in the field of image processing, which aims to reconstruct clear and accurate high-resolution images from degraded low-resolution images. In other words, SR technology can recover a high-resolution (HR) image from one or more low-resolution (LR) input images [1]. In recent years, image SR has attracted great attention from academia and industry. e quality of the reconstructed SR image also greatly affects the accuracy of other computer vision tasks (such as image classification, image segmentation, and target detection). At present, SR technology has made great progress; it has been widely used in high-definition digital television, remote sensing monitoring, medical image reconstruction, image restoration, military reconnaissance, and other fields [2]. However, as an ill-posed problem, the existing image SR technology cannot always obtain satisfactory reconstruction results. erefore, the study of SR is still a challenging but significant research topic.
According to the number of input LR images, the SR can be classified into multi-image super-resolution (MISR) [3][4][5] and single-image super-resolution (SISR) [6][7][8]. MISR uses the complementary information provided by multiple images of the same scene to enhance resolution. However, in a practical application, it is much difficult to obtain image sequences of the same scene and achieve precise subpixel registration. Meanwhile, SISR is suitable for more scenarios, which enhances the resolution of one input image based on some prior information [9,10]. To date, SISR methods can be broadly classified into interpolation-based methods, reconstruction-based methods, and learning-based methods [11].
Interpolation methods like nearest-neighbor interpolation, bilinear interpolation, and bicubic Interpolation [12,13] estimate the intensity at a point using the information of adjacent pixels. While these methods are easier calculated and very speedy, they tend to generate excessive smoothing and jagged artifacts. Reconstruction-based SISR methods [14][15][16] generally utilize image priors to restrict the possible solution space with an advantage of recovering sharp details. However, these methods are usually timeconsuming, and their performance suffers a rapid degradation when the amplification factor increases.
Learning-based SISR methods include example-based methods, neighbor-embedding-based methods, sparse-representation-based methods, and deep-learning-based methods [17,18]. e basic idea of these approaches is to obtain the correspondence between the LR image and the HR image by a priori knowledge training and then reconstruct an HR image from the input LR image using the correspondence that has been learned. Freeman et al. [19] presented an example-based method with a Markov network. Chang et al. [20] proposed a neighbor-embedding-based method that needs less training set. Yang et al. [21] proposed a sparse-representation-based algorithm to obtain more favorable results and further improve their work by jointly training coupled dictionaries for LR and HR image patch pairs [22]. Very recently, deep-learning-based SISR methods, using convolutional network [23,24], residual network [25,26], generative adversarial network [27], attention network [28,29], and so on, have become popular and demonstrated good performance with the rapid development of deep learning technology. Pan-sharpening [30] is a multispectral image fusion and super-resolution technique; the panchromatic (PAN) image has higher spatial resolution than the corresponding multispectral image but only has a single band. Pansharpening fuses PAN image and multispectral image together to enhance both spatial and spectral resolutions of the data. Motivated by the idea of pan-sharpening, we present a novel variational approach for single-image super-resolution. In this paper, we assume that there is a PAN image of the input LR image. By utilizing the structure tensor [31] of input image, we can construct the gradient of the PAN image; then we build a variational model combined with low-rank and sparse representation of similar patches to effectively fuse the constructed PAN information with LR image and obtain HR image. e proposed approach addresses the super-resolution problem of both single nature image and the multispectral data without panchromatic image. e remainder of this paper is organized as follows. In Section 2, we present the proposed SISR model in detail. e numerical algorithm for solving our model is given in Section 3. e experimental results and analysis are discussed in Section 4. Finally, conclusion is presented in Section 5. In addition, the descriptions of the acronyms used in this paper are listed in Table 1.

Variational Model
Drawing on the idea of pan-sharpening, we construct a variational fusion model to realize single-image super-resolution. As panchromatic image contains the edges and details needed for resolution enhancement, we first construct the information of PAN image which does not actually exist. According to [31], the key information of the assumed PAN image can be exploited from the structure tensor of the input image. en the gradient of PAN image can be fused into the LR image to enhance the HR details. e similar image patch pairs extracted from LR and HR images are constrained by low-rank and sparse regularization to improve the SR performance.

Image Degradation Model.
We consider that the observed low-resolution image H is a downsampled version of the high-resolution image F, and the degradation model [32] can be written as where D represents a downsampling operator and ϵ is the additive Gaussian white noise.

Construct the Gradient of PAN Image.
We construct the gradient of the PAN image corresponding to input LR image using the knowledge of structure tensor [31]. Firstly, we apply quadratic interpolation to the given LR image H to obtain the enlarged image H * in desired scale. en, matrix G known as the structure tensor [31] can be expressed as where N denotes the number of spectral bands; s i denotes the weight and it can be determined by . As a symmetric matrix, G can be decomposed as where Q is an orthogonal matrix and the columns of which are the eigenvectors of G; Λ is a diagonal matrix and the diagonal elements of which are the eigenvalues of G. Equation (3) can also be written as where λ 1 is the maximum eigenvalue and it gives the maximum rate of change of H * ; λ 2 is the minimum Iterative-back-projection-based SC method eigenvalue and it gives the minimum rate of change. e corresponding eigenvectors θ 1 and θ 2 give the directions of change. e PAN image P can theoretically capture the basic geometry of the image H * . Consider image P with λ max � |∇P| 2 and λ min � 0, whose structure tensor ∇P(∇P) T should approximately be equal to G.
us, we have the following equation: In order to solve ∇P, equation (5) can be rewritten as follows: To specify the sign of the eigenvectors, the target gradient of the constructed PAN image P is thus obtained as

Variational Model Based on Gradient Consistency.
In this subsection, we build a variational model to fuse the lowresolution image H and the gradient of its corresponding PAN image P; then the model is optimized together with the nonlocal sparse and low-rank regularization to obtain the high-resolution image F.

Consistency of Gradient.
Bands of an ideal HR image F often have similar structural information with the PAN image P; namely, bands of HR image closely approximate the PAN image in gradient information. erefore, we can use constraint about the consistency of gradient to establish a relationship model. We use F i to denote the i-th band of image F and then set up the gradient consistency regularization term as follows: where α i denotes the weight of the i-th band in F, and we generally set α i � 1/N in the absence of special requirements.

Sparse and Low-Rank Decomposition.
Motivated by the effective work of optical flow estimation with nonlocal sparse and low-rank regularization in [33], we adopt a joint sparse and low-rank decomposition term formulated as follows: where R k denotes the operator that extracts the k-th patch of image and groups its nonlocal similar patches together. en the grouped similar patches are decomposed to the low-rank component L k and the sparse component S k , and ε denotes the Gaussian white noise.

Total Energy.
Combining the equality constraint in (1), (8), and (9), the proposed variational SISR model can be formulated as the following optimization objective: where M is the number of image patches that are extracted from the i-th band of image F, N is the number of spectral bands, μ 1 , μ 2 , μ 3 , μ 4 are parameters, rank(L k ) denotes the rank of low-rank component, and ℓ 0 norm denotes the sparse regularization. It is well known that the rank-solution and ℓ 0 norm are discrete combinatorial optimization problems; they are both NP-hard. erefore, in order to make the minimization tractable, we resort to their convex relaxed modification, and model (10) can be reformulated as the following convex optimization problem: where ‖·‖ * denotes the nuclear norm, which is the convex envelope (tightest convex surrogate) of rank, and ℓ 1 norm is the nearest convex norm to ℓ 0 norm.

Numerical Algorithm
In this section, the numerical procedure of the proposed model will be implemented. We adopt the alternating direction method of multipliers (ADMM) [34] to solve the minimization problem (11). e flexibility of the ADMM lies in the fact that it splits the initial optimization problem into several subproblems.

Solution of the Sparse Matrix S k .
For each band, we fix F and the low-rank matrix L k ; the optimization problem of sparse matrix S k is as follows: We can obtain the solution of equation (12) by soft threshold algorithm as Mathematical Problems in Engineering 3 where soft(a, τ) � max ‖a‖ 2 2 − τ, 0 sign(a).

Solution of the Low-Rank Matrix L k .
For each band, we fix F and the sparse matrix S k ; the optimization problem of low-rank matrix L k is as follows: We can obtain the solution of equation (14) by singular value threshold algorithm as where (σ i − τ) + � max σ i − τ, 0 , σ i are singular values of matrix R k F − S k , U and V are left and right singular vectors of σ i , and τ � μ 3 /μ 2 .

Solution of the SR Image F.
In order to separately solve the fidelity term and two regularization terms in (11), we introduce two intermediate variables U and V which are close to F. us, we need to iteratively solve the three following subproblems to obtain the solution of SR image.
(1) We use the intermediate variable U to replace F in (1/2)‖DF − H‖ 2 2 . e optimization problem for U is as follows: where 〈·, ·〉 denotes the inner product, W 1 denotes the Lagrange multiplier matrix, and ρ is a positive penalty parameter. We set the derivative of (16) with respect to U to 0: which yields e Lagrange multiplier W 1 is updated as follows: (2) We introduce an intermediate variable V i to replace F i in the gradient consistency term. e optimization problem for V i is as follows: where W 2 denotes the Lagrange multiplier matrix. We set the derivative of (20) with respect to V i to 0: Using the fast Fourier transform (FFT), we can obtain the closed-form solution of V i from (21) directly: where F and F − 1 denote the FFT and inverse FFT, respectively. e Lagrange multiplier W 2 is updated as follows: (3) After U and V are substituted into equation (11), the optimization subproblem of F becomes We set the derivative of (24) with respect to F i to 0, which yields Overall, taking all above analyses into account, we can summarize the complete numerical procedure for the proposed method. e detailed descriptions are shown in Algorithm 1.

Experimental Results and Analysis
In this section, we first evaluate the sensitivity of the proposed algorithm with respect to the main parameters, that is, μ 1 , μ 2 , μ 3 , and μ 4 . en, to demonstrate the SR effectiveness on both natural images and remote sensing images with different magnifying factors (2 and 4), we compare our method with the Bicubic interpolation method, sparsecoding-based method (SC) [22], and its back-projection enhanced version (SCBP) [35]. We evaluate the outcome of various methods by using quantitative indexes: peak signalto-noise ratio (PSNR), root mean square error (RMSE), and structure similarity index (SSIM) [36]; they are calculated on RGB channels for nature image and on each band for multispectral image. e higher PSNR and SSIM values and lower RMSE value represent the better SR performance.

Adjustment of Parameters.
e parameter μ 1 is the weight to the consistency of gradient, μ 2 controls the lowrank and sparse regularization of the similar patches, and μ 3 and μ 4 control the contribution of low-rank and sparse regularization terms. We adjust the parameter values in proper ranges and demonstrate the average quantitative metrics as functions of the parameters on the experimental dataset; then the optimal parameter combination is found with extensive numerical tests. Figure 1 shows the average PSNR, RMSE, and SSIM results when the parameter μ 1 varies from 0 to 0.5 with fixed μ 2 � 1.6, μ 3 � 0.2, and μ 4 � 0.2. From Figure 1, we can see that when μ 1 � 0, all quantitative metrics get the worst values, which illustrate the effectiveness of the gradient consistency constraint. Meanwhile, Figure 1 shows that our method achieves the best performance when μ 1 � 0.4. Figure 2 shows the average PSNR, RMSE, and SSIM results under various combinations of μ 2 and μ 3 with fixed μ 1 � 0.4 and μ 4 � 0.2. e analysis of μ 2 and μ 3 shows that the trend of the overall performance gets better and then worse. When the maximal PSNR is selected, we have μ 2 � 1.8 and μ 3 � 0.3; when the maximal SSIM is selected, we have μ 2 � 1.4 and μ 3 � 0.1. In order to synthesize the performance of our algorithm in these two indexes, we normalize PSNR and SSIM by introducing E � 0.8 * PSNR + 0.2 * SSIM, and E gets its maximum when μ 2 � 1.6 and μ 3 � 0.2. Besides, we find that the performance of the proposed method is insensitive to the parameter μ 4 , and we simply set μ 4 � 0.2.

Natural Images Super-Resolution.
To evaluate the performance of the proposed super-resolution algorithm, we use bicubic interpolation and SC and SCBP algorithms for comparison. Experimental results of the images "Flower," "Leaf," and "Lena" with magnification factor of 2 are shown in Table 2 and Figure 3, while the results with magnification factor of 4 are shown in Table 3 and Figure 4.
As we can see from the visual comparison, our method could obtain desired results; it effectively reduces zigzag and distortion of the recovered SR images so that it better preserves the image details compared to other competitors. From the close-ups of the regions in red boxes, we can see clearly that the SR images obtained by bicubic interpolation and SC methods are more blurred than ours, and SCBP method could restore clear edges but generate some obvious artifacts in SR images. In summary, our algorithm is superior to the contrast methods in terms of reconstructing clearer SR images with sharper edges and without artifacts. e above observations can be quantitatively confirmed by the Tables 2  and 3 which record  e SR image F * .

Initialization:
Compute an initial SR image F 0 via bicubic interpolation; Calculate the target gradient ∇P via equation (7).

Multispectral Images Super-Resolution
. e adopted multispectral images contain four bands (red, blue, green, and infrared bands), and our algorithm is executed at each band. Same as the nature image experiments, we compared the performance of our method with bicubic interpolation, SC, and SCBP. Experimental results for the multispectral datasets "Field," "Tree," and "Beach" with magnification factor of 2 are shown in Figure 5 and Table 4, while the results for magnification factor of 4 are shown in Table 5. e observation of visual comparison is close to that of natural image experiment. Bicubic interpolation and SC algorithms tend to reconstruct overly smooth images that lack detailed information and have fuzzy image quality, while the SR images obtained by SCBP method are overenhanced and have obvious artifacts. Comparatively speaking, our algorithm can better restore the detailed information of images and make a better texture effect. By comparing the numerical results in Tables 4 and 5, our algorithm also gets the best results for all indicators. It is worth noting that, because of the lack of training set, the results of SCBP algorithm are worse than bicubic    interpolation. In conclusion, our proposed algorithm is superior to other algorithms compared in the subjective visual impression and objective indexes.

Conclusion
is paper presented an efficient variational method for single-image super-resolution. Motivated by image fusion and pan-sharpening, we exploit the gradient of imaginary PAN image from the input image via structure tensor and then construct a gradient consistency constraint which is supposed to fuse HR edges and details into the optimization target image. Meanwhile, we adopt the regularization of sparse and low-rank decomposition for similar image patches to further improve the SR accuracy. e variables of our variational model are iteratively optimized by ADMMbased algorithm. Extensive experiment results demonstrate the edge recovery and artifacts suppression capabilities of our method.

Data Availability
e datasets used to support the findings of this study are public available and included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.