# A Variational Bayesian Superresolution Approach Using Adaptive Image Prior Model

Authors: Shengrong Zhao; Renchao Jin; Xiangyang Xu; Enmin Song; Chih-Cheng Hung
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The objective of superresolution is to reconstruct a high-resolution image by using the information of a set of low-resolution images. Recently, the variational Bayesian superresolution approach has been widely used. However, these methods cannot preserve edges well while removing noises. For this reason, we propose a new image prior model and establish a Bayesian superresolution reconstruction algorithm. In the proposed prior model, the degree of interaction between pixels is adjusted adaptively by an adaptive norm, which is derived based on the local image features. Moreover, in this paper, a monotonically decreasing function is used to calculate and update the single parameter, which is used to control the severity of penalizing image gradients in the proposed prior model. Thus, the proposed prior model is adaptive to the local image features thoroughly. With the proposed prior model, the edge details are preserved and noises are reduced simultaneously. A variational Bayesian inference is employed in this paper, and the formulas for calculating all the variables including the HR image, motion parameters, and hyperparameters are derived. These variables are refined progressively in an iterative manner. Experimental results show that the proposed SR approach is very efficient when compared to existing approaches.


---

## Body

## 1. Introduction

Superresolution (SR) technique [1–5] is an important branch in image fusion technology which targets the reconstruction of a high-resolution (HR) image from a set of degraded low-resolution (LR) images. The LR images are usually affected by warping, blurring, downsampling, and noising. Due to advancement of the technology, superresolution has been widely used in a broad range of applications such as computer vision, medical imaging, public safety, and military reconnaissance.SR is a typically ill-posed inverse problem that is hard to solve without the introduction of some prior image information [6–8]. Thus, a number of regularization-based SR approaches have been proposed [9–11] by incorporating the prior knowledge of the unknown high-resolution image in the regularization strategy [12]. The authors in [9] proposed a regularization-based SR approach based on the Tikhonov regularizer using L
2 norm. It can remove noises effectively; however, it also blurs the image edges. Farsiu et al. [10] proposed the bilateral total variation (BTV) regularizer, which penalized gradient magnitudes measured by L
1 norm, with the intention of preserving edges in the reconstruction. However, the BTV regularizer often produces artifacts in the smoothed regions due to the use of L
1 norm [11]. In [11], the authors improved BTV with the norm ρ
x
,
a
=
a
a
2
+
x
2
-
a
2, where x is the gradient value and a is a positive scale parameter, and proposed the bilateral edge-preserving (BEP) regularizer. This new norm can adaptively use L
1 and L
2 norms, and the transition from L
1 to L
2 can be controlled by modifying the positive scale parameter a. Unfortunately, the BEP regularizer employed a fixed scale parameter for the whole image and ignored the local image features. With this drawback, it cannot well control the level of punishment to local gradients.It is well known that accurate registration which refers to the estimation of motion information is a very important factor in the success of the SR reconstruction method [13]. The above SR methods only registered the observations once before reconstructing the HR image. This is not a robust method. The error caused by an inaccurate registration would produce poor effects in the reconstruction process [13]. Other than the motion parameters, the hyperparameters related to the image prior model and the noise model could intimately affect the image reconstructing quality [12]. Recently, the variational Bayesian method has been used for SR [14–16]. This method approximates the joint posterior distribution of the unknowns, including the HR image, the motion parameters, and the hyperparameters, using a tractable distribution by minimizing the Kullback-Leibler (KL) divergence and by estimating the unknowns simultaneously in each iteration.The Bayesian approach uses a prior model to estimate a priori knowledge of the unknown HR image. For edge-preservation, the variational Bayesian approaches based on TV prior model [14] and l
1 prior model [16] have been proposed. In these two prior models, the gradient magnitudes are measured by the L
1 norm, and consequently both of them can preserve edges. Compared with the TV prior, the l
1 prior contains two model parameters, which constrain certain preferred edge directions [15, 16]. However, the L
1 norm may not always distinguish the real image features from the effects of errors caused by inaccurate registration and additive Gaussian noise. Therefore, some artifacts will be produced in the smoothed image regions. Villena et al. made progress by combining the TV prior model and simultaneous autoregressive (SAR) (i.e., based on L
2 norm) prior model [15]. A similar approach is also applied to the l
1 prior model. The problem is that it is difficult to choose a proper parameter to balance the priors in the combination.In this paper, a variational Bayesian estimate of the HR image, the motion parameters, and the hyperparameters are proposed and derived and then given a set of LR images. In order to overcome the difficulty of aforementioned SR methods, an adaptive image prior model is proposed, in which the degree of interaction between pixels is adjusted adaptively by using the normρ
(
x
,
a
) in order to preserve edges and remove noises. It is important to note that we propose a method for automatically estimating the scale parameters (i.e., a) to the proposed prior model. Information needed to determine these scale parameters is updated iteratively based on the available estimated HR image. Experiment results show the effectiveness of the proposed SR method.The rest of this paper is organized as follows: the observation model and a Bayesian framework for image SR are described in Section2. The adaptive image prior model is sketched in Section 3. The proposed SR approach is presented in Section 4. Experimental results are demonstrated in Section 5. Finally, the conclusion is given in Section 6. In addition, Appendix is presented at the end of the paper showing the solving process in detail.


## 2. The Problem Formulation

We formulate the problem to be solved in this section. The formulation is divided into the observation model and hierarchical Bayesian framework which are discussed below.


### 2.1. The Observation Model

Formulating an observation model that relates the original HR image to the LR observations is the first step to fully analyze the SR reconstruction problem.Letu denote the original HR image of size P
1
N
1
×
P
2
N
2, and each observed LR image is of size N
1
×
N
2. P
1 and P
2 denote the downsampling factors in the horizontal and vertical directions, respectively. The HR image and LR images are written in lexicographical notations as the vectors of sizes P
1
N
1
P
2
N
2
×
1 and N
1
N
2
×
1, respectively. In this work, the image acquisition process is modeled by geometric transformation, blurring, downsampling, and adding with white Gaussian noise. Thus, the following popular matrix notation [15] is adopted to describe the acquisition process:(1)
v
k
=
A
H
k
C
s
k
u
+
ε
k
1
≤
k
≤
L
,where v
k is one of a set (L) of LR images, A (the dimension of A is N
1
N
2
×
P
1
N
1
P
2
N
2) is the downsampling matrix, H
k (the dimension of H
k is P
1
N
1
P
2
N
2
×
P
1
N
1
P
2
N
2) is the blur matrix, C
(
s
k
) (the dimension of C
(
s
k
) is P
1
N
1
P
2
N
2
×
P
1
N
1
P
2
N
2) is the warp matrix, and ε
k (the dimension of ε
k is N
1
N
2
×
1) is the additive white Gaussian noise.The motion that occurs during the image acquisition is represented by warp matrixC
(
s
k
). In this paper, we assume that the motion includes global rotation and translation; that is, s
k
=
θ
k
,
x
k
,
y
k
T, where θ
k, x
k, and y
k are the rotation angle and the displacement in horizontal direction and vertical direction of the kth HR image, respectively, relative to the reference frame. In this study, we assume that the blur matrices H
k are determined by sensors. The sensor PSF is usually modeled as a spatial averaging operator, and the characteristics of this blur are usually assumed to be known. Under these assumptions, the matrices H
k and C
(
s
k
) have a block circulant with circulant block structure. The downsampling matrix A is determined by the downsampling factors P
1 and P
2, and it generates LR images from the warped and blurred HR image.For convenience, we defineB
(
s
k
)
=
A
H
k
C
(
s
k
) and v
=
{
v
k
}, where · is a set.


