Image Magniﬁcation Based on Bicubic Approximation with Edge as Constraint

: Image magniﬁcation can be reduced to construct an approximation surface with data points in the image while keeping image details and edge features. In this paper, a new image magniﬁcation method is proposed by constructing piecewise bicubic polynomial surfaces constrained by edge features. The main innovation includes three parts. First, on the small adjacent area of each pixel, the new method constructs a quadratic polynomial sampling patch to approximate pixels on the small neighborhood with edge features as constraints. All quadric polynomial sampling patches are weighted to generate piecewise whole bicubic polynomial sampling surface. Second, a technique for calculating the error image is proposed: the error image is used to construct a correction surface to improve the accuracy and visual effect of the magniﬁed image. Finally, in order to improve the accuracy of the approximation surface, a technology of balancing polynomial coefﬁcients is put forward. Experimental results show that, compared with other methods, the proposed method makes better use of the local feature information of the image, which not only improves the PSNR/SSIM numerical accuracy of the magniﬁed image but also improves the visual effect of the magniﬁed image.


Introduction
Image magnification is one of the important issues in the fields of image processing, computational medicine, computer graphics, computer vision, and virtual reality. It has a wide range of applications in medical image-aided diagnosis, product defect inspection, and visualization. The purpose of image magnification is to increase the number of pixels in the image so as to increase the details of the image and to make its features clearer, helping doctors better understand the details of the lesions. Image magnification methods can be divided into three categories: fitting-based methods, learning-based methods, and other methods.
In a fitting-based method, a surface is constructed from a low-resolution image. Then, more sampling points are obtained by the surface to generate a high-resolution image. Early fitting methods include bilinear interpolation [1] and bicubic interpolation [2,3]. In these methods, the magnified images produced by bicubic interpolation have higher accuracy, while the least square [4] is similar to the two types of interpolation methods. Sometimes, there are obvious sawtooth and ringing, and serious distortion in the magnified image [1][2][3][4]. The common feature of the magnified images by these methods is that the effect is better in flat areas, but it is not ideal in non-flat area (edge and texture). In order to reduce the sawtooth and blur in the magnified image, Zhang et al. proposed a soft the source of the training patch, this method obtains competitive performance by using LR image patch information. Based on error minimization, iterative back projection can significantly improve the quality of the image. Reference [24] proposed an image restoration and amplification method combining nonlocal self-similarity and global structure sparsity. A group of similar image blocks is reconstructed by adaptive regularization techniques based on weighted kernel norms, and a new strategy is used to maintain the global structure. The strategy decomposes the image into smooth components and sparse residuals. The latter uses L1 norm to regularization, which provides a new method for image restoration and magnification. Irani and Peleg [25] continue to degrade the high resolution, to calculate errors, and to optimize the error until the error is less than a certain threshold. The quality of the reconstructed image is significantly enhanced, but there are obvious sawtooth and ringing in the edge area. In [26,27], a nonlocal iterative back projection method was proposed. By searching for similar blocks in a certain area, the average image of the magnified image is filtered and the sawtooth and ringing phenomena can be corrected to improve the correctness of the error image. This method improves the shortcomings of the method [25] and improves the image quality to a certain extent. If the features of an image can be extracted first and the different edges of the image are segmented [28,29], then the surface is fitted based on the features and the edge structure; the accuracy and the visual effect of the magnified image can be significantly improved.
As the image data points are complex, this paper uses a piecewise polynomial surface to fit the given image. A sampling quadratic polynomial patch is constructed on a small adjacent area of each pixel. All the sampling quadratic polynomial patches are weighted averagely to generate the whole sampling surface of the magnified image. In order to improve the quality of the magnified image, this paper generates a magnified image by the whole surface with the correcting surface. The correction surface is constructed by the error image, which is computed by the magnified image. To make the magnified image with higher quality, the constructed quadratic polynomial sampling surface patch is required to satisfy two conditions: high precision and the shape described by the corresponding image block. For satisfying the two conditions, the quadratic polynomial patch is constructed with the edge feature as constraints and it has quadratic polynomial approximation accuracy. A model describing constraint based on the edge feature is proposed. It is inevitable that there are errors in the approximation by the quadratic polynomial patches, especially the approximation effect is not ideal at the complex edge. The correcting surface is constructed to reduce the errors in the magnified image. A method for calculating the error image is proposed. The image magnified twice is compressed twice to obtain a compressed image, and the error image is calculated from the difference between the compressed image and the given image, which provides a new technique for constructing the correcting surface.
The following content of this paper is arranged as follows. The second section is the basic idea and description of the method. The third and fourth sections describe the construction process of our image magnification method. The fifth section compares the results of the magnified image by different methods. The sixth section is the summary of the method and the future work.

