Single image super-resolution based on approximated Heaviside functions and iterative refinement

One method of solving the single-image super-resolution problem is to use Heaviside functions. This has been done previously by making a binary classification of image components as “smooth” and “non-smooth”, describing these with approximated Heaviside functions (AHFs), and iteration including l1 regularization. We now introduce a new method in which the binary classification of image components is extended to different degrees of smoothness and non-smoothness, these components being represented by various classes of AHFs. Taking into account the sparsity of the non-smooth components, their coefficients are l1 regularized. In addition, to pick up more image details, the new method uses an iterative refinement for the residuals between the original low-resolution input and the downsampled resulting image. Experimental results showed that the new method is superior to the original AHF method and to four other published methods.


Introduction
Image super-resolution (SR) is to generate or recover high-resolution (HR) images from one or multiple low-resolution (LR) images. If we generate/recover the HR image from only one LR image, we call it single-frame SR. Otherwise, we call it multiple-frame SR. The multipleframe SR methods are available in that multiple LR images are of the same scene with different sub-pixel shifts taken as input. It has direct applications in video SR problems (see [1]). Singleframe SR methods are quite popular and challenging when only one LR image is avaliable. In particular, we focus on single-frame SR problems in this paper.
To produce a HR image, the simplest and effective way is to interpolate, e.g., bicubic and nearest interpolations. Recently, more interpolation-based methods have been proposed(see [2][3][4][5][6][7][8]). For instance, Youngjoon et al. [8] utilize a generalized curvature source term estimated from the LR image to construct a HR image. In particular, the resulting HR image has a reliable curvature profile which minimizes ringing artifacts. In addition, they propose an iterative application of the curvature interpolation method [9]. The method utilizes the gradientweighted curvature measured from the LR image, being an interpolator to suppress texture oversmoothing. In [10], Wang et al. present a fast image upsampling method to preserve the represented by multiple classes of AHFs with sharp edges. In particular, due to the sparsity of non-smooth components, we give l 1 regularization model and solve the model via block-wise alternating direction method of multipliers (ADMM) [26]. Furthermore, we design a novel iterative refinement algorithm to pick up more image details. Finally, the proposed method has been numerically proved competitively to some state-of-the-art methods.
The paper is organized as follows. A brief review of the AHFM algorithm [7] is introduced in Section 2. In Section 3, we give the proposed model and its corresponding algorithms. In Section 4, we present the numerical results of different methods. It demonstrates that the proposed algorithms are more efficient. Finally, we conclude the paper in Section 5.

Some preliminaries
This section gives general remarks on approximated Heaviside functions. In addition, a brief review of the method based on approximated Heaviside functions can be found from [7].

General remarks on Heaviside function
Heaviside function or Heaviside step function, is a discontinuous function whose value is zero for negative argument and one for positive argument. Heaviside function could be defined as following alternative form of the unit step, as a function of a discrete variable x (see Fig 1(a)), The definition of ϕ(0) = 0 is significant. In the practical applications, some logistic functions to the Heaviside functions are often used for smooth approximations, called approximated Heaviside functions (AHFs), such as As illustrated in Fig 1(b), a smaller ξ corresponds to a sharper transition at x = 0. In our work, we employ Eq (3) where ψ is an approximated functions. H d is the set of characteristic functions of closed halfspace of R d . Theorem 1 [27] For any positive integer d, Theorem 2 [27] For any positive integer m,d and every p 2 [1, 1), span m H d is approximately a compact sub-set of (L p ([0, 1] d ), k Á k p ).
Consequently, we can use span m H d for a finite m in practical computing.

