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 multiple-frame 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]). Single-frame 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-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 gradient-weighted 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 sharpness via two-scale sharpness preserving operations. On the one hand, the low-frequency of image is recovered based on a well-constructed displacement field. On the other hand, the local high-frequency structures are reconstructed via a sharpness preserving reconstruction algorithm. However, due to images zooming with the interpolation-based methods to be solved are frequently estimated by information of unknown locations without other priors. The interpolation-based methods usually introduce jagged artifacts or blur effect.

Learning-based methods (see [11-14]) have attracted attention in image processing. Recently, they are also performed impressively in image SR problems. A key issue for the learning-based methods is to learn high frequency correspondences from a database generated by LR and HR image patches pairs. And then we could apply the correspondences to the LR image input to obtain its HR output. Purkait et al. [11] develop fuzzy rules to find different possible HR patches and combine them according to different rule strength to obtain the estimated HR patches. The rule parameters are learned from LR-HR patch pairs and then they use the Takagi-Sugeno (TS) model [15] with the rule parameters expressed as a linear combination of the different input possible HR patches. In addition, Yang et al. [12] apply the theory of sparse coding to SR problems effectively. This method jointly trains two dictionaries for LR and HR image patches, and then they could use the LR dictionary to generate sparse representations of the LR input to obtain the corresponding HR output. And Dong and Loy [13] first learn a mapping between low-resolution and high-resolution images, which is represented by a deep convolutional neural network (CNN). And then they take the LR image as input via the CNN to generate the HR image output. However, these methods typically reply on the similarity between test images and the database. Consequencely, they also involve expensive computation. In recent years, there has been tremendous interest in developing statistics-based methods [14, 16, 17], such as the popular tool Maximum a Posteriori (MAP) and Maximum Likelihood estimator (MLE). As proposed in [14], Peleg and Elad assume that prediction of high resolution patches can be obtained by MMSE estimation and the resulting scheme has the useful interpretation of a feedforward neural network.

Apart from the methods mentioned above, a hybrid method [18], reconstruction methods [19-21], and other methods [22-25] have been used. And these methods are not completely independent with each other. For instance, Peleg and Elad propose a SR method via combining a statistics method and a learning-based method.

In this paper, we study an effective single-frame SR approach, which is an improvement of the so called approximated Heaviside functions method (AHFM) proposed in [7]. It shows that the underlying image, can be viewed as a instensity function, which can be approximately represented by two classes of AHFs. Deng et al. cast the image super-resolution problem as an intensity function estimation problem. Defined on a continuous domain, the underlying intensity image, which belongs to a space with redundant basis, could be approximately represented via two classes of AHFs. Using only one LR image input, we could compute the representation coefficients by the proposed iterative AHF method. Then the high-resolution image is generated by combining the representation coefficients with two classes of AHFs. Here, the two classes of AHFs are corresponding to the high-resolution image.

However, there are only two classes of AHFs, which may not describe/represent the whole information of one image. It may generate some oversharp information for the final HR image. Forced by that, we tend to extend the AHFM algorithm to general form to suppress the resulted image texture oversharp and improve super-resolution images qualities. We consider that an image is normally consisted of different smooth components and non-smooth components. In particular, there are some details, such as edges and corners, which have different sharpness. Based on this, we propose that the smooth components should be represented by multiple classes of AHFs with smooth edges, and the non-smooth components are also represented by multiple classes of AHFs with sharp edges. In particular, due to the sparsity of non-smooth components, we give l1 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].

2.1 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 or As illustrated in Fig 1(b), a smaller ξ corresponds to a sharper transition at x = 0. In our work, we employ Eq (3) to approximate Heaviside functions.

Fig 1

(a) Heaviside function. (b) two approximated Heaviside functions with ξ = 1 (blue solid line) and ξ = 0.01 (red dash line).

From [27], Kainen et al. propose that any functions in L([0, 1]), p ∈ [1, ∞) could be well approximated by linear combinations of m characteristic functions of half-space, and m is any positive integer. Let H be a set of functions on [0, 1] defined as: where ψ is an approximated functions. H is the set of characteristic functions of closed half-space of .

Theorem 1 [27] For any positive integer d, define spanH as , where and , and , then it is known that UspanH is dense in (L([0, 1]), ∥ ⋅ ∥), p ∈ [1, ∞).

Theorem 2 [27] For any positive integer m,d and every p ∈ [1, ∞), spanH is approximately a compact sub-set of (L([0, 1]), ∥ ⋅ ∥).

Consequently, we can use span H for a finite m in practical computing.

2.2 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 ∈ L([0, 1]2) with p ∈ [1, ∞). According to the theorems stated in section 2.1, f can be approximated by the following equation: where , v = {(cosθ, sinθ)′, t = 1, 2, …, p} denote p different directions, and is to denote discrete positions, m = pq, z = (x, y)′, where q is the total number of pixels of the input image. Consequently, the function is called a class of AHF with a specific ξ. For an image , 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, , with n = n1n2, m = pq. We compute coefficient ω and get the high resolution image with equation , where s is an upscaling factor and with size N = s2n1n2.

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 θ at the locations c with two ξ1, ξ2. One is big parameter ξ1 to represent smooth components (forming Ψ1), another one is the smaller ξ2 to represent non-smooth components (forming Ψ2). Thus, the vector-form image L can be approximated by the following discrete formula: Since non-smooth components are sparse, l1 regularization could be given on the coefficient β2. The optimal model is expressed as: The Eq (6) is solved via ADMM [28-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, which can be obtained from by the Eq (7), and the ring artifacts could be reduced significantly.

where where G 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: , s: upscaling factor. τ: maximum number of iteration.

Output: high-resolution image

According to Eq (4), construct matrices on coarse grids, on fine grids construct , where N = s2n.

Initialization: L(1) = L.

for k = 1:

Compute the coefficients:

Update the high-resolution image:

Downsampling H( to coarse grid: .

Compute residual:

end

Assemble the high-resolution outputs:

Compute the final high-resolution image:

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 non-smooth components.

3.1 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 according to Eq (4). Ψ(i = 1, 2, ⋯, l) represents smooth components formed by larger ξ. Φ(j = 1, 2, ⋯, k) represents non-smooth components formed by smaller ξ. And are the corresponding representation coefficients with Ψ(i = 1, 2, ⋯, l) and Φ(j = 1, 2, ⋯, k). After computing the representation coefficients, we apply them into following the equation to obtain the high-resolution image: where 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 l1 regularization on β(j = 1, 2, ⋯, k) to characterize this feature. Thus, the optimization model can be written as following: where are regularization parameters.

To simplify Eq (11), we set , Ψ = (Ψ1, Ψ2, ⋯, Ψ), Φ = (Φ1, Φ2, ⋯, Φ), u = (α1, α2, ⋯, α, β1, β2, ⋯, β), Thus, Eq (11) could be rewritten as following:

As l1 term is not differentiable, we make a set of variable substitutions for Bu(j = 1, 2, ⋯, k) and then rewrite Eq (12) as: where p(j = 1, 2, ⋯, k) are the substitution variable. In particular, the Eq (13) is separable, w.r.t (u, p)(j = 1, 2, ⋯, k). Many methods are used to solve the l1 problem (see [28-33]). The work in [7] solve the l1 problem via alternating direction method of multipliers(ADMM) [31]. Here, we solve Eq (13) via block-wise ADMM [26].

The augmented Lagrangian of Eq (13) is where b(j = 1, 2, ⋯, k) are Lagrangian multipliers with proper size. The problem of minimizing could be solved by iteratively and alternatively via the following subproblems:

Hence, the solution to the problem (11) is solved by block-wise ADMM, which is shown in Algorithm 2:

Algorithm 2

Input: Given the low-resolution image , , λ(i = 1, 2, ⋯, l), γ > 0, τ > 0(j = 1, 2, ⋯, k)

Output: coefficient u

Initialization:

t ← t + 1

u( ← solve subproblem (15) for

← solve subproblem (16) for

.

The u-subproblem is a smooth quadratic problem, we can solve it by least squares method: where , The p-subproblem(j = 1, 2, ⋯, k) has a closed form solution for each (p) (see [33]), where shrink(a, b) = sign(a) max(|a| − b, 0) and 0.(0/0) = 0 is assumed.

By Algorithm 2, we have computed the representation coefficients α, β(i = 1, 2, ⋯, l; j = 1, 2, ⋯, k). We can get the high resolution image via Eq (10).

3.2 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 , where 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 (forming ) to depict the residual image. Since the non-smooth components are also sparse in the residual. We apply l1 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 l1 norm, we solve it by ADMM scheme. We make two variable substitutions for , and rewrite Eq (19) as: We rewrite Eq (20) as: where , A = (I, 0), and B = (0, I). And 0 is zero matrix. I is identity matrix. Since the optimization Eq (21) is separable, w.r.t (β*, u1, u2). The augmented Lagrangian of Eq (21) is where 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: where , and The subproblem u1 and u2 have a closed form solution as mentioned in Algorithm 2. Thus, the solutions to the subproblems u1 and u2 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.

Algorithm 3 (Iterative refinement)

Input: residual image R, μ1, μ2, , ,

Output: coefficient β

Initialization:

t ← t + 1

β*( ← solve subproblem (23) for ,

← solve subproblem (24) for

← solve subproblem (25) for

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

Take different behaviours of Eqs (9) and (19) into consideration, we get the modified single super-resolution based on approximated Heaviside functions and iterative refinement, which is shown in Algorithm 4. The details of the proposed method could be found from the flow chart shown in the Fig 2.

Fig 2

The flow chart of the proposed work.

Algorithm 4

Input: Given the low-resolution image L (a vector form), λ, γ > 0(i = 1, 2, ⋯, l;j = 1, 2, ⋯, k), μ1 > 0, μ2 > 0, , . s: the upscaling factor. t: maximum iterations of the iterative refinement.

Output: high-resolution H

According to Eq (4), construct matrices , on coarse grids, construct on corresponding fine grids where N = s2n.

Compute the coefficients (α, β) according to Algorithm 2:

Get high-resolution image :

Initialization: .

for

Compute the residual image R(: R( = L − DH(.

Compute the coefficients according to Algorithm 3:

Update the high-resolution image:

end

Get non-smooth components:

Compute the final high-resolution image: where conv represents a convolution operator, p is a Gaussian kernel with a small size.

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 in Eq (12). For example, if a LR image is 512 × 512 and we choose 10 different directions, the size of matrix will be equal to (5122 × 10 × 5122). 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 set , , M = (I, 0, 0), , , where 0 is zero matrix, I is identity matrix. The Eq (29) can be rewritten as following.

Since l1 term is not differentiable, we make two variable substitutions p1, p2 for to rewrite the Eq (30) as: The Eq (31) is separable, w.r.t . Here, we solve this problem by block-wise ADMM.

The augmented Lagrangian of Eq (31) is where b1, b2 are Lagrangian multipliers with proper size. The problem of minimizing could be solved by iteratively and alternatively by following subproblems: The is a smooth quadratic problem. We solve it by least squares method as following: where The p-subproblem(j = 1, 2) has a closed form solution Eqs (39) and (40) for each (p)(j = 1, 2): where shrink(a, b) = sign(a) max(|a| − 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)

Input: Given the low-resolution image L (vector form), λ1, γ1, γ2 > 0, μ1 > 0, μ2 > 0, , , s: the upscaling factor. t: maximum iterations of the iterative refinement.

Output: high-resolution H

According to Eq (4), construct matrices on coarse grids, construct on corresponding fine grids where N = s2n.

Compute the coefficients according to Algorithm 2:

Get high-resolution image :

Initialization: .

for

Compute the residual image R(: R( = L − DH(.

Compute the coefficients according to Algorithm 3:

Update the high-resolution image:

end

Get non-smooth components:

Compute the final high-resolution image: where conv represents a convolution operator, p is a Gaussian kernel with a small size.

Experiments

4.1 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.

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.

Table 1

Comparison with AHFM algorithm and MAHFM algorithm(factor = 2).

Algorithms	PSNR/RMSE	woodpecker	workman	flowers	butterfly	children	blackbutterfly
AHFM	RMSE	3.3077	5.5733	10.3057	11.7942	3.7469	6.4146
	PSNR	38.4155	33.2085	27.8692	26.6974	036.6573	31.9874
MAHFM	RMSE	3.0950	5.4714	6.7661	11.2689	3.5336	6.2337
	PSNR	38.6273	33.3689	31.5241	27.0932	37.1665	32.2359

Table 2

Quantitative results by different super-resolution algorithms (factor = 2).

Images	PSNR/RMSE	bicubic	07’TIP [34]	08’TOG [22]	10’TIP [12]	AHFM [7]	MAHFM
leaves	RMSE	15.6427	16.5913	16.1366	15.1104	14.8658	14.8190
	PSNR	24.2446	23.7332	23.9746	24.5433	24.6871	24.7145
father	RMSE	6.1870	8.5823	6.8589	5.6152	5.6024	5.5610
	PSNR	32.3012	29.4587	31.4057	33.1435	33.1633	33.2354
comic	RMSE	12.6995	16.4580	13.9571	11.1757	11.5335	11.1000
	PSNR	26.0551	23.8033	25.2349	27.1653	26.8916	27.2243
babyface	RMSE	4.6689	5.9649	5.1108	4.4681	4.4029	4.3713
	PSNR	34.7466	32.6188	33.9610	35.1284	35.2560	35.3185
lions	RMSE	7.3928	10.4452	8.2584	6.7658	6.4465	6.2582
	PSNR	30.7546	27.7524	29.7929	31.5244	31.9444	32.2018
tree	RMSE	13.2687	15.8486	13.8727	12.7048	12.4237	12.3883
	PSNR	25.6742	24.1310	25.2876	26.0515	26.2458	26.2705
peppers	RMSE	5.6703	7.1575	5.9815	5.1910	4.9485	4.6897
	PSNR	33.0587	31.0356	32.5945	33.8258	34.2413	34.7080

Table 3

Quantitative results by different super-resolution algorithms (factor = 3).

Images	PSNR/RMSE	07’TIP [34]	08’TOG [22]	10’TIP [12]	14’TIP [10]	AHFM [7]	MAHFM
barbara	RMSE	15.7546	11.9587	12.2959	22.1094	12.0127	11.8880
	PSNR	24.1826	26.5771	26.3356	21.2393	26.1975	26.6286
view	RMSE	14.4619	10.9190	11.2941	21.9687	11.3169	11.2920
	PSNR	24.9263	27.3672	27.0731	21.2947	27.0562	27.0754
baboon	RMSE	20.3930	17.8625	17.8047	25.4575	18.1387	17.7188
	PSNR	21.9412	23.0920	23.1201	20.0145	23.1621	23.1621
fish	RMSE	18.4114	17.7298	17.5758	24.0806	17.3458	17.2849
	PSNR	22.8291	23.1567	23.2325	20.4975	23.3469	23.3775
tool	RMSE	8.5852	5.1631	5.2033	18.1161	4.8748	4.8210
	PSNR	29.4558	33.8726	33.8053	22.9695	34.3716	34.4681

4.2 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.

Fig 3

Results of “baboon” with the upscaling factor s = 3.

(a) LR image, (b) bicubic, (c) 11’IPOL [2], (d) 14’TIP [10], (e) 10’TIP [12], (f) AHFM [7], (g) MAHFM.

Fig 4

Results of “peppers” with the upscaling factor s = 2.

(a) LR image, (b) ground truth, (c) bicubic, (d) 11’IPOL [2], (e) 14’TIP [10], (f)10’TIP [12], (g) AHFM [7], (h) MAHFM.

4.3 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).

Fig 5

Comparisons between MAHFM and AHFM with factor s = 2.

(a) MAHFM (RMSE = 8.7501, PSNR = 29.2905 dB). (b) AHFM [7] (RMSE = 9.0810, PSNR = 28.9681 dB). (c) residual (For better visualization, we add 0.5 to the intensities of the residual).

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 is the true high resolution image, and ∥ ⋅ ∥ 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.

Fig 6

Comparisons between two iterative refinement methods.

(a) Reerror. (b) Reerror1.

4.4 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 ∈ [0.8, 1], ξ2 ∈ [10−2, 10−1] and ξ3 ∈ [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.

Fig 7

(a) The RMSE trend of the MAHFM algorithm as the different parameter varies when other parameters are roughly estimated (To simplify the experiment, e1, e2, e3 represents ξ1, ξ2, ξ3, respectively in this figure). (b) The RMSE trend of the MAHFM algorithm as the parameter ξ2 varies, where ξ1 = 0.8, ξ3 = 10−4.

Table 4

Parameters of MAHFM algorithm.

ξ1	ξ2	ξ3	λ1	γ1	γ2	τ1
1	6 × 10−2	1 × 10−4	5 × 10−4	1 × 10−3	1 × 10−5	3 × 10−3
τ2	ξ~_1	ξ~_2	μ1	μ2	τ~_1	τ~_2
1 × 10−4	1 × 10−4	1 × 10−5	1 × 10−4	5 × 10−3	1 × 10−5	1 × 10−6

4.5 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.

Fig 8

(a) Computation time of Algorithm 2 (l = 1, k = 2) and Algorithm 3 vs. upscaling factor for the low-resolution image “butterfly” with LR size 60 × 60. (b) Computation time of Algorithm 2 (l = 1, k = 2) and Algorithm 3 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, to reduce the instability of Matlab, the computation time is the average of 10 runs.

Fig 9

(a) Computation time of MAHFM (i.e., Algorithm 5) vs. upscaling factors for the low-resolution image “butterfly” with LR size 60 × 60. (b) 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.

4.6 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 non-smooth 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 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.

Table 5

Comparison with super resolution based on the basis of wavelets and Heaviside function.

Algorithms	RMSE/PSNR	peppers (s = 2)	babyface (s = 2)	barbara (s = 3)	girl (s = 3)	baboon (s = 3)
wavelets	RMSE	6.1932	5.1663	13.0254	6.2360	18.5015
	PSNR	32.2925	33.9518	25.8350	32.2327	22.7867
Heaviside function	RMSE	4.6897	4.3713	11.8880	4.8210	17.7188
	PSNR	34.7080	35.3185	26.6286	34.4681	23.1621

Fig 10

Results of “peppers” with the upscaling factor s = 2 via the representation of Haar wavelets.

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. In particular, we discussed how to choose parameters. Extensive examples have been provided in the experiments to show the effectiveness of the MAHFM algorithm both quantitatively and visually.