Basic Idea and Description of the Method
For a given image P, composed of n × n pixels P i,j , i, j = 1, 2, 3 · · · , n (shown by black dots in Figure 1), it is obtained by sampling from an original scene [17,18] according to the image generation principle. These pixels P i,j , i, j = 1, 2, 3 · · · , n can be regarded as the points (x, y, P) in a 3D coordinate system with coordinates (i, j, P i,j ). For the convenience of the following description, the points (i, j, P i,j ), i, j = 1, 2, 3 · · · , n are regarded as being sampled from a surface P(x, y); then, the sampling formula P s (x, y) of P(x, y) can be written as follows: where 2h is the side length of the sampling square area.
Without loss of generality, the side length of the sampling square region of a given pixel P i,j can be taken as 1, so h = 0.5. Then, there is The essence of image magnification is that more points on the surface P(x, y) are obtained. If the surface P(x, y) can be reconstructed, the image magnification becomes a problem of calculating the point on the surface P(x, y). Because of the complexity of the scene and sampling mechanism, it is impossible to reconstruct P(x, y) accurately. Therefore, only an approximating surface representing P(x, y) can be constructed. That is, the approximating surface is inevitably subject to errors or noise. Because the original scene is usually complex, the corresponding surface P(x, y) is complex. Construction of a whole polynomial surface to fit P(x, y) will not only produce large errors but also generate a large amount of calculation. According to approximation theory, if any given complex function is piecewise continuous and differentiable, it can be approximated with arbitrary precision by piecewise surfaces. Therefore, we construct a sampling patch in a small neighborhood of each pixel, and all patches are weighted to form a whole approximation surface F(x, y). F(x, y) approximates the sampling surface P s (x, y) of P(x, y). Because polynomials has a good mathematical basis, is simple, and is easy to calculate, the polynomial function is widely used in numerical fitting. This paper uses polynomial surfaces to construct surface patches. Based on the above discussion, the process of image magnification in this paper can be described as the following three parts: (a) For each pixel P i,j , i, j = 1, 2, 3 · · · , n, its coordinate in 3D coordinate system is (i, j, P i,j ). In order to simplify the complexity of the problem and to reduce the amount of calculation, the approximate surface patch corresponding to P i,j is constructed with a quadratic polynomial patch on the small area where the pixel P i,j is located. Another reason for constructing approximate patches using quadratic polynomials is that, although biquadratic polynomials can be constructed uniquely from P i,j and the 8 adjacent points, our computation results show that the biquadratic polynomials surface interpolating 9 pixels are often difficult to have the shape represented by the 9 data points. The approximate quadratic polynomial surface patch P i,j (x, y) of the original scene corresponding to P i,j can be written as follows: where u = x − i, v = y − j. a 0 , a 1 , a 2 , a 3 , a 4 , and a 5 are the undetermined coefficients.
(b) Note that the sampling surface patch of P i,j (x, y) on the square area with side length 2h is approximated by a quadratic polynomial patch f i,j (x, y), i, j = 1, 2, 3 · · · , n. The whole sampling surface F(x, y) is generated by the weighted average of f i,j (x, y). Since F(x, y) has an approximation error, correcting surface patches e i,j (x, y), i, j = 1, 2, 3 · · · , n are constructed by the approximation errors. The sampling surface patch f i,j (x, y) is corrected by e i,j (x, y) so as to improve the accuracy and shape of F(x, y).
(c) The magnified image is obtained by calculating the sampling points of F(x, y).
The discussion above showed that the accuracy and shape of the whole sampling surface F(x, y) is determined by the quadratic polynomial patch f i,j (x, y), i, j = 1, 2, 3 · · · , n; hence, the key point is to construct f i,j (x, y), i, j = 1, 2, 3 · · · , n. Following, we will discuss how to construct the quadratic polynomial sampling patches f i,j (x, y), i, j = 1, 2, 3 · · · , n.