Single image super-resolution via iterative AHF method (AHFM)
The single image super-resolution via iterative AHF method (AHFM) proposed in [7] for image super-resolution gives a selection of sharp-related terms, which are measured from the LR image input and apply them to fine grids to generate the HR image. They assume the underlying image intensity function f is defined on [0, 1] 2 , then f 2 L p ([0, 1] 2 ) with p 2 [1, 1).
According to the theorems stated in section 2.1, f can be approximated by the following equation: where o j 2 R; v j 2 R 2 , v = {(cosθ t , sinθ t ) 0 , t = 1, 2, . . ., p} denote p different directions, and c j ¼ f 1 q ; 2 q ; 3 q . . . ; 1g is to denote discrete positions, m = pq, z = (x, y) 0 , where q is the total number of pixels of the input image. Consequently, the function fcðv j Á z þ c j Þg m j¼1 is called a class of AHF with a specific ξ. For an image L 2 R n 1 n 2 , we assume it is a discretization of intensity function f on [0, 1] 2 , i.e., Therefore, Eq (4) could be rewritten as matrix-vector form, L % Co; f 2 R n ; o 2 R m , with n = n 1 n 2 , m = pq. We compute coefficient ω and get the high resolution image with equatioñ Co, where s is an upscaling factor andC 2 R NÂm with size N = s 2 n 1 n 2 . By this strategy, based on the observation that an image consists of smooth components and non-smooth components. We use two classes of AHFs to depict an image. Different components of an image may be described by different orientations θ t at the locations c j with two ξ 1 , ξ 2 . One is big parameter ξ 1 to represent smooth components (forming C 1 ), another one is the smaller ξ 2 to represent non-smooth components (forming C 2 ). Thus, the vector-form image L can be approximated by the following discrete formula: Since non-smooth components are sparse, l 1 regularization could be given on the coefficient β 2 . The optimal model is expressed as: The Eq (6) is solved via ADMM [28][29][30][31] in article [7]. The single image super-resolution via iterative AHF method (AHFM) proposed in [7] is outlined as shown in Algorithm 1.
We close this section with the following remarks.
• In order to pick up more details, such as edges, they design an iterative strategy to conduct an iterative refinement. They consider the difference (L − DH) as a new low-resolution input of Eq (6) to recompute a residual high-resolution image.
• The AHFM algorithm performs well for natural images. However, the images with smooth backgrounds usually appear ring artifacts along the large scale edges, which mainly come from the added non-smooth components (see step (b) in Algorithm 1). Aiming to discard the ring artifacts of non-smooth components E, they use bicubic interpolation as the intermediate method to make a mask. Actually, in the step (b), only the E needs to update by E new , which can be obtained from by the Eq (7), and the ring artifacts could be reduced significantly.
where G i, j is a vector-form of gradient at location (i, j) of image B. The image B is generated via the bicubic interpolation. Notation . Ã stands for dot product. t is a threshold value and t = 0.05 is in the experiments. Algorithm 1 (Single image super-resolution via iterative AHF method (AHFM) [7]) Input: one vector-form low-resolution image: L 2 R nÂ1 ; l 1 > 0; l 2 > 0, s: upscaling factor. τ: maximum number of iteration. Output: high-resolution imageĤ 2 R NÂ1 1. According to Eq (4), construct matrices C 1 ; C 2 2 R nÂm on coarse grids, on fine grids constructC 1 ;C 2 2 R NÂm , where N = s 2 n.
2. Initialization: L (1) = L. for k = 1:τ 1. Compute the coefficients: where conv represents a convolution operator, and p is a Gaussian kernel with a small size.

Modified AHFM and iterative refinement
This section presents a new image super-resolution algorithm, which is an extension of the AHFM algorithm. Note that only two classes of AHFs are not enough for the whole image components. Hence, taking into account the varying sharpness of the whole image, we will consider a modified AHFM with the different sharpness components and further present its algorithm for implementation in next section. Besides, the modified AHFM algorithm contains a new iterative refinement strategy, aiming to pick up more information about nonsmooth components.

Modified AHFM with different sharpness components
The AHFM has its respective advantages, which is completely a single image super-resolution method without extra training data. The AHFM algorithm has only used two classes of AHFs to represent the whole image components. The sharper components or the smoother components are pivotal but are not represented well. Motivated by the point and the works proposed in [7], we propose a modified AHFM algorithm to comprise the whole information of one image as well as possible.
First, we assume that a low-resolution image L is consisted of smooth components and non-smooth components. We exploit different ξ to form different AHFs to describe smooth components and non-smooth components. Hence, the low-resolution image L could be approximated by the following discrete formula: where l, k represent numbers of smooth components and non-smooth components, respectively. We could obtain kÞ are the corresponding representation coefficients with C i (i = 1, 2, Á Á Á, l) and F j (j = 1, 2, Á Á Á, k). After computing the representation coefficients, we apply them into following the equation to obtain the high-resolution image: whereC i ;F j 2 R NÂm ði ¼ 1; 2; Á Á Á ; l; j ¼ 1; 2; Á Á Á ; k; N ¼ s 2 nÞ are obtained on fine grids. s is the upscaling factor. Since the non-smooth components, such as edges and corner details, are sparse in generic images. Hence, we apply l 1 regularization on β j (j = 1, 2, Á Á Á, k) to characterize this feature.

