Proceedings of 2010 IEEE 17th International Conference on Image Processing

September 26-29, 2010, Hong Kong

EDGE-ADAPTIVE IMAGE INTERPOLATION USING CONSTRAINED LEAST SQUARES
Kazu Mishiba, Taizo Suzuki and Masaaki Ikehara
Department of Electronics and Electrical Engineering, Keio University
Yokohama, Kanagawa 223-8522, Japan
ABSTRACT
Some adaptive image interpolation methods have been proposed to create higher visual quality images than traditional
interpolation methods such as bicubic interpolation. Thesemethods, however, often suffer from high computational costs
and unnatural texture interpolation. This paper proposes anovel edge-adaptive image interpolation method using anedge-directed smoothness filter. Our approach estimates theenlarged image from the original image based on an observation model. The estimated image is constrained to have manyedge-directed smooth pixels which are measured by using
the edge-directed smoothness filter introduced in this paper.
Simulation results show that the proposal method produces
images with higher visual quality, higher PSNRs and fastercomputational times than the conventional methods.
Index Terms— image processing, adaptive image interpolation, edge-directed smoothness.
1. INTRODUCTION
Image interpolation is useful in many applications, for example, video upconversion is well-known as one of the mostcommon ones. In recent years, conversion from SDTV signals to HDTV signals is especially performed with the spreadof high definition devices. Several image interpolation methods have been developed because such upconversion requires
higher quality image with faster implementation and lower
cost. For example, there are convolution-based interpolation, polynomial interpolation and spline interpolation. The
biggest advantage of these methods is low computational
complexity. However, produced images using these methods
sometimes suffer from visual artifacts such as unsmooth andblurred edge.
To improve these artifacts, adaptive interpolation methods
have been proposed in recent years. They are classified into
explicit and implicit methods. Implicit methods tend to create
higher visual quality images compared with explicit methods.
In [1], missing HR pixels are interpolated based on HR local
covariance estimated from LR local covariance. To interpolate missing HR pixels to fit the local pixel structures, Zhang
and Wu use a 2-D piecewise autoregressive model whose parameters are estimated in a moving window in the LR image

978-1-4244-7993-1/10/$26.00 ©2010 IEEE

2837

(a)

(b)

(c)

Fig. 1. Filtering results: (a)original image, (b)Laplacian filter,(c)edge-directed smoothness filter with σ = 0.5, s = 3 andβ = 0.01.
[2]. Dong et al.[3] use the iterative back-projection techniquewhich is regularized by a non-local regularizer to minimize
the reconstruction error. Although these methods create highvisual quality images, they require high computational costand large memory to calculate. Also in some cases, textures
which differ from ones in original images are produced.
This paper proposes a new edge-adaptive interpolation
method to improve above problems. To estimate the HR
image from its LR counterpart, our proposal method solves
a constrained least squares problem based on an observation model and edge-directed smoothness measured by theedge-directed smoothness filter introduced in this paper. In
addition, to obtain the HR image with faster implementationand lower cost, the problem is solved as a linear least squares
problem.
2. EDGE-DIRECTED SMOOTHNESS FILTERIn this section, we introduce a new filter to quantify the degreeof edge-directed smoothness in a local area. When the filter isapplied in areas where are not smooth along edge directions,the filtering result has large energy. To obtain an interpolated
image with edge-directed smoothness, we constrain the enlarged image applied the filter to have low energies.
Local smoothness in an image is evaluated by measuringhigh-frequency energy in the local area. Such energy measurement methods have been researched widely. As a representative example, Laplacian filter is simple and useful operator to measure it. Fig. 1(b) shows the energies of the image
applied Laplacian filter to Fig. 1(a). As shown in Fig. 1(b),ICIP 2010

edge and texture areas have large high-frequency energies andsmooth areas have small energies.
The second spatial derivative of intensity values in edge
direction is one of the simplest measures of edge-directedsmoothness. This is calculated by the difference between apixel and the average of its neighboring pixels in edge direction. This calculation, however, is highly sensitive to noise.
Additionally, it is not work well in an area with slight edge
intensities, whose smoothness should be measured by calculating the difference between a pixel and the average of
its neighboring pixels like Laplacian filter. We introduce a
new edge-directed smoothness filter which solves the problem. The filter calculates the difference between a pixel andthe weighted average of its neighboring pixels whose weights
are determined by a 2D Gaussian function and vary with thecondition of an edge. In an area with larger edge intensities,
the weights in an edge direction become larger than in the
direction perpendicular to the edge. In an area without edge,
the weight distribution becomes isotropic like Laplacian filter.
This weight control is achieved by changing a spread parameter σ of a Gaussian function in the direction perpendicular to
the edge depending on edge intensities.
Let us denote the relation between a spatial filter c and theimage d applied c to an image h byd(u) =c(u,v) h(u + v).
(1)
v∈Ts

where Ts is a set of the spatial positions with their origin at thecenter of the filter size s, u = [ux , uy ]T and v = [vx , vy ]Tindicate spatial positions in vector notation. In the caseof s = 3, for example, Ts = {[1, 1]T , [1, 0]T , [1, −1]T ,[0, 1]T , [0, 0]T , [0, −1]T , [−1, 1]T , [−1, 0]T , [−1, −1]T }. The
edge-directed smoothness filter c(u,v) is given as−1v = [0, 0]Tc(u,v)∈Ts =(2)α(u,v) otherwisewherep(u, v)2),
2σ 2p(u, v) = G(u) R(θ(u) ) v,




1
0
cos θ
sin θ
, Rθ =
G(u) =, (3)
0 E(g(u) )− sin θ cos θν is a normalization parameter satisfyingc = 0, σ is aspread parameter, g(u) and θ(u) are an edge intensity andan edge direction at u, respectively, and E(g(u) ) is a monotonically increasing function satisfying E(0) ≥ 1. We use
E(g(u) ) = βg(u) + 1 where β(> 0) is an edge intensity parameter. Scaling matrix G(u) makes pixel weights in an edge
direction larger than ones in a direction perpendicular to theedge direction.
Here we show the calculation process of the coefficientsof the edge-directed smoothness filter. First, the edge intensities and directions in local areas are calculated. These canα(u,v) = ν exp(−

2838

(a)

(b)

(c)

(d)

Fig. 2. Zoomed images of part of the original image andthe image applied the edge-directed smoothness filter. (a)A
complex area in the original image. (b)The applied image of(a). (c)An edge-directed smooth area in the original image.
(d)The applied image of (c).
be obtained by using the Sobel operator as follows. Let hxand hy denote the results obtained by applying the Sobel operator in the horizontal and vertical directions to an image h,respectively. The edge directions θ and intensities g are given
by
 
h
θ = tan−1 − x , g = hx 2 + hy 2 .
(4)hy
Next, substitution of (4) in (3) gives the coefficients of theedge-directed smoothness filter. Fig. 1(c) shows the resultof applying the edge-directed smoothness filter to Fig. 1(a)with σ = 0.5, β = 0.01 and s = 3. The zoomed images ofthe original image and the applied image are shown in Fig.
2. The applied image has low energies in an edge-directedsmooth area and high energies in a complex area.
3. INTERPOLATION ALGORITHM
In this section, we describe our new image interpolation algorithm using the edge-directed smoothness filter introduced in
section 2. Our interpolation algorithm is based on an observation model, which formulates the relation between an HR
image and an LR image. We use the following model: an
LR image f of size N1 × N2 results from warping, blurring,subsampled on an HR image h of size L1 N1 × L2 N2 , andadded noise [4]. L1 and L2 represent subsampling factors in
the horizontal and vertical directions, respectively. From theview point of image interpolation, L1 and L2 are considered
as enlargement factors. For simplicity, we consider the casethat L1 and L2 are integers satisfying L1 = L2 = L ≥ 2.
Let f and h denote lexicographically ordered vectors of sizeN1 N2 × 1 and L2 N1 N2 × 1, respectively. Here f and h arecomposed of pixel values of the image f and h, respectively.
The relation between f and h can be expressed as
f = Wh + n(5)where W is a warping, blurring and subsampling matrix ofsize N1 N2 × L2 N1 N2 and n is a lexicographically ordered
noise vector.
In image interpolation, only subsampling and blurring
caused by the point spread function (PSF) of the sensor areconsidered since there is no need to consider a camera motion
between the HR image and the LR image. This blurring andsubsampling are accomplished by averaging a square blockof HR pixels. Here let Bfi denote a set of HR pixels which
are contained in the HR square block corresponding to theposition of the LR pixel fi . The relation between LR pixels
f (x, y) and HR pixels h(x, y) is represented by
f (u) =1 h(Lu + v).
L2


(6)

v∈TL

where v ∈ TL is a set of the spatial positions satisfying 0 ≤vx , vy < L. From (5) and (6), the interpolated image h isexpressed in matrix notation asĥ = arg min f − W h2h(7)

where entry wi,j of W satisfy wi,j = 1/L2 if hj ∈ Bfi
and otherwise wi,j = 0. However, the problem of solving(7) is an ill-posed. To make the problem well-posed, (7) isregularized by using the prior knowledge which is represented
by the edge-directed smoothness constraint. We regularize (7)asĥ = arg min f − W h2 + λCh2(8)hwhere λ is a regularization parameter and entry ci,j of C satisfy ci,j = c(u,v) if u ∈ Ts and otherwise ci,j = 0, whereh(u) equals to the position hi and h(u + v) to hj . This regularization is to obtain an interpolated image which has manyedge-directed smooth pixels. In (8), λ controls the balance
between the fidelity to the LR image intensity and the smoothness in edge directions.
The HR edge information to calculate C in (8) is derived
not from an estimated HR image h but from the LR image f .
There are two reasons to do so. One is to preserve the LR edge
structure in the HR image. Another is to solve (8) as a linear
least squares problem. The process of obtaining the HR edge
information is as follows. First, the LR horizontal edge fx andvertical edge fy are obtained by applying the Sobel operator
to the original image f . Next, the HR horizontal edge hx andthe vertical edge hy are estimated by interpolating fx and fy
to the HR image size by using bilinear interpolation method.
Last, substitution of hx and hy in (4) gives the estimated HRedge information. Using the obtained C, (8) can be rewritten
as the linear least squares problem
Ah ≈ b

(9)

where


W
T
, b = fT 0 · · · 0 .
A= √λCSolving (9) gives an interpolated image.

2839

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 3. Test images: (a)Lena, (b)Barbara, (c)Room, (d)Tank,(e)Chart, (f)Airplane.
Table 1. PSNRs[dB] and average computation times[sec].
Name
Lena
Barbara
Room
Tank
Chart
Airplanetime[sec]

Bicubic
34.07
25.40
28.40
31.18
18.28
31.26
0.012

NEDI[1]
33.63
24.67
29.04
30.48
18.05
29.64
20.28

SAI[2]
34.65
25.34
29.78
31.12
17.38
31.13
∗
2.44

NLBP[3]
31.46
24.51
26.80
29.46
19.44
28.88
∗
5.74

Proposal
35.46
25.66
30.92
31.80
19.34
32.72
2.03

*Computation times of SAI and NLBP are obtained with executable
files. The others are implemented with Matlab.

4. SIMULATION RESULTS
In this section, we compare the proposal method with bicubic method, NEDI method [1], SAI method [2] and NLBP
method [3] to validate our proposal algorithm. Fig. 3 lists
six example images in our test set. We obtained simulated
LR images by using the function imresize in MATLAB Image Processing Toolbox with a downsampling factor of 2. All
algorithms run using default parameters. In our method, we
found that σ = 0.5, β = 0.001, s = 3 and λ = 0.1 produce
sufficiently good results. We use these values in all our experimental results. To solve the linear least squares problem (9),we use the conjugate gradient method. All the experiments
were taken on an Intel Core i7 3.2GHz with 6GB of memory.
The PSNR and the average computational time for eachinterpolation algorithm are shown in Table 1. The PSNR results of the proposal method show better performance thanother interpolation methods for most test images. The proposal method has the advantage of computational speed over
the other adaptive interpolation methods. The bicubic interpolation method is faster than other methods. However, theproduced images have unsmooth and blurred edges more thanother methods as shown in Fig. 4(b). The NEDI method creates steep edges but suffers from noisy interpolation artifact
in areas which have multiple edges as shown in Fig. 4(c)and Fig. 5(c). Additionally, texture regions are smoothed
and these details are lost as shown in Fig. 6(c). The SAI
method creates the smoothest edges in all methods. However, this method produces inadequate pixel structures which
do not exist in original images as shown in Fig. 4-6. In Fig.
4(d), some unconnected stripes in the original image are connected due to over-smoothing. Loss of details can be seen in
Fig. 6(d). The NLBP method over-enhances edges as shownin Fig. 4(e) and Fig. 6(e). This leads to some artifacts and results in lower PSNRs. Only in Chart image, which has sharpedges and many flat regions, the method produces an image

(a)
(b)
(c)
(d)
(e)
(f)
Fig. 4. Comparison on Room image: (a)original image, (b)bicubic, (c)NEDI, (d)SAI, (e)NLBP, (f)proposal method.

(a)

(b)

(c)

(a)

(b)

(c)

(d)

(e)

(f)

(d)

(e)
(f)
Fig. 5. Comparison on Chart image: (a)original image,(b)bicubic, (c)NEDI, (d)SAI, (e)NLBP, (f)proposal method.
with the highest PSNR. The proposal method produces sharpedges and edge-directed smoothness pixels as shown in Fig.
4(f). In Fig. 5(f), production of inadequate interpolated pixels is improved compared with other adaptive interpolation
methods. Keeping details can be seen in Fig. 6(f).

5. CONCLUSION
We propose an adaptive image interpolation method using theedge-directed smoothness filter. Our method constrains theinterpolated image to have edge-directed smoothness and thefidelity to the original image data based on the observation
model. Thereby our proposal method produces images withhigh visual quality and improves some artifacts. Furthermore,
the proposed method also performs well on PSNRs and computational times compared with other adaptive interpolation
methods.

2840

Fig. 6. Comparison on Tank image: (a)original image,(b)bicubic, (c)NEDI, (d)SAI, (e)NLBP, (f)proposal method.
Acknowledgement This work was supported by Global
COE Program ”High-Level Global Cooperation for LeadingEdge Platform on Access Spaces (C12)”.
6. REFERENCES[1] Xin Li and Michael T. Orchard, “New edge-directed interpolation,” IEEE Transactions on Image Processing,vol. 10, pp. 1521–1527, 2001.
[2] Xiangjun Zhang and Xiaolin Wu, “Image interpolation by
adaptive 2-d autoregressive modeling and soft-decision
estimation,” Image Processing, IEEE Transactions on,
vol. 17, no. 6, pp. 887–896, June 2008.
[3] W.S. Dong, L. Zhang, G.M. Shi, and X.L. Wu, “Nonlocal
back-projection for adaptive image enlargement,” in ICIP
2009, 2009, pp. 349–352.
[4] Sung C. Park, Min K. Park, and Moon G. Kang, “Superresolution image reconstruction: a technical overview,”
Signal Processing Magazine, IEEE, vol. 20, no. 3, pp.
21–36, 2003.