Constructing Quadratic Sampling Patches
Firstly, we construct a quadratic polynomial sampling patch f i,j (x, y) to approximate the pixels on the small neighborhood with the edge features as constraints. Considering the symmetry, we construct P i,j (x, y) (Equation (3)) with P i,j and 8 pixel points adjacent to P i,j (as shown in Figure 1, where P i,j is the center of 9 points) for f i,j (x, y). That is, f i,j (x, y) is determined by 9 pixels P i+l,j+k , l, k = −1, 0, 1. The constructed f i,j (x, y) should have the quadratic polynomial precision and the shape described by 9 pixels. Ideally, f i,j (x, y) passes P i,j and 8 adjacent points, i.e., passes (i + l, j + k, P i+l,j+k ), l, k = −1, 0, 1. Replace P(x, y) in Equation (1) with P i,j (x, y) in Equation (3) and integrate to get the sampling surface patch f i,j (x, y) on the square area with side length 2h as follows: where u = x − i, v = y − j and where a 0 , a 1 , a 2 , a 3 , a 4 , and a 5 are the undetermined coefficients. In Equation (4), taking h = 0.5 and (x, y, f i,j (x, y)) = (i + l, j + k, P i+l,j+k ), l, k = −1, 0, 1 gets the following 9 equations: P i+l,j+k = a 0 (l 2 + 1/12) + a 1 lk + a 2 (k 2 + 1/12) + a 3 l + a 4 k + a 5 , l, k = −1, 0, 1.
Equation (5) is the integral of Equation (2) over the 9 unit squares as shown in Figure 1. The accuracy and shape of the sampling surface F(x, y) is completely determined by the sampling patch f i,j (x, y), i, j = 1, 2, 3, · · · , n (4). In order to improve the accuracy of F(x, y), the technology of balancing polynomial coefficients is put forward. Next, we will discuss how to calculate the 6 undetermined coefficients a 0 , a 1 , a 2 , a 3 , a 4 , and a 5 using the 9 equations in Eqauation (5).
In Equation (5), there are 9 equations and 6 unknowns. The common method is to determine 6 unknowns by the least square method. Because 9 points in the neighborhood of P i,j may belong to different surfaces in the original scene, the role of 9 points should be different in the construction of f i,j (x, y). To make f i,j (x, y) closer to 9 pixels, we use the weighted least squares to determine 6 unknowns. Different weights are used to make each pixel have different effects on determining 6 unknowns. In order to determine the 6 unknowns better, we first analyze the coefficient influence of f i,j (x, y) on the shape of the surface. Because f i,j (x, y), i, j = 1, 2, 3 · · · , n are local surface patches, they form a whole sampling surface F(x, y) by the weighted average. f i,j (x, y) only plays a large role on the F(x, y) in its adjacent area. According to the characteristics of Taylor expansion, in the local area, the influence of constant term, line term, and quadratic term of f i,j (x, y) on the shape of F(x, y) decrease in turn. Therefore, we will use different strategies to determine these unknowns. The six coefficients are divided into three groups according to the constant term, the first term, and the quadratic term. The calculation of the three groups of coefficients will be discussed below.
Let us firstly discuss the calculation of the coefficients a 3 and a 4 ; then discuss the calculation of a 0 , a 1 , and a 2 ; and finally calculate a 5 . The process of calculating a 0 , a 1 , a 2 , a 3 , a 4 , and a 5 is divided into three steps, which not only reduces the amount of calculation but also improves the accuracy of the patches. The following discusses how to calculate the coefficients a 3 and a 4 .