Modified AHFM with iterative refinement
The Eq (9) takes different smooth components and non-smooth components into consideration. A natural remedy for extracting more details is to find the structure of the difference between LR image and last updated resulted image. We find that residual R ¼ L À DĤ , wherê H is the last updated HR image and D is the downsampled operator. Fortunately, we find the residual is mostly non-smooth components. To pick up more edge details to make image less blurry, we design a new iterative refinement model based on the special structure. In the model, we use two smaller to depict the residual image. Since the nonsmooth components are also sparse in the residual. We apply l 1 regularization to the corresponding coefficients. The iterative refinement model could be written as following: are the coefficients of the non-smooth components. μ 1 , μ 2 are regularization parameters. Since the Eq (19) is the form of l 1 norm, we solve it by ADMM scheme. We make two variable substitutions for b Ã 3 ; b Ã 4 , and rewrite Eq (19) as: We rewrite Eq (20) as: , 0), and B = (0, I). And 0 is zero matrix. I is identity matrix. Since the optimization Eq (21) is separable, w.r.t (β Ã , u 1 , u 2 ). The augmented Lagrangian of Eq (21) is are the proper size Lagrangian multipliers. By ADMM, the minimizing of Eq (22) is solved by following three subproblems: We could use least squares method to solve β Ã -subproblem: The subproblem u 1 and u 2 have a closed form solution as mentioned in Algorithm 2. Thus, the solutions to the subproblems u 1 and u 2 are shown as following: The main algorithm for iterative refinement is shown in Algorithm 3. We will illustrate the necessary of the Algorithm 3 in section 4.

Single image super-resolution based on approximated Heaviside functions and iterative refinement
Take different behaviours of Eqs (9) and (19)   Single image SR based on approximated Heaviside functions and iterative refinement
We use a Gaussian kernel p with a small size to make a convolution to avoid the oversharp information on the non-edge parts. And D is the bicubic downsampling operator.
However, the computation of the Algorithm 4 is very expensive, due to the large scale and non-sparse matrix C in Eq (12). For example, if a LR image is 512 × 512 and we choose 10 different directions, the size of matrix C will be equal to (512 2 × 10 × 512 2 ). It is obviously large. Thus, it is observed that if we increase the number of smooth components and non-smooth components. The quality of an image usually gets improved at the cost of increased computation time and memory requirement. However, for large number of smooth components and non-smooth components, it is difficult to do SR on a limited hardware. We have experimentally seen that the use of more than l = 1, k = 2 in Algorithm 4 does not improve quality of the images significantly. In addition, l = 1, k = 2 work well on a regular desktop. Hence, we choose l = 1, k = 2 in our work. The optimization model can be simplified as: To simplify Eq (29), we setC ¼ ðC 1 ; F 1 ; F 2 Þ,ũ ¼ ða 1 ; b 1 ; b 2 Þ, M = (I, 0, 0), B 0 1 ¼ ð0; I; 0Þ, B 0 2 ¼ ð0; 0; IÞ, where 0 is zero matrix, I is identity matrix. The Eq (29) can be rewritten as following.
Since l 1 term is not differentiable, we make two variable substitutions p 1 , p 2 for B 0 1ũ ; B 0 2ũ to rewrite the Eq (30) as: The Eq (31) is separable, w.r.t ðũ; p 1 ; p 2 Þ. Here, we solve this problem by block-wise ADMM. The augmented Lagrangian of Eq (31) is where b 1 , b 2 are Lagrangian multipliers with proper size. The problem of minimizing Lðũ; p 1 ; p 2 ; b 1 ; b 2 Þ could be solved by iteratively and alternatively by following subproblems: Theũ-subproblem is a smooth quadratic problem. We solve it by least squares method as following:ũ where The p j -subproblem(j = 1, 2) has a closed form solution Eqs (39) and (40) for each (p j ) i (j = 1, 2): where shrink(a, b) = sign(a) max(jaj − b, 0). We summarize l = 1, k = 2 in Algorithm 4 as modified single image based on approximated Heaviside functions and iterative refinement, which we call this algorithm as MAHFM, shown in Algorithm 5.
Algorithm 5 (MAHFM algorithm) end 5. Get non-smooth components: where conv represents a convolution operator, p is a Gaussian kernel with a small size.

Numerical results
We have implemented our algorithm to some benchmark images whose high-resolution versions are available. For a quantitative comparison, we downscale actual HR images to their LR versions via bicubic interpolation and then generate HR images by MAHFM algorithm and other methods. For gray images, we perform the proposed algorithm to them directly. While working for color images, the input image is first converted from RGB to YCbCr. We only perform the MAHFM algorithm on the illuminance channel because the human are more sensitive to the brightness information. The other two channels contain chromaticity information and they could be upsampled by bicubic interpolation. Finally, we convert three channels back to RGB to get the estimated color HR image. We make the quantitative comparisons between the recovered HR images and the actual HR images. And we use root mean square error (RMSE) and peak signal to noise ratio (PSNR) on the illuminance channel to evaluate numerical performance. ð41Þ where h,ĥ are the vector-form of ground-truth image and the resulted high-resolution image, respectively. N represents testing times for one image. IJ is the original image consisting of (I × J) pixels. n is the number of bits per sample. The smaller RMSE is, the better performances of the method usually are. But PSNR has the opposite property. As mentioned in [7], if we apply MAHFM algorithm to one image, it would involve expensive computation. It follows that we could utilize image patches to reduce computation and storage significantly. In our work, we set patch size to be 6 × 6 and overlap to be 3. First, we compare the MAHFM algorithm with AHFM algorithm [7], and the numerical results are shown in Table 1. The performance of MAHFM algorithm is compared with that of bicubic interpolation, a kernel regression method (denoted as "07'TIP" [34]), a fast upsampling method (denoted as "08'TOG" [22]), one state-of-the art learning-based method (denoted as "10'TIP" [12]), a fast image upsampling method via the displacement field (denoted as "14'TIP" [10]) and AHFM [7]. The numerical results by these methods are displayed in Table 2. Furthermore, we have also carried out our experiments for the upscaling factor s = 3. The quantitative measurements for them are shown in Table 3. As can be seen, these results reveal that the MAHFM algorithm is effective.

Visual comparisons
The test images (LR images) in Figs 3 and 4 are shown in Figs 3(a) and 4(a) of every figures, respectively. The image "baboon" and "peppers" are downloaded from http://decsai.ugr.es/ cvg/dbimagenes/c512.php. The bottom of the web page has shown that these images are Copyright free. The version of the images in our paper are other special formats which are converted by Photoshop from the sources above. We have shown the results on these images with different upscaling factors. For example, we have set the upscaling factor s = 3 in the Fig  3. The proposed algorithm is compared with bicubic interpolation, a learning-based method ("10'TIP" [12]), one state-of-the-art interpolation-based methods ("11'IPOL" [2]), one deep learning method("14'TIP" [10]) and the AHFM [7]. As can be seen, blur effect is generated by bicubic interpolation method. The learning-based method has comparable vision while we have to need extra training data to learn the relation between the test images and the training images. In addition, the interpolate-based methods ("11'IPOL") also show impressive performance. In particular, "14'TIP" shows superior vision on image edges and runs quite fast. But as shown in figures, not only does it lose small-scale texture but also introduces minimal jagged artificial and smooths non-edge regions. By contrast, the recovered HR images by the MAHFM algorithm are preserved sharpness on the edges, such as the texture of the mustache without ringing and blurring artifacts. Compared with the ground truth image, the MAHFM alorithm generates superior effect than others not only visually but also quantitatively. Single image SR based on approximated Heaviside functions and iterative refinement

Comparisons between the AHFM and the proposed algorithm (MAHFM)
In this section we will illustrate the comparisons between the AHFM algorithm [7] and MAHFM algorithm. First, it is necessary to illustrate the difference between Algorithm 2 (let l = 1, k = 2) and AHFM algorithm. Although the only difference between the two methods is that the MAHFM algorithm is consisted of another term, which characters more sharp components. These could be seen much more clearly in Fig 5, which shows that the resulted HR image of Algorithm 2 (let l = 1, k = 2) is more sharp than that in the AHFM algorithm, especially at edges in Fig 5(c). Alternatively, we also compare the performance of the iterative refinement generated with that of the iterative AHFM algorithm [7]. To measure the effect of iterative refinement, we employ two classes of relative error defined as following: and where L is a low resolution image. D is the downsampling operator. H is the last updated resulted image. H true is the true high resolution image, and k Á k F is the Frobenius norm. After computing the relative errors of all the test images, we take the average of all the relative errors as the final results. The final results compared with the iterative AHFM algorithm in [7] are plotted in Fig 6, illustrating the Algorithm 3 could produce more details and less error than the original iterative refinement method. Besides, the precision of the proposed method is higher than that in original method. From the Fig 6, we can also get the iterative steps of the iterative refinement method. Only two iterations could be achieved the goal of refinement, which would also reduce the operation time.