### 2.2. The Hierarchical Bayesian Framework

We use a hierarchical Bayesian framework [17] to model the acquisition process, the HR image u, the motion parameters {
s
k
}, and the hyperparameters α and {
β
k
}. Parameters α and {
β
k
} are the model parameters of our proposed prior model and the noise model, respectively. Thus, we can obtain the following joint posterior distribution of u, {
s
k
}, {
β
k
}, and α given v by using the Bayes rule:(2)
p
u
,
s
k
,
β
k
,
α
∣
v
=
p
v
∣
u
,
s
k
,
β
k
p
u
∣
α
p
s
k
p
β
k
p
α
p
v
,where p
v
∣
u
,
s
k
,
β
k represents the conditional distribution of the LR images v and p
u
∣
α, p
(
{
s
k
}
), p
(
{
β
k
}
), and p
(
α
) are prior distributions for the unknown image u, the motion parameters {
s
k
}, and the hyperparameters {
β
k
} and α, respectively.Suppose that the noiseε
k in the LR image v
k (k
=
1
,
2
,
…
,
L) is the white Gaussian noise with a zero-mean distribution ε
k
~
N
(
0
,
β
k
-
1
), and we can write(3)
p
ε
k
∣
β
k
∝
β
k
N
1
N
2
/
2
exp
⁡
-
β
k
2
ε
k
2
2
,where ·
2 denotes the L
2 norm of a given vector ·.Due toε
k
=
v
k
-
B
(
s
k
)
u, we can obtain the conditional distribution of the LR image v
k,(4)
p
v
k
∣
u
,
s
k
,
β
k
∝
β
k
N
1
N
2
/
2
exp
⁡
-
β
k
2
v
k
-
B
s
k
u
2
2
.Moreover, if we assume the statistical independence of noises among the LR images, we can obtain the conditional distribution of the LR images v
(5)
p
v
∣
u
,
s
k
,
β
k
=
∏
k
=
1
L
p
v
k
∣
u
,
s
k
,
β
k
.By using the available registration algorithm, the initial values of motion parameters can be obtained. However, the obtained initial values are often inaccurate. In this paper, we refine the values of motion parameters by modeling them as stochastic variables following Gaussian distributions(
s
k
~
N
(
s
k
0
,
δ
k
)
), which is similar to the model used in [15].The prior distributions for the hyperparameters are defined to be Gamma distributions (t
~
Γ
(
a
t
,
b
t
), t
∈
{
α
,
{
β
k
}
}). Gamma distributions for the hyperparameters were selected because they are conjugate for the variance of the Gaussian distribution; therefore, the posteriors will also have Gamma distributions in the Bayesian formulation [18].Our proposed adaptive image prior model, denoted byp
u
∣
α, will be presented in the next section.


## 2.1. The Observation Model

Formulating an observation model that relates the original HR image to the LR observations is the first step to fully analyze the SR reconstruction problem.Letu denote the original HR image of size P
1
N
1
×
P
2
N
2, and each observed LR image is of size N
1
×
N
2. P
1 and P
2 denote the downsampling factors in the horizontal and vertical directions, respectively. The HR image and LR images are written in lexicographical notations as the vectors of sizes P
1
N
1
P
2
N
2
×
1 and N
1
N
2
×
1, respectively. In this work, the image acquisition process is modeled by geometric transformation, blurring, downsampling, and adding with white Gaussian noise. Thus, the following popular matrix notation [15] is adopted to describe the acquisition process:(1)
v
k
=
A
H
k
C
s
k
u
+
ε
k
1
≤
k
≤
L
,where v
k is one of a set (L) of LR images, A (the dimension of A is N
1
N
2
×
P
1
N
1
P
2
N
2) is the downsampling matrix, H
k (the dimension of H
k is P
1
N
1
P
2
N
2
×
P
1
N
1
P
2
N
2) is the blur matrix, C
(
s
k
) (the dimension of C
(
s
k
) is P
1
N
1
P
2
N
2
×
P
1
N
1
P
2
N
2) is the warp matrix, and ε
k (the dimension of ε
k is N
1
N
2
×
1) is the additive white Gaussian noise.The motion that occurs during the image acquisition is represented by warp matrixC
(
s
k
). In this paper, we assume that the motion includes global rotation and translation; that is, s
k
=
θ
k
,
x
k
,
y
k
T, where θ
k, x
k, and y
k are the rotation angle and the displacement in horizontal direction and vertical direction of the kth HR image, respectively, relative to the reference frame. In this study, we assume that the blur matrices H
k are determined by sensors. The sensor PSF is usually modeled as a spatial averaging operator, and the characteristics of this blur are usually assumed to be known. Under these assumptions, the matrices H
k and C
(
s
k
) have a block circulant with circulant block structure. The downsampling matrix A is determined by the downsampling factors P
1 and P
2, and it generates LR images from the warped and blurred HR image.For convenience, we defineB
(
s
k
)
=
A
H
k
C
(
s
k
) and v
=
{
v
k
}, where · is a set.


## 2.2. The Hierarchical Bayesian Framework

We use a hierarchical Bayesian framework [17] to model the acquisition process, the HR image u, the motion parameters {
s
k
}, and the hyperparameters α and {
β
k
}. Parameters α and {
β
k
} are the model parameters of our proposed prior model and the noise model, respectively. Thus, we can obtain the following joint posterior distribution of u, {
s
k
}, {
β
k
}, and α given v by using the Bayes rule:(2)
p
u
,
s
k
,
β
k
,
α
∣
v
=
p
v
∣
u
,
s
k
,
β
k
p
u
∣
α
p
s
k
p
β
k
p
α
p
v
,where p
v
∣
u
,
s
k
,
β
k represents the conditional distribution of the LR images v and p
u
∣
α, p
(
{
s
k
}
), p
(
{
β
k
}
), and p
(
α
) are prior distributions for the unknown image u, the motion parameters {
s
k
}, and the hyperparameters {
β
k
} and α, respectively.Suppose that the noiseε
k in the LR image v
k (k
=
1
,
2
,
…
,
L) is the white Gaussian noise with a zero-mean distribution ε
k
~
N
(
0
,
β
k
-
1
), and we can write(3)
p
ε
k
∣
β
k
∝
β
k
N
1
N
2
/
2
exp
⁡
-
β
k
2
ε
k
2
2
,where ·
2 denotes the L
2 norm of a given vector ·.Due toε
k
=
v
k
-
B
(
s
k
)
u, we can obtain the conditional distribution of the LR image v
k,(4)
p
v
k
∣
u
,
s
k
,
β
k
∝
β
k
N
1
N
2
/
2
exp
⁡
-
β
k
2
v
k
-
B
s
k
u
2
2
.Moreover, if we assume the statistical independence of noises among the LR images, we can obtain the conditional distribution of the LR images v
(5)
p
v
∣
u
,
s
k
,
β
k
=
∏
k
=
1
L
p
v
k
∣
u
,
s
k
,
β
k
.By using the available registration algorithm, the initial values of motion parameters can be obtained. However, the obtained initial values are often inaccurate. In this paper, we refine the values of motion parameters by modeling them as stochastic variables following Gaussian distributions(
s
k
~
N
(
s
k
0
,
δ
k
)
), which is similar to the model used in [15].The prior distributions for the hyperparameters are defined to be Gamma distributions (t
~
Γ
(
a
t
,
b
t
), t
∈
{
α
,
{
β
k
}
}). Gamma distributions for the hyperparameters were selected because they are conjugate for the variance of the Gaussian distribution; therefore, the posteriors will also have Gamma distributions in the Bayesian formulation [18].Our proposed adaptive image prior model, denoted byp
u
∣
α, will be presented in the next section.


## 3. The Proposed Adaptive Image Prior Model

In this section, we propose an adaptive image prior model and then present a method to calculate the scale parametera used in the proposed prior model.


### 3.1. The Adaptive Image Prior Model

The adaptive normρ
x
,
a
=
a
a
2
+
x
2
-
a
2 is used to measure the horizontal and vertical gradients in which x is the gradient value and a is a positive scale parameter. These parameters are to determine the severity of penalizing image gradients. A new adaptive image prior model is proposed and defined as follows:(6)
p
u
∣
α
∝
α
1
α
2
N
/
4
exp
⁡
-
∑
i
=
1
N
α
1
ρ
∇
i
x
u
,
a
i
,
1
+
α
2
ρ
∇
i
y
u
,
a
i
,
2
,where N
=
P
1
N
1
×
P
2
N
2 is the number of pixels in the HR image. Symbols ∇
i
x
u and ∇
i
y
u represent the horizontal and vertical gradient components, respectively, for the pixel i, α
=
{
α
1
,
α
2
} is the hyperparameter of this prior controlling the degree of regularization, and a
i
,
1 and a
i
,
2 are the scale parameters of the adaptive norms measuring the horizontal and vertical gradient component, respectively, of pixel i.If we take the logarithm, (6) becomes(7)
log
⁡
p
u
∣
α
∝
-
∑
i
=
1
N
α
1
ρ
∇
i
x
u
,
a
i
,
1
+
α
2
ρ
∇
i
y
u
,
a
i
,
2
,where log
⁡
p
u
∣
α is proportional to a linear combination of ρ
(
(
∇
i
x
u
)
,
a
i
,
1
) and ρ
(
(
∇
i
y
u
)
,
a
i
,
2
), i
=
1
,
…
,
N. Equation (7) is strictly a convex function because ρ
(
x
,
a
) is strictly convex [11].In addition,log
⁡
p
u
∣
α is adaptive, because the following approximations of ρ
(
x
,
a
) can be obtained with the assumption that parameter a is fixed:(8)
ρ
x
,
a
≈
1
2
x
2
when
x
⟶
0
,
(9)
ρ
x
,
a
≈
a
x
-
a
2
when
x
⟶
∞
.Thus, ρ
(
(
∇
i
x
u
)
,
a
i
,
1
) or ρ
(
(
∇
i
y
u
)
,
a
i
,
2
), i
=
1
,
…
,
N, in (7) can use either L
1 norm (i.e., (9)) or L
2 norm (i.e., (8)) adaptively based on the evaluation of the local image features.Thus, the expressionF
=
α
1
ρ
(
∇
i
x
u
,
a
i
,
1
)
+
α
2
ρ
(
∇
i
y
u
,
a
i
,
2
) can be approximated as(10)
F
≈
α
1
2
∇
i
x
u
2
+
α
2
a
i
,
2
∇
i
y
u
-
a
i
,
2
2
,
when
∇
i
x
u
⟶
0
,
∇
i
y
u
⟶
∞
;
(11)
F
≈
α
1
a
i
,
1
∇
i
x
u
-
a
i
,
1
2
+
α
2
2
∇
i
y
u
2
,
when
∇
i
x
u
⟶
∞
,
∇
i
y
u
⟶
0
;
(12)
F
≈
α
1
a
i
,
1
∇
i
x
u
-
a
i
,
1
2
+
α
2
a
i
,
2
∇
i
y
u
-
a
i
,
2
2
,
when
∇
i
x
u
⟶
∞
,
∇
i
y
u
⟶
∞
;
(13)
F
≈
α
1
2
∇
i
x
u
2
+
α
2
2
∇
i
y
u
2
,
when
∇
i
x
u
⟶
0
,
∇
i
y
u
⟶
0
.For example, if there exists an edge along the horizontal direction in an LR image,|
∇
i
y
u
| of the pixel at the edge is relatively large while |
∇
i
x
u
| is small. Thus, a
i
,
2
(
a
i
,
2
)
2
+
|
∇
i
y
u
|
2
-
(
a
i
,
2
)
2 approximates to a
i
,
2
|
∇
i
y
u
|
-
(
a
i
,
2
)
2 and a
i
,
1
(
a
i
,
1
)
2
+
|
∇
i
x
u
|
2
-
(
a
i
,
1
)
2 approximates to 1
/
2
|
∇
i
x
u
|
2. Therefore, F approximates to (10). In the SR process based on (10), when the LR pixel splits in both horizontal and vertical directions into a block of HR pixels, the large gradient magnitude in the vertical direction will be preserved due to the L
1 norm, and the small gradient magnitude in the horizontal direction will be smoothed due to the L
2 norm. Thus whenever there exists an edge along horizontal or vertical direction, a
a
2
+
|
·
|
2
-
a
2 can act as the L
1 norm in the direction perpendicular to the edge, preserving large gradients, and, at the same time, a
a
2
+
|
·
|
2
-
a
2 can act as the L
2 norm in the direction along the edge, obtaining relatively ideal smoothing effects. In (12), when |
∇
i
x
u
|
→
∞ and |
∇
i
y
u
|
→
∞, the L
1 norm acts on both horizontal and vertical directions, to preserve the edges. In (13), when |
∇
i
x
u
|
→
0 and |
∇
i
y
u
|
→
0, the noise will be effectively removed.


### 3.2. Adaptive Calculation of Parametera

It is known that the gradients produced by edges should be preserved, while the gradients produced by errors, which suffer from inaccurate registration and additive Gaussian noise, need to be smoothed. Thus, the severity of penalizing image gradients should be determined according to the local image features. Inρ
(
x
,
a
), the scale parameter, a
, can control the severity of penalizing local gradients. For the gradients produced by errors, L
2 norm is a good choice. In this case, the scale parameter should be set as a large value. When the value of the scale parameter increases, ρ
(
x
,
a
) accepts a larger range of errors and handles them like L
2 norm [11]. Otherwise, it should be set as a small value. Therefore, the property of scale parameter should be constrained with the following conditions: (1) its value is determined according to the local image features, (2) the value is inversely proportional to gradient magnitudes, and (3) the value is larger than zero.We use the well-known monotonically decreasing function [19] to calculate a
i
,
1, and a
i
,
2 in (6),(14)
a
i
,
1
=
1
1
+
∇
i
x
u
2
,
a
i
,
2
=
1
1
+
∇
i
y
u
2
,
1
≤
i
≤
N
,where ∇
i
x
u and ∇
i
y
u represent the magnitudes of horizontal and vertical gradient components, respectively, for the pixel i.For convenience, we use the symbolR to represent a
i
,
1, and a
i
,
2; that is, R
=
{
a
i
,
1
,
a
i
,
2
,
1
≤
i
≤
N
}.


## 3.1. The Adaptive Image Prior Model