Calculating the Coefficients of First Term
Based on Equation (5), the following four equations can be obtained [18]: where, Equation (6) has four equations with two unknowns a 3 and a 4 . Using the weighted least squares method, a 3 and a 4 are determined by minimizing the objective function G(a 3 , a 4 ) defined by where w k , k = 1, 2, 3, 4 is the weight function. Let us consider how to define the weight function in Equation (7). If 9 pixels P i+l,j+k , l, k = −1, 0, 1 are sampling points on the same object, we could set w 1 = w 3 = 1, w 2 = w 4 = 1/ √ 2, considering the relative positions of the other 8 pixels and P i,j . w 2 and w 4 take smaller weights because their corresponding first-order differences e 2 and e 4 are on the diagonal (as shown in Figure 1). The points relatively far away from P i,j should be assigned a smaller weight. If the 9 pixels are sampling points on the different objects, there are edges in the 9 pixels; that is, at least three pixels are on an edge. To simplify the complexity of the problem, we assume that there are up to four edges in a block of 9 pixels, which are on the horizontal, vertical, and diagonal lines of Figure 1. For ensuring the edge features of the magnified image, the pixels on an edge should play a relatively large role in the determination of a 3 and a 4 , which is equivalent to assigning a relatively large value to the corresponding weight in Equation (7). For example, if P i−1,j−1 , P i,j , and P i+1,j+1 are on the diagonal edge, w 2 should give a relatively large weight. Note δ k , k = 1, 2, 3, 4 is a second-order difference along four directions.
The feature of P i−1,j−1 , P i,j , and P i+1,j+1 on the edge is that |δ 2 | in Equation (8) has a relatively small value. Therefore, assigning a relatively large weight to w 2 is equivalent to w 2 being inversely proportional to |δ 2 |. If P i−1,j−1 , P i,j , and P i+1,j+1 are on the edge, the corresponding |δ 4 | in Equation (8) of the edge's normal direction is usually relatively large, and then, giving a relatively large weight to w 2 is equivalent to w 2 being proportional to |δ 4 |. Based on the above discussion, the weight w k , k = 1, 2, 3, 4 can be defined as follows: where δ k , k = 1, 2, 3, 4 are defined in Equation (8). Experimental results show that β = 0.2, α = 0.5 is preferable in Equation (9).