Parameters
As mentioned in section 3, there are many parameters in the MAHFM algorithm, e.g. ξ 1 , ξ 2 , ξ 3 , λ 1 , λ 2 , λ 3 and others. Although the parameters are so many, they are easy to select. In our work, we select these parameters mainly according to the experience. In particular, the parameters ξ 1 , ξ 2 , ξ 3 are especially important, because the sharpness for different components of an image are decided by them. First, we test ξ 1 by fixing ξ 2 , ξ 3 and then tune ξ 1 . In this way, we obtain a rough estimation of the three parameters (see Fig 7(a)). Note that there are some special points, which we can roughly estimate ξ 1 2 [0.8, 1], ξ 2 2 [10 −2 , 10 −1 ] and ξ 3 2 [10 −4 , 10 −2 ] .  Fig 7(a) as well as reveals ξ 1 and ξ 3 are not sensitive to the selection of parameters. Thus we can select ξ 1 = 0.8 and ξ 3 = 10 −4 at first, and then tune ξ 2 from 10 −2 to 10 −1 by 10 −2 time change, which is plotted in Fig 7(b). We could find that ξ 2 = 6 × 10 −2 is the best fit. We conduct the same method to select other parameters. The final parameters of MAHFM algorithm are shown in Table 4.

Computation time
From Figs 8 and 9, we can see the proposed method is a little slower in operation time. The main reason is that our algorithm is mainly to compute the coefficients u, however, due to many components in the proposed algorithm have to be character thereafter there are a lot of loop programs used Matlab. Hereto, to speed up the computation time, there is a lot of room. We believe that the limitations could be improved by using cmex in Matlab involving a lot of loops to speed up the code. From Fig 8, we can also find that Algorithm 3 runs more time than Algorithm 2, since it involves more matrix-vector representation. In addition, in our work, we also use the theory of patch to speed up. Based on this, we can also consider utilizing the parallel computing.

The differences between Heaviside function and wavelets
In the image processing field, there exist some other basis functions such as wavelets [35]. Therefore, we discuss the difference between the Heaviside function and wavelets. To test the difference, we write a program about the representation of super resolution based on the haar wavelet. The quantitative comparisons are shown in Table 5. Compared Fig 10 with Fig 4, we can see that the algorithm based on Heaviside function has got more details about the nonsmooth components than the wavelets do. The algorithm based on wavelets mainly splits one image into different sub bands. The first layer scale coefficients and wavelet coefficients of the image are extracted from the structure of wavelet decomposition. The different sub bands focus on the different ways, such as the information about horizontal, vertical and other   Computation time of MAHFM vs. the size of low-resolution image "butterfly", the size of low-resolution image is increased from 30 × 30 to 110 × 110 and the upscaling factor is always set to be 4. Note that, at the common point with Fig 8 (i.e., upscaling factor 4 and 60 × 60 LR size), it has slightly different computation due to the instability of Matlab, thus we present the computation time by the average of 10 runs to reduce the gap.
https://doi.org/10.1371/journal.pone.0182240.g009 directions. Every sub band persists the information only about the corresponding direction. If some operators, such as bicubic interpolation, are applied, error would be generated. And then the error would make pixels overlap, causing Gibbs effect. Differently, the basis is generated by the corresponding Heaviside functions and the representation coefficients are extracted from the low resolution image. The MAHFM algorithm could represent different sharp components via different classes of Heaviside functions as well as possible, which makes the information minimize the loss.

Conclusion
In this paper, we have proposed a general framework of approximated Heaviside functions model that can describe different smooth components and non-smooth components in an image. The different components of an image could be approximately represented by different classes of approximated Heaviside functions (AHFs). With only one LR image input, we could compute the representation coefficients of AHFs, and then utilized these representation coefficients to obtain the high-resolution images. In addition, to pick up more image details, we employed an iterative refinement algorithm according to the residuals between the LR input and the downsampled LR image. To reduce computation and storage size, we applied the MAHFM algorithm to image patches.