Image Interpolation via Gradient Correlation-Based Edge Direction Estimation

,


Introduction
Image interpolation is used to produce a high-resolution image from a low-resolution one. Image interpolation is used in many applications and fields of image processing [1][2][3]. e increase in the resolution of modern-day display screens demands artifact-free upscaling of videos and images [4]. e TV industry has and is still evolving through the introduction of advanced and high-resolution technologies. e evolution of the TV industry from analog to digital and then from HDTV to SUHDTV increases the significance of image/video interpolation. Video interpolation is needed to produce a pleasant experience of viewing a low-resolution video on high-resolution displays. A video essentially consists of a sequence of frames, and each frame is represented by an image. Upscaling a video sequence typically amounts to upscaling a series of frames. An image interpolation method should interpolate a low-resolution image in such a way that it produces no artifacts or few artifacts to provide a pleasant viewing experience.
Image interpolation is also used in medical image processing [5,6] to enhance magnetic resonance images (MRI), picture archiving and communication systems (PACS), computer-assisted surgery (CAS), digital subtraction angiography (DSA), and computer-aided diagnosis. After the acquisition of medical images, the zooming and rotation of images are required for the purpose of diagnosis and treatment, and interpolation methods are required to perform such operations.
An image contains two regions: uniform regions and edges. Interpolating an edge correctly is the most important requirement of an interpolating algorithm since human eyes are very sensitive to image edges. e production of artifacts along an edge results in an unpleasant viewing experience. Most of the interpolation algorithms that are prominent in the field [6][7][8][9][10][11][12][13][14][15][16][17][18] have focused on the interpolation of edges. ese algorithms range from simple linear algorithms to complex edge-preserving methods. Linear algorithms fail to preserve most of the edges, whereas complex algorithms preserve edges very well compared to linear algorithms. However, these complex algorithms fail to detect certain edge patterns such as thin edges and sometimes produce artifacts along edges.
Nearest neighbor, bilinear, and bicubic [6] are the most prominent simple methods, providing time-efficient interpolations at the cost of blurring, jagging, and smoothing along the edges. Paik et al. [10] used cross correlations using horizontal and vertical direction vectors that fail in restoring edges that are not correlated with predefined horizontal and vertical vectors. Allebach and Wong [11] interpolated smooth regions using bilinear interpolation and applied data correction along edges. e problem with EDI is its sensitivity to the texture that can spuriously produce artifacts along the edges. Dong et al. [12] combined a nonlocal aggressive model with a sparse representation model. e sparse representation-based method provides reasonable results, but they are always time-inefficient. Sparse-based methods are dictionary-based, and they change local properties such as image intensity and edge properties, as well as introduce artifacts along the edges. Shi and Shan [14] proposed a method that correlates neighbors of the unknown pixel with six modes and interpolates unknown pixels along the mode that has the closest correlation. e problem with this method is that it uses sixteen pixels for each mode, and the major variation of only one pixel along an appropriate mode might result in the rejection of that mode.
Some researchers have provided particularly promising contributions in the field of image interpolation. Li and Orchard [15] proposed the interpolation method new edgedirected interpolation (NEDI). NEDI calculates local covariance coefficients from a low-resolution input image and then uses these coefficients to a produce high-resolution image on the basis of geometric duality between the coefficients of both images. However, NEDI fails to preserve thin edges or high-frequency regions.
Giachetti and Asuni [17] proposed the iterative curvaturebased interpolation (ICBI) method that provides better results compared to other conventional algorithms and preserves edges well. ICBI is a two-step interpolation that uses a secondorder derivative in the first step and minimization of energy functions in the second step. Minimization of the energy function refines edges, which can result in oversmoothing of edges. Yang et al. [13] and Timofte et al. [18] worked on a dictionary learning-based interpolation algorithm. Timofte et al. proposed an anchored neighborhood regression (ANR) approach that provides a much faster dictionary learningbased interpolation compared to the method proposed by Yang et al. e problem with a dictionary-based learning approach is that information in an image is replaced using a limited size dictionary that sometimes changes image properties such as intensity and edge shapes.
In this paper, a new interpolation algorithm, named directional gradient-based edge interpolation (DGEI), is proposed. e proposed method categorizes the input image into obvious edges, transitional regions, and uniform regions. Uniform regions are interpolated using simple line averaging, whereas DGEI is applied to obvious edge pixels. DGEI predicts the direction of an edge and interpolates unknown edge pixels in the direction with a high gradient correlation. Transitional edges are interpolated in the direction of neighboring pixels. e proposed method preserves edges with no or minimum artifacts using all possible predefined edge patterns. Simulation results show that the proposed algorithm provides better results compared to bicubic, NEDI, ICBI, and ANR. e rest of this paper is organized as follows. Section 2 describes some of the conventional methods. Section 3 describes the proposed algorithm. Section 4 describes the performance analysis. Section 5 is the conclusion.

Conventional Methods
Before describing the proposed algorithm, NEDI, ICBI, and ANR are briefly discussed since they are state-of-the-art methods and provide good results.
NEDI assumes that the input low-resolution image of size H × W directly comes from high-resolution output image of size c with pixel to pixel relation of I 2i,2j � X i,j . Hence, unknown pixel can be calculated using the following equation: Generalized equation for ∝ is given as where R � [R kl ] and r → � [r k ] are high-resolution covariance. Both R kl and r k can be obtained from Figure 1. e problem with NEDI is that is assuming that lowresolution images and high-resolution images are symmetrical.
is assumption works well for nonedges and smooth edges, but it fails to interpolate when there is an abrupt change in edge, thin edges, or region with texture.
ICBI is a two-step algorithm that first interpolates image in the main diagonal or auxiliary diagonal direction based on the minimum response of second-order derivative. It then refines the initially interpolated values by minimizing three energy terms U c , U e , and U i which are edge sharpness, artifact removal, and preserving isolevel curves. e energy function that requires minimization is given as 2 Scientific Programming U c , U e , and U i are given in the following equations, respectively.
ANR [18] is a patch-based dictionary learning approach based on minimization of the following function: where N l is set that consists of K nearest neighbors of the input feature image and c. allows to alleviate the singularity of ill-posed problems and stabilize the coefficient β.
e problem with ANR is that it completely changes some properties of an image such as intensity, production of artifacts in the transitional region, and unnecessary sharpness of edges.

The Proposed Method
e objective of this work is to propose a method that provides artifact-free upscaling and preserves the natural properties of the image. Figure 2 shows the flowchart of the algorithm. e proposed algorithm initially calculates the edge map of the low-resolution image by first calculating gradient G low of the input image of size M × N and then applies hysteresis thresholding to G low to obtain a lowresolution edge map E low . E low contains both obvious and

Edge Detection.
Hysteresis thresholding-based edge detection is used for edge detection. Hysteresis thresholding [19][20][21][22][23][24] involves the calculation of two thresholds T low and T high and uses edge-map linking for better true detection and noise reduction compared to single thresholding [25,26]. e method proposed by Medina-Carnicer et al. [21] was used for the calculation of T low and T high as it provides reasonable pair of low and high thresholds. e input of the Medina-Carnicer method is the gradient image G low obtained by applying horizontal and vertical Sobel operators [27].
After calculation of T L and T H , low-and high-thresholded edge maps E L and E H are obtained. E L contains all the obvious edges along with noisy edges, whereas E H contains incomplete true edges. For each incomplete edge in E H , the corresponding connected edge pixels from E L were obtained and added to the final low-resolution edge map E low . In that way, all the true edges were taken from E L while ignoring noisy edges, as E H does not have the representation of noisy edges.

Isolating Transitions and Obvious Edges.
Edges of the edge map E low are thick, as transition regions also have some reasonable gradient values and are detected in the low threshold edge map. To separate transition edges from obvious edges, nonmaximal suppression is applied to gradients at detected edge locations. Nonmaximal suppression is a technique used by Canny [22]. is method is used to thin edges by selecting the edge that corresponds to a maximum gradient value. Canny applied nonmaximal suppression to the gradient image before thresholding. e proposed method, however, applies nonmaximal suppression only after the application of edge detection. Nonmaximal suppression is applied to the image G low × E low to obtain an edge map E edge low that has edges with a thickness of only one pixel. ese edges are considered as obvious edges because they have a maximum gradient value in a thick detected edge. e transition map that corresponds to the transition regions can be calculated using Figure 3 shows the results of the isolation of the edge map and transition map. It can be seen from Figure 3 that the edge map obtained using nonmaximal suppression has edges of a thickness of one pixel. e direction obtained using E edge low contains a small chance of error because the direction of an edge can be predicted more precisely if the thickness of the edge is one pixel. E trans low , on the other hand, contains thick transition regions, and the probability of producing an error is high if the direction is obtained separately for transition regions. However, it is clear from Figure 3 that transitional edges in E trans low have the same direction as those of E edge low . erefore, if an accurate direction for an edge pixel of E edge low is predicted, an interpolation in the same direction of its neighboring transition pixels in E trans low can be performed to produce artifact-free image upscaling.

Interpolating Edge/Transition Maps. Both E trans
low and E edge low are subjected to simple edge map interpolation. e proposed algorithm upscales an edge/transition map by first interpolating along the row direction and then along the column direction. Both row and column interpolation follow the same procedure. e edge/transition map of size M × N is used as the input edge/transition map for performing edge/transition interpolation along the row direction in order to obtain the row-interpolated edge/transition map E Row of size 2M × N. e transpose of E Row is then used as the input image for column interpolation to obtain an interpolated image of size 2N × 2M. e transpose of interpolated 2N × 2M image is needed to obtain a final interpolated image E edge high (or E trans high in the case of transition map interpolation) of size 2M × 2N.  A 3 × 3 window is moved across every unknown row. e pixel that has edge/transition pixels in any of the three directions (main diagonal, vertical, and auxiliary diagonal) shown in Figure 4 is selected as an edge/transition pixel. Figure 5 shows the results of the application of edge interpolation on E edge low .

Gradient-Based Direction Estimation.
All the unknown pixels that are not detected as edge pixels or transitional regions in E edge high and E trans high are interpolated using simple line averaging. Information for known edge pixels in E edge high is used to estimate the direction of unknown edge pixels. Unknown edge pixels are then interpolated along the estimated direction. Unknown transitional pixels in E trans high are interpolated in the estimated direction of their neighboring edge pixel in E edge high . is section provides the direction calculation methodology for the unknown edges pixels in E edge high . Figure 6 shows all possible directions that an edge might have.
ese directions include horizontal, vertical, main diagonal, auxiliary diagonal, left concave, right concave, upward concave, and downward concave. e direction for every unknown edge pixel is determined by using eight available directions.
Every unknown pixel has some cases that show the positions of unknown and known pixels around it. For evaluating the direction, only the properties of the known pixels can be used. Given that known pixels belong to an odd row and odd column, three cases show the positions of known and unknown pixels.
(1) Unknown pixel is at an even row and even column (2) Unknown pixel is at an odd row and even column (3) Unknown pixel is at an even row odd column Figure 7 shows the three cases for the positions for known and unknown pixels. e red circle in each case shows the center unknown pixel, and blue circles show the positions of the known pixels, whereas black pixels show the positions of the surrounding unknown pixels. e case of even rows and even columns does not have known pixels along with adjacent vertical and horizontal directions, and the case of the odd row and even columns does not have known pixels along the vertical direction, whereas the case of even rows and odd columns do not have known pixels along the horizontal direction. erefore, the number of possible directions for the first case is six, whereas in the second and third cases, the number of possible directions is seven. e proposed direction estimation-based method uses edge map E edge high to see if a certain direction is a candidate direction or not and then uses gradient correlation to determine the pixels that can be used for interpolation. A certain direction cannot be selected as a candidate direction if it does not follow the crossing condition, for example, if one wants to see whether or not the main diagonal is a candidate for an unknown edge pixel that is on an even row and even column. e pixels C31 and C51 should be edge pixels in order to satisfy the crossing condition. If any of these two pixels are not edge pixels, then the edge is certainly not moving along the main diagonal. Also, for gradient correlation, a reference gradient value is needed because gradient correlation takes the absolute difference between the reference gradient with respective gradients of possible pixels in that direction. ese differences are finally compared with a threshold value to see whether the correlation is in the wanted range or not. Every direction has its own possible reference gradient position, its own "must be edges" condition to satisfy the crossing condition, and its own positions of possible directional members. Table 1 shows all the information needed to estimate the direction for an unknown pixel that has a position of an even row and even column. Steps for estimating the direction of the unknown edge pixel are (1) Take a 9 × 9 patch around the unknown edge pixel.
(2) For all directions, if possible references contain an edge pixel, move to step 3 and select the gradient value of that possible reference as a reference gradient. If none of the possible references for a direction has an edge, ignore that direction as the candidate direction. (3) If the direction follows the "must be edges" condition given in Table 1, select this direction as the candidate direction and move to step 4. If the "must be edges"  condition is not justified, ignore that direction as the candidate direction. (4) Take the absolute difference of the gradient of "possible directional members" with the reference gradient to calculate the gradient correlation. Select all those pixels as interpolating pixels that have a gradient correlation less than T � 21. (5) Among all candidate directions, select the direction that has more pixels with a gradient correlation value of less than T. Store their gradient correlation in an array D G . D G will then be used in the coefficient calculation for interpolation.
Directions estimated for unknown edge pixels that have positions of odd rows and even columns or even rows and odd columns can be calculated in a similar way using Tables 2 and 3, respectively.

Interpolation Using the Estimated
Direction. An edge pixel can be interpolated using where I GC r , C GC G(r) , and C GC D(r) are the image intensity, the gradient coefficient, and the distance coefficient of the rth pixel with a gradient correlation less than T. C G and C D can be calculated using Table 1: Information for evaluating edge directions for the case of even rows and even columns.
ere is no need to dynamically calculate C D for every pixel. is is because each pixel in a 9 × 9 patch is located at a fixed distance from the center pixel. erefore, a 9 × 9 matrix of coefficients C D can be calculated only once and used for all pixels. C G , however, requires dynamic calculation for every pixel because the gradient correlation is different for all pixels. Before applying C D and C G for interpolation, both coefficients should be normalized in such a way that it follows the given condition: After interpolating an edge pixel, its neighboring region needs to be scanned for the interpolation of transitional pixels. If the edge is moving in the vertical direction with zero slope or some slope, all of the transitional pixels that are connected to the left and right are scanned and interpolated in the same direction as the edge pixel. Unknown transitional pixels on the left or right direction of the interpolated edge pixel are checked for whether they are on an even row and even column, even row and odd column, or odd row and even column and are interpolated in that direction using information available in the corresponding Tables 1-3. A transitional pixel is interpolated by selecting a reference gradient and then following steps 4 and 5 of the previous subsection by calculating the gradient correlation D G , using the gradient of the selected transitional pixel as the reference. After calculating D G and indices for members of the array D G , interpolation can be applied to the transitional pixel using equations (10)- (14). However, if an edge is moving in the horizontal direction, transitional pixels connected to its upward and downward directions are interpolated similarly. A few exceptional cases might occur when the direction estimated for an edge pixel is vertical or horizontal; however, some neighboring transitional pixels are located, where they do not have known vertical or horizontal pixels. In that case, the transitional pixel is interpolated by using the mean of the two closest pixels in that direction using line averaging.

Subjective Evaluation.
e proposed method is compared with the bicubic, NEDI, ICBI, and ANR approaches by applying algorithms on grayscale images provided in [28]. Both subjective and objective evaluations of all algorithms are performed in order to quantitatively compare the quality of the images created by different methods and the related computational costs. e tic toc function in Matlab 2018b was used to calculate the execution time. Bicubic interpolation was applied using the MATLAB built-in function named interp2, whereas the source codes for NEDI, ICBI, and ANR were kindly provided by Prof. Xin Li, Prof. Andrea Giacheti, and Dr. Yang, respectively. Figure 8 shows some cropped regions of some grayscale test images that are upscaled by a factor of two. e input image was first downscaled by a factor of two, and then all algorithms were applied to the downscaled images to see whether or not they preserve distorted edges. It can be seen that the bicubic and NEDI algorithms failed to preserve edges. ANR provided some sharp edges, but changed the edge structures and intensity of those patches. ICBI preserved edges well; however, it oversmoothed the edges with the postprocessing step of refinement, as seen in Figure 8(j). e proposed method preserved edges well compared to other algorithms and also effectively preserved the contrast and shapes of the edges. e error index (EI) defines an algorithm's cumulative magnitude of error for a patch and is calculated for all patches of size M × N, given as Error indices calculated for patches in Figure 9 are given in Table 4. Error indices for the ICBI method and the proposed method are the lowest with the proposed method showing superior results. ANR and bicubic have the highest EI, whereas NEDI provides EI that is better than ANR and Bicubic.

Objective Evaluation.
e proposed algorithm outperformed conventional algorithms based on not only subjective evaluations but also objective evaluations. Table 5 shows the peak signal to noise ratio (PSNR) comparison of all algorithms for 20 images. On average, the proposed algorithm provides a PSNR that is 2.5181 dB higher than bicubic, 0.959 dB higher than NEDI, 0.594 dB higher than ICBI, and 9.13 dB higher than ANR.   Table 6 shows the CPU elapsed time for all algorithms. CPU elapsed time for bicubic is the lowest compared to other methods because it is a very basic interpolation method and the function interp2 interpolates an image using the mex compiler of MATLAB. On average, the proposed method is almost 9 times faster than NEDI, 2.61 times faster than ICBI, and 1.2 times faster than ANR.

Conclusions
e preservation of edges is considered to be the most important parameter in the field of image interpolation. Linear interpolation methods produce blurring effects while complex algorithms result in oversmoothed edges or add artifacts to the image. ese artifacts are a direct consequence of postprocessing conducted on the interpolated images to either enhance or refine edges. In this paper, a gradient-based edgepreserving interpolation method is proposed. e proposed method uses a combination of low-complexity and efficient hysteresis thresholding based edge detection, the isolation of edges and transitional regions, simple edge map interpolation, direction estimation using the edge map and gradient map, followed by gradient and distance coefficient based interpolation. Simulation results show that the proposed method provides better performance with regard to both subjective and objective evaluations, while accurately preserving edges and image details when compared to leading alternative methods.

Abbreviations
ICBI: Iterative curvature-based interpolation ANR: Anchored neighborhood regression NEDI: New edge-directed interpolation EI: Error index PSNR: Peak signal to noise ratio.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

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