Calculating the Coefficients of Quadratic Term
Since a 3 and a 4 are computed, there are 9 equations with four unknowns in Equation (5): a 0 , a 1 , a 2 , and a 5 . The third group of coefficients for a 0 , a 1 , and a 2 can be determined based on 9 equations. In order to improve the accuracy of interpolation, two sets of values of a 0 , a 1 , and a 2 are determined.
The final values of a 0 , a 1 , and a 2 are calculated using the weighted averages of the two sets of values.
The following discusses how to determine the first set of values for a 0 , a 1 , and a 2 . Using the 9 equations in Equation (5), we have the following 8 difference equations: In these 8 equations, there are 3 unknowns: a 0 , a 1 , and a 2 . Using the weighted least squares method, the 3 unknowns are determined by minimizing the objective function G(a 0 , a 1 , a 2 ) defined by the following: How to determine the weight function in Equation (10) is discussed below. Edge features and distance play an important role in determining a 3 and a 4 and should also play an important role in determining a 0 , a 1 , a 2 , and a 5 . Therefore, based on the discussion of Equation (9), the weight function can be defined as follows: where w i , i = 1, 2, 3, 4 are defined by Equation (9). The 9 pixels shown in Figure 1 are in horizontal, vertical, and two diagonal directions. For each pixel, it is unreasonable to simply assign the same weight to the two pixels on an edge in the same direction. For example, for horizontal P i−1,j , P i,j , and P i+1,j , if P i−1,j and P i,j are sampling points on the same object and P i+1,j is sampling point on another object, it is reasonable to assign relatively large weights to P i−1,j and P i,j and to assign a relatively small weight to P i+1,j , while according to Equation (11), P i−1,j and P i+1,j are assigned the same weight. For this reason, the new weights are defined as follows: where w l,k , l, k = −1, 0, 1 on the right side of Equation (12) is defined by Equation (11). For the convenience of the following discussion, the first set of values for a 0 , a 1 , and a 2 computed by Equations (10)-(12) are denoted as a A 0 , a A 1 , and a A 2 . In the following, we discuss how to determine the second set of values for a 0 , a 1 , and a 2 . From Equation (5), one gets the following three equations: These three equations correspond to the four pixels, P i,j , P i+1,j , P i,j+1 , and P i+1,j+1 , which form a square in Figure 1. When a 3 and a 4 are determined, a 0 , a 1 , a 2 can be calculated from Equation (13) and denoted as a 0 0 , a 0 1 , a 0 2 . Similarly, the other three sets of values for a 0 , a 1 , and a 2 can be calculated from the pixels forming the other three square in Figure 1, which are denoted as a 1 0 , a 1 1 , a 1 2 ; a 2 0 , a 2 1 , a 2 2 ; and a 3 0 , a 3 1 , a 3 2 , respectively. Let us consider how to determine the second set of values for a 0 , a 1 , a 2 by a k 0 , a k 1 , a k 2 , k = 0, 1, 2, 3. In order to make the quadratic polynomial as simple as possible, |a 0 |, |a 1 |, and |a 2 | should be as small as possible. According to the approximation theory, the interpolation function should be as simple as possible on the premise of satisfying the approximation accuracy. Therefore, a 0 is determined by weighted average of a k 0 , k = 0, 1, 2, 3, and the weight corresponding to a k 0 should be inversely proportional to |a k 0 |. Now, the second set of values for a 0 , a 1 , a 2 , denoted as a B 0 , a B 1 , and a B 2 , is defined by where w l,k = 1/(|a k l | + λ), l, k = 0, 1, 2 are weight functions, λ is a very small positive value. We do not simply use Equation (15) to compute the second set of values for a 0 , a 1 , a 2 . Equation (15) is a simple weighted average. The experiment results and the following discussion also showed that the magnified image is not good using Equation (15) to compute the second set of values for a 0 , a 1 , and a 2 .
Experimental results also show that the quadratic polynomial defined by Equation (14) does integrate the advantages of the four quadratic polynomials, and its surface shape is closer to the shape represented by the 9 pixels in Figure 1. From a k 0 , a k 1 , and a k 2 , k = 0, 1, 2, 3, 4 and Equation (14), we can define 6 patches, respectively, denoted as f k (x, y), k = 0, 1, 2, 3, 4 and f 5 (x, y), which are all approximation to the 9 pixels in Figure 1. The approximation errors of six patches to 9 pixels are denoted as E k i,j , k = 0, 1, 2, 3, 4, 5, i, j = 2, 3 · · · , n − 1, respectively, and hence, 6 error images E k , k = 0, 1, 2, 3, 4, 5 can be obtained by E k i,j , k = 0, 1, 2, 3, 4, 5, i, j = 2, 3 · · · , n − 1. Using Baboon and Kod in the 12 typical images in Section 5, two error images can be generated with the 6 patches. Figure 2 is a histogram of two error images, where the top-down histograms correspond to Baboon and Kod images, respectively. In Figure 2, e k (n) denotes the number of pixels with value n in the error image E k , k = 0, 1, 2, 3, 4, 5. For the sake of clarity, each histogram is divided into two parts. Since the number of pixels with value 0 is too small compared with the number of pixels with value 1, the former is not displayed to make the histogram visually clear. In Figure 2, the left figures are the histograms of pixels with error values of 1-20 and the right figures are the histograms of the rest error pixels. It is obvious from the curve e 5 (n) in Figure 2 that, in the error image E 5 generated using a 0 , a 1 , a 2 by Equation (14), the pixels with small error value are the most and the pixels with large error value are the least, while it is also obvious from the curve e 4 (n) that, compared with the error images E 5 , the pixels with large error value in the error image E 4 generated using a 0 , a 1 , a 2 by Equation (15) are relatively larger and the pixels with small error value are small.
From Equations (10) and (12), the first set of values a A 0 , a A 1 , and a A 2 are computed. From Equations (13) and (14), the second set of values a B 0 , a B 1 , and a B 2 are computed. The final values of a 0 , a 1 , and a 2 , are computed by the weighted average of the two sets of values a A 0 , a A 1 ,a A 2 , and a B 0 , a B 1 ,a B 2 . Similar to the discussion of Equation (14), the coefficient with small absolute value is assigned a large weight. That is, the coefficients a 0 , a 1 , and a 2 are defined by where w l,A = (a B l + λ)/(|a A l | + |a B l | + 2λ) is weight function, l = 0, 1, 2, λ is a very small positive value.

Calculating the Constant Term
Now, how to compute a 5 is discussed. a 0 , a 1 , a 2 , a 3 , and a 4 are obtained above, so that only one unknown a 5 is left in Equation (5). a 5 is determined by the weighted least squares method from the 9 equations in Equation (5). That is, a 5 is determined by minimizing the objective function G(a 5 ) defined by where g(l, k) = a 0 (l 2 + 1/12) + a 1 lk + a 2 (k 2 + 1/12) + a 3 l + a 4 k.

