![Fig. 11. Interpolation results of the image Butterfly. (a) Original image, interpolated image by (b) the cubic convolution, (c) the method in [8], (d) the method in [9], (e) the proposed LMMSE_INTR_cubic, and (f) the proposed OW_INTR_cubic.](/figures/fig-11-interpolation-results-of-the-image-butterfly-a-2ywk05fc.webp)
![Fig. 7. Interpolation results of the image Lena. (a) Original image, interpolated image by (b) the cubic convolution, (c) the method in [8], (d) the method in [9], (e) the proposed LMMSE_INTR_cubic, and (f) the proposed OW_INTR_cubic.](/figures/fig-7-interpolation-results-of-the-image-lena-a-original-axhcsta2.webp)
![Fig. 10. Interpolation results of the image Cloth. (a) Original image, interpolated image by (b) the cubic convolution, (c) the method in [8], (d) the method in [9], (e) the proposed LMMSE_INTR_cubic, and (f) the proposed OW_INTR_cubic.](/figures/fig-10-interpolation-results-of-the-image-cloth-a-original-1e7sdx1r.webp)
![Fig. 8. Interpolation results of the image Splash. (a) Original image, interpolated image by (b) the cubic convolution, (c) the method in [8], (d) the method in [9], (e) the proposed LMMSE_INTR_cubic, and (f) the proposed OW_INTR_cubic.](/figures/fig-8-interpolation-results-of-the-image-splash-a-original-206oy4oz.webp)
![Fig. 12. Interpolation results of the image Butterfly when the LR image is generated by low-pass filtering and downsampling the HR image. (a) Original image, interpolated image by (b) the cubic convolution, (c) the method in [8], (d) the method in [9], (e) the proposed LMMSE_INTR_cubic, and (f) the proposed OW_INTR_cubic.](/figures/fig-12-interpolation-results-of-the-image-butterfly-when-the-1736yput.webp)
![Fig. 9. Interpolation results of the image Peppers. (a) Original image, interpolated image by (b) the cubic convolution, (c) the method in [8], (d) the method in [9], (e) the proposed LMMSE_INTR_cubic, and (f) the proposed OW_INTR_cubic.](/figures/fig-9-interpolation-results-of-the-image-peppers-a-original-13nwue7h.webp)
![Fig. 6. Interpolation results of the image Holes. (a) Original image; interpolated image by (b) the cubic convolution; (c) the method in [8]; (d) the method in [9]; (e) the proposed LMMSE_INTR_cubic; (f) the proposed OW_INTR_cubic.](/figures/fig-6-interpolation-results-of-the-image-holes-a-original-3agsr3t6.webp)
2226 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 15, NO. 8, AUGUST 2006
An Edge-Guided Image Interpolation
Algorithm via Directional Filtering
and Data Fusion
Lei Zhang, Member, IEEE, and Xiaolin Wu, Senior Member, IEEE
Abstract—Preserving edge structures is a challenge to image
interpolation algorithms that reconstruct a high-resolution image
from a low-resolution counterpart. We propose a new edge-guided
nonlinear interpolation technique through directional filtering
and data fusion. For a pixel to be interpolated, two observation sets
are defined in two orthogonal directions, and each set produces an
estimate of the pixel value. These directional estimates, modeled
as different noisy measurements of the missing pixel are fused
by the linear minimum mean square-error estimation (LMMSE)
technique into a more robust estimate, using the statistics of the
two observation sets. We also present a simplified version of the
LMMSE-based interpolation algorithm to reduce computational
cost without sacrificing much the interpolation performance.
Experiments show that the new interpolation techniques can
preserve edge sharpness and reduce ringing artifacts.
Index Terms—Data fusion, edge preservation, image interpola-
tion, linear minimum mean square-error estimation (LMMSE).
I. INTRODUCTION
M
ANY users of digital images desire to improve the
native resolution offered by imaging hardware. Image
interpolation aims to reconstruct a higher resolution (HR)
image from the associated low-resolution (LR) capture. It has
applications in medical imaging, remote sensing and digital
photographs [3]–[5], etc. A number of image interpolation
methods have been developed [1], [2], [5], [6], [8]–[16]. While
the commonly used linear methods, such as pixel duplication,
bilinear interpolation, and bicubic convolution interpolation,
have advantages in simplicity and fast implementation [7],
they suffer from some inherent defects, including block effects,
blurred details and ringing artifacts around edges. With the
prevalence of inexpensive and relatively LR digital imaging
devices and the ever increasing computing power, interests in
and demands for high-quality image interpolation algorithms
have also increased.
The human visual systems are highly sensitive to edge struc-
tures, which convey much of the image semantics, so a key re-
Manuscript received December 31, 2004; revised October 13, 2005. X.
Wu was supported in part by the Natural Sciences and Engineering Research
Council of Canada under Grants IRCPJ 283011-01 and RGP45978-2000. The
associate editor coordinating the review of this manuscript and approving it for
publication was Prof. Vicent Caselles.
L. Zhang is with the Department of Computing, The Hong Kong Polytechnic
University, Kowloon, Hong Kong (e-mail: cslzhang@comp.polyu.edu.hk).
X. Wu is with the Department of Electrical and Computer Engi-
neering, McMaster University Hamilton, ON L8S 4K1 Canada (e-mail:
xwu@mail.ece.mcmaster.ca).
Digital Object Identifier 10.1109/TIP.2006.877407
quirement for image interpolation algorithms is to faithfully re-
construct the edges in the original scene. The traditional linear
interpolation methods [1]–[3], [5], [6] do not work very well
under the edge preserving criterion. Some nonlinear interpola-
tion techniques [8]–[15] were proposed in recent years to main-
tain edge sharpness. The interpolation scheme of Jensen and
Anastassiou [8] detects edges and fits them by some templates
to improve the visual perception of enlarged images. Li and
Orchard [9] used the covariance of the LR image to estimate
the HR image covariance, which represents the edge direction
information to some extent, and proposed a Wiener-filtering
like interpolation scheme. Since this method needs a relatively
large window to compute the covariance matrix for each missing
sample, it may introduce some artifacts in local structures due to
sample statistics change and, hence, the incorrect estimation of
covariance. The image interpolator by Carrato and Tenze [10]
first replicates the pixels and then corrects them by using some
preset 3
3 edge patterns and optimizing the parameters in the
operator. Muresan [15] detected the edge in diagonal and nondi-
agonal directions and then recovered the missing samples along
the detected direction by using one-dimensional (1-D) polyno-
mial interpolation.
Some nonlinear interpolation methods try to enlarge an
image by predicting the fine structures in the HR image from
its LR counterpart. To do so, a multiresolution representation
of the image is needed. Takahashi and Taguchi [11] represented
an image by Laplacian pyramid, and with two empirically
determined parameters, they estimated the unknown high-fre-
quency components from the LR Laplacian detail signal. In the
past two decades, wavelet transform (WT) theory [17] has been
well developed and it endows a good multiresolution frame-
work for signal representation. WT decomposes a signal into
several scales, along which the signal sharp edges have some
correlation. Carey,
et al. [12] exploited the Lipschitz property
of sharp edges in wavelet scales. They used the modulus
maxima information at coarse scales to predict the unknown
wavelet coefficients at the finest scale. Then, the HR image is
constructed by inverse WT. Muresan and Parks [14] extended
this strategy by using the entire cone influence of a sharp edge
in wavelet scale space, instead of only the modulus maxima, to
estimate the finest scale coefficients through optimal recovery
theory. The wavelet-based interpolation method by Zhu, et
al. [13] uses a parametric discrete time model to characterize
important edges. With this model, in wavelet domain the
lost edge information at the finest scale is recovered through
linear minimum mean square-error estimation (LMMSE). The
1057-7149/$20.00 © 2006 IEEE
ZHANG AND WU: EDGE-GUIDED IMAGE INTERPOLATION ALGORITHM 2227
above schemes employ, implicitly or explicitly, an isolated
sharp edge model, such as an ideal or smoothed step edge,
in the algorithm development. For real images, however, the
wavelet coefficients of a sharp edge may be interfered by the
neighboring edges. Generally, nonlinear interpolation methods
are better at edge preservation than linear methods. In [16],
Malgouyres and Guichard analyzed some linear and nonlinear
image enlargement methods theoretically and experimentally.
Compared with the discontinuities in 1-D signals, edges in
two-dimensional (2-D) images have one more property: the
direction. In the linear interpolation methods, 1-D filtering is
done alternatively in horizontal and vertical directions without
heeding the local edge structures. In the presence of a sharp
edge if a missing sample is interpolated across instead of along
the edge direction, large and visually disturbing artifacts will
be introduced. A conservative strategy to avoid most severe
artifacts is to use an isotropic 2-D filter. This, however, reduces
the edge sharpness. A more “assertive” approach is to interpo-
late in an estimated edge direction. The problem with the latter
is that the penalty to image quality is high if the estimated edge
direction is wrong, which can happen due to the difficulty in
determining the edge direction from insufficient data provided
by the LR image.
This paper proposes a new balanced approach to the problem.
A missing sample is interpolated in not one but two mutually or-
thogonal directions. The two interpolation results are treated as
two estimates of the sample and adaptively fused using the sta-
tistics of a local window. Specifically, we partition the neigh-
borhood of each missing sample into two oriented subsets in or-
thogonal directions. The hope is that the two observation sets
will exhibit different statistics, since the missing sample has
higher correlation with its neighbors in the edge direction. Each
oriented subset yields an estimate of the missing pixel. The
pixel is finally interpolated by combining the two directional
estimates in the principle of LMMSE. This process can dis-
criminate the two subsets in terms of their coherence to the
missing sample, and make the subset perpendicular to the edge
direction contribute less to the LMMSE estimate of the missing
sample. The above new approach performs significantly better
than linear interpolation methods in preserving edge sharpness
while suppressing artifacts, by adapting interpolation to local
image gradient. A drawback of the proposed interpolation ap-
proach is its relatively high computational complexity. We also
develop a simplified interpolation algorithm of greatly reduced
computation requirement but without significant degradation in
performance.
The paper is organized as follows. Section II describes our
edge-guided LMMSE-type image interpolation algorithm.
Section III presents a simplified version of the algorithm,
striving for fast, practical implementations. Section IV reports
the experimental results. Section V concludes.
II. E
DGE-GUIDED LMMSE-BASED INTERPOLATION
As in many previous papers, we assume an LR image is
directly downsampled from an associated HR image
through
, , . Re-
ferring to Fig. 1, the black dots represent the available samples
of
and the white dots represent the missing samples of .
Fig. 1. Formation of an LR image from an HR image by directly down sam-
pling. The black dots represent the LR image pixels and the white dots represent
the missing HR samples.
The interpolation problem is to estimate the missing samples in
HR image
, whose size is , from the samples in LR
image
, whose size is .
The central issue of image interpolation is how to infer and
utilize the information on the missing sample that is hidden
in the neighboring pixels. If the downsampled signal of the
LR image exceeds the Nyquist sampling limit, the convolu-
tion-based interpolation methods will suffer from the aliasing
problem in reconstructing the HR image. This is the cause
of artifacts such as ringing effects in the interpolated images
which are common to linear interpolation methods. Given
that the human visual system is highly sensitive to the edges,
especially in their spatial locations, it is crucial to suppress the
interpolation artifacts while retaining the edge sharpness and
geometry.
The edge direction is the most important information for the
interpolation process. To extract and use this information, we
partition the neighboring pixels of each missing sample into two
directional subsets that are orthogonal to each other. From each
subset, a directional interpolation is made, and then the two in-
terpolated values are fused to arrive at an LMMSE estimate of
the missing sample. We recover the HR image
in two steps.
First, those missing samples
at the center locations
surrounded by four LR samples are interpolated. Second, the
other missing samples
and
are interpolated with the help of the already recovered samples
.
A. Interpolation of Samples
Referring to Fig. 2, we can interpolate the missing HR sample
along two orthogonal directions: 45 diagonal and
135
diagonal. Denote by and the
two directional interpolation results by some linear methods,
such as bilinear interpolation, bicubic convolution, or spline in-
terpolation [1]–[5]. Consider the directional interpolation out-
puts as the noisy measurements of the missing HR sample
(2-1)
where the random noise variables
and represent the
interpolation errors in the corresponding direction.
2228 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 15, NO. 8, AUGUST 2006
Fig. 2. Interpolation of the HR samples
I
(2
n;
2
m
)
. Two estimates of
I
(2
n;
2
m
)
are made in the 45 and 135 directions as two noisy measurements of
I
(2
n;
2
m
)
.
To fuse the two directional measurements and into a
more robust estimate, we rewrite (2-1) into matrix form
(2-2)
where
and
Now, the interpolation problem is to estimate the unknown
sample
from the noisy observation . This estimation can
be optimized in minimum mean square-error sense. To obtain
the minimum mean square-error estimation (MMSE) of
,
i.e.,
, we need to know the
probability density function
. In practice, however, it
is very hard to get this prior information or
cannot be
estimated at all. Thus, in real applications, LMMSE is often
employed instead of MMSE. To implement LMMSE, only the
first and second order statistics of
and are needed, which
may be estimated adaptively.
From (2-2), the LMMSE of
can be calculated as [18]
(2-3)
where
,
is the co-variance operator, and we abbreviate as
, the variance operator. Through the LMMSE operation,
fuses the information provided by directional measurements
and .
Let
, . Through intensive ex-
periments on 129 images, including outdoor and indoor images,
portraits, MRI medical images, and SAR images, etc., we found
that
and . Thus, noise vector can be con-
sidered to be zero mean. Denote by
and the normalized
correlation coefficients of
and with
Our experiments also show that the values of and are very
small. Thus, we consider
and and, consequently, to
be nearly uncorrelated with
. With the assumption that is
zero mean and uncorrelated with
, it can be derived from (2-3)
that
(2-4)
where
and . To implement the
above LMMSE scheme for
, parameters , , and need
to be estimated for each sample
in a local window.
First, let us consider the estimation of
and . Again, re-
ferring to Fig. 2, the available LR samples around
are used to estimate the mean and variance of . De-
note by
a window that centers at and contains
the LR samples in the neighborhood of
. For esti-
mation accuracy, we should use a sufficiently large window as
long as the statistics is stationary in
. However, in a locality
of edges, the image exhibits strong transient behavior. In this
case, drawing samples from a large window will be counter-
productive. To balance the conflicting requirements of sample
size and sample consistency, we propose a Gaussian weighting
in the sample window
to account for the fact that the corre-
lation between
and its neighbors decays rapidly in
the distance between them. The further an LR sample is from
, the less it should contribute to the mean value of
. We compute as
(2-5)
where
is a 2-D
Gaussian filter with scale
. The variance of is
computed as
(2-6)
Next, we discuss the estimation of
, the co-variance ma-
trix of
. Using (2-1) and the assumption that and are
zero mean and uncorrelated with
, it can be easily derived that
(2-7)
Since
has been estimated by (2-6), we need to estimate
and in a local window to arrive at
ZHANG AND WU: EDGE-GUIDED IMAGE INTERPOLATION ALGORITHM 2229
and . For this, we associate with a set
of its neighbors in the 45
diagonal direction. Denote by
the vector that centers at
(2-8)
Set
encompasses and its neighbors, i.e., the original
samples and the directional (45
diagonal) interpolated samples.
Symmetrically, we define the sample set
for associ-
ated with interpolated results in the 135
diagonal
The estimates of and are computed as
and (2-9)
where
is a 1-D Gaussian
filter with scale
.
Now,
and can be computed by (2-8),
and finally the co-variance matrix
can be estimated as
(2-10)
where
is the normalized correlation coefficient of with
Although and are nearly uncorrelated with , they are
somewhat correlated to each other because
and have
some similarities due to the high local correlation. We found
that the values of
are between 0.4 and 0.6 for most of the test
images. The correlation between
and varies, from rela-
tively strong in smooth areas to weak in active areas. In the areas
where sharp edges appear, which is the situation of our concern
and interests, the values of
are sufficiently low, and we can as-
sume that
and are uncorrelated with each other without
materially affecting the performance of the proposed interpola-
tion algorithm in practice. In practical implementation
, the
correlation coefficient between
and , can be set as 0.5
or even 0 for most natural images. Our experiments reveal that
the interpolation results are insensitive to
. Varying from 0
to 0.6 hardly changes the PSNR value and visual quality of the
interpolated image.
If a sharp edge presents in
in or near one of the two di-
rections (the 45
diagonal or the 135 diagonals), the corre-
sponding noise variances
and will differ
significantly from each other. By the adjustment of
in (2-4),
Fig. 3. Interpolations of the missing HR samples (a)
I
(2
n
0
1
;
2
m
)
and
(b)
I
(2
n;
2
m
0
1)
. The symbols “
” represent the already recovered samples
I
(2
n;
2
m
)
. The two estimates of
I
(2
n
0
1
;
2
m
)
or
I
(2
n;
2
m
0
1)
are made
in horizontal and vertical directions.
I
(2
n
0
1
;
2
m
)
and
I
(2
n;
2
m
0
1)
are
estimated similarly to
I
(2
n;
2
m
)
.
the interpolation value
or , whichever is in the direction
perpendicular to the edge, will contribute far less to the final es-
timation result
. The presented technique removes much of
the ringing artifacts around the edges, which often appear in the
interpolated images by cubic convolution and cubic spline in-
terpolation methods.
B. Interpolation of Samples
and
After the missing HR samples are estimated, the
other missing samples
and can
be estimated similarly, but now with the aid of the just estimated
HR samples. Referring to Fig. 3(a) and (b), the LR image pixels
are represented by black dots “ ,” while the estimated
samples
by symbols “ .” The samples that are to
be estimated are represented by white dots “
.” As illustrated in
Fig. 3, the missing sample
or
can be estimated in one direction by the original pixels of the LR
image, and in the other direction by the already interpolated HR
samples. Similar to (2-2), the two directional approximations of
the missing sample are considered as the noisy measurements of
or , and then the LMMSE of the
missing sample can be computed in a similar way as described in
the previous section. Finally, the whole HR
is reconstructed
by the proposed edge-guided LMMSE interpolation technique.
III. S
IMPLIFIED LMMSE INTERPOLATION ALGORITHM
In interpolating the HR samples, the LMMSE technique of
(2-4) needs to estimate
, , , and compute the inverse
of a 2
2 matrix. This may amount to too heavy a computa-
tion burden for some applications that need high throughput.
Specifically, if we set
be the average of the four nearest
LR neighbors of
to reduce computation, then computing
needs three additions and one division and computing needs
seven additions, four multiplications, and one division. By set-
ting the size of vector
and as 5 and setting
, i.e., , in (2-10) to reduce
the computational cost, we still need 20 additions and 20 mul-
tiplications to compute
. The remaining operations in (2-4)
include nine additions, eight multiplications, and one division.
In total, the algorithm needs 39 additions, 32 multiplications,
and three divisions to compute a
with (2-4).
2230 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 15, NO. 8, AUGUST 2006
One way to reduce the computational complexity is to ju-
diciously invoke the LMMSE algorithm only for pixels where
high local activities are detected, and use a simple linear inter-
polation method in smooth regions. Since edge pixels represent
minority in the total sample population, this will lead to signif-
icant savings in computations. In addition, we present a simpli-
fied version of the LMMSE-based interpolation algorithm while
only slightly decreasing the performance.
We can see that the LMMSE estimate of HR sample
is
actually a linear combination of
, and . Referring to
(2-4) and let
, then is a 2-D
vector and we rewrite (2-4) as
(3-1)
where
and are the first and second elements of .We
empirically observed that
is close to zero, and, hence,
has a light effect on . In this view, can be simplified to
a weighted average of
and , while the weights depend
largely on the noise covariance matrix
.
Instead of computing the LMMSE estimate of
, we deter-
mine an optimal pair of weights to make
a good estimate of
. The strategy of weighted average leads to significant reduc-
tion in complexity over the exact LMMSE method. Let
(3-2)
where
. The weights and are de-
termined to minimize the mean square-error (MSE) of
:
.
Although the measurement noises of and , i.e.,
and , are correlated to some extent, their correlation is suf-
ficiently low to consider
and as being approximately
uncorrelated. This assumption holds better in the areas of edges
that are critical to the human visual system and of interests to
us. In fact, if
and are highly correlated, that is to say,
the two estimates
and are close to each other, then
varies little in and anyway. With the assumption that
and are approximately uncorrelated, we can show that
the optimal weights are
(3-3)
It is quite intuitive why the weighting scheme works. For in-
stance, for an edge in or near the 135
diagonal direction, the
variance
is higher than . From (3-3), we see
that
will be less than and, consequently, has less
influence on
than , and vice versa. To compute
and as described in Section II, however, we still need
30 additions, 24 multiplications, and two divisions. In order to
TABLE I
O
PERATIONS
NEEDED FOR THE
LMMSE ALGORITHM
AND THE
SIMPLIFIED
ALGORITHM
further simplify and speed up the computation of and ,
we use the following approximations:
(3-4)
where “
” means nearly equivalent to. With the above sim-
plification, we only need 23 additions, two multiplications, and
two divisions to obtain
and . Finally, with (3-2), we
only need 24 additions, four multiplications, and two divisions
to get
. This yields a significant computational savings com-
pared with (2-4), which needs 39 additions, 32 multiplications,
and three divisions. Table I lists the operation counts of the sim-
plified algorithm and the LMMSE algorithm.
As we will see in the next section, the simplified algorithm
only has slightly lower PSNR result than the LMMSE algo-
rithm, while its output images are virtually indistinguishable
from those of the latter.
IV. E
XPERIMENTAL RESULTS
The proposed image interpolation algorithms were imple-
mented and tested, and their performance was compared with
some existing methods. We downsampled some HR images
to get the corresponding LR images, from which the original
HR images were reconstructed by the proposed and competing
methods. Since the original HR images are known in the sim-
ulation, we can compare the interpolated results with the true
images, and measure the PSNR of those interpolated images.
The presented LMMSE-based interpolator was compared
with the bicubic convolution interpolator, the bicubic spline
interpolator, the subpixel edge detection-based interpolator
of Jensen and Anastassiou [8], and the Wiener filtering-like
interpolator of Li and Orchard [9]. To assess the sensitivity
of the proposed interpolation algorithms to different initial
directional estimates prior to fusing, they were tested when
coupled with bicubic and bilinear convolution interpolators re-
spectively. In the table and figure legends, the LMMSE method
developed in Section II is labeled LMMSE_INTR_cubic or
LMMSE_INTR_linear, depending on if bicubic or bilinear
convolution interpolator is used to obtain initial directional
estimates. Likewise, the simplified method in Section III is
labeled OW_INTR_cubic (OW stands for optimal weighting)
or OW_INTR_linear. In the experiments, we set the scale
of
2-D Gaussian filter
[referring to (2-5)] around 1 and scale
of 1-D Gaussian filter [referring to (2-9)] around 1.5.
The PSNR results of the eight algorithms for the six test im-
ages are listed in Table II. The proposed LMMSE_INTR_cubic