The adaptive normρ
x
,
a
=
a
a
2
+
x
2
-
a
2 is used to measure the horizontal and vertical gradients in which x is the gradient value and a is a positive scale parameter. These parameters are to determine the severity of penalizing image gradients. A new adaptive image prior model is proposed and defined as follows:(6)
p
u
∣
α
∝
α
1
α
2
N
/
4
exp
⁡
-
∑
i
=
1
N
α
1
ρ
∇
i
x
u
,
a
i
,
1
+
α
2
ρ
∇
i
y
u
,
a
i
,
2
,where N
=
P
1
N
1
×
P
2
N
2 is the number of pixels in the HR image. Symbols ∇
i
x
u and ∇
i
y
u represent the horizontal and vertical gradient components, respectively, for the pixel i, α
=
{
α
1
,
α
2
} is the hyperparameter of this prior controlling the degree of regularization, and a
i
,
1 and a
i
,
2 are the scale parameters of the adaptive norms measuring the horizontal and vertical gradient component, respectively, of pixel i.If we take the logarithm, (6) becomes(7)
log
⁡
p
u
∣
α
∝
-
∑
i
=
1
N
α
1
ρ
∇
i
x
u
,
a
i
,
1
+
α
2
ρ
∇
i
y
u
,
a
i
,
2
,where log
⁡
p
u
∣
α is proportional to a linear combination of ρ
(
(
∇
i
x
u
)
,
a
i
,
1
) and ρ
(
(
∇
i
y
u
)
,
a
i
,
2
), i
=
1
,
…
,
N. Equation (7) is strictly a convex function because ρ
(
x
,
a
) is strictly convex [11].In addition,log
⁡
p
u
∣
α is adaptive, because the following approximations of ρ
(
x
,
a
) can be obtained with the assumption that parameter a is fixed:(8)
ρ
x
,
a
≈
1
2
x
2
when
x
⟶
0
,
(9)
ρ
x
,
a
≈
a
x
-
a
2
when
x
⟶
∞
.Thus, ρ
(
(
∇
i
x
u
)
,
a
i
,
1
) or ρ
(
(
∇
i
y
u
)
,
a
i
,
2
), i
=
1
,
…
,
N, in (7) can use either L
1 norm (i.e., (9)) or L
2 norm (i.e., (8)) adaptively based on the evaluation of the local image features.Thus, the expressionF
=
α
1
ρ
(
∇
i
x
u
,
a
i
,
1
)
+
α
2
ρ
(
∇
i
y
u
,
a
i
,
2
) can be approximated as(10)
F
≈
α
1
2
∇
i
x
u
2
+
α
2
a
i
,
2
∇
i
y
u
-
a
i
,
2
2
,
when
∇
i
x
u
⟶
0
,
∇
i
y
u
⟶
∞
;
(11)
F
≈
α
1
a
i
,
1
∇
i
x
u
-
a
i
,
1
2
+
α
2
2
∇
i
y
u
2
,
when
∇
i
x
u
⟶
∞
,
∇
i
y
u
⟶
0
;
(12)
F
≈
α
1
a
i
,
1
∇
i
x
u
-
a
i
,
1
2
+
α
2
a
i
,
2
∇
i
y
u
-
a
i
,
2
2
,
when
∇
i
x
u
⟶
∞
,
∇
i
y
u
⟶
∞
;
(13)
F
≈
α
1
2
∇
i
x
u
2
+
α
2
2
∇
i
y
u
2
,
when
∇
i
x
u
⟶
0
,
∇
i
y
u
⟶
0
.For example, if there exists an edge along the horizontal direction in an LR image,|
∇
i
y
u
| of the pixel at the edge is relatively large while |
∇
i
x
u
| is small. Thus, a
i
,
2
(
a
i
,
2
)
2
+
|
∇
i
y
u
|
2
-
(
a
i
,
2
)
2 approximates to a
i
,
2
|
∇
i
y
u
|
-
(
a
i
,
2
)
2 and a
i
,
1
(
a
i
,
1
)
2
+
|
∇
i
x
u
|
2
-
(
a
i
,
1
)
2 approximates to 1
/
2
|
∇
i
x
u
|
2. Therefore, F approximates to (10). In the SR process based on (10), when the LR pixel splits in both horizontal and vertical directions into a block of HR pixels, the large gradient magnitude in the vertical direction will be preserved due to the L
1 norm, and the small gradient magnitude in the horizontal direction will be smoothed due to the L
2 norm. Thus whenever there exists an edge along horizontal or vertical direction, a
a
2
+
|
·
|
2
-
a
2 can act as the L
1 norm in the direction perpendicular to the edge, preserving large gradients, and, at the same time, a
a
2
+
|
·
|
2
-
a
2 can act as the L
2 norm in the direction along the edge, obtaining relatively ideal smoothing effects. In (12), when |
∇
i
x
u
|
→
∞ and |
∇
i
y
u
|
→
∞, the L
1 norm acts on both horizontal and vertical directions, to preserve the edges. In (13), when |
∇
i
x
u
|
→
0 and |
∇
i
y
u
|
→
0, the noise will be effectively removed.


## 3.2. Adaptive Calculation of Parametera

It is known that the gradients produced by edges should be preserved, while the gradients produced by errors, which suffer from inaccurate registration and additive Gaussian noise, need to be smoothed. Thus, the severity of penalizing image gradients should be determined according to the local image features. Inρ
(
x
,
a
), the scale parameter, a
, can control the severity of penalizing local gradients. For the gradients produced by errors, L
2 norm is a good choice. In this case, the scale parameter should be set as a large value. When the value of the scale parameter increases, ρ
(
x
,
a
) accepts a larger range of errors and handles them like L
2 norm [11]. Otherwise, it should be set as a small value. Therefore, the property of scale parameter should be constrained with the following conditions: (1) its value is determined according to the local image features, (2) the value is inversely proportional to gradient magnitudes, and (3) the value is larger than zero.We use the well-known monotonically decreasing function [19] to calculate a
i
,
1, and a
i
,
2 in (6),(14)
a
i
,
1
=
1
1
+
∇
i
x
u
2
,
a
i
,
2
=
1
1
+
∇
i
y
u
2
,
1
≤
i
≤
N
,where ∇
i
x
u and ∇
i
y
u represent the magnitudes of horizontal and vertical gradient components, respectively, for the pixel i.For convenience, we use the symbolR to represent a
i
,
1, and a
i
,
2; that is, R
=
{
a
i
,
1
,
a
i
,
2
,
1
≤
i
≤
N
}.


## 4. Variational Bayesian Superresolution

The variational Bayesian inference [14] is utilized to estimate the unknowns. This variational Bayesian technique approximates the true posterior distribution p
(
u
,
{
s
k
}
,
α
,
{
β
k
}
∣
v
) analytically by a tractable distribution q
(
u
,
{
s
k
}
,
α
,
{
β
k
}
) that minimizes the KL divergence, which can measure the difference between the two distributions p
(
u
,
{
s
k
}
,
α
,
{
β
k
}
∣
v
) and q
(
u
,
{
s
k
}
,
α
,
{
β
k
}
). That is to say, by minimizing this KL divergence, we can obtain the optimal approximation distribution. Thus, we can obtain the following expression:(15)
q
Ω
=
arg
⁡
min
q
Ω
⁡
C
KL
q
Ω
∥
p
Ω
∣
v
=
∫
q
Ω
log
⁡
q
Ω
p
Ω
∣
v
d
Ω
=
∫
q
Ω
log
⁡
q
Ω
p
Ω
,
v
d
Ω
+
const
=
E
q
Ω
,
v
+
const
, where Ω
=
{
u
,
{
s
k
}
,
{
β
k
}
,
α
} and p
(
Ω
,
v
)
=
p
v
∣
u
,
s
k
,
β
k
p
u
∣
α
p
α
p
β
k
p
s
k.Thus, (11) can be rewritten as(16)
q
Ω
=
arg
⁡
min
q
Ω
⁡
E
q
Ω
,
v
.Because (16) cannot be evaluated due to the form of (6), the inequality e
≤
e
+
g
/
2
g with e
≥
0, g
>
0 is employed in (6). This inequality was used in [14, 15]. Then, we can obtain(17)
p
u
∣
α
≥
α
1
α
2
N
/
4
exp
⁡
-
∑
i
=
1
N
α
1
a
i
,
1
∇
i
x
u
2
+
a
i
,
1
2
+
g
i
x
2
g
i
x
-
a
i
,
1
2
+
α
2
a
i
,
2
∇
i
y
u
2
+
a
i
,
2
2
+
g
i
y
2
g
i
y
-
a
i
,
2
2
=
Z
α
,
u
,
g
, where g
=
g
i
x
,
g
i
y
∣
i
=
1
,
…
,
N, the auxiliary variables g
i
x and g
i
y are related to the unknown HR image, and they need to be estimated (see (22)). Thus, we can obtain the upper bound of E
(
q
(
Ω
)
,
v
)
(18)
E
q
Ω
,
v
≤
∫
q
Ω
log
⁡
q
Ω
K
Ω
,
v
,
g
=
E
^
q
Ω
,
v
,
g
,where K
(
Ω
,
v
,
g
)
=
p
v
∣
u
,
s
k
,
β
k
Z
(
α
,
u
,
g
)
p
(
α
)
p
(
{
β
k
}
)
p
(
{
s
k
}
).As shown in [20], the minimization of E
(
q
(
Ω
)
,
v
) can be replaced by the minimization of its upper bound E
^
(
q
(
Ω
)
,
v
,
g
). Thus, we can obtain the distributions q
(
u
), q
(
{
s
k
}
), q
(
{
β
k
}
), and q
(
α
) and g by minimizing E
^
(
q
(
Ω
)
,
v
,
g
); that is, (19)
q
Ω
=
arg
⁡
min
q
Ω
⁡
E
^
q
Ω
,
v
.The details of equations derivation are given in Appendix at the end of this paper. We obtain for the posterior distributionq
(
u
) the multivariate Gaussian with mean value(20)
u
=
Σ
u
∑
k
β
k
B
s
k
T
v
kand covariance (21)
Σ
u
=
∑
k
β
k
B
s
k
T
B
s
k
+
∑
k
∑
i
∑
j
ξ
k
i
j
O
k
i
s
k
T
O
k
j
s
k
+
α
1
∇
x
T
W
g
x
∇
x
+
α
2
∇
y
T
W
g
y
∇
y
-
1
,where · denotes the expected value of ·, T denotes the transposed operator, ξ
k
i
j is a 3 × 3 region in the covariance matrix Λ
k in (24), O
k
r
s
k
=
A
H
k
N
k, r
=
1,2
,
3, where N
1, N
2, and N
3 denote partial derivatives of C
(
s
k
) for θ
k, x
k, and y
k, respectively, and W
g, ∀
g
∈
R
+
N is a diagonal matrix, in which the element on the diagonal is W
g
i
i
=
1
/
g
i.Then, the following expressions are obtained for:(22)
g
i
x
=
∇
i
x
u
2
+
a
2
q
u
,
g
i
y
=
∇
i
y
u
2
+
a
2
q
u
,where ·
q
∗ denotes the expected value of · using q
∗.We also obtain for the posterior distributionq
(
{
s
k
}
) the multivariate Gaussian with mean value(23)
s
k
=
Λ
k
δ
k
-
1
s
k
0
+
β
k
Γ
k
s
k
+
Ψ
k
s
k
+
Q
k
-
Φ
kand covariance(24)
Λ
k
=
δ
k
-
1
+
β
k
Γ
k
+
Ψ
k
-
1
.In (23) and (24), Φ
k
i
=
t
r
a
c
e
B
s
k
T
O
k
j
s
k
Σ
u,  Ψ
k
i
j
=
t
r
a
c
e
O
k
i
s
k
T
O
k
j
s
k
Σ
u,  Γ
k
i
j = Υ
k
i
T
Υ
k
j with Υ
k
=
O
k
1
s
k
u
,
O
k
2
s
k
u
,
O
k
3
s
k
u, and Q
k
i
=
v
k
-
B
s
k
T
O
k
i
s
k
u, for i
,
j
=
1
,
2
,
3.Finally, we obtain the mean values of the hyperparameter distributions, which are used to estimate hyperparameters,(25)
β
k
=
N
1
N
2
/
2
+
a
β
k
1
/
2
v
k
-
B
s
k
u
2
2
q
u
+
b
β
k
,
(26)
α
1
=
N
+
a
α
1
2
∑
i
a
g
i
x
-
a
2
+
b
α
1
,
(27)
α
2
=
N
+
a
α
2
2
∑
i
a
g
i
y
-
a
2
+
b
α
2
.At each iteration step, after obtaining the HR image, the scale parametersR are updated using (10). We now conclude our algorithm (Algorithm 1).Algorithm 1:Our iterative SR algorithm.
Input:
The observationsv
k; The motion parameters s
k
0; The hyper-parameters β
k
0,
α
0; The scale parameters  R
0; The parameters g
0.
Output:
The HR imageu;
(1) Use the Bicubic interpolation of the first LR imagev
1 as initial value u
0 for the HR image.
(2)n
=
1;
(3)repeat
(4)  Computeu
n by using (20) and (21);
(5)  Computeg
n  by using (22);
(6)  Computes
k
n by using (23) and (24);
(7)  Computeβ
k
n by using (25);
(8)  Computeα
n by using (26) and (27);
(9)  ComputeR
n by using (14);
(10)n
=
n
+
1;
(11)until convergence criterion is met
(12) Setu
=
u
n.


## 5. Experimental Results

We will give our experimental results in detail in this section. The experimental environment and parameter settings are presented in Section5.1. The evaluation standards and simulated experimental results obtained by four popularly used SR approaches and our proposed SR approach are illustrated in Section 5.2. Some of the results on real dataset are given in Section 5.3.


### 5.1. General Description

We compare the performance of our proposed approach with the bicubic approach, the BEP model [11] based approach, the l
1 model [16] based approach, the TV model [14] based approach, and the SAR model [15] based approach in terms of Peak Signal-to-Noise Ratio (PSNR), Structural Similarity (SSIM) values [21], and thorough visual inspection. The images presented in Figure 1 were chosen for the simulation test. We used MATLAB R2009a on a 3.0 GHz Pentium Dual core computer with 4.0 GB RAM.Figure 1
The 128 × 128 test images: (a) “zebra” and (b) “car.”
(a)
(b)We need to derive a set of LR images from original images in the experiments. These LR images should have been generated with subpixel motion, rotation, blurring, downsampling, and noise addition. In the experiments, we letθ
k
∈ (−3°, 3°), x
k
∈ (−1, 1), y
k
∈ (−1, 1), k = 1, 2, 3, 4, 5, and the transformation of LR images be variant from each other. 3 × 3 average and Gaussian blurring operators with deviation equal to 1 were used in the simulation. The downsampling factors were P
1 = 2 and P
2 = 2. Finally, additive zero-mean white Gaussian noises were added to the LR images with Signal-to-Noise Ratio (SNR) levels of 1 dB, 5 dB, 10 dB, and 30 dB, respectively. Thus, in each level, five LR images were obtained from each of the original images. The initial parameters used in the experiments were set as follows: our proposed method: g
i
x
0
=
∇
i
x
u
0
2
+
a
i
,
1
0
2, g
i
y
0
=
∇
i
y
u
0
2
+
a
i
,
2
0
2, a
i
,
1
0
=
1
/
1
+
|
∇
i
x
u
0
|
2, a
i
,
2
0
=
1
/
1
+
|
∇
i
y
u
0
|
2, α
1
0
=
N
/
4
/
∑
i
(
a
i
,
1
0
g
i
x
0
-
(
a
i
,
1
0
)
2
), and α
2
0
=
N
/
4
/
∑
i
(
a
i
,
2
0
g
i
y
0
-
(
a
i
,
2
0
)
2
); the l
1 method: g
i
x
0
=
∇
i
x
u
0
2, g
i
y
0
=
∇
i
y
u
0
2, α
1
0
=
N
/
4
/
∑
i
g
i
x
0, and α
2
0
=
N
/
4
/
∑
i
g
i
y
0; the TV method: g
i
0
=
∇
i
x
u
0
2
+
∇
i
y
u
0
2, α
0
=
N
/
4
/
∑
i
g
i
0; the SAR method: α
0
=
N
/
2
/
ħ
u
0
2
2 with Laplacian operator ħ; the BEP method: g
i
x
0
=
u
i
0
-
s
x
l
s
y
m
u
i
0
2
+
a
i
,
1
0
2, where a is the scale parameter, s
x
l and s
y
m shift u
i
0 by l and m pixels in the horizontal and vertical directions, respectively, and these parameters are selected to obtain the best reconstructions; consider α
0
=
N
/
4
/
∑
i
g
i
0. For all the methods, a bicubic interpolation of v
1 is used as initial value u
0, β
k
0
=
N
/
v
k
-
B
(
s
k
)
u
0
2
2, the initial motion parameters were estimated using the method proposed in [22], δ
k
=
0, and a
t
=
b
t
=
0, t
∈
α
1
,
α
2
,
β
k. In our work, the convergence of the algorithms was set with the same criterion u
n
+
1
-
u
n
2
2
/
u
n
2
2
≤
1
0
-
5.


### 5.2. Simulation Experiments

We used LR images derived from artificial HR images to test our proposed approach. The PSNR and SSIM are used to evaluate the reconstruction quality of different approaches. The first set of experiments was with the high level of noise (i.e., SNR = 1 dB and SNR = 5 dB). For the second set of experiments, white Gaussian noises with SNR = 10 dB and SNR = 30 dB were used and the motion parameters were estimated in each iteration.In the first set of experiments, we will show the improvement of our proposed approach in heavy noise (i.e., SNR = 1 dB and SNR = 5 dB). Tables1 and 2 show the PSNR values of all the estimated HR images in heavy noise with average blur and Gaussian blur. Based on the PSNR values, our proposed method shows a good performance. Moreover, Tables 1 and 2 give the corresponding SSIM values that also demonstrate the efficiency of our proposed method.Table 1
Comparisons of PSNR and SSIM with average blur and white Gaussian noise with SNR = 1 dB and SNR = 5 dB.
SNR
Bicubic
SAR
l1
TV
BEP
Ours
Zebra
1 dB
PSNR
18.51
21.15
21.62
20.24
21.72
22.27
SSIM
0.5645
0.6894
0.6780
0.6410
0.6953
0.7086
5 dB
PSNR
18.63
22.39
22.63
22.29
22.70
24.06
SSIM
0.5645
0.7858
0.7889
0.7840
0.7905
0.8020
Car
1 dB
PSNR
26.20
30.33
28.27
28.48
29.93
30.73
SSIM
0.7369
0.8365
0.7601
0.7697
0.8025
0.8563
5 dB
PSNR
26.20
32.27
30.35
32.57
31.57
33.03
SSIM
0.7370
0.8938
0.8640
0.8957
0.8475
0.9121Table 2
Comparisons of PSNR and SSIM with Gaussian blur and white Gaussian noise with SNR = 1 dB and SNR = 5 dB.
SNR
Bicubic
SAR
l1
TV
BEP
Ours
Zebra
1 dB
PSNR
18.97
21.44
21.00
21.06
21.46
22.04
SSIM
0.6072
0.7053
0.6733
0.6742
0.7108
0.7255
5 dB
PSNR
19.08
22.90
23.12
23.95
22.74
24.41
SSIM
0.6072
0.8012
0.7977
0.8173
0.7197
0.8236
Car
1 dB
PSNR
26.56
30.27
29.40
29.43
30.41
31.45
SSIM
0.7511
0.8417
0.7771
0.7956
0.8573
0.8660
5 dB
PSNR
26.77
32.86
33.06
33.61
32.16
34.64
SSIM
0.7512
0.8930
0.8869
0.9021
0.6576
0.9288For the second set of experiments, the motion parameters were estimated. In order to further illustrate the effectiveness of our approach, we will present the reconstruction results for “zebra” and “car” images in noise with SNR = 10 dB and SNR = 30 dB, respectively. Tables3 and 4 show the comparisons of PSNR values of bicubic approach, BEP approach, l
1 approach, SAR approach, TV approach, and our approach. Tables 3 and 4 also give the SSIM results of these approaches with average blur and Gaussian blur, respectively. From Tables 3 and 4, we can see the improvement of our proposed approach over other tested methods. Among all the approaches tested, the proposed method achieves the highest PSNR and SSIM values on both test images.Table 3
Comparisons of PSNR and SSIM with average blur and white Gaussian noise with SNR = 10 dB and SNR = 30 dB.
SNR
Bicubic
SAR
l1
TV
BEP
Ours
Zebra
10 dB
PSNR
19.22
23.73
26.41
25.62
25.81
26.90
SSIM
0.6745
0.8247
0.8420
0.8360
0.8102
0.8766
30 dB
PSNR
19.23
37.09
42.88
41.17
43.34
43.87
SSIM
0.6768
0.9207
0.9252
0.8744
0.9261
0.9413
Car
10 dB
PSNR
26.45
34.51
35.95
35.96
34.98
36.32
SSIM
0.7690
0.8898
0.8857
0.8866
0.8816
0.9197
30 dB
PSNR
26.45
38.73
40.29
39.81
44.78
48.74
SSIM
0.7710
0.9036
0.9160
0.9169
0.9288
0.9461Table 4
Comparisons of PSNR and SSIM with Gaussian blur and white Gaussian noise with SNR = 10 dB and SNR = 30 dB.
SNR
Bicubic
SAR
l1
TV
BEP
Ours
Zebra
10 dB
PSNR
19.39
24.20
26.74
25.33
25.80
27.72
SSIM
0.6868
0.8364
0.8537
0.8510
0.8756
0.9167
30 dB
PSNR
19.39
37.49
40.34
39.01
40.97
42.32
SSIM
0.6874
0.9128
0.9208
0.9194
0.9265
0.9491
Car
10 dB
PSNR
26.77
34.98
36.09
36.01
36.26
36.55
SSIM
0.7793
0.8969
0.8956
0.8951
0.8325
0.9246
30 dB
PSNR
26.78
44.11
46.14
46.17
47.58
49.01
SSIM
0.7812
0.9238
0.9246
0.9242
0.8624
0.9468For visual quality comparison, the error images, that is, the difference between the estimated HR image and the original image, are shown in Figures2 and 3. Our proposed approach can preserve edge details well. Figures 2(g)–2(l) and 3(g)–3(l) show the corresponding error images to the reconstructed images obtained with different approaches. Brighter pixels represent a large error. From these error images, the difference between different SR approaches is clearly observed. In general, our proposed approach can obtain the highest PSNR and SSIM values, with the best visual quality.Figure 2
Results obtained by applying different approaches to LR zebra images corrupted with Gaussian blur and white Gaussian noise with SNR = 10 dB and their corresponding error images. Brighter pixels represent a large error in error image.
(a)
Bicubic
(b)
SAR
(c)
l
1
(d)
TV
(e)
BEP
(f)
Ours
(g)
Error image of bicubic interpolation
(h)
Error image of SAR
(i)
Error image ofl
1
(j)
Error image of TV
(k)
Error image of BEP
(l)
Error image of oursFigure 3
Results obtained by applying different approaches to LR car images corrupted with average blur and white Gaussian noise with SNR = 30 dB and their corresponding error images. Brighter pixels represent a large error in error image.
(a)
Bicubic
(b)
SAR
(c)
l
1
(d)
TV
(e)
BEP
(f)
Ours
(g)
Error image of bicubic interpolation
(h)
Error image of SAR
(i)
Error image ofl
1
(j)
Error image of TV
(k)
Error image of BEP
(l)
Error image of ours


### 5.3. Experiments on Real Data

The performance of our approach is tested with real dataset. The datasets are those popular used video sequences, downloaded from Milanfar’s websitehttps://users.soe.ucsc.edu/~milanfar/software/sr-datasets.html. The first ten LR images were used to reconstruct the HR image.Figures4 and 5 present the obtained HR images for two of the image sequences downloaded. In Figures 4 and 5, the HR images obtained by the bicubic and SAR approaches are still blurred. In addition, there exist artifacts in the HR images obtained by the BEP, l
1, and TV approaches. Our approach is superior to the approaches compared.Figure 4
The reconstructed images for “Adyoron” sequence obtained by different approaches. A subimage is enlarged and displayed at the right bottom corner in the image.
(a)
Bicubic
(b)
SAR
(c)
l
1
(d)
TV
(e)
BEP
(f)
OursFigure 5
The reconstructed images for “EIA” sequence obtained by different approaches. A subimage is enlarged and displayed at the right bottom corner in the image.
(a)
Bicubic
(b)
SAR
(c)
l
1
(d)
TV
(e)
BEP
(f)
Ours


## 5.1. General Description

We compare the performance of our proposed approach with the bicubic approach, the BEP model [11] based approach, the l
1 model [16] based approach, the TV model [14] based approach, and the SAR model [15] based approach in terms of Peak Signal-to-Noise Ratio (PSNR), Structural Similarity (SSIM) values [21], and thorough visual inspection. The images presented in Figure 1 were chosen for the simulation test. We used MATLAB R2009a on a 3.0 GHz Pentium Dual core computer with 4.0 GB RAM.Figure 1
The 128 × 128 test images: (a) “zebra” and (b) “car.”
(a)
(b)We need to derive a set of LR images from original images in the experiments. These LR images should have been generated with subpixel motion, rotation, blurring, downsampling, and noise addition. In the experiments, we letθ
k
∈ (−3°, 3°), x
k
∈ (−1, 1), y
k
∈ (−1, 1), k = 1, 2, 3, 4, 5, and the transformation of LR images be variant from each other. 3 × 3 average and Gaussian blurring operators with deviation equal to 1 were used in the simulation. The downsampling factors were P
1 = 2 and P
2 = 2. Finally, additive zero-mean white Gaussian noises were added to the LR images with Signal-to-Noise Ratio (SNR) levels of 1 dB, 5 dB, 10 dB, and 30 dB, respectively. Thus, in each level, five LR images were obtained from each of the original images. The initial parameters used in the experiments were set as follows: our proposed method: g
i
x
0
=
∇
i
x
u
0
2
+
a
i
,
1
0
2, g
i
y
0
=
∇
i
y
u
0
2
+
a
i
,
2
0
2, a
i
,
1
0
=
1
/
1
+
|
∇
i
x
u
0
|
2, a
i
,
2
0
=
1
/
1
+
|
∇
i
y
u
0
|
2, α
1
0
=
N
/
4
/
∑
i
(
a
i
,
1
0
g
i
x
0
-
(
a
i
,
1
0
)
2
), and α
2
0
=
N
/
4
/
∑
i
(
a
i
,
2
0
g
i
y
0
-
(
a
i
,
2
0
)
2
); the l
1 method: g
i
x
0
=
∇
i
x
u
0
2, g
i
y
0
=
∇
i
y
u
0
2, α
1
0
=
N
/
4
/
∑
i
g
i
x
0, and α
2
0
=
N
/
4
/
∑
i
g
i
y
0; the TV method: g
i
0
=
∇
i
x
u
0
2
+
∇
i
y
u
0
2, α
0
=
N
/
4
/
∑
i
g
i
0; the SAR method: α
0
=
N
/
2
/
ħ
u
0
2
2 with Laplacian operator ħ; the BEP method: g
i
x
0
=
u
i
0
-
s
x
l
s
y
m
u
i
0
2
+
a
i
,
1
0
2, where a is the scale parameter, s
x
l and s
y
m shift u
i
0 by l and m pixels in the horizontal and vertical directions, respectively, and these parameters are selected to obtain the best reconstructions; consider α
0
=
N
/
4
/
∑
i
g
i
0. For all the methods, a bicubic interpolation of v
1 is used as initial value u
0, β
k
0
=
N
/
v
k
-
B
(
s
k
)
u
0
2
2, the initial motion parameters were estimated using the method proposed in [22], δ
k
=
0, and a
t
=
b
t
=
0, t
∈
α
1
,
α
2
,
β
k. In our work, the convergence of the algorithms was set with the same criterion u
n
+
1
-
u
n
2
2
/
u
n
2
2
≤
1
0
-
5.


## 5.2. Simulation Experiments