Constructing the Whole Sampling Surface
In this section, we first discuss the construction of the initial whole sampling surface F(x, y) and then discuss how to construct a correcting surface to improve the accuracy of F(x, y). Finally, we perform an error analysis on F(x, y).
In general, the weight functions can be defined by the following bilinear basis functions: In the unit square area [i, i + 1] × [j, j + 1], if the approximation accuracy of the four patches, f i,j (x, y), f i+1,j (x, y), f i+1,j+1 (x, y), and f i,j+1 (x, y) is about the same, the ideal result can be obtained with the weight function defined by the above equation. If the approximation accuracy of the four patches is significantly different, it is more reasonable to define the weight functions with the following rational functions.
where α l,k , l, k = 0, 1 is the parameter, which determines the importance of the corresponding weight w i+l,j+k (u, v).

Constructing Correcting Surface
The magnified image obtained from the initial whole sampling surface F(x, y) (Equation (19)) is inevitably subject to errors. This section discusses constructing a correcting surface to improve the accuracy of F(x, y) (Equation (19)). Firstly, the error image is computed, and the correcting surface is constructed by the error image.
The error image is calculated as follows. First, the initial whole sampling surface F(x, y) is sampled to obtain the doubled image p, and the pixels in p are denoted as p i,j , i, j = 1, 2, 3 · · · , 2n. Let Ideally, for a given pixels P i,j , i, j = 1, 2, 3, · · · , n, P i,j and given pixels P i,j should satisfy the relationship P i,j = P i,j . Therefore, the error pixels are defined by the difference between P i,j and P i,j as follows: From the error image E i,j , i, j = 1, 2, 3, · · · , n, a quadratic correcting sampling patch e i,j (x, y) similar to Equation (4) can be constructed in the neighborhood of each pixel P i,j , i, j = 2, 3, 4, · · · , n − 1, using the method in Section 3 above: Then, the two quadratic polynomials defined by Equations (4) and (23), respectively, are merged; that is, a k + e k → a k , k = 0, 1, 2, 3, 4, 5. The new quadric sampling surface is still denoted as Equation (4). We still denote the new whole sampling surface as F(x, y). The new error image can be obtained by sampling the new F(x, y), and a new correcting surface patch of the form in Equation (23) can be obtained. Then, Equations (4) and (23) are merged to get a new F(x, y) again. In this way, the surface F(x, y) (Equation (19)) includes bicubic sampling surface patches and correcting surface patches.
Experimental results show that, if Equation (4) corrects more than four times, the correction effect will not be improved much and the amount of calculation will be increased. Therefore, we only correct the quadratic polynomial in Equation (4) four times under comprehensive consideration.

Error Analysis
The Taylor polynomial expansion function is an effective tool for analyzing the approximation error. In this paper, the Taylor polynomial expansion is used for error analysis. For the error of the method in this paper, we have the following theorem.

Theorem 1. The bicubic surface F(x, y) defined by Equation (19) has quadratic polynomial precision.
Proof: If n × n pixels P i,j , i, j = 1, 2, 3, · · · , n of the image are points on the quadratic polynomial surface, f i,j (x, y) (Equation (4)) determined by Equations (4)-(23) has quadratic polynomial precision due to the uniqueness of the constructing result. Since all f i,j (x, y), i, j = 1, 2, 3, · · · , n have quadratic polynomial precision, the surface patches F i,j (x, y) by their weighted average of f i,j (x, y) also have quadratic polynomial precision. Thus, the bicubic surface F(x, y) formed by the combination of the surface patches F i,j (x, y) also has quadratic polynomial precision. This completes the proof.

Experimental Results
The experiment consists of three parts: numerical accuracy (PSNR (Peak Signal-to-Noise Ratio), the most common and widely used image objective evaluation index based on error-sensitive image-quality evaluation, and SSIM (Structural Similarity), a measure of the similarity of two images measuring image similarity from three aspects: brightness, contrast, and structure), visual effects, and runtime. The low-resolution image (LR) is obtained by averaging downsampling of the high-resolution image (HR), and the magnified images of the LR images are compared with the HR images. The images used for comparison are 12 typical images, BSD200, Urban100, and the aliasing image set (AIS) [30]. These 12 typical images are commonly used in the experiment of image magnification. In order, they are Barbara, Baboon, Boat, Chest, Couple, Crowd, Dollar, Goldhill, Kod, Lake, Lenna, and Peppers, as shown in Figure 3. The size of each typical image is 512 × 512 pixels (HR). BSD200 uses only 100 images for learning, denoted as BSD100. AIS is from the website with http://www.cgl.uwaterloo.ca/csk/projects/alias/. The 8 methods used for comparison are bicubic, DCCI [7], ICBI [8], SelfxSR [15], NARM [13], SISRRFI [23], LPRGSS [24], and the method (New).

Numerical Accuracy (PSNR and SSIM)
First, PSNR and SSIM are used as criteria to compare the magnified images. The 8 methods are used to magnify LR images twice, and the PSNR and SSIM values of the magnified images are calculated using HR images. The PSNR and SSIM values of the magnified images of 12 typical images are shown in Table 1, and the average PSNR and SSIM values of the images in BSD100, UrBan100, and AIS are shown in Table 2, where the boldface indicates the maximum PSNR and SSIM values. The last row in Table 1 is the average of the PSNR and SSIM values for each method. When PSNR is used to measure image quality, the larger the PSNR value, the better the image quality, while when the image quality is measured in terms of SSIM, the larger the SSIM value, the higher the similarity between the magnified image and the HR image. The results in Table 1 show that, in the 12 magnified images, except for Barbara, the PSNR and SSIM values of the new method are the largest and the average PSNR and SSIM value of the new method is also the largest. The results in Table 2 show that, for BSD100, UrBan100, and AIS image sets, the average PSNR and SSIM values of the new method are the largest except for the SSIM value for AIS. Therefore, for the four sets of comparison data, the image quality and similarity of the magnified images by the new method are the best among the 8 methods used for comparison.

Visual Effect
Then In order to facilitate the comparison of visual effects, we have marked some visually distorted parts on the magnified image, as shown in the red box in Figure 4-8. Figures 4 and 5 show the results of the 8 methods for magnifying the two typical grayscale images. In Figure 4, it can be seen that the images (f) and (i) obtained by NARM and new method, respectively, are closest to the HR image (a); the images (c) and (d) obtained by DCCI and ICBI methods are relatively smooth at the edges and producing obvious distortion of lung tissue information; and the images (b), (g), and (h) obtained by the BiCubic, SISRRFI, and LPRGSS, respectively, have sawtooth problems at the edge. In Figure 5, the images (c), (d), and (f) obtained by DCCI, ICBI, and NARM methods, respectively, are relatively smooth at the edges, resulting in reduced detail information. The images (b), (e), (g), and (h) obtained by the bicubic, SelfxSR, SISRRFI, and LPRGSS methods, respectively, have different degrees of sawtooth problems at the edges and the edges are very rough. Overall, the image by the new method (i) is closest to the HR image (a). Figures 6 and 7 show the results of 8 methods for magnifying two color images in BSD100. For the images in Figure 6, there are small differences between these images, the images (c), (d), and (f) obtained by DCCI, ICBI, and NARM methods, respectively, are relatively smooth, and the texture details of the branches are missing; the images (b), (e), (g) and (h) obtained by the bicubic, SelfxSR, SISRRFI, and LPRGSS methods, respectively, are slightly distorted, such as some details of the branches on the far right side of the tail are missing; and by comparison, the image (i) obtained by the new method is closest to the HR image (a). In Figure 7, the images (c), (d), and (f) obtained by the DCCI, ICBI, and NARM methods, respectively, have shallow wrinkles in the corners of the eye of the elderly and loss of texture on the left side of the eyelids; the image (e) obtained by SelfxSR has the defect of deepening wrinkles in the elderly and does not look real; and in the images (b), (g) and (h) obtained by the bicubic, SISRRFI, and LPRGSS methods, respectively, there are severe jagged edges in the ears and eyelids of the elderly. Overall, the image (i) by the new method is closest to the HR image (a). Figure 8 shows the results of 8 methods for magnifying one aliasing image in AIS. Experimental results show that, for the aliasing images in AIS, the magnified images by the 8 methods have very small visual differences. In Figure 8, relatively speaking, the images (h) and (i) obtained by LPRGSS and the new method are closest to the HR image (a); the images (b), (e), (f), and (g) by bicubic, SelfxSR, NARM, and SISRRFI, respectively, are the second; and the images (c) and (d) by DCCI and ICBI are third. Obviously, in Figure 8, all eight images have the problem of losing details and distortion. Table 2, Figure 8 shows that, for the images in AIS, these 8 methods generally can not get good results; one of the main reasons is that the images obtained by downsampling may destroy the details of the original aliasing image. For example, Figure 9 is the downsampling of the aliasing image (a) in Figure 8, and these two images are obviously different.

Runtime
We compare the runtime of the 8 methods, and the comparison results are shown in Table 3, where the time unit is seconds. In Table 3, the runtime of the new method is basically the same as the bicubic method. The time in Table 3 of the new method is the time used to magnify the image, but the new method takes 0.207 s to construct the approximation surface. The approximation surface only needs to be constructed once. After the surface is constructed, the image can be magnified by any multiple, and the time taken is basically the same as that of the bicubic method.

Discussion
The results in Tables 1-3 and Figures 4-8 show that, among the 8 methods in the comparison, (1) overall, the quality of the magnified images by the new method is the best, followed by LPRGSS method; (2) the bicubic method has the shortest running time, followed by the new method, these two methods are suitable for real-time image magnification; (3) for areas without edges and textures in the image, there is little visual difference between the magnified images by the 8 methods, the 9 images in Figure 6 are examples. In fact, for areas without edges and textures in the image, bicubic method, as the most classical basic image magnification method, can produce the same PSNR and SSIM values as the new method, while for areas with edges and textures in the image, bicubic method generally produces image blocks with poor visual effect. The images (b) in Figures 4, 5, and 7 are examples. The main reason is that the bicubic method constructs each surface patch by passing 16 image pixels, which means that, when bicubic method is used, there is no degree of freedom in constructing bicubic surface; the constructed surface is unique. If the pixels are sampled from different objects, the surface patch constructed by bicubic method will produce larger oscillations, resulting in jagged or blurred pixels at the edges of the image. For magnifying the image without edges and textures, the new method has the advantages of bicubic method. Because the new method uses edge features as constraints to construct quadratic polynomial surface, for magnifying the image with edges and textures, the new method gets better magnified images than bicubic method.
We also tested the images in data sets set5 and set14 and database BSD300. The test results were basically the same as those of the 12 typical images, the images in BSD200, Urban100, and AIS.

Conclusions
It is an effective method to construct a surface to approximate the image and then sampling the surface to generate the magnified image, especially when magnifying medical images. The key problem is how to construct the approximation surface. Based on the complexity of image shape, it is an effective strategy to construct local surfaces and to generate a whole sampling surface by the weighted combination of the local surfaces. In this paper, the second-order difference quotient is used to describe the characteristics of the pixels at the edge, which better distinguishes the characteristics difference between the pixels at the edge and the pixels at the non-edge. In the process of local surface construction, we use the edge features of the image as constraints to construct a quadratic polynomial surface patch on the adjacent area of each pixel. Four adjacent quadratic polynomial patches are weighted to form a bicubic surface. All the bicubic surfaces are put together to form the whole sampling surface. The constructed bicubic surface has quadratic polynomial precision. It is inevitable that there are errors in the magnified image obtained from the whole sampling surface. The correction surface can be formed by the approximation error, so as to correct the errors and to improve the approximation accuracy of the whole sampling surface. Due to the use of edge feature constraints to construct local surfaces, each local surface can better maintain the edge features or texture features near the pixels. Experiment results showed that it is an effective method to magnify the image by constructing local surface patches with edge feature constraints, especially for medical image.
Although the magnified image keeps information such as the edge texture of the image well, there is still a lack of edge and texture information. Moreover, the bicubic surface should have bicubic polynomial precision, but the method in this paper has only quadratic polynomial precision. Therefore, our future research work is to improve the retention of edge and texture information and the interpolation accuracy, so that the magnified image has higher image quality and better visual effect.
In addition, Table 2 and Figure 8 show that, for the aliasing images in AIS, the eight methods are not ideal for enlarging the images. It is necessary to study a new image magnification method according to the characteristics of the aliasing images. Because in the aliasing images, there are many edges and the difference of pixel values at the edge is large. One idea to consider is to magnify the aliasing image by the combination of segmentation and fitting. Firstly, the image is segmented, and then, the piecewise fitting patches are constructed based on the segmentation results so as to reduce the oscillation of the patchs as much as possible. This is our second research work we will do next.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: