# Local and Nonlocal Regularization to Image Interpolation

Authors: Yi Zhan; Sheng Jie Li; Meng Li
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper presents an image interpolation model with local and nonlocal regularization. A nonlocal bounded variation (BV) regularizer is formulated by an exponential function including gradient. It acts as the Perona-Malik equation. Thus our nonlocal BV regularizer possesses the properties of the anisotropic diffusion equation and nonlocal functional. The local total variation (TV) regularizer dissipates image energy along the orthogonal direction to the gradient to avoid blurring image edges. The derived model efficiently reconstructs the real image, leading to a natural interpolation which reduces blurring and staircase artifacts. We present experimental results that prove the potential and efficacy of the method.


---

## Body

## 1. Introduction

Image interpolation has been an active research topic in both academia and industry due to its huge number of applications. It aims at producing a high-resolution (HR) imageu from its low-resolution (LR) counterpart u
0. The relationship between HR and LR images can be described as
(1)
u
0
=
H
u
,
where the matrix H combines both filtering and downsampling processes. Equation (1) is an ill-posed inverse problem since the null space of H is nonzero as stated in [1]. This ill-posed inverse problem is generally approached in a regularization-based framework, which was formulated as an energy functional,
(2)
E
(
u
)
=
J
d
(
u
,
u
0
)
+
λ
J
r
(
u
)
,
where λ is a regularization parameter that controls the tradeoff between the data fidelity term J
d and the regularization term J
r. The data fidelity function J
d generally is formulated in the classical least-squares sense as J
d
(
u
,
u
0
)
=
(
1
/
2
)
|
H
u
-
u
0
|
2.Many regularizersJ
r have been developed to specify (2). Generally speaking, variational regularization approaches to image processing can be divided into two categories: the local regularization and the nonlocal regularization. In the local regularization approaches, the p-Laplacian functional and variations of it have been extensively studied in the literature. It possesses the following form:
(3)
J
p
(
u
)
=
∫
Ω
×
Ω
‍
|
∇
u
|
p
d
x
d
y
.
TV-based regularization, p
=
1, as first proposed by Rudin et al. [2] does an excellent job at preserving edges in image denoising, image interpolation [1], image restoration [3–7], and image inpainting [8]. Choosing p
=
2 results in isotropic diffusion which solves the staircasing problem. Together with other constraints, it formulates the regularization for image segmentation [9]. Different values of 1
<
p
<
2 result in anisotropic diffusion between TV-based and isotropic smoothing.These local regularization methods which penalize derivatives involve only the values and derivatives ofu at the same point [10]. Moreover, the classical gradients which are used to extract the image feature direction are local operators [11]. The information contained in the gradient is limited to a point and its immediate neighbors, while the edge curve is not a local event. As pointed out in the literatures, for example, [12, 13], this sometimes causes “blocky" solutions, staircasing effects, and “false edges," which can mislead a human or computer into identifying erroneous features not present in the true image.In order to better respect edges in images, the nonlocal regularization schemes have been proposed. Kindermann et al.  [10] were the first to interpret nonlocal means and neighborhood filters as regularization based on nonlocal functionals which possesses the following general form:
(4)
J
nl
(
u
)
=
∫
Ω
×
Ω
‍
g
(
|
u
(
x
)
-
u
(
y
)
|
2
h
2
)
w
(
x
-
y
)
d
x
d
y
with an appropriate positive weight function w, a differentiable filter function g
:
ℝ
+
→
R, and a parameter h. Choosing g
(
s
)
=
(
h
p
/
2
p
)
s
p
/
2, the nonlocal p-Laplacian functional is obtained as follows:
(5)
J
nl
p
(
u
)
=
1
2
p
∫
Ω
×
Ω
‍
|
u
(
x
)
-
u
(
y
)
|
p
w
(
x
-
y
)
d
x
d
y
,
which shares many properties of the classical p-Laplacian functional (3). The nonlocal p-Laplacian regularization of (5) has been used recently in the study of saliency detection [14] and image denoising and segmentation [15, 16].Another type of nonlocal regularization is based on thegraph gradient operator. The graph gradient operator is defined by [17–19]
(6)
∀
x
,
∇
x
w
u
=
(
w
(
x
,
y
)
(
u
(
y
)
-
u
(
x
)
)
)
y
∈
ℝ
n
,
|
∇
w
u
(
x
)
|
=
∫
Ω
‍
(
u
(
y
)
-
u
(
x
)
)
2
w
(
x
,
y
)
d
y
.
Thus, this type of nonlocal regularization is defined by
(7)
J
nl
g
(
u
)
=
∫
Ω
‍
g
(
|
∇
w
u
|
2
)
d
x
.
Taking g
(
s
)
=
s in the functional (7), Gilboa and Osher [17] focus on the nonlocal TV regularization
(8)
J
nl
_
TV
=
∫
Ω
‍
∫
Ω
‍
(
u
(
y
)
-
u
(
x
)
)
2
w
(
x
,
y
)
d
y
d
x
and use it to detect and remove irregularities from textures and image inpainting. Elmoataz et al. [18] present nonlocal discrete p-Laplacian regularization on image and manifold by the choice g
(
s
)
=
(
1
/
p
)
s
p
/
2 in (7) which leads to the corresponding nonlocal regularization
(9)
J
nl
_
g
p
=
1
p
∫
Ω
‍
|
∇
w
u
(
x
)
|
p
d
x
.
The above nonlocal functionals do not take into account the gradient of the image u. As pointed out in [10], if the problem is ill-posed a functional involving gradients will do better.In order to make use of the advantages of local and nonlocal functional, local and nonlocal regularization are presented in this paper. It is a linear combination of local TV and nonlocal BV. The nonlocal BV function is an exponential function including gradient, which shares the properties of the classical Perona-Malik equation. This makes our nonlocal BV functional possess the property of anisotropic diffusion, which preserves image edges better. The rest of the paper is organized as follows. In Section2, the motivation of our nonlocal BV regularization was introduced. The local and nonlocal regularization appear in Section 3. In Section 4, we demonstrate the experimental results to verify the effectiveness of our method, and the last section is for conclusion.


## 2. Background

Here we briefly review some related nonlocal functional presented in [10]. Considering the functional (4), the authors chose g
(
x
)
=
1
-
e
-
x, which gives the functional
(10)
J
nl
1
(
u
)
=
∫
Ω
×
Ω
‍
(
1
-
e
-
|
u
(
x
)
-
u
(
y
)
|
2
/
h
2
)
w
(
|
x
-
y
|
)
d
x
d
y
.
The corresponding Euler-Lagrange descent flow of (10) is
(11)
u
t
(
x
)
=
∫
Ω
‍
e
|
u
(
x
)
-
u
(
y
)
|
2
/
h
2
(
u
(
x
)
-
u
(
y
)
)
w
(
|
x
-
y
|
)
d
y
.As stated in [10], the above functionals might suffice for denoising purposes, but not for the deblurring or interpolation case which is ill-posed in the space L
2, since adding a Prokhorov term will in general not regularize the problem in a reasonable space. Hence, a functional involving gradients is beneficial for ill-posed problem. So for image deblurring or image interpolation it is a good choice to use nonlocal functionals derived from the space of BV. The more general functional of the this type is
(12)
J
nl
b
v
(
u
)
=
∫
Ω
‍
g
(
|
∇
u
(
x
)
-
∇
u
(
y
)
|
2
)
w
(
x
-
y
)
d
y
.
In [10], the following nonlocal regularization term was selected (g
(
s
)
=
s
+
ɛ
2) to image denoising:
(13)
NLBV
(
u
)
:
=
∫
Ω
×
Ω
‍
(
∇
u
(
x
)
-
∇
u
(
y
)
)
2
+
ɛ
2
d
x
d
y
,
where u
∈
W
1,1
(
Ω
).Forɛ
=
0, the corresponding Euler-Lagrange descent flow of (13) is
(14)
u
t
(
x
)
=
-
∫
Ω
‍
div
⁡
(
∇
u
(
x
)
-
∇
u
(
y
)
|
∇
u
(
x
)
-
∇
u
(
y
)
|
w
(
x
-
y
)
)
d
y
.
Note that w
(
x
-
y
) does not include image information. Equation (14) indicates that (13) bears some resemblance to the TV regularization of which the corresponding Euler-Lagrange descent flow is u
t
=
div
⁡
(
∇
u
/
|
∇
u
|
). This kind of dissipation process diffuses u along the image contours. It will preserve the location and the intensity transitions of the contours to avoid blurring image edges. In order to enhance image edges it is important to consider the diffusion along the direction orthogonal to the image contours, similar to the Perona-Malik anisotropic diffusion equation [20],
(15)
u
t
=
div
⁡
(
c
(
x
,
y
,
t
)
|
∇
u
|
)
.
The general selection of c
(
x
,
y
,
t
) is c
(
x
,
y
,
t
)
=
e
-
(
|
∇
u
|
/
k
)
2 or c
(
x
,
y
,
t
)
=
1
/
(
1
+
(
|
∇
u
|
/
k
)
2
), where k is a constant. It is difficult to obtain the Euler-Lagrange equation as above from a functional, but we can derive its nonlocal form similar to (15) from a nonlocal functional motivated by (10) and (11). In the next section we consider this kind of nonlocal functional.


## 3. Local and Nonlocal Regularization to Image Interpolation

In this section we formulate a nonlocal BV functional by incorporating (10) and (12), which has the form
(16)
J
BV
(
u
)
=
∫
Ω
×
Ω
‍
(
1
-
e
-
|
∇
u
(
x
)
-
∇
u
(
y
)
|
2
/
h
2
)
w
(
x
-
y
)
d
x
d
y
.
The directional derivative of J
BV
(
u
) is given by
(17)
J
BV
′
(
u
)
v
=
lim
⁡
α
→
0
⁡
1
α
[
J
BV
(
u
+
α
v
)
-
J
BV
(
u
)
]
=
2
h
2
∫
Ω
×
Ω
‍
e
-
|
∇
u
(
x
)
-
∇
u
(
y
)
|
2
/
h
2
(
∇
u
(
x
)
-
∇
u
(
y
)
)
x
x
x
x
x
×
(
∇
v
(
x
)
-
∇
v
(
y
)
)
w
(
x
-
y
)
d
x
d
y
.
We can split the integral into a term involving ∇
v
(
x
) and ∇
v
(
y
). In the latter a change of variables (
x
;
y
)
→
(
y
;
x
) gives the same integral as that with ∇
v
(
x
); hence, we end up with
(18)
J
BV
′
(
u
)
v
=
4
h
2
∫
Ω
×
Ω
‍
e
-
|
∇
u
(
x
)
-
∇
u
(
y
)
|
2
/
h
2
(
∇
u
(
x
)
-
∇
u
(
y
)
)
x
x
x
x
x
×
w
(
x
-
y
)
d
x
∇
v
(
y
)
d
y
.
Integrated by parts, the Fréchet derivative of J
BV as a functional from L
2
(
Ω
) to ℝ is given by
(19)
J
BV
′
(
u
)
=
4
h
2
∫
Ω
‍
div
⁡
x
x
x
x
x
×
(
e
-
|
∇
u
(
x
)
-
∇
u
(
y
)
|
2
/
h
2
(
∇
u
(
x
)
-
∇
u
(
y
)
)
w
(
x
-
y
)
)
d
y
.
The corresponding Euler-Lagrange descent flow is
(20)
u
t
=
-
4
h
2
∫
Ω
‍
div
⁡
(
e
-
|
∇
u
(
x
)
-
∇
u
(
y
)
|
2
/
h
2
x
x
x
x
x
x
x
x
x
x
x
×
(
∇
u
(
x
)
-
∇
u
(
y
)
)
w
(
x
-
y
)


e
-
|
∇
u
(
x
)
-
∇
u
(
y
)
|
2
/
h
2
)
d
y
.
As can be seen, (20) is similar to the Perona-Malik equation (15) with c
(
x
,
y
,
t
)
=
e
-
(
|
∇
u
|
/
k
)
2. At the points lying in flat area or gentle slope, |
∇
u
(
x
)
-
∇
u
(
y
)
| is small and (20) acts as a Laplacian evolution equation, which is used in image denoising and so forth, while, for the points on image edges, |
∇
u
(
x
)
-
∇
u
(
y
)
| is large and (20) diffuses image energy backward like the Perona-Malik equation. This process enhances image edges, which represents our desirable requirement. To eliminate oscillatory image contour, the TV regularizer will suffice.Based on the above, our proposed energy functional for regularized image interpolation is given by(21)
E
(
u
)
=
α
∫
Ω
‍
|
∇
u
|
d
x
+
β
∫
Ω
×
Ω
‍
(
1
-
e
-
|
∇
u
(
x
)
-
∇
u
(
y
)
|
2
/
h
2
)
w
(
x
-
y
)
d
x
d
y
+
1
2
|
H
u
-
u
0
|
2
.
The gradient flow associated with the functional E
(
u
) is
(22)
u
t
(
x
,
t
)
=
α
div
⁡
(
∇
u
|
∇
u
|
)
+
β
J
(
u
)
-
H
T
(
H
u
-
u
0
)
,
u
(
x
,
0
)
=
H
T
u
0
,
where
(23)
J
(
u
)
=
∫
Ω
‍
div
⁡
×
(
e
-
|
∇
u
(
x
)
-
∇
u
(
y
)
|
2
/
h
2
(
∇
u
(
x
)
-
∇
u
(
y
)
)
w
(
x
-
y
)
)
d
y
.
The kernel w
:
Ω
→
R is assumed to be nonnegative, bounded continuous radial function, with sup
⁡
(
w
)
⊂
B
(
0
,
d
) and ∫
Ω
‍
w
(
z
)
d
z
=
1.We develop a fully discrete numerical method to approximate problem (22). We recall first the notations in the finite differences scheme used in our paper. Let τ and Δ
t be the space and time steps, respectively, and let (
x
1
i
;
x
2
j
)
=
(
i
τ
;
j
τ
) be the grid points. Let u
n
(
i
;
j
) be an approximation of the function u
(
n
Δ
t
;
x
1
i
;
x
2
j
), with n
≥
0. Equation (22) can be discretized as follows:
(24)
u
n
+
1
(
i
,
j
)
-
u
n
(
i
,
j
)
Δ
t
=
(
α
div
⁡
(
∇
u
|
∇
u
|
)
-
H
T
(
H
u
-
u
0
)
)
i
j
+
β
∑
(
k
,
l
)
∈
Ω
‍
div
⁡
(
exp
⁡
(
-
|
∇
u
(
i
,
j
)
-
∇
u
(
k
,
l
)
|
2
h
2
)
x
x
x
x
x
x
x
x
x
x
x
x
x
×

(
∇
u
(
i
,
j
)
-
∇
u
(
k
,
l
)
)
w
(
(
i
,
j
)
-
(
k
,
l
)
)


(
-
|
∇
u
(
i
,
j
)
-
∇
u
(
k
,
l
)
|
2
h
2
)
)

.
In all numerical experiments, we choose the following kernel function:
(25)
w
(
x
)
:
=
{
C
exp
⁡
(
1
|
x
|
2
-
d
2
)
,
if
x
<
d
;
0
,
if
x
≥
d
;
the constant C
>
0 is selected such that ∫
Ω
‍
w
(
x
)
d
x
=
1 and d
=
2.


## 4. Experimental Results

In this section, we tested the proposed interpolation method (24) on a variety of images. The nine natural images of them are shown in Figure 1. The image sizes are 512
×
512 (Barbara), 350
×
324 (race), 360
×
364 (fingerprint), 400
×
320 (eye), 354
×
532 (ball), and 256
×
256 (hat, parrot, flower, and boat). The images were downsampled without low-pass filtering and upsampled two times, such that the peak signal-to-noise ratio (PSNR) and the mean structural similarity (MSSIM) index [21] were measured, as shown in Tables 1 and 2. The algorithms for comparison are the NEDI [22], BSAI [23], adaptive gradient magnitude self-interpolation (AGMS) (using the default parameters) [24], the LMMSE estimation [25], and nonlocal BV [10]. The MATLAB codes of the NEDI, LMMSE, and AGMS and the C codes of the BSAI were obtained from the authors' corresponding websites, respectively. The MATLAB code of AGMS is available at http://www.datatang.com/data/44058. Some of the results are shown in Figures 2–6. The choice of the parameters in (24) is based on subjective quality of the results assessed informally by our personal preference as human viewers in terms of edge sharpness, contour crispness, no ringing in smooth regions, and no ringing near edges. We use the parameters α
=
0.2, β
=
2000, h
=
0.1, and time step Δ
t
=
0.15 for the proposed interpolation method. There is no visible improvement on subjective or objective quality of the results when the parameters are not badly changed. The iteration is terminated, when |
u
n
+
1
-
u
n
|
2
<
1
0
-
6, normally within a few decade iterations.Table 1
Comparison of different interpolation algorithms using PSNR values.
Image
NEDI
LMMSE
BSAI
AGMS
NLBV
Proposed
Hat
28.6606
28.6328
28.8784
30.0459
30.9239
31.3346
Parrot
27.6830
27.6966
28.0056
28.1065
31.3969
31.6225
Boat
25.6305
25.6840
25.8219
27.3422
28.4654
28.5086
Race
23.6013
23.6792
23.8557
24.3376
26.0261
26.1027
Eye
29.7578
29.8821
29.9564
31.0001
32.7364
32.6488
Ball
31.9757
32.1733
32.2449
34.6244
36.3935
36.5681
Flower
26.7856
26.7789
26.9631
28.6568
30.2598
30.5425
Barbara
23.8678
24.3268
24.3028
25.1123
25.8749
25.5891
Fingerprint
23.7430
24.2690
24.3663
25.5444
26.3480
26.8346Table 2
Comparison of different interpolation algorithms using MSSIM values.
Image
NEDI
LMMSE
BSAI
AGMS
NLBV
Proposed
Hat
0.8575
0.8548
0.8590
0.8807
0.8949
0.9037
Parrot
0.9018
0.9037
0.9066
0.9182
0.9378
0.9433
Boat
0.7786
0.7823
0.7862
0.8309
0.8653
0.8634
Race
0.7286
0.7367
0.7445
0.7732
0.8412
0.8344
Eye
0.8834
0.8879
0.8883
0.9226
0.9443
0.9430
Ball
0.9093
0.9147
0.9149
0.9468
0.9595
0.9630
Flower
0.8311
0.8335
0.8401
0.8741
0.9073
0.9131
Barbara
0.9046
0.9011
0.9030
0.9391
0.9811
0.9650
Fingerprint
0.8911
0.8942
0.8973
0.9364
0.9420
0.9480Figure 1
Eight test images: (left to right and top to bottom) Barbara, race, hat, parrot, ball, eye, fingerprint, flower, and boat.Figure 2
Portions of eye image. From left to right and top to bottom, they are the interpolated images using NEDI, LMMSE, BSAI, AGMS, NLBV, and the proposed method.Figure 3
Portions of parrot image. From left to right and top to bottom, they are the interpolated images using NEDI, LMMSE, BSAI, AGMS, NLBV, and the proposed method.Figure 4
Portions of race image. From left to right and top to bottom, they are the original HR image and interpolated images using NEDI, LMMSE, BSAI, AGMS, NLBV, and the proposed method.Figure 5
Portions of flower image. From left to right and top to bottom, they are the original HR image and interpolated images using NEDI, LMMSE, BSAI, AGMS, NLBV, and the proposed method.Figure 6
Portions of fingerprint image. From left to right and top to bottom, they are the original HR image and interpolated images using NEDI, LMMSE, BSAI, AGMS, NLBV, and the proposed method.Figure2 shows a portion of the eye image. It is obvious that the NEDI suffers from zigzagging around the fence in the pupil. It is observable that the fence in the eye is blurring in LMMSE and BSAI interpolation and is deformed and distorted in AGMS interpolation. The artificial effects are clear in NEDI and NLBV. The proposed method produces crisp and smooth image edges, which contributes to 2.6 dB improvement in PSNR in the eye image. In Figure 3, the interpolation artifacts around black plumage strip appear in LMMSE and NLBV result. The deformation on black plumage strip appears in AGMS result. The shape of black plumage strip is natural and clear in our result, since local TV regularizer smoothes image edges and avoids staircase effects while nonlocal BV regularizer acts as anisotropic diffusion to enhance image edges in our method. But the blurring is visible in NEDI and BSAI. The same effects arise in Figures 4 and 5. In Figure 4, the teeth in tire are clear in our result. Although the teeth in tire are the clearest in NLBV, it is a false appearance produced by zigzagging effect. In Figure 5, the NEDI and our method obtain clear pistils, but our result possesses natural petal. Figure 6 shows that our method works well for texture (fingerprint) image.We use two measures, the classic PSNR and MSSIM, to characterize the difference between the reference images and the outputs of a method. IfI and K are two m
×
n images, the PSNR is defined as PSNR
=
10
·
log
⁡
10
(
Max
⁡
I
2
/
MSE
), where MSE
=
(
1
/
m
n
)
∑
i
=
0
m
-
1
‍
∑
j
=
0
n
-
1
‍
|
I
(
i
,
j
)
-
K
(
i
,
j
)
|
2. The detailed definition of MSSIM refers to [21]. MSSIM index takes values in [
0
,
1
] and increases as the quality increases. It is calculated by the code available at http://www.cns.nyu.edu/~lcv/ssim/, using the default parameters. From the two tables, the proposed method is about 2 dB PSNR and 0.05 MSSIM higher than the other algorithms on average.


## 5. Conclusion

In this paper, a local and nonlocal image interpolation model based on local TV and nonlocal BV regularization is proposed. It combines the advantages of TV regularizer and nonlocal BV regularizer, that is, allowing discontinuities and preserving 1D image structures and the diffusion of the grey values of images along image feature direction. The direction of anisotropic diffusion is indicated by the information of image feature in a larger neighborhood. This results in minimal smoothing in the directions across the image features preserving sharp edges and maximal smoothing in the directions along the image features reducing zigzagging artifacts and oscillatory. We have shown improvement over local and nonlocal regularization on a subjective scale and in many cases with an improvement in PSNR and MSSIM.



---

*Source: 230348-2014-03-12.xml*