We used LR images derived from artificial HR images to test our proposed approach. The PSNR and SSIM are used to evaluate the reconstruction quality of different approaches. The first set of experiments was with the high level of noise (i.e., SNR = 1 dB and SNR = 5 dB). For the second set of experiments, white Gaussian noises with SNR = 10 dB and SNR = 30 dB were used and the motion parameters were estimated in each iteration.In the first set of experiments, we will show the improvement of our proposed approach in heavy noise (i.e., SNR = 1 dB and SNR = 5 dB). Tables1 and 2 show the PSNR values of all the estimated HR images in heavy noise with average blur and Gaussian blur. Based on the PSNR values, our proposed method shows a good performance. Moreover, Tables 1 and 2 give the corresponding SSIM values that also demonstrate the efficiency of our proposed method.Table 1
Comparisons of PSNR and SSIM with average blur and white Gaussian noise with SNR = 1 dB and SNR = 5 dB.
SNR
Bicubic
SAR
l1
TV
BEP
Ours
Zebra
1 dB
PSNR
18.51
21.15
21.62
20.24
21.72
22.27
SSIM
0.5645
0.6894
0.6780
0.6410
0.6953
0.7086
5 dB
PSNR
18.63
22.39
22.63
22.29
22.70
24.06
SSIM
0.5645
0.7858
0.7889
0.7840
0.7905
0.8020
Car
1 dB
PSNR
26.20
30.33
28.27
28.48
29.93
30.73
SSIM
0.7369
0.8365
0.7601
0.7697
0.8025
0.8563
5 dB
PSNR
26.20
32.27
30.35
32.57
31.57
33.03
SSIM
0.7370
0.8938
0.8640
0.8957
0.8475
0.9121Table 2
Comparisons of PSNR and SSIM with Gaussian blur and white Gaussian noise with SNR = 1 dB and SNR = 5 dB.
SNR
Bicubic
SAR
l1
TV
BEP
Ours
Zebra
1 dB
PSNR
18.97
21.44
21.00
21.06
21.46
22.04
SSIM
0.6072
0.7053
0.6733
0.6742
0.7108
0.7255
5 dB
PSNR
19.08
22.90
23.12
23.95
22.74
24.41
SSIM
0.6072
0.8012
0.7977
0.8173
0.7197
0.8236
Car
1 dB
PSNR
26.56
30.27
29.40
29.43
30.41
31.45
SSIM
0.7511
0.8417
0.7771
0.7956
0.8573
0.8660
5 dB
PSNR
26.77
32.86
33.06
33.61
32.16
34.64
SSIM
0.7512
0.8930
0.8869
0.9021
0.6576
0.9288For the second set of experiments, the motion parameters were estimated. In order to further illustrate the effectiveness of our approach, we will present the reconstruction results for “zebra” and “car” images in noise with SNR = 10 dB and SNR = 30 dB, respectively. Tables3 and 4 show the comparisons of PSNR values of bicubic approach, BEP approach, l
1 approach, SAR approach, TV approach, and our approach. Tables 3 and 4 also give the SSIM results of these approaches with average blur and Gaussian blur, respectively. From Tables 3 and 4, we can see the improvement of our proposed approach over other tested methods. Among all the approaches tested, the proposed method achieves the highest PSNR and SSIM values on both test images.Table 3
Comparisons of PSNR and SSIM with average blur and white Gaussian noise with SNR = 10 dB and SNR = 30 dB.
SNR
Bicubic
SAR
l1
TV
BEP
Ours
Zebra
10 dB
PSNR
19.22
23.73
26.41
25.62
25.81
26.90
SSIM
0.6745
0.8247
0.8420
0.8360
0.8102
0.8766
30 dB
PSNR
19.23
37.09
42.88
41.17
43.34
43.87
SSIM
0.6768
0.9207
0.9252
0.8744
0.9261
0.9413
Car
10 dB
PSNR
26.45
34.51
35.95
35.96
34.98
36.32
SSIM
0.7690
0.8898
0.8857
0.8866
0.8816
0.9197
30 dB
PSNR
26.45
38.73
40.29
39.81
44.78
48.74
SSIM
0.7710
0.9036
0.9160
0.9169
0.9288
0.9461Table 4
Comparisons of PSNR and SSIM with Gaussian blur and white Gaussian noise with SNR = 10 dB and SNR = 30 dB.
SNR
Bicubic
SAR
l1
TV
BEP
Ours
Zebra
10 dB
PSNR
19.39
24.20
26.74
25.33
25.80
27.72
SSIM
0.6868
0.8364
0.8537
0.8510
0.8756
0.9167
30 dB
PSNR
19.39
37.49
40.34
39.01
40.97
42.32
SSIM
0.6874
0.9128
0.9208
0.9194
0.9265
0.9491
Car
10 dB
PSNR
26.77
34.98
36.09
36.01
36.26
36.55
SSIM
0.7793
0.8969
0.8956
0.8951
0.8325
0.9246
30 dB
PSNR
26.78
44.11
46.14
46.17
47.58
49.01
SSIM
0.7812
0.9238
0.9246
0.9242
0.8624
0.9468For visual quality comparison, the error images, that is, the difference between the estimated HR image and the original image, are shown in Figures2 and 3. Our proposed approach can preserve edge details well. Figures 2(g)–2(l) and 3(g)–3(l) show the corresponding error images to the reconstructed images obtained with different approaches. Brighter pixels represent a large error. From these error images, the difference between different SR approaches is clearly observed. In general, our proposed approach can obtain the highest PSNR and SSIM values, with the best visual quality.Figure 2
Results obtained by applying different approaches to LR zebra images corrupted with Gaussian blur and white Gaussian noise with SNR = 10 dB and their corresponding error images. Brighter pixels represent a large error in error image.
(a)
Bicubic
(b)
SAR
(c)
l
1
(d)
TV
(e)
BEP
(f)
Ours
(g)
Error image of bicubic interpolation
(h)
Error image of SAR
(i)
Error image ofl
1
(j)
Error image of TV
(k)
Error image of BEP
(l)
Error image of oursFigure 3
Results obtained by applying different approaches to LR car images corrupted with average blur and white Gaussian noise with SNR = 30 dB and their corresponding error images. Brighter pixels represent a large error in error image.
(a)
Bicubic
(b)
SAR
(c)
l
1
(d)
TV
(e)
BEP
(f)
Ours
(g)
Error image of bicubic interpolation
(h)
Error image of SAR
(i)
Error image ofl
1
(j)
Error image of TV
(k)
Error image of BEP
(l)
Error image of ours


## 5.3. Experiments on Real Data

The performance of our approach is tested with real dataset. The datasets are those popular used video sequences, downloaded from Milanfar’s websitehttps://users.soe.ucsc.edu/~milanfar/software/sr-datasets.html. The first ten LR images were used to reconstruct the HR image.Figures4 and 5 present the obtained HR images for two of the image sequences downloaded. In Figures 4 and 5, the HR images obtained by the bicubic and SAR approaches are still blurred. In addition, there exist artifacts in the HR images obtained by the BEP, l
1, and TV approaches. Our approach is superior to the approaches compared.Figure 4
The reconstructed images for “Adyoron” sequence obtained by different approaches. A subimage is enlarged and displayed at the right bottom corner in the image.
(a)
Bicubic
(b)
SAR
(c)
l
1
(d)
TV
(e)
BEP
(f)
OursFigure 5
The reconstructed images for “EIA” sequence obtained by different approaches. A subimage is enlarged and displayed at the right bottom corner in the image.
(a)
Bicubic
(b)
SAR
(c)
l
1
(d)
TV
(e)
BEP
(f)
Ours


## 6. Conclusions

The errors caused by inaccurate registration and noises in the traditional regularization-based SR methods often produce unsatisfactory results in the reconstruction. Thus, the variational Bayesian method, which can simultaneously estimate the HR image, the motion parameters, and the hyperparameters, has been used to improve the reconstruction quality. However, the existing variational Bayesian approaches cannot adapt to local image features. Therefore, in this paper, a Bayesian SR approach is proposed by designing a new adaptive image prior model, based on an adaptive norm.Our adaptive image prior model can be adjusted automatically based on the evaluation of the local image features; therefore, this new model not only preserves edge details but also avoids artifacts in the smoothed regions. We also propose a method for automatically estimating the scale parameters for the proposed adaptive image prior model. Information needed to determine these scale parameters is updated in each iteration step based on the available estimated HR image, and they are calculated by using a monotonically decreasing function. In our approach, the acquisition process, the HR image, the motion parameters, and the hyperparameters are modeled in a stochastic sense by using a hierarchical Bayesian framework. And all unknowns are estimated by employing the variational Bayesian inference.The experimental results show that the HR images obtained with our SR approach are better than those previously tested. In the proposed adaptive image prior model, we only calculate the gradients of neighboring pixels at the top and left of the center pixel for fast computation. We will investigate further for the prior on the selection of 4-neighbourhood or 8-neighbourhood in order to improve the proposed algorithm in our future work.



---

*Source: 469859-2015-10-05.xml*