Adaptive Interpolation Scheme for Image Magnification Based on Local Fractal Analysis

In this paper, we constructed a new rational fractal interpolation model with scale factors and shape parameters. The model obtains different expressions by changing the scale factor, according with the diversity of image features. On the basis of the constructed model, this paper presents a novel adaptive rational fractal magnification (ARFM) algorithm based on local fractal feature analysis. Firstly, the image is divided into different regions according to an adaptive threshold determined by a calculated fractal dimension, and the scaling factor is calculated based on the relationship between global fractal dimension (GFD) and local fractal dimension (LFD). Secondly, for different regions, different interpolation forms are selected according to regional characteristics. Thirdly, parameter optimization and sub-block selection are studied to further enhance the quality of the magnified image. The experimental results illustrate that the performance of the proposed ARFM algorithm is very competitive compared with the latest magnification algorithms.


I. INTRODUCTION
In recent years, the high resolution (HR) image has been an urgent demand as a response to new and various applications. Without changing the resolution of the image sensor, the magnification technique can obtain a HR image with the visual contents of the original image preserved as much as possible. For example, standard-definition videos can be upscaled to match the high definition screen format with the spatial interpolator in a high-definition television (HDTV) [1].
Interpolation is widely employed in image magnification. Many magnification methods have been introduced. These methods can be roughly divided into linear filtering interpolation (LFI) and edge directional interpolation (EDI). LFI approaches [2]- [6] design a particular interpolation kernel applied to the entire image. In particular, these methods allow The associate editor coordinating the review of this manuscript and approving it for publication was Feng Shao . the image to be resized arbitrarily, which is a key for image upscaling applications. The methods in [7]- [9] can preserve edge structures, but do not preserve image detail very well. Contrarily, the methods in [10]- [12] can preserve texture areas, but suffer from blocking and ringing artifacts in the edges. The Human Vision System (HVS) pays great attention to both edge structures and texture areas. Therefore, it is considered a challenge of the interpolation algorithm to retain the edge structure while keeping the details of the magnified image.
To improve the performance of the aforementioned methods, several adaptive methods have been proposed. Li and Orchard [7] proposed a new edge oriented interpolation (NEDI) method, which estimates local covariance features at low resolution and uses them to guide high resolution interpolation based on resolution invariance in the edge direction. Takeda et al. [8] presented an adaptive steering kernel regression (KR) method. In KR, the dominant direction is estimated by the initial estimate of the image gradient, which is then applied to adaptively guide the local edge structure according to the local edge structure. Zhang and Wu [12] developed a directional filtering and data fusion (DFDF) method where the two orthogonal directions are used for the interpolation of the missing pixel and a linear minimum mean square error estimation is used to fuse the results. Zhang and Wu [13] proposed a soft-decision interpolation method where missing pixels can be estimated in groups rather than one. Based on a 2-D piecewise autoregressive model, this technique is capable of learning diverse scene structures with the ability of estimating parameters in a moving window. Considering that the conventional sparse representation model (SRM) became less effective due to the lack of structural constraint of the data fidelity term on the missing pixels, Dong et al. [14] developed an image interpolation method by nonlocal autoregressive modeling (NARM) which is further embedded in the SRM. In NARM, the coherence between sampling matrix and sparse representation dictionary is reduced, which improves the performance of SRM in image interpolation. Recently, Dong et al. [15] proposed a deep learning method (SRCNN) based on image super-resolution reconstruction, which can directly learn the end-to-end mapping between LR and HR images. Compared with almost current classical magnification methods, it improves the visual effect well. However, from the perspective of human vision, there is still space for improvement in respect to preserve the image details and edge structures. Furthermore, it is the incapability to resize images with arbitrary ratios, which meets the needs of image magnification applications well.
Presently, the rational function has been applied as a new method to image interpolation, which has the ability to approximate the ideal kernel [16]. The reconstructed images based on rational function generally show good visual result with no blocky artifacts and less blurring and ringing artifacts. Based on an adaptive osculatory rational interpolation kernel function, Hu and Tan [16] presented a method for preserving edges. Compared with linear polynomial interpolation kernel functions, this function can more accurately approximate the ideal interpolation because it is established on the approximation to the ideal interpolation kernel function. Carrato and Tenze [17] proposed an interpolator with good interpolation performance not only on synthetic images but also on real-world images, where the interpolation of edge sensitive data is implemented by the rational operator to obtain images without artifacts. Liu et al. [18] developed an adaptive interpolation function with weight whose coefficients can keep the edge attributions through adaptive calculation of distance, gradient, and difference quotient based on point sampling. Although rational interpolation function can eliminate ringing artifact and other artificial marks efficiently and has capabilities of efficiently preserving texture features, it shows powerless in dealing with edge structure.
Fractal analysis has been a powerful tool in image processing due to that the fractal can effectively describe texture details. They have found application in classification [19]- [21], segmentation [22], [23], synthesis [24] and other important texture problems [25], [26]. At present, the fractal is widely used in image super resolution. The fractal method is one of the techniques with the most potential, which can achieve a high-resolution enhancement for the sharpness across edges. Generally, global fractal analysis can fully describe the texture features because the global self-similarity shown in the texture of the image. However, the image which has a complex geometric structure usually contains texture and non texture region, it means that image measurement has local adaptation. For the representation and analysis of the image, the local fractal is more suitable than the global fractal. Thus, as for image magnification, according to local fractal feature, various interpolation forms which coincide with image features should be chosen to estimate the unknown point. The magnification method based on local fractal analysis is an effective means to improve the quality of the interpolated image.
To preserve the edge structure and the details of the image and resize the image with arbitrary ration simultaneously, we propose a novel rational fractal interpolation model on the basis of previous researches on rational spline [27]- [29], which can be uniquely identified by the values of scaling factor and shape parameters. The expression of the model varies with the scale factors, which accords with diverse image features. Rational fractal form and rational form are used for fractal feature salient region (FFSR) and non-fractal feature salient region (NFFSR) of the image, respectively. Shape parameters can be optimized to further enhance the quality of the magnified image. Hence, we present a novel adaptive rational fractal magnification (ARFM) algorithm based on local fractal feature analysis. First, based on the calculated fractal dimension, an adaptive threshold is determined, which can be used to divide the image into different regions. Second, different interpolation forms are selected in different regions according to regional characteristics. Third, the scaling factor is determined through the calculation of fractal dimension. Finally, the quality of the interpolation image is improved by a parameter optimization technique. The schematic diagram of the proposed ARFM algorithm is shown in Fig. 1.
The major contributions of this work are summarized as follows: (1) We construct a new type of rational fractal interpolation model, which is organic unity that integrate rational fractal model with rational spline model. (2) According to image features, a new adaptive threshold selection method is employed to regional division, which is based on the distribution characteristics of local fractal dimension. (3) Based on the relationship between scaling factor and fractal dimension, a new method for accurately calculating scaling factor is proposed.
The remainder of this paper is organized as follows: Section II introduces the basic knowledge of fractal and rational functions briefly. In section III, based on the classical binary rational fractal interpolation method, a novel binary rational fractal interpolation method is proposed. In section IV, a hybrid interpolation algorithm with fractal VOLUME 8, 2020 dimension and local adaptive threshold is proposed. Finally, the quality of the interpolation image is further improved by an optimization technology. In section V, the validity of the algorithm is discussed.

A. REVIEW OF FRACTAL
A typical fractal can be constructed with the following strategies. Firstly, we decompose a geometry G into N similar copies, each of which is scaled down by a multiple s. Similarly, the decomposition process can be applied to each of the N copies. Then, through repeating the similar decomposition infinitely, a fractal F can be obtained. Fractal geometry contains the feature that the Hausdorff-Besicovitch dimension exceeds the topological dimension strictly. Such characteristics determine the fractals with infinite complexity. Besides, fractals are self-similar, that is, every single part of the object is similar to the whole.
Recently, some definitions of fractal dimensions have been proposed. Among these, what we know and cite is the Hausdorff-Besicovitch dimension, the box-counting dimension, etc. All these methods have a common point that they are measured for a different scale. The fractal dimension must express the behavior of the measure as ε → 0. In fact, a measure M ε (F) of a set F must generally obey the power law: where c is a constant and D represents the fractal dimension of F. The value of D can be acquired from: ], as follows: where This interpolation is a rational cubic interpolation based on function value and derivative values, which satisfies the following conditions, It is obvious that the interpolating function ; y j , y j+1 ] is defined by using the x-direction function P * i,j (x) as follows [30]: is a bivariate rational interpolation based on function values and partial derivative values, which satisfies the following conditions.
Obviously, this form of the interpolating function

III. CONSTRUCTION OF RATIONAL FRACTAL INTERPOLATION FUNCTION
Rational fractal interpolation functions are the basis of the proposed blending algorithm. This section proposes a new rational fractal interpolation based on the classical binary interpolation function.
], as follows: where ; y j , y j+1 ] is defined by using the x-direction function P * i,j (x) as follows: (4) VOLUME 8, 2020 Denote ,j , s i,j is treated as vertical scaling factor of the iterated function system and |s i,j | < 1.
Thus, the iterated function system {F; (φ i (x), ϕ j (y), F i,j (x, y, z)) : i ∈ I , j ∈ J } defined above admits a unique graph of continuous functions ψ(x, y) which is called attractor G. The bivariate rational fractal interpolation function (BRFIF) is defined by Eq. (3), which has the following form: Remark : The interpolation model is uniquely determined by scaling factor s i,j and shape parameters α i,j , β i,j . When there is scaling factor s i,j = 0 for all i ∈ I , j ∈ J , the Fractal Interpolation Function (FIF) ψ(x, y) is consistent with the bivariate rational interpolation function P(x, y). When there is scaling factor s i,j = 0 for all i ∈ I , j ∈ J , the FIF ψ(x, y) is rational fractal interpolation function. It means that the model has advantages over the current interpolation schemes in terms of flexibility and diversity.

IV. RATIONAL FRACTAL INTERPOLATION FOR IMAGE MAGNIFICATION A. LOCAL FRACTAL ANALYSIS
Fractal dimension (FD), also known as a global fractal dimension (GFD), is an important concept in fractal theory. FD and partial shape dimension (LFD) are powerful tools to describe the texture, which are closely related to human's recognition of image roughness. FD represents the complexity of the entire image, while LFD only describes the complexity of a single block [31]. ε− blanket method is widely used in dimension calculation. The fractal dimension calculated by ε− blanket method can cover the whole range of dimensions. Next, we introduce the calculation of FD and LFD using ε− blanket method.

1) CALCULATION OF FRACTAL DIMENSION
The image can be thought as function H (i, j), the covering blanket is defined by the upper surface T ε and the lower surface B ε . Initially, the function H (i, j) is given.
where d(i, j, m, n) is distance of two points. Expansion and corrosion technology are applied in image processing. Choosing size of structure elements, the corresponding surface area can be obtained by the following formula: where T ε represents the area of the upper surface, B ε represents the area of the lower surface lower. Surface area is expressed as: where A(ε) is the area with thickness of ε. When plotting A(ε) versus ε on a log-log scale, the straight line of slope as follows: 2) THRESHOLD DETERMINATION As a simple texture detection method, threshold method can detect the sharp change of gray value, which is the most obvious features of the texture. Obtaining an accurate threshold is the key of region division, which directly affects the performance of the interpolated algorithm. The algorithm treats the non-fractal feature salient region as a fractal feature salient region when the value of the threshold is small. Conversely, the fractal feature salient region is considered to be a nonfractal feature salient region when the value of the threshold is large. The automatic threshold based on image features is an automatic calculation of the threshold and has been widely used. The equation for calculating the automatic threshold is as follows [32], where n represents the number of points, and x denotes the gray scale of each point. Average gray level of images is given as follows, This method ignores the structural information of the image itself, and threshold value based on global gray value is not accurate. Therefore, on the basis of the local fractal dimension, we introduce an automatic threshold method to describe complexity of texture. It is found that the local fractal dimension of image is approximate compliance with normal distribution. Normal distribution, also known as the gaussian distribution, played a vital role in fields of mathematics, physics, and engineering, etc. Here given the definition of normal distribution: as for a random variable x, which obeys a position parameter for µ, scale parameter for the probability distribution of σ and probability density function is given.
is denoted as x ∼ N (µ, σ 2 ). Gaussian model is used to quantify things with Gaussian probability density function in image processing, and has been widely used. The principle and procedure of establishing the gaussian model for texture area of fractal image are as follows: the LFD histogram reflects the estimate of the probability density of the fractal dimension of the image. If the difference between the fractal texture area and the non-texture area is relatively large, the histogram of the LFD will display a Bell-shaped curve. Hence, the whole image can be easily divided into two parts: fractal feature salient region, LFD > µ non fractal feature salient region, LFD < µ In our interpolation algorithm, ε -Blanket method is used to estimate the local fractal dimension in each block of the size 5×5. As mentioned above, µ can be treated as a threshold to distinguish textures from images. That is, the block is regarded as the texture region the value of LFD is greater than µ. Otherwise, it is regarded as the non-textured region. Fig. 2 is an original image that consists of fractal feature salient region and non fractal feature salient region. The LFD values are ranged from 2.0 to 3.0. LFD shows a greater value in fractal feature salient region. The distribution and density functions of LFD values are respectively shown in Fig. 3 and Fig. 4. Threshold images generated by fractal dimension are shown in Fig. 5.

3) SCALING FACTOR DETERMINATION
Apparently, it is meaningless to give the scaling factor in random. Hence, it is imperative that we should calculate the value of the scale factor with the given information. The more suitable the value of s i,j is, the more accurate the fractal functions are. Scaling factor is usually given a range of values or as a free parameter [33], [34]. In this paper, the scaling factor is calculated by fractal dimension, since the fractal dimension   has a strong correlation with the scaling factor. Scaling factor is obtained from the following formula.
where N=2, FD and LFD are the fractal dimensions of the entire image and a block with the size of 5 × 5 respectively. Considering the effect of the scaling factor, we keep all free shape parameters on the fractal surface and change the scaling factors s i,j . It can be found from Fig. 6 that the interpolation is very sensitive to scaling factor s i,j . Furthermore, we find VOLUME 8, 2020  out that the appropriate conditions of s i,j can ensure that the interpolation function is monotonic and convex.

B. ADAPTIVE IMAGE MAGNIFICATION
In this section, we discuss how to get a high-resolution (HR) images from a low-resolution images (LR) using our proposed method. First, we divide the image into FFSR and NFFSR by calculating the fractal dimension of the image. Second, according to regional characteristics, interpolation models are adopted in different regions. FFSR and NFFSR adopt rational fractal interpolation and rational interpolation, respectively. Third, the scaling factor is determined by the fractal dimension. Finally, the quality of the interpolation image is improved by an optimization technique.
The proposed model is a hybrid interpolation model that is effective for FFSR and NFFSR interpolation by adjusting parameters. In NFFSR, the bivariate rational interpolation is used, and in FFSR, the rational fractal interpolation is used. The proposed method based on rational fractal function and rational spline function is showed in Fig. 7.

1) FRACTAL FEATURE SALIENT REGION
For fractal feature salient region, interpolation function is written as follows: where α and β are optimized shape parameters. They can be settled in Section V. Except for α and β, another parameter named scaling factor needs to be settled in fractal interpolation formula.

2) FRACTAL FEATURE NON-SALIENT REGION
For fractal feature non-salient region, interpolation function is written as follows: where

A. PARAMETERS OPTIMIZATION
Image quality assessment (IQA) [35] aims to use computational models to measure the image quality consistently with subjective evaluations. The feature-similarity (FSIM) index is proposed according to the fact that HVS understands an image from its low-level features, which brings IQA from pixel based stage to structure based stage. The FSIM index is calculated in two stages: the local similarity map is first calculated, which is then pooled into a single similarity score. There are two components, phase congruency (PC) and gradient magnitude (GM) need the measurement of feature similarity. For the PC values PC 1 (x) and PC 2 (x), the similarity measure is defined as: where T 1 represents a positive constant used to improve the stability of S PC . Similarly, the similarity measure of G 1 (x) and G 2 (x), the values of GM, is defined as: where T 2 represents a positive constant depended on the dynamic range of GM values. To make the FSIM more convenient to use, T 1 and T 2 are fixed to all databases. Define The FSIM index between f 1 and f 2 is given as follows: Substitute Eq. (9) into Eq. (15), denote as U (α, β), the local optimal parameters, α and β are obtained by minimizing U (α, β) as follows: ∂U (α, β) ∂α = 0, Obviously, the method of obtaining the optimal parameters has a high degree of complexity. Here, a numerical method for optimizing the parameters will be given. The detail of the method is shown in Algorithm 1.
For assigned interpolation function, the interpolation quality is improved by adjusting parameters. Shape parameters are optimized by FSIM as the objective function. Optimal parameters can be obtained by the experiment, which is the intersection of two parameters optimized respectively.
As shown in Fig. 8(a), FSIM varies with the parameter values. When α ∈ [4,5], β ∈ [14,16], FSIM reached the peak value. In Fig. 8(c), when α takes a fixed value, FSIM achieved peak when β ∈ [13,15]. Similarly, when β takes a fixed value, FSIM also reached maximal value as α ∈ [3,5]. The intersection of α, β are the optimal value. The constructed model is highly consistent with image features. Shape parameters contained in the model pave the way for improving the accuracy of interpolation to some degrees. The quality of interpolated images can be further improved by optimizing shape parameters.

B. PERFORMANCE ANALYSIS
In this subsection, the proposed method is compared with state-of-the-art image interpolation methods, containing directional cubic convolution interpolation (DCCI) [9], KR [8], NEDI [7], iteractive curvature based interpolation (ICBI) [10], SRCNN [15], NARM [14]. In our experiments, the six color images listed in Fig. 9 are applied to test, including Snow, River, Night, Build, Hall and Farm. For color images, since the human visual system is more sensitive to changes in brightness, we only exploit the method proposed in this paper to the luminance channel and apply the bicubic interpolation to the color channel. We down-sampled HR images to obtain the corresponding LR images, from which the original HR images are reconstructed through the presented and comparing methods. In the following subsections, we present some cropped portions of HR images reconstructed by the competing methods, and conduct qualitative and quantitative analysis of experimental results.

1) QUANTITATIVE ANALYSIS
The quantitative results of the test images (magnified by 2 times) are displayed in Table 1, with the best results in bold type. For color images, we compare the PSNR and SSIM metrics for the luminance channel. From these results, it is apparent that the proposed method performs the approximate optimal results in the comparison algorithms of the two quantitative evaluations. For the image of Night, the proposed algorithm has achieved the best PSNR and SSIM values compared with other methods in Table 1. For the Farm image,   the performance of NARM and the proposed method in terms of quantitative measures are roughly equivalent, which is consistent with qualitative analysis. The results acquired by the presented method are nearly the best among all the methods of comparison.

2) QUALITATIVE ANALYSIS
We tested natural images to verify the visual improvements of the proposed algorithm. These images were visually examined to demonstrate that the proposed method can reduce common known artifacts (containing blurring, ringing, and blocking) and successfully reconstruct image detail information. Fig. 10 gives the portions of the enlagred image Snow. The proposed algorithm has clear advantages compared with other algorithms, which shows a strong image texture preservation ability. The KR image suffers from blurred artifacts and loses shape, as shown in Fig. 10(c). It can be observed from Fig. 10(d) that the resulting image of the NEDI method has directional artifacts. In the SRCNN magnified image, fewer artifacts appear. Although the visual quality of NARM is similar to the presented algorithm, NARM is time consuming.
As shown in Fig. 11, the proposed method performs better on both edge preservation and texture preservation due to that the algorithm can appropriately interpolate different structures in LR images. The magnified KR image indicates that the algorithm is not good enough for preserving the edge structure in the River cropped image as shown in Fig. 11(c).  In the magnified images of KR, DCCI and NEDI, the blocking artifacts around the tree are evident and the texture twist heavily. The image magnified by ICBI suffers from blocking artifacts and ringing artifacts as shown in Fig. 11(e). The texture in the magnified image of SRCNN does not get well preserved as shown in Fig. 11(f). It can be observed from Fig. 11(g) that the texture regions of the tree are relatively preserved by the NARM. Another image is presented to further evaluate the proposed algorithm. The comparisons are given in Fig. 12.
As shown in Fig. 13 and Fig. 14, smooth edges can be observed by the proposed algorithm, and noise artifacts in the magnified image can be greatly reduced. However, the models used in DCCI, KR, NEDI, ICBI, SRCNN, and NARM are not effective in adapting to changing scene structures in the edge region. Further, Fig. 15 is used to verify the ability of algorithm to preserve both edge structures and details. As shown in Fig. 15(c), the edge structure is not efficiently preserved in KR. In the images magnified by DCCI and NEDI, artifacts appear and the texture twist heavily as shown in Fig. 15(b) and Fig. 15(d). Similarly, the ICBI and SRCNN magnified images suffer from texture distortion at the wall as shown in Fig. 15(e) and Fig. 15(f). The magnified NARM image has almost the same visual effects as our method, but our method has an advantage in edge preservation and spends less time than NARM. VOLUME 8, 2020

3) CLASSICAL DATABASE ANALYSIS
To further evaluate the performance of the proposed algorithm, we tested the Set 5 and Set 14 database. Table 2 present the objective quality of the test images magnified by all methods.As shown in Fig. 16 and Fig. 17, this method can  restore more details compared to other comparison methods. The images generated by DCCI, NEDI, ICBI are fuzzy and distorted. In KR and SRCNN magnified images, artifacts appear and the texture is severely distorted. The images magnified by NARM produce overly smooth edges and lose some texture details. In summary, our method is superior to other algorithms in terms of the objective quality of the magnified images. Furthermore, from the presented visual results, the proposed algorithm outperforms other methods in retaining edge structures and keeping texture structures. In view of this application, the proposed algorithm can be used to efficiently magnify LR images and applied to different applications.

VI. CONCLUSION
This work proposes an novel image magnification algorithm based on interpolation. Firstly, based on a rational fractal interpolation method, a novel method of surface modeling is proposed, which integrates rational fractal interpolation with rational spline. The model can coincide with diversity of image features due to its different expression with different values of scaling factor and shape parameters, and thus it shows good performance in describing complex geometric structure of image. Further, the quality of image can be improved by adjusting shape parameters in respect of image quality assessment; scaling factor reveals the complexity of texture; the model is stable for window selection. Next, a new adaptive rational fractal magnification (ARFM) algorithm on the basis of surface model is proposed. Basically, ARFM algorithm has advantages from fractal interpolation and rational spline, and it can preserve the edge structures and the details of the image efficiently. The experimental results illustrate that the proposed method achieves competitive performance both subjectively and objectively.