1. Introduction

Resolution of an image has been always an important issue in many image- and video-processing applications. Image interpolation can be used for image resolution enhancement and many interpolation techniques have been developed to increase the quality of this task. Image and video codings, such as spatial scalability and transcoding, rely on image interpolation methods. Hence, these require fast and accurate interpolation methods.

Recently, various interpolation methods are proposed to improve the performance of image resolution enhancement, such as bilinear [1] and bicubic [2] interpolation methods. Bilinear is simple to be implemented and the computation complexity is low. However, this method fails to capture the value changes near the edge and blurs the image details. Bicubic interpolation is another popular approach used for resolution enhancement. It has high accuracy but its high computation complexity results in heavy time consumption. Besides, the performance of bicubic convolution is poor in edge regions of images. In [3] , Guo et al. combined bilinear interpolation and discrete cosine transform (DCT) to improve qualities of images, but it still suffers from high computation time.

Another approach for image interpolation is wavelet-based method which can mitigate the problems of blurring and contrast abating. These algorithms predict the high-pass filtered coefficients, and the reconstructed image can be obtained by these estimated coefficients. Hence, the main problem of wavelet-based interpolation is how to estimate the high-pass coefficients correctly. In [4] , detail coefficients in the strong edge area are estimated using the regularity of wavelet transform in coarser subbands. However, it only estimates the coefficients on the region with strong edges, so that its improvement is limited. Temizel and Vlachos [5] used a method that exploits wavelet coefficient correlation in a local neighborhood sense and employs linear least-squares regression to estimate the unknown detail coefficients for image-resolution enhancement. Piao et al. [6] proposed an image-resolution enhancement approach using intersubband correlation in the wavelet domain. However, linear least-squares regression cannot obtain the optimal solutions. The methods based on the probability model usually employ the hidden Markov tree (HMT) model. The HMT scheme can effectively capture the interscale transition characteristics and intrascale statistics. This scheme reflects magnitude of coefficients, but none of the signs of coefficients is reflected. Besides, the accuracy of the prediction of the signs of coefficients is limited [7] [8] .

Kim et al. [9] used the multilayer perceptron (MLP) neural networks to train a mapping from the coarser scale to the finer scale for each specified image. The weakness of this method is the limited scaling factor. Lee et al. [10] proposed a different Haar-transform-based lifting filter to predict high-frequency subbands from the LL-band for wavelet-based interpolation scheme to enhance resolution of medical images. Guo et al. [11] proposed novel 5/3 and 9/7 wavelet-based image interpolations in lifting structures for image resolution enhancement. Demirel and Anbarjafari [12] proposed the methods by using difference between low resolution image and interpolated low-frequency subband as an intermediate stage for estimating the high-frequency subbands to achieve a sharper satellite image. Chavez-Roman et al. [13] designed the technique which is based on additional edge preservation procedure and mutual interpolation between the high frequency subband images performed via the discrete wavelet transform (DWT) and the input low resolution image. In [14] , a wavelet-domain approach based on the dual-tree complex wavelet transform (DT-CWT) and nonlocal means (NLM) is proposed for resolution enhancement of the satellite images. In [15] [16] , an image resolution enhancement technique based on interpolation of the high frequency subband images is obtained by DWT and the difference image. The edges are enhanced by using stationary wavelet transform (SWT). Then all these subbands are combined to generate a new high resolution image by using inverse DWT (IDWT). Furthermore, in [16] , generalized histogram equalization (GHE) is used to obtain high resolution and high contrast satellite image.

Invertible wavelet transforms that map integers to integers [17] are important in the applications of lossless coding [18] [19] . In this paper, we propose floating-point wavelet-based and integer wavelet-based image interpolations in lifting structures. It uses the low-band wavelet coefficients to estimate high-frequency wavelet coefficients of the original image, so that the original image can be reconstructed from the low-band wavelet coefficients by using the proposed prediction algorithms. Furthermore, the proposed method utilizes linear relations between original high-frequency coefficients and estimated high-frequency coefficients to further improve the accuracy of the predicted high-frequency coefficients. The coefficients of predicted high-frequency subbands can be formed in floating-point or in integers. And then, the inverse wavelet transform is performed for synthesis of an interpolated image.

The paper is organized as follows: The standard interpolation method and a brief overview of wavelet transform and reversible integer wavelet transforms are introduced in Section 2. Section 3 derives the proposed methods. Experiment results for the proposed methods are given in Section 4. Finally, Section 5 provides conclusions and future works.

2. Background

This section first introduces the common interpolation concept and its main drawback in Section 2.1. The detailed description of wavelet transform is in [20] . Section 2.2 simply presents the concepts of the wavelet transform. Then, the wavelet-based interpolation framework with lifting structure [21] is described in Section 2.3. Reversible integer wavelet transforms and their advantages in compression systems are introduced in the following subsection.

2.1. Image Interpolation

The definition of interpolation is to determine the parameters of a continuous image representation from a set of discrete points. The resolution enhancement process can be conceptually regarded as a two-step operation. Initially, the discrete data is interpolated into a continuous curve. Second, for the additional samples, we need to stuff those added points with values to be determined. Here, the simplest method―bilinear interpolation―is selected to demonstrate the process of interpolation. Let the original image be denoted by f and the interpolated image be $\tilde{f}$ . In the following example, the interpolation ratio is assumed to be 2. In order to simplify the process, the one-dimensional linear interpolation is separately applied in vertical and horizontal directions of an image to achieve the two-dimensional bilinear interpolation. The bilinear interpolation with an interpolation ratio 2 can be formulated as

$\tilde{f}\left(2x\right)=f\left(x\right)$ (1a)

$\tilde{f}\left(2x+1\right)=f\left(x\right)+f\left(x+1\right)/2$ (1b)

where $x=0,1,2,\cdots$ . After performing the bilinear interpolation, the image can be zoomed in or out.

Wavelet-based interpolation is to predict high-frequency subbands, LH-band, HL-band, and HH-band, from the low-frequency subband. Figure 1 shows the forward two-dimensional wavelet transform, and the right side is the backward wavelet transform. By assuming that the input is a low-pass filtered image (LL-band), the unknown high-frequency subbands can be predicted from the low-frequency subband. The easiest prediction algorithm is shown in Figure 1. is to pad zero values to high-frequency subbands as LH-band, HL-band, and HH-band. However, the image interpolation model in wavelet transform is desired to present a better prediction algorithm for effectively predicting the high-frequency subbands.

Because bilinear interpolation assumes the original data are first-derivative continuous, the result is usually blurred when it interpolates the points at edges. Nevertheless, the wavelet-based interpolation can avoid this artifact by its good approximation property.

2.2. Wavelet Transform

Wavelet transform is a valuable tool for image compression. It provides efficient time-frequency localization and multiresolution analysis; thus, it is suitable for image interpolation. Wavelet transformation decomposes data into different subbands hierarchically and each high-frequency subband can locate the regions of edges and details in the original image. Figure 2 shows single stage of a two-channel analysis and synthesis filter bank. In Figure 2, h₀ and h₁ are analysis filters, and g₀ and g₁ are synthesis filters. In the analysis step an input signal x is filtered with h₀ and h₁ and downsampled to generate the low-pass band s and the high-pass band d. In the synthesis step s and d are upsampled and filtered with g₀ and g₁. The sum of the filter outputs results in the reconstructed signal $\hat{x}$ .

Two-dimensional wavelet transform can be implemented using a one-dimensional filter on an image in each column vertically and in each row horizontally, which induces four subbands (i.e., LL-, LH-, HL-, and HH-band). The LL-band is the approximation of the original image, the LH-band represents vertical information, the HL-band represents horizontal information, and the HH-band represents diagonal information shown in Figure 3. The LH-, HL-, and HH-subbands will be estimated from the LL-band in our framework for image resolution enhancement.

2.3. Lifting Scheme

The wavelet transform uses two filters L and H to conduct the convolution operation. However, such convolution operation suffers from high computation cost and requires more memory for storage. Thus, an improved approach called lifting scheme for computing the discrete wavelet transform was developed. Any discrete wavelet transform can be computed with this scheme, and almost all these transforms have reduced computational complexity compared with the standard filtering algorithm. In this scheme a trivial wavelet transform, called lazy transform, is computed. This transform splits the input signal into even- and odd-indexed sequences.

$s_0\left[n\right]=x\left[2n\right]$ (2a)

$d_0\left[n\right]=x\left[2n+1\right]$ (2b)

Next, dual lifting and lifting steps are applied to obtain

$d_i\left[n\right]=d_{i-1}\left[n\right]-\underset{k}{{\displaystyle \sum }}{p}_i\left[k\right]s_{i-1}\left[n-k\right]$ (3a)

$s_i\left[n\right]=s_{i-1}\left[n\right]-\underset{k}{{\displaystyle \sum }}{u}_i\left[k\right]d_i\left[n-k\right]$ (3b)

Figure 4 illustrates the above process with M pairs of dual and lifting steps. The samples $s_M\left[n\right]$ become the low-pass coefficients $s\left[n\right]$ and the samples $d_M\left[n\right]$ become the high-pass coefficients $d\left[n\right]$ when scaled with a factor K:

$s\left[n\right]=\frac{s_M\left[n\right]}{K}$ (4a)

$d\left[n\right]=K{d}_M\left[n\right]$ (4b)

We can find the inverse transform by reversing the operations and flipping the signs. The inverse transform is illustrated in Figure 5.

$s_M\left[n\right]=Ks\left[n\right]$ (5a)

$d_M\left[n\right]=\frac{d\left[n\right]}{K}$ (5b)

Then undo the M lifting steps and dual lifting steps to obtain

$s_{i-1}\left[n\right]=s_i\left[n\right]+\underset{k}{{\displaystyle \sum }}{u}_i\left[k\right]d_i\left[n-k\right]$ (6a)

$d_{i-1}\left[n\right]=d_i\left[n\right]+\underset{k}{{\displaystyle \sum }}{p}_i\left[k\right]s_{i-1}\left[n-k\right]$ (6b)

Finally, the even and odd samples is retrieved

$x\left[2n\right]=s_0\left[n\right]$ (7a)

$x\left[2n+1\right]=d_0\left[n\right]$ (7b)

2.4. Integer Wavelet Transform

The wavelet transform produces floating-point coefficients in most cases. Although this allows perfect reconstruction of the original image, the use of finite-precision

arithmetic and quantization results in lossy compression. For lossless compression, integer transform are needed. Traditionally, integer wavelet transform are difficulty to construct. However, the construction becomes very simple with lifting scheme. Since we can write every wavelet transform using lifting, it was shown that an integer version of every floating-point wavelet transform can be built by use of the lifting scheme. Integer wavelet transforms can be created by rounding-off the result of each dual lifting and lifting steps before adding and subtracting.

The dual lifting and the lifting step thus becomes

$d_i\left[n\right]=d_{i-1}\left[n\right]-\lfloor {\displaystyle \underset{k}{\sum }{p}_i\left[k\right]s_{i-1}\left[n-k\right]}+\frac{1}{2}\rfloor$ (8a)

$s_i\left[n\right]=s_{i-1}\left[n\right]-\lfloor {\displaystyle \underset{k}{\sum }{u}_i\left[k\right]d_i\left[n-k\right]}+\frac{1}{2}\rfloor$ (8b)

It is invertible and the inverse can be obtained by reversing the lifting and the dual lifting steps and flipping signs:

$s_{i-1}\left[n\right]=s_i\left[n\right]+\lfloor \underset{k}{{\displaystyle \sum }}{u}_i\left[k\right]d_i\left[n-k\right]+\frac{1}{2}\rfloor$ (9a)

$d_{i-1}\left[n\right]=d_i\left[n\right]+\lfloor \underset{k}{{\displaystyle \sum }}{p}_i\left[k\right]s_{i-1}\left[n-k\right]+\frac{1}{2}\rfloor$ (9b)

This obviously results in an integer to integer transform. But the coefficients $p_i\left[k\right]$ and $u_i\left[k\right]$ are not necessarily integers. Thus computing the integer transform coefficients still requires floating-point operations. However, all floating-point operations can be avoided when rational coefficients are power of two denominators in the transform. Here, we use the following integer wavelet transforms of the form $\left(N,\tilde{N}\right)$ , where N is the number of vanishing moments of the analyzing high pass filter, while $\tilde{N}$ is the number of vanishing moments of the synthesizing high pass filter. The S is for sequence, and the P is for prediction in $S+P$ transform. $\left(2+2,2\right)$ transform is inspired by the $S+P$ transform, using one extra lifting step to build the earlier $\left(2,2\right)$ into a transform with four vanishing moments of the high pass analyzing filter. The resulting transform is different from the earlier $\left(4,2\right)$ transform and therefore is called the $\left(2+2,2\right)$ transform. The analysis low pass filter of $9/7-F$ has nine coefficients, while the analysis high pass filter has seven coefficients. Both analysis and synthesis high pass filters have four vanishing moments. 2/6 transform has two coefficients of the analysis low pass filter, while the analysis high pass filter has six coefficients. 2/6 transform is a version of the $\left(3,1\right)$ transform.

$S+P\{\begin{array}{l}{d}_1\left[n\right]=d_0\left[n\right]-s_0\left[n\right]\\ s\left[n\right]=s_0\left[n\right]+\lfloor \frac{1}{2}{d}_1\left[n\right]\rfloor \\ d\left[n\right]=d_1\left[n\right]+\lfloor \frac{1}{4}\left(s\left[n-1\right]-s\left[n\right]\right)+\frac{3}{8}\left(s\left[n\right]-s\left[n+1\right]\right)+\frac{2}{8}{d}_1\left[n+1\right]+\frac{1}{2}\rfloor \end{array}$ (10)

$\left(2+2,2\right)\{\begin{array}{l}{d}_1\left[n\right]=d_0\left[n\right]-\lfloor \frac{1}{2}\left(s_0\left[n\right]+s_0\left[n+1\right]\right)+\frac{1}{2}\rfloor \\ s\left[n\right]=s_0\left[n\right]+\lfloor \frac{1}{4}\left(d_1\left[n-1\right]+d_1\left[n\right]\right)+\frac{1}{2}\rfloor \\ d\left[n\right]=d_1\left[n\right]-\lfloor \frac{1}{16}\left(-s\left[n-1\right]+s\left[n\right]+s\left[n+1\right]-s\left[n+2\right]\right)+\frac{1}{2}\rfloor \end{array}$ (11)

$\left(2,2\right)\{\begin{array}{l}d\left[n\right]=d_0\left[n\right]-\lfloor \frac{1}{2}\left(s_0\left[n\right]+s_0\left[n+1\right]\right)+\frac{1}{2}\rfloor \\ s\left[n\right]=s_0\left[n\right]+\lfloor \frac{1}{4}\left(d\left[n-1\right]+d\left[n\right]\right)+\frac{1}{2}\rfloor \end{array}$ (12)

$\left(4,2\right)\{\begin{array}{l}d\left[n\right]=d_0\left[n\right]-\lfloor \frac{9}{16}\left(s_0\left[n\right]+s_0\left[n+1\right]\right)-\frac{1}{16}\left(s_0\left[n-1\right]+s_0\left[n+2\right]\right)+\frac{1}{2}\rfloor \\ s\left[n\right]=s_0\left[n\right]+\lfloor \frac{1}{4}\left(d\left[n-1\right]+d\left[n\right]\right)+\frac{1}{2}\rfloor \end{array}$ (13)

$\left(2,4\right)\{\begin{array}{l}d\left[n\right]=d_0\left[n\right]-\lfloor \frac{1}{2}\left(s_0\left[n\right]+s_0\left[n+1\right]\right)+\frac{1}{2}\rfloor \\ s\left[n\right]=s_0\left[n\right]+\lfloor \frac{19}{64}\left(d\left[n-1\right]+d\left[n\right]\right)-\frac{3}{64}\left(d\left[n-2\right]+d\left[n+1\right]\right)+\frac{1}{2}\rfloor \end{array}$ (14)

$\left(6,2\right)\{\begin{array}{l}d\left[n\right]=d_0\left[n\right]-\lfloor \frac{75}{128}\left(s_0\left[n\right]+s_0\left[n+1\right]\right)-\frac{25}{256}\left(s_0\left[n-1\right]+s_0\left[n+2\right]\right)+\frac{3}{256}\left(s_0\left[n-2\right]+s_0\left[n+3\right]\right)+\frac{1}{2}\rfloor \\ s\left[n\right]=s_0\left[n\right]+\lfloor \frac{1}{4}\left(d\left[n-1\right]+d\left[n\right]\right)+\frac{1}{2}\rfloor \end{array}$ (15)

$\left(4,4\right)\{\begin{array}{l}d\left[n\right]=d_0\left[n\right]-\lfloor \frac{9}{16}\left(s_0\left[n\right]+s_0\left[n+1\right]\right)-\frac{1}{16}\left(s_0\left[n-1\right]+s_0\left[n+2\right]\right)+\frac{1}{2}\rfloor \\ s\left[n\right]=s_0\left[n\right]+\lfloor \frac{9}{32}\left(d\left[n-1\right]+d\left[n\right]\right)-\frac{1}{32}\left(d\left[n-2\right]+d\left[n+1\right]\right)+\frac{1}{2}\rfloor \end{array}$ (16)

$2/6\{\begin{array}{l}{d}_1\left[n\right]=d_0\left[n\right]-s_0\left[n\right]\\ s\left[n\right]=s_0\left[n\right]+\lfloor \frac{1}{2}{d}_1\left[n\right]\rfloor \\ d\left[n\right]=d_1\left[n\right]+\lfloor \frac{1}{4}\left(s\left[n-1\right]-s\left[n+1\right]\right)+\frac{1}{2}\rfloor \end{array}$ (17)

$9/7-F\{\begin{array}{l}{d}_1\left[n\right]=d_0\left[n\right]-\lfloor \frac{203}{128}\left(s_0\left[n+1\right]+s_0\left[n\right]\right)+\frac{1}{2}\rfloor \\ s_1\left[n\right]=s_0\left[n\right]-\lfloor \frac{217}{4096}\left(d_1\left[n\right]+d_1\left[n-1\right]\right)+\frac{1}{2}\rfloor \\ d\left[n\right]=d_1\left[n\right]+\lfloor \frac{113}{128}\left(s_1\left[n+1\right]+s_1\left[n\right]\right)+\frac{1}{2}\rfloor \\ s\left[n\right]=s_1\left[n\right]+\lfloor \frac{1817}{4096}\left(d_1\left[n\right]+d_1\left[n-1\right]\right)+\frac{1}{2}\rfloor \end{array}$ (18)

3. Proposed Method

This section illustrates the proposed prediction algorithm using the low-frequency coefficients to estimate the unknown high-frequency coefficients in order to enhance resolution of images to obtain a well-reconstructed image.

3.1. Prediction

Padding zero values into high-frequency subbands is an easy way to reconstruct the image. However, it causes blurred image. The proposed prediction method can solve this problem. The proposed prediction algorithm can estimate high-frequency subbands by

$\widetilde{LH}=J_{Row}\ast I_{Col}\ast LL$ (19a)

$\widetilde{HL}=I_{Row}\ast J_{Col}\ast LL$ (19b)

$\widetilde{HH}=J_{Row}\ast J_{Col}\ast LL$ (19c)

where I is the proposed low-pass prediction filter and J is the proposed high-pass prediction filter to each row and column of an image, notation * represents a convolution operator, and LL is the input low-pass filtered image. The next subsection derives $\left\{I_{\left(4,2\right)},J_{\left(4,2\right)}\right\}$ according to $\left(4,2\right)$ filters as an example to explain the algorithm.

3.2. (4,2) Filter-Based Prediction Algorithm

The lifting of $\left(4,2\right)$ filter is shown in Figure 6, which is equal to (13).

$\begin{array}{l}{d}_i=\frac{1}{16}\cdot x_{2i-2}+0\cdot x_{2i-1}+\left(-\frac{9}{16}\right)\cdot x_{2i}\\ +1\cdot x_{2i+1}+\left(-\frac{9}{16}\right)\cdot x_{2i+2}+0\cdot x_{2i+3}+\frac{1}{16}\cdot x_{2i+4}\end{array}$ (20a)

$\begin{array}{l}{s}_i=\frac{1}{64}\cdot x_{2i-4}+0\cdot x_{2i-3}+\left(-\frac{8}{64}\right)\cdot x_{2i-2}+\frac{1}{4}\cdot x_{2i-1}+\frac{46}{64}\cdot x_{2i}\\ +\frac{1}{4}\cdot x_{2i+1}+\left(-\frac{8}{64}\right)\cdot x_{2i+2}+0\cdot x_{2i+3}+\frac{1}{64}\cdot x_{2i+4}\end{array}$ (20b)

where $i=0,1,2,\cdots ,\left(M/2-1\right)$ (M is the input data size). The boundary conditions of an input image are extended symmetrically by $x_i=x_{-i}$ , for $i<0$ and $x_i=x_{2M-i}$ , for $i>M-1$ .

$d_i$ and $s_i$ induce the corresponding low-pass filter $\left[\frac{1}{64},0,\left(-\frac{8}{64}\right),\frac{1}{4},\frac{46}{64},\frac{1}{4},\left(-\frac{8}{64}\right),0,\frac{1}{64}\right]$ and the high-pass filter $\left[\frac{1}{16},0,\left(-\frac{9}{16}\right),1,\left(-\frac{9}{16}\right),0,\frac{1}{16}\right]$ .

In the prediction problem, we need to find out the relationships between the input and the low-band images. In (20b), $s_i$ shows the input image has nine sequential pixels connected to one wavelet coefficient in the low-band image by dotted lines in Figure 7. Here, the input image is downsampled to the low-band image in the even sequence, where $s_i$ is dominated by $x_{2i}$ for $i=0,1,2,\cdots ,\left(M/2-1\right)$ . In order to produce the low-band image from the input image, split the input image into even sequence ( $x_i$ for $i=0,2,4,\cdots ,M-1$ ) for $s_i$ and odd sequence ( $x_i$ for $i=1,3,5,\cdots ,M-2$ ) for $d_i$ . The relations from the low-band image to the even and odd sequences of the input image are given by Figure 8.

${\tilde{x}}_{2i}=a_0\cdot s_{i-2}+a_1\cdot s_{i-1}+a_2\cdot s_i+a_3\cdot s_{i+1}+a_4\cdot s_{i+2}$ (21a)

and

${\tilde{x}}_{2i+1}=b_0\cdot s_{i-1}+b_1\cdot s_i+b_2\cdot s_{i+1}+b_3\cdot s_{i+2}$ (21b)

respectively, where $i=0,1,2,\cdots ,\left(M/2-1\right)$ and $a_0,a_1,a_2,a_3,a_4$ and $b_0,b_1,b_2,b_3$ have constraints as below

$a_0+a_1+a_2+a_3+a_4=1$ (22a)

and

$b_0+b_1+b_2+b_3=1$ (22b)

Then, substituting ${\tilde{x}}_{2i}$ , ${\tilde{x}}_{2i+1}$ in (21a) and (21b) to $x_{2i}$ , $x_{2i+1}$ in (20a) and (20b) yields

$\begin{array}{c}{\tilde{s}}_i=\left(\frac{1}{64}\right)\left(a_0\cdot s_{i-4}+a_1\cdot s_{i-3}+a_2\cdot s_{i-2}+a_3\cdot s_{i-1}+a_4\cdot s_i\right)\\ +\left(0\right)\left(b_0\cdot s_{i-3}+b_1\cdot s_{i-2}+b_2\cdot s_{i-1}+b_3\cdot s_i\right)\\ +\left(-\frac{8}{64}\right)\left(a_0\cdot s_{i-3}+a_1\cdot s_{i-2}+a_2\cdot s_{i-1}+a_3\cdot s_i+a_4\cdot s_{i+1}\right)\\ +\left(\frac{1}{4}\right)\left(b_0\cdot s_{i-2}+b_1\cdot s_{i-1}+b_2\cdot s_i+b_3\cdot s_{i+1}\right)\\ +\left(\frac{46}{64}\right)\left(a_0\cdot s_{i-2}+a_1\cdot s_{i-1}+a_2\cdot s_i+a_3\cdot s_{i+1}+a_4\cdot s_{i+2}\right)\end{array}$

$\begin{array}{l}+\left(\frac{1}{4}\right)\left(b_0\cdot s_{i-1}+b_1\cdot s_i+b_2\cdot s_{i+1}+b_3\cdot s_{i+2}\right)\\ +\left(-\frac{8}{64}\right)\left(a_0\cdot s_{i-1}+a_1\cdot s_i+a_2\cdot s_{i+1}+a_3\cdot s_{i+2}+a_4\cdot s_{i+3}\right)\\ +\left(0\right)\left(b_0\cdot s_i+b_1\cdot s_{i+1}+b_2\cdot s_{i+2}+b_3\cdot s_{i+3}\right)\\ +\left(\frac{1}{64}\right)\left(a_0\cdot s_i+a_1\cdot s_{i+1}+a_2\cdot s_{i+2}+a_3\cdot s_{i+3}+a_4\cdot s_{i+4}\right)\end{array}$

$\begin{array}{l}=\left(\frac{a_0}{64}\right)\cdot s_{i-4}+\left(-\frac{8a_0}{64}+\frac{a_1}{64}\right)\cdot s_{i-3}+\left(\frac{46a_0}{64}-\frac{8a_1}{64}+\frac{a_2}{64}+\frac{b_0}{4}\right)\cdot s_{i-2}\\ +\left(-\frac{8a_0}{64}+\frac{46a_1}{64}-\frac{8a_2}{64}+\frac{a_3}{64}+\frac{b_0}{4}+\frac{b_1}{4}\right)\cdot s_{i-1}\\ +\left(\frac{a_0}{64}-\frac{8a_1}{64}+\frac{46a_2}{64}-\frac{8a_3}{64}+\frac{a_4}{64}+\frac{b_1}{4}+\frac{b_2}{4}\right)\cdot s_i\\ +\left(\frac{a_1}{64}-\frac{8a_2}{64}+\frac{46a_3}{64}+\frac{b_2}{4}+\frac{b_3}{4}\right)\cdot s_{i+1}\\ +\left(\frac{a_2}{64}-\frac{8a_3}{64}+\frac{46a_4}{64}+\frac{b_3}{4}\right)\cdot s_{i+2}+\left(\frac{a_3}{64}-\frac{8a_4}{64}\right)\cdot s_{i+3}+\left(\frac{a_4}{64}\right)\cdot s_{i+4}\end{array}$ (23a)

$\begin{array}{c}{\tilde{d}}_i=\left(\frac{1}{16}\right)\left(a_0\cdot s_{i-3}+a_1\cdot s_{i-2}+a_2\cdot s_{i-1}+a_3\cdot s_i+a_4\cdot s_{i+1}\right)\\ +\left(0\right)\left(b_0\cdot s_{i-2}+b_1\cdot s_{i-1}+b_2\cdot s_i+b_3\cdot s_{i+1}\right)\\ +\left(-\frac{9}{16}\right)\left(a_0\cdot s_{i-2}+a_1\cdot s_{i-1}+a_2\cdot s_i+a_3\cdot s_{i+1}+a_4\cdot s_{i+2}\right)\\ +\left(1\right)\left(b_0\cdot s_{i-1}+b_1\cdot s_i+b_2\cdot s_{i+1}+b_3\cdot s_{i+2}\right)\\ +\left(-\frac{9}{16}\right)\left(a_0\cdot s_{i-1}+a_1\cdot s_i+a_2\cdot s_{i+1}+a_3\cdot s_{i+2}+a_4\cdot s_{i+3}\right)\end{array}$

$\begin{array}{l}+\left(0\right)\left(b_0\cdot s_i+b_1\cdot s_{i+1}+b_2\cdot s_{i+2}+b_3\cdot s_{i+3}\right)\\ +\left(\frac{1}{16}\right)\left(a_0\cdot s_i+a_1\cdot s_{i+1}+a_2\cdot s_{i+2}+a_3\cdot s_{i+3}+a_4\cdot s_{i+4}\right)\\ =\left(\frac{a_0}{16}\right)\cdot s_{i-3}+\left(-\frac{9a_0}{16}+\frac{a_1}{16}\right)\cdot s_{i-2}+\left(-\frac{9a_0}{16}-\frac{9a_1}{16}+\frac{a_2}{16}+b_0\right)\cdot s_{i-1}\\ +\left(\frac{a_0}{16}-\frac{9a_1}{16}-\frac{9a_2}{16}+\frac{a_3}{16}+b_1\right)\cdot s_i+\left(\frac{a_1}{16}-\frac{9a_2}{16}-\frac{9a_3}{16}+\frac{a_4}{16}+b_2\right)\cdot s_{i+1}\\ +\left(\frac{a_2}{16}-\frac{9a_3}{16}-\frac{9a_4}{16}+b_3\right)\cdot s_{i+2}+\left(\frac{a_3}{16}-\frac{9a_4}{16}\right)\cdot s_{i+3}+\left(\frac{a_4}{16}\right)\cdot s_{i+4}\end{array}$ (23b)

which implies the low-pass prediction filter $I_{\left(4,2\right)}$ and the high-pass prediction filter $J_{\left(4,2\right)}$ are respectively given by

$\begin{array}{c}{I}_{\left(4,2\right)}=[\left(\frac{a_0}{64}\right),\left(-\frac{8a_0}{64}+\frac{a_1}{64}\right),\left(\frac{46a_0}{64}-\frac{8a_1}{64}+\frac{a_2}{64}+\frac{b_0}{4}\right),\\ \left(-\frac{8a_0}{64}+\frac{46a_1}{64}-\frac{8a_2}{64}+\frac{a_3}{64}+\frac{b_0}{4}+\frac{b_1}{4}\right),\\ \left(\frac{a_0}{64}-\frac{8a_1}{64}+\frac{46a_2}{64}-\frac{8a_3}{64}+\frac{a_4}{64}+\frac{b_1}{4}+\frac{b_2}{4}\right),\\ \left(\frac{a_1}{64}-\frac{8a_2}{64}+\frac{46a_3}{64}+\frac{b_2}{4}+\frac{b_3}{4}\right),\\ \left(\frac{a_2}{64}-\frac{8a_3}{64}+\frac{46a_4}{64}+\frac{b_3}{4}\right),\left(\frac{a_3}{64}-\frac{8a_4}{64}\right),\left(\frac{a_4}{64}\right)]\end{array}$ (24a)

$\begin{array}{c}{J}_{\left(4,2\right)}=[\left(\frac{a_0}{16}\right),\left(-\frac{9a_0}{16}+\frac{a_1}{16}\right),\left(-\frac{9a_0}{16}-\frac{9a_1}{16}+\frac{a_2}{16}+b_0\right),\\ \left(\frac{a_0}{16}-\frac{9a_1}{16}-\frac{9a_2}{16}+\frac{a_3}{16}+b_1\right),\left(\frac{a_1}{16}-\frac{9a_2}{16}-\frac{9a_3}{16}+\frac{a_4}{16}+b_2\right),\\ \left(\frac{a_2}{16}-\frac{9a_3}{16}-\frac{9a_4}{16}+b_3\right),\left(\frac{a_3}{16}-\frac{9a_4}{16}\right),\left(\frac{a_4}{16}\right)]\end{array}$ (24b)

The prediction filters $I_{\left(4,2\right)}$ and $J_{\left(4,2\right)}$ should be symmetric and the low-pass wavelet filter of $\left(4,2\right)$ is dominated by the intermediate coefficient (46/64 in (20b)), thus they have extra constraints as following

$a_0=a_4$ (25a)

$a_1=a_3$ (25b)

$b_0=b_3$ (25c)

$b_1=b_2$ (25d)

$a_2>a_1>a_0$ (25e)

$b_1>b_0$ (25f)

By setting $a_0=0$ , $a_1=1/4$ , $a_2=1/2$ , $a_3=1/4$ , $a_4=0$ and $b_0=0$ , $b_1=1/2$ , $b_2=1/2$ , $b_3=0$ empirically as one solution for general images, one has

$I_{\left(4,2\right)}=\left[0,\frac{1}{256},\left(-\frac{3}{128}\right),\frac{63}{256},\frac{35}{64},\frac{63}{256},\left(-\frac{3}{128}\right),\frac{1}{256},0\right]$ (26a)

$J_{\left(4,2\right)}=\left[0,\frac{1}{64},\left(-\frac{7}{64}\right),\frac{3}{32},\frac{3}{32},\left(-\frac{7}{64}\right),\frac{1}{64},0\right]$ (26b)

As a result, the unknown subbands ( $\widetilde{LH}$ , $\widetilde{HL}$ , $\widetilde{HH}$ ) can be constructed by substituting $I_{\left(4,2\right)}$ to I and $J_{\left(4,2\right)}$ to J in (19a)-(19c).

If we know the actual high-frequency coefficients of the reconstructed image, then we can build the relationship between actual high-frequency coefficients and estimated high-frequency coefficients. The relationship improves prediction algorithm which is based on polynomial curve fitting. In Section 4, we denote this situation as one layer. If we do not know the actual high-frequency coefficients of the reconstructed image, we can decompose the original low-frequency image into four subbands as layer two in wavelet transform of the reconstructed image. And we use LL in layer two to predict the other three high-frequency subbands in layer two. Then, we can build the relationship between actual high-frequency coefficients and estimated high-frequency coefficients in layer two to improve prediction algorithm on polynomial curve fitting. In Section 4, we denote this situation as two layers. By denoting the exact coefficients at position $\left(m,n\right)$ in LH as $y_{LH\left(m,n\right)}$ and the estimated coefficients at position $\left(m,n\right)$ in LH as $x_{LH\left(m,n\right)}$ . Applying to all coefficients of LH, we can obtain the weights $a_{LH}$ and $b_{LH}$ by

$y_{LH\left(m,n\right)}=a_{LH}\times x_{LH\left(m,n\right)}+b_{LH}$ (27)

These weights are subsequently used to gain accuracy of the coefficients of the estimated LH subband. Then, we round-off the result to get the integer version of wavelet coefficients of LH.

${y^{\prime }}_{LH\left(m,n\right)}=\lfloor a_{LH}\times x_{LH\left(m,n\right)}+b_{LH}\rfloor$ , (28)

where ${y^{\prime }}_{LH\left(m,n\right)}$ is the improved integer coefficient at position $\left(m,n\right)$ .

A similar process is possible for HL and HH subbands. After coefficient improvement, ${y^{\prime }}_{LH}$ , ${y^{\prime }}_{HL}$ and ${y^{\prime }}_{HH}$ are carried out as above, then the high-resolution image is obtained by applying the inverse integer wavelet transform. Table 1

shows prediction filters of all integer wavelet filters used in this paper.

4. Experimental Results

In this paper, we evaluate the quality of the proposed interpolation method. The test images are shown in Figure 9. Lena has various types of image components, Baboon has many detailed components, Barbara has mainly diagonal edges, Boat has mainly horizontal edges, and Peppers has strongly vertical edges. Medical 1 and Medical 2 are x-ray images of the breast and the tooth, respectively. These $512\times 512$ images (Lena, Baboon, Barbara, Boat, Peppers) are from the USC image database [22] . Color images are transformed to gray images by Matlab function rgb2gray. All the experiments are done on a PC with AMD Phenom(tm) II × 4955 Processor 3.20 GHz with 4.00 GB ram, and implemented using MATLAB R2013a codes. To evaluate the objective quality of the reconstructed high-resolution image based on the proposed method, the peak signal noise ratio (PSNR) measurement is manipulated. The LL-band of images is totally stored, and it is used to recover to the original size of images.

Tables 2-8 show that the proposed interpolation method reconstructed for test images in each integer wavelet transform for a zooming ratio 1/2. These test images are reduced by the different integer wavelet transforms and then the reduced images are enlarged with the proposed method. One layer in the table represents that the original input images are decomposed into one layer and the LL-band in layer one is used to predict the other high frequency subbands. Then build the linear relations between predicted high frequency coefficients in layer one and actual high frequency coefficients in layer one to reconstruct the original images. Two layers in the table represent that the original input images are decomposed into two layers and the LL-band in layer two is used to predict the other high frequency subbands. Then build the linear relations between predicted high frequency coefficients in layer two and actual high frequency coefficients in layer two to reconstruct the original images. Tables 2-8 show that the linear relationship built from one layer decomposition and two layers decomposition has similar effects.

Tables 9-15 show that reconstructed images using proposed method in different wavelet transforms and different integer wavelet transforms which are compared with the original image in PSNR. Tables 9-15 present that the estimated high-frequency coefficients in floating-point are better than the estimated high-frequency coefficients in integer, but the advantages in compression for integer wavelet transform can let us tolerate the distortion in integer.

The proposed method is compared with other interpolation methods, such as, zero padding, [11] without evolutionary programming (EP), [11] with EP, DT-CWT-NLM [14] and SWT-DWT [15] . The methodologies in [11] are only for $\left(2,2\right)$ and $9/7-F$ analysis filters, so we compare the results with $\left(2,2\right)$ and $9/7-F$ analysis filters between the methodologies in [11] and the proposed algorithm. First, we show PSNR values of different prediction algorithms

for $\left(2,2\right)$ wavelet filter with zooming ratio 1/2 in Table 16. These test images are reduced by the $\left(2,2\right)$ wavelet transform in one layer and then the reduced images are enlarged with corresponding methods respectively. Then, we compare different prediction algorithms with $9/7-F$ wavelet filter for zooming ratio 1/2 in Table 17. These test images are reduced by the $9/7-F$ wavelet transform in one layer and then the reduced images are enlarged with corresponding methods respectively.

The results for $\left(2,2\right)$ wavelet filter in Table 16 show that the proposed method outperforms other methods for Lena, Baboon, Barbara, Boat and Medical 1. Because the parameters of the proposed method derived in Section 3 assume the

input images are continuous, so it is suitable to predict images which are smooth and soft, and not so well with images with a lot of edges and high frequency components. Peppers and Medical 2 have worse results because of the continuous assumption. By applying the linear relation to tune the coefficients in the proposed method, it reconstructs better results. Literature [11] with EP has better performance than the proposed method for Peppers and Medical 2. The results for $9/7-F$ wavelet filter in Table 17 show that the proposed method outperforms other methods for Lena, Baboon, Barbara and Medical 1. Boat, Peppers, Medical 2 have worse results because of the continuous assumption. But, the methodology [11] with EP needs more computational time to obtain the better results and it is used to verify the method in [11] without EP closes to optimal solution. The approach in [14] is used in satellite images; DTCWT can preserve more high-frequency components of original images. The interpolation algorithm used in this method does not estimate the high-frequency information, so it cannot obtain a better performance.

One can find which filter is the best in different condition from Table 18. Using high-order wavelet analysis filter can lead to a better performance than low-order wavelet analysis filter, because high-order wavelet analysis filter has higher vanishing moments. Consequently, the reconstructed high-order filtered signal will produce a better approximation and get a better reconstruction result.

Figure 10 and Figure 11 present the subjective image quality of the proposed

algorithm. Figure 10 and Figure 11 show the proposed method is sharper than zero padding method and show the bilinear and bicubic methods blur the images.

4. Conclusion

In this paper, wavelet-based image interpolation for high performance image resolution enhancement is proposed to predict the detail coefficients in the original image from the low-pass filtered image by observing the wavelet transform in lifting scheme. The proposed method can predict the vertical and horizontal subbands and the diagonal subband, so the proposed method is quite suitable to reconstruct the original image. Then, we utilize polynomial curve fitting to build the linear relationships between the actual high frequency subbands and the predicted high frequency subbands. These relationships improve the accuracy of prediction algorithm. Reconstructed images by the proposed method are also in the form of integers. Experimental results compare the performances of different integer wavelet transforms with the proposed method in some well-known natural images. The proposed methods have been shown better performance than bilinear, bicubic, and other wavelet-based methods.

Acknowledgements

This work was supported by the Ministry of Science and Technology of Republic of China under contacts, MOST 107-2221-E-006-203-MY2 and MOST 105-2221-E-006-102-MY3.