# Al1 Norm Based Image Prior Combination in Multiframe Superresolution

Authors: Lei Min; Ping Yang; Lizhi Dong; Wenjin Liu; Shuai Wang; Bing Xu; Yong Liu
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

We address the multiframe superresolution problem using the variational Bayesian method in this paper. In the variational Bayesian framework, the prior is crucial in transferring the ill-posed reconstruction problem to a well-posed one. We propose a prior combination method based on filter bank andl1 norm. Multiple filters are used in our prior model, and the corresponding combination coefficient vector can be estimated by the characteristics of the filtered image and noise. Furthermore, the local adaptive coefficients of every filter are more effective in removing noise and preserving image edges. Extensive experiments demonstrate the advantages of the proposed method.


---

## Body

## 1. Introduction

Spatial resolution is a crucial parameter to indicate the image quality. Images with high spatial resolution have more distinguishable details and are usually desirable in real applications. In this paper, we address the multiframe superresolution (SR) problem. It uses a sequence of undersampled, shifted, rotated, and noisy low resolution (LR) observation images of the same scene to reconstruct a high resolution (HR) image by the algorithm. It has been widely used in remote sensing [1], surveillance [2], medical imaging [3], and astronomy imaging [4].The superresolution method has been developed for more than three decades [5]. Although the literature is rich (see [6, 7] for reviews), it is still an open and wildly investigated topic. The superresolution concept was first proposed by Tsai and Huang in the frequency domain [5]. However, these frequency approaches [8] have difficulty in including spatial prior knowledge. Therefore, methods in spatial domain have become more popular in recent years due to their convenient incorporation of prior knowledge. These methods fall into two categories: deterministic approach, such as projection onto convex sets (POCS) [9] and iterative back projection (IBP) [10]; stochastic approach, including maximum likelihood (ML) [11], maximum a posterior (MAP) [12–14], and Bayesian method [15–17]. Compared with the deterministic approach, the stochastic approach has more flexibility in modeling the imaging formation process and the prior knowledge [5], and using the hierarchical Bayesian method, both the high resolution image and model parameters can be estimated simultaneously and automatically. So in this paper, we use the Bayesian framework to address the superresolution problem.Using the Bayesian theorem, the prior knowledge of the desired HR image and the LR observation data are combined to get the posterior of the HR image and the parameters. Using the MAP principle, the HR image and the parameters maximize the posterior. However, for many prior models, the posterior is intractable mathematically. This can be efficiently solved by the variational Bayesian method, which was first introduced into superresolution by Babacan et al. [16]. To suppress the noise, the prior models are designed to constraint the high spatial frequency component in the HR image. However, important edges and textures are also high frequency signals. So the prior model is very critical to keep a balance between removing noise and preserving edges and textures. How to choose the appropriate prior model is still an open question. The simultaneous autoregressive model (SAR) prior [18] and the total variation (TV) prior [19] were first used in the variational Bayesian framework [16]. Later the l1 prior [20] was introduced and proven to have better performance. In [21], the combination of the sparse prior and nonsparse SAR prior was proposed, but the combination coefficient had to be determined by the hand empirically. In [17], the author suggested a nonstationary image prior combination method based on a general combination of spatially adaptive image filters, and the combination coefficient can be estimated automatically. However, this method favors oversmooth image estimates because of the heavy penalty of the l2 norm in the prior model. Also the reported PSNR values are low in lots of filter combinations. Based on the above problem, we proposed a prior combination method based on filter bank and l1 prior. The high resolution image is reconstructed according to the characteristics of the filtered image, and the locally adaptive property is guaranteed by the majorization-minimization (MM) approximation of the l1 prior. Furthermore, our work is useful in exploring the superresolution limit performance with respect to the prior model in the variational Bayesian framework. Extensive experiments demonstrate the advantages of the proposed method.The paper is organized as follows. In Section2, we formulate the underlying mathematical framework for variational Bayesian superresolution. In Section 3, we present the variational Bayesian SR algorithms using prior combination based on filter bank and l1 norm. Experimental results are given in Section 4. Finally, Section 5 concludes the paper.


## 2. Problem Formulation

### 2.1. Imaging Model

The imaging process is based on a mathematical model that describes the physical process of obtainingL LR images yk,k=1,…,L, from the latent HR image x we desire. We adopt matrix-vector notation such that each yk and x are arranged as N×1 and PN×1 vectors in lexicographic order, where the integer P>1 is the factor of increase in resolution, and P is the magnification factor. The imaging process (Figure 1) consists of warping, blurring, and downsampling, which is modeled as [22, 23](1)yk=DHkMkskx+nk=Bkskx+nk,where nk is the additive observation noise, D is the N×PN downsampling matrix, Hk is the PN×PN blurring matrix, and Mk is the PN×PN warping matrix, encoding the motion of yk with respect to reference coordinate defined by x. In this work, we assume Hk is known and decided by the linear space invariant point spread function (PSF), and Mksk is generated by the motion vector sk=sk1,sk2,sk3=(θk,ck,dk)T, where θk is the rotation angle and ck and dk are horizontal and vertical translation pixels, respectively. We denote the motion vectors by s=skk=1L. The effects of downsampling, blurring, and warping matrix are combined into a N×PN degradation matrix Bk, describing the generation of LR frames from the HR image. The joint imaging model for all LR frames is given by(2)y=Bx+n,where y=y1T,…,yLTT, B=B1T,…,BLTT, n=n1T,…,nLTT, and (·)T is the transposition operator.Figure 1
The imaging process.In this work, we assume that the additive noisenk of LR frame is zero-mean white Gaussian noise with inverse variance (precision) βk, and the conditional distribution of yk is(3)pyk∣x,βk,sk∝βkN/2exp⁡-βk2yk-Bkskx22,where ·2 is the l2 norm of the vector, for a vector a=a1,a2,…,aMT, a2=∑i=1Mai2. Assuming statistical independence of observation noise among LR frames, the conditional probability of LR images y given x, β=(β1,…,βL) and s=skk=1L, is(4)py∣x,β,s∝∏k=1Lpyk∣x,βk,sk.


### 2.2. Image Prior Model

In this paper, we propose the following image prior:(5)px∣α=∏j=1dpxj∣αj=∏j=1d1Zαjexp⁡-αjxj1=∏j=1d1Zαjexp⁡-αj∑i=1PNxji,where xj=Fjx, xj(i) is the ith component of the vector xj, Fj is the PN×PN convolution filter matrix, which is decided by the convolution kernel, and ·1 is the l1 norm of the vector. For a vector a=a1,a2,…,aMT, a1=∑i=1M|ai|. Z(αj) is the partition function that we approximate as(6)Zαj∝αj-PN/2d,and α=(α1,α2,…,αd) is the prior combination coefficient vector which makes this model more adaptable than TV prior and existing l1 norm prior. Here d filters are used in our prior model, and in a hierarchical Bayesian framework, the corresponding combination coefficient vector can be determined by the image and noise characteristics of filtered image. It is fully data-driven and cumbersome parameter-tuning by hand is no longer necessary. Both our proposed prior method and the one in [17] use filter band which contains a group of high-pass filter, such as Laplacian filter and some edge detection filters. Using the filter band, high spatial frequency signals in the reconstructed image are extracted and penalized by the l2 norm for prior model in [17] and by l1 norm for our proposed one. As is well known, for the large amplitude high spatial frequency component, the l2 norm penalizes it heavily than the l1 norm, which means steep edges and textures are oversmoothed by the l2 norm. Also the image prior model can be considered as regularization term and, as is well known, the l1 norm favors a more sparse solution than the l2 norm. And the sparsity of the high spatial frequency signals means crisp edges and local smoothness. The nonsparse property of the solution induced by l2 norm usually means blur edges, especially for low number of LR images. Experiments demonstrate that good performance can be achieved using appropriate filter combinations.


### 2.3. Hyperprior of the Hyperparameters

We call the parameter in the imaging model and image prior model as hyperparameter, and they form the hyperparameter setΘ=β1,…,βL,α1,…,αd. Hyperprior of the hyperparameter gives the prior knowledge on the different hyperparameter θ∈Θ, which will be modeled using(7)pΘ=∏θ∈Θpθ,where pθ is the Gamma hyperprior(8)pθ=Γθ∣aθ,bθ=bθaθΓaθθaθ-1exp⁡-bθθ,where θ>0, aθ>0, bθ>0, Γ· is the Gamma function, and the expectation E(θ)=aθ/bθ. The Gamma hyperpriors are used because they are the conjugate priors [24] for the variance of the Gaussian distribution and convenient to effectively incorporate the prior knowledge.


### 2.4. Prior of the Motion Vectors

In the preprocessing step, the initial motion vectors are useful to accelerate the convergence rate of the SR algorithm. For example, using the algorithms reported in [25], we can acquire an initialization of sk, denoted by s-kp, from the LR images. However, this initialization is usually inaccurate. As in [16], we use the Gaussian distribution(9)psk=Nsk∣s-kp,Λ-kp,ps=∏k=1Lpskto model the prior of this initialization, where Λ-kp is the covariance matrix to describe the uncertainty of the initialization.


## 2.1. Imaging Model

The imaging process is based on a mathematical model that describes the physical process of obtainingL LR images yk,k=1,…,L, from the latent HR image x we desire. We adopt matrix-vector notation such that each yk and x are arranged as N×1 and PN×1 vectors in lexicographic order, where the integer P>1 is the factor of increase in resolution, and P is the magnification factor. The imaging process (Figure 1) consists of warping, blurring, and downsampling, which is modeled as [22, 23](1)yk=DHkMkskx+nk=Bkskx+nk,where nk is the additive observation noise, D is the N×PN downsampling matrix, Hk is the PN×PN blurring matrix, and Mk is the PN×PN warping matrix, encoding the motion of yk with respect to reference coordinate defined by x. In this work, we assume Hk is known and decided by the linear space invariant point spread function (PSF), and Mksk is generated by the motion vector sk=sk1,sk2,sk3=(θk,ck,dk)T, where θk is the rotation angle and ck and dk are horizontal and vertical translation pixels, respectively. We denote the motion vectors by s=skk=1L. The effects of downsampling, blurring, and warping matrix are combined into a N×PN degradation matrix Bk, describing the generation of LR frames from the HR image. The joint imaging model for all LR frames is given by(2)y=Bx+n,where y=y1T,…,yLTT, B=B1T,…,BLTT, n=n1T,…,nLTT, and (·)T is the transposition operator.Figure 1
The imaging process.In this work, we assume that the additive noisenk of LR frame is zero-mean white Gaussian noise with inverse variance (precision) βk, and the conditional distribution of yk is(3)pyk∣x,βk,sk∝βkN/2exp⁡-βk2yk-Bkskx22,where ·2 is the l2 norm of the vector, for a vector a=a1,a2,…,aMT, a2=∑i=1Mai2. Assuming statistical independence of observation noise among LR frames, the conditional probability of LR images y given x, β=(β1,…,βL) and s=skk=1L, is(4)py∣x,β,s∝∏k=1Lpyk∣x,βk,sk.


## 2.2. Image Prior Model

In this paper, we propose the following image prior:(5)px∣α=∏j=1dpxj∣αj=∏j=1d1Zαjexp⁡-αjxj1=∏j=1d1Zαjexp⁡-αj∑i=1PNxji,where xj=Fjx, xj(i) is the ith component of the vector xj, Fj is the PN×PN convolution filter matrix, which is decided by the convolution kernel, and ·1 is the l1 norm of the vector. For a vector a=a1,a2,…,aMT, a1=∑i=1M|ai|. Z(αj) is the partition function that we approximate as(6)Zαj∝αj-PN/2d,and α=(α1,α2,…,αd) is the prior combination coefficient vector which makes this model more adaptable than TV prior and existing l1 norm prior. Here d filters are used in our prior model, and in a hierarchical Bayesian framework, the corresponding combination coefficient vector can be determined by the image and noise characteristics of filtered image. It is fully data-driven and cumbersome parameter-tuning by hand is no longer necessary. Both our proposed prior method and the one in [17] use filter band which contains a group of high-pass filter, such as Laplacian filter and some edge detection filters. Using the filter band, high spatial frequency signals in the reconstructed image are extracted and penalized by the l2 norm for prior model in [17] and by l1 norm for our proposed one. As is well known, for the large amplitude high spatial frequency component, the l2 norm penalizes it heavily than the l1 norm, which means steep edges and textures are oversmoothed by the l2 norm. Also the image prior model can be considered as regularization term and, as is well known, the l1 norm favors a more sparse solution than the l2 norm. And the sparsity of the high spatial frequency signals means crisp edges and local smoothness. The nonsparse property of the solution induced by l2 norm usually means blur edges, especially for low number of LR images. Experiments demonstrate that good performance can be achieved using appropriate filter combinations.


## 2.3. Hyperprior of the Hyperparameters

We call the parameter in the imaging model and image prior model as hyperparameter, and they form the hyperparameter setΘ=β1,…,βL,α1,…,αd. Hyperprior of the hyperparameter gives the prior knowledge on the different hyperparameter θ∈Θ, which will be modeled using(7)pΘ=∏θ∈Θpθ,where pθ is the Gamma hyperprior(8)pθ=Γθ∣aθ,bθ=bθaθΓaθθaθ-1exp⁡-bθθ,where θ>0, aθ>0, bθ>0, Γ· is the Gamma function, and the expectation E(θ)=aθ/bθ. The Gamma hyperpriors are used because they are the conjugate priors [24] for the variance of the Gaussian distribution and convenient to effectively incorporate the prior knowledge.


## 2.4. Prior of the Motion Vectors

In the preprocessing step, the initial motion vectors are useful to accelerate the convergence rate of the SR algorithm. For example, using the algorithms reported in [25], we can acquire an initialization of sk, denoted by s-kp, from the LR images. However, this initialization is usually inaccurate. As in [16], we use the Gaussian distribution(9)psk=Nsk∣s-kp,Λ-kp,ps=∏k=1Lpskto model the prior of this initialization, where Λ-kp is the covariance matrix to describe the uncertainty of the initialization.


## 3. Variational Bayesian Inference

In this section, we use the variational Bayesian method in [16, 24] to infer the HR image x, hyperparameters Θ, and the motion vector s. The Bayesian inference will be based on posterior probability distribution(10)px,Θ,s∣y=px,Θ,s,ypy,where py has fixed value for given LR observation images y.Express the joint probability distributionpx,Θ,s,y in (10) in terms of the hyperprior pΘ, the prior of the motion vectors ps, the image prior px∣α, and the imaging model py∣x,β as(11)px,Θ,s,y=pΘpspx∣Θpy∣x,Θ,s.Since px,Θ,s∣y is intractable due to the fixed unknown py and the l1 norm in the image prior, as in [16], we apply variational methods to approximate px,Θ,s∣y by the simpler distribution qx,Θ,s by minimizing the Kullback-Leibler (KL) divergence, which is given by(12)CKLqx,Θ,s∥px,Θ,s∣y=∫qx,Θ,slog⁡qx,Θ,spx,Θ,s∣ydxdΘds.The KL divergence is nonnegative and equals zero only when qx,Θ,s=px,Θ,s∣y. Substitute px,Θ,s∣y in (12) with (10); we have(13)CKLqx,Θ,s∥px,Θ,s∣y=∫qx,Θ,slog⁡qx,Θ,spypx,Θ,s,ydxdΘds=∫qx,Θ,slog⁡qx,Θ,spx,Θ,s,ydxdΘds+∫qx,Θ,slog⁡pydxdΘds=∫qx,Θ,slog⁡qx,Θ,spx,Θ,s,ydxdΘds+log⁡py.To further make the computation of CKL(qx,Θ,s∥px,Θ,s∣y) tractable, we assume that qx,Θ,s can be factorized as(14)qx,Θ,s=qx∏j=1dqαj∏k=1Lqβk∏k=1Lqsk.Due to the l1 norm in the image prior, we still cannot calculate the KL divergence. We resort to the majorization-minimization (MM) approach [26] to overcome this difficulty. Utilizing the inequality which states that, for real numbers a≥0 and b>0,(15)a≤a2+b22b,we can make the minimization of (13) tractable. For the image prior model in (5), let a=xj(i), and b=wji; using (15), then(16)px∣α=∏j=1d1Zαjexp⁡-αj∑i=1PNxji≥∏j=1d1Zαjexp⁡-αj∑i=1PNxji2+wjiwji=∏j=1dSxj,αj,wj,where wj is a PN×1 auxiliary vector with wji>0 as its ith component. wj is the local adaptive coefficient for the jth filtered image. Using (16), a lower bound of the joint distribution px,Θ,s,y can be found(17)px,Θ,s,y≥pΘps∏j=1dSxj,αj,wjpy∣x,Θ,s=Fx,Θ,s,y,wj.Substitute px,Θ,s,y in (13) with Fx,Θ,s,y,wj and use (17); we have(18)CKLqx,Θ,s∥px,Θ,s∣y≤∫qx,Θ,slog⁡qx,Θ,sFx,Θ,s,y,wjdxdΘds+log⁡py=CKLqx,Θ,s∥Fx,Θ,s,y,wj+log⁡py.Using (18), we have(19)minqx,Θ,sCKLqx,Θ,s∥px,Θ,s∣y≤minqx,Θ,s⁡CKLqx,Θ,s∥Fx,Θ,s,y,wj+log⁡py=minqx,Θ,s⁡CKLqx,Θ,s∥Fx,Θ,s,y,wj,where the last equality holds because the solution of the minimization problem (19) is not effected by the fixed valued log⁡py even if we do not know the exact value of it.Since inequality (19) holds for any distribution qx,Θ,s and auxiliary variable {wj}, we can solve the problem minqx,Θ,sCKL(qx,Θ,s∥Fx,Θ,s,y,wj) to get an approximation minimum of CKL(qx,Θ,s∥px,Θ,s∣y). We circularly update {wj} and the factors of qx,Θ,s in (14) given by the minimum of CKL(qx,Θ,s∥Fx,Θ,s,y,{wj}) with other distributions held constant. For fixed {wj}, to solve the problem(20)minqx,Θ,sCKLqx,Θ,s∥Fx,Θ,s,y,wj,the standard solutions of variational Bayesian methods [24] can be used to estimate the unknown distribution by(21)qξ=const∗exp⁡E\ξlog⁡Fx,Θ,s,y,wj,where E\ξ[·] denotes expected value with respect to all stochastic variable (i.e., x,Θ,s) with ξ∈{x,Θ,s} removed. For fixed qx,Θ,s, we differentiate CKL(qx,Θ,s∥Fx,Θ,s,y,{wj}) with respect to {wj} to find the minimum to tighten the upper bound in (19). By this procedure, we get a sequence of distribution qnx,Θ,s and auxiliary variable {wjn}, which satisfy(22)CKLqnx,Θ,s∥Fx,Θ,s,y,wjn≥CKLqn+1x,Θ,s∥Fx,Θ,s,y,wjn≥CKLqn+1x,Θ,s∥Fx,Θ,s,y,wjn+1,where the first inequality holds because qn+1x,Θ,s is the solution of the problem minqx,Θ,sCKL(qx,Θ,s∥Fx,Θ,s,y,{wjn}), and the second one holds because {wjn+1}=min{wj}CKL(qn+1x,Θ,s∥Fx,Θ,s,y,{wj}).The proximity of the estimated posterior distributionqx,Θ,s=limn→∞qnx,Θ,s to the true posterior px,Θ,s∣y is still an open question. However, the experiments demonstrate that it is a good approximation. The following (A), (B), (C), and (D) give the estimation of the unknown distributions and auxiliary variable.(A) Estimation of the HR Image Distribution. From (21), the distribution of HR image x is(23)qx∝exp⁡-12∑j=1dαj∑i=1PNxji2wji-12∑k=1Lβkyk-Bkskx22sk,which is Gaussian distribution(24)qx=Nx∣x,Σx,where the mean and inverse variance are(25)x=Σx∑k=1LβkBkskTskyk,(26)Σx-1=∑j=1dαjFjTWjFj+∑k=1LβkBkskTBksksk,where(27)Wj=diag1wjii=1PN.In this paper, for a random variableξ, ξ is the mean value of it, and f(ξ)ξ is the mean value of f(ξ) with respect to ξ; for a vector a, diag(a) generates a square diagonal matrix with the elements of vector a on the main diagonal; for a matrix A, diag(A) generates a column vector using the main diagonal elements of A.(B) Estimation of Motion Vector Distributions. From (21), the distribution of the motion vector sk is(28)qsk∝exp⁡-12βkyk-Bkskx22x-12sk-s-kpTΛ-kp-1sk-s-kp.After simple algebra manipulation, qsk can be expressed as multivariate Gaussian distribution(29)qsk=Nsk∣sk,Λk.(C) Update the Auxiliary Variable
{
w
j
}
. We can update {wj} by differentiating CKL(qx,Θ,s∥Fx,Θ,s,y,{wj}) with respect to {wj} and set the result equal to zero vector, that is, by(30)∂CKLqx,Θ,s∥Fx,Θ,s,y,wj∂wj=0,j=1,2,…,d.Equation (30) gives the unique solution(31)wj=diagFjxxTFjTx=diagFjΣxFjT+FjxxTFjT,j=1,2,…,d.(D) Estimation of the Hyperparameter Distributions. Using (21), we find that the distributions qαj and qβk are still Gamma distributions, expressed as(32)qαj∝αjaαj-1+PN/2exp⁡-αjbαj+∑i=1PNxji2x+wji2wji,(33)qβk∝βkaβk-1+Nexp⁡-βkbβk+yk-Bkskx22x,sk2,where aαj, bαj, aβk, and bβk are the parameters of the hyperpriors. Substitute the first equality of (31) into (32); we get(34)qαj∝αjaαj-1+PN/2exp⁡-αjbαj+∑i=1PNwji,and, as a result,(35)αj=aαj+PN/2bαj+∑i=1PNwji,βk=2aβk+N2bβk+yk-Bkskx22x,sk.In (A), (B), and (D), the explicit form of distribution qx, qsk, qαj, and qβk depends on expectation values BkskTsk, BkskTBksksk, yk-Bkskx22x, and yk-Bkskx22x,sk. Since Bksk is nonlinear with respect to sk, we use the first-order Taylor approximation at sk as in [16]. For the above estimated Gaussian distribution qx=Nx∣x,Σx of high resolution image x, the mean value x maximizes the posterior and is the optimal estimation of x. The algorithm is summarized in Algorithm 1.Algorithm 1:Variational Bayesian superresolution.
Given the valuesaαj, bαj, aβk, bβk, s-kp, Λ-kp, initial high resolution image x^0, initial
auxiliary variablewj0. Set the initial distribution q0x=δ(x-x^0),
q
0
s
k
=
N
s
k
∣
s
-
k
p
,
Λ
-
k
p and n=1.
while convergence criterion is not met do
(1) Given qn-1x and qn-1sk, obtain wjn, qnαj, qnβk using Eq. (31),
Eq. (32), Eq. (33), respectively.
(2) Given qn-1x and qnβk, obtain qnsk using Eq. (28).
(3) Given qnαj, qnβk, wjn and qnsk, obtain qnx using Eq. (23).


## 4. Experiments and Results

### 4.1. Experimental Setup

This section presents experimental results of our proposed method on simulated and real images and compares it with(1) bicubic interpolation (BBC), that is, bicubic interpolation of the LR reference frame, (2) the variational SR method using simultaneous autoregressive model (SAR) prior [16] (denoted by SAR), (3) the variational SR method using TV prior model [16] (denoted by TV), and (4) the nonstationary image prior combination method [17] with the filter combination NF3, because it has the best performance in the combination selection (denoted by NS_NF3). It is easy to see that the l1 prior method in [20] is a special case of our proposed method. We have not compared our method with the one in [21], because how to determine the combination coefficients is not given in [21].For the sake of comparison, the filters in Section2.2 and the filter combinations are mainly designed as the ones in [17] except the Haar filters. The image filters used in this paper are the first-order difference (f.o.d.) filters f1, f2, f3, and f4, the second-order Laplacian filter f5, the second-order difference filters f6, f7, f8, and f9, the Prewitt filters f10 and f11, the Sobel filters f12 and f13, and the Haar filters f14, f15, and f16. All the corresponding convolution kernels are listed below:(36)f1=000-110000,f2=0-10010000,f3=00-1010000,f4=00001000-1,f5=0-10-14-10-10,f6=000-12-1000,f7=0-100200-10,f8=00-1020-100,f9=-10002000-1,f10=-1-1-1000111,f11=-101-101-101,f12=-1-2-1000121,f13=-101-202-101,f14=-110-110000,f15=-1-10110000,f16=1-10-110000.In our method, we use the combinations: NF2=f1,f2, NF3=f1,f2,f5, and NF4=f14,f15,f16, and our method with these combinations is denoted by NF2, NF3, and NF4, respectively. Our proposed method with NF2 is equivalent to the l1 prior method in [20].The initial values of the algorithm are chosen as follows: a noninformative hyperprior model assumption, that is,aθ=bθ=0,θ∈Θ. For the motion vectors, the initial values are estimated using the LK method [25], assuming that (Λ-kp)-1 is equal to zero matrix. The HR image estimate is initialized using bicubic interpolation of the reference frame. With the initialized HR image, the parameters wj, αj, and βk are initialized using the formula in (31) and (35) with the expectation operator removed.We usexn-xn-12/xn-12<10-5 as convergence criterion, where xn and xn-1 are the image estimates at iterations n and n-1, respectively.


### 4.2. Assessment of SR Algorithms

For experiments with simulated images, the objective methods have been used to assess the results of SR algorithms. With ground truth HR images, peak signal-to-noise ratio (PSNR) and structural similarity measure (SSIM) are used to assess the reconstructions of the SR algorithms. PSNR is defined as(37)PSNR=10log10NPx^-x2.SSIM is defined as(38)SSIM=2μxμx^+C12σxx^+C2μx2+μx^2+C1σx2+σx^2+C2,where x is the ground truth grayscale HR image, from which the sequence of LR images is generated, x^ is the reconstructed HR image using the SR algorithms, and C1,C2 are constant. In this paper, we set C1=C2=0.01. μx and μx^ are the mean of the component of vectors x and x^, respectively. σx and σx^ are the standard deviations of the associated images, and σxx^ is defined as(39)σxx^=1PN-1∑i=1PNxi-μxx^i-μx^.In this paper, high resolution images are all reconstructed iteratively for superresolution methods mentioned in this paper except the BBC. And for experiments with real images, the ground truth images are not available. So we choose the reconstruction images empirically such that they produce the visually most appealing results in the iteration process.


### 4.3. Experiments with Simulated Images

The simulated images are generated using the grayscale high resolution images in Figure2. For every image, synthetic sequences of 5 LR images are generated according to (1), that is, the following step: (1) warp, including transition and rotation, the transitions are(40)00,00.5,0.50,10,01pixels, respectively, and the rotation angles are(41)0∘,3∘,-3∘,5∘,-5∘;(2) blur using isotropic 3×3 Gaussian PSF with standard deviation 1; (3) downsample the row and column of the image by factor of P=2; (4) add independent identically distributed (i.i.d.) Gaussian noise of SNR (signal noise ratio): 10 dB, 15 dB, 20 dB, 25 dB, and 30 dB; SNR is defined as(42)SNR=10log10σx2σn2,where σx and σn are the standard derivation of the HR image and noise, respectively. We conduct the simulations with 10 different noise realizations at each SNR, and the average and variance of the performance are reported.Figure 2
Images used in the synthetic experiments.
(a)
(b)
(c)
(d)For our proposed method with different filter combinations, the performance comparisons in terms of PSNR and SSIM of the reconstructed image are shown in Figures3, 4, and 5. Here we take Figure 2(b) as an example. The filter combinations in Figure 3 are coupled by the distance between the positive and negative values of the filter, distance 1: (1) f1,f2, (2) f14,f15,f16; distance2: (3) f3,f4; distance 2: (4) f10,f11, and (5) f12,f13. From Figure 3, we can observe that the performance degrades as the distance increases. Figure 4 shows the reconstruction performance for the second-order filters combinations, distance 1: (1) f6,f7, distance 2: (2) f8,f9 and (3) f5. The performance degrades as the distance increases, but f5 has better performance than f8,f9. Figure 5 shows the performance comparison when combining the first- and second-order filters, the performance is mainly decided by the first-order filters, and the addition of the second-order filters slightly degrades the performance.Figure 3
(a) PSNR and (b) SSIM of the reconstructed images from the image in Figure2(b), for our proposed method, using first-order combinations: case 1 f1,f2, case 2 f14,f15,f16, case 3 f3,f4, case 4 f10,f11, and case 5 f12,f13.
(a)
(b)Figure 4
(a) PSNR and (b) SSIM of the reconstructed images from the image in Figure2(b), for our proposed method, using second-order combinations: case 1 f6,f7, case 2 f5, and case 3 f8,f9.
(a)
(b)Figure 5
(a) PSNR and (b) SSIM of the reconstructed images from the image in Figure2(b), for our proposed method, using first- and second-order combinations: case 1 f1,f2,f6,f7, case 2 f1,f2,f8,f9, case 3 f1,f2,f5, case 4 f14,f15,f16,f6,f7, case 5 f14,f15,f16,f8,f9, case 6 f14,f15,f16,f5, case 7 f3,f4,f6,f7, case 8 f3,f4,f8,f9, and case 9 f3,f4,f5.
(a)
(b)Figure6 shows the performance comparison in terms of PSNR and SSIM for the images presented in Figure 2 with different methods. It is observed that the proposed method has the highest PSNR and SSIM for all images with all SNR and this performance advantage is more obvious for low SNR. For our method, the best performance is achieved by filter combination of NF2 and NF4. In the following comparison, only NF2 and NF4 are used for our method.Figure 6
PSNR and SSIM for the reconstructed images for the images in Figure2 with different method: (a)-(b) for Figure 2(a), (c)-(d) for Figure 2(b), (e)-(f) for Figure 2(c), and (g)-(h) for Figure 2(d).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)To further compare the performance of different method, we set up another experiment. Here, we use the synthesized sequenceEIA from the MDSP benchmark dataset [27]. Detailed description about this synthesized sequence is available in [14]. We approximate the PSF using isotropic 3×3 Gaussian function with standard deviation 1, first 15 frames are used to reconstruct the HR images, and the first frame in the sequence is set as the reference frame. Figure 7 shows the HR reconstructions for different method. It is observed that our method has best overall performance. Ringing artifacts can be found for the SAR. The NS_NF3 has the best noise removing performance, but there are some unnatural artifacts around the edges, for example, in the red circle domain of the Figure 7(e). Compared to TV, our proposed method has better noise removing performance. So our proposed method has better trade-off between preserving the edges and removing noise and artifacts.Figure 7
Superresolution reconstructions on theEIA sequence (L = 15 frames, magnification factor 2). (a) Four low resolution frames (frames 1, 4, 7, and 11). Reconstruction using superresolution method: (b) BBC, (c) SAR, (d) TV, (e) NS_NF3, (f) NF2, and (g) NF4.
(a)
(b)
(c)
(d)
(e)
(f)
(g)


### 4.4. Experiments with Real Images

In this section, we demonstrate the performance comparison of our proposed method with real images. The color image sequencecar is downloaded from the website [28]. From Figures 8(a) and 9(a), we can observe that it is a challenging sequence because the LR car images are severely degraded by the blur and unknown noise model. We reconstruct each color channel separately using the superresolution algorithm. We approximate the PSF using isotropic 3×3 Gaussian function with standard deviation 1, and the first frame in the sequence is set as the reference frame.Figure 8
Superresolution reconstruction on thecar sequence (L = 5 frames, magnification factor 2). (a) The first four low resolution frames. Reconstruction using superresolution method: (b) BBC, (c) SAR, (d) TV, (e) NS_NF3, (f) NF2, and (g) NF4.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Figure 9
Superresolution reconstruction on thecar sequence (L = 15 frames, magnification factor 2). (a) Four low resolution frames (frames 1, 4, 7, and 11). Reconstruction using superresolution method: (b) BBC, (c) SAR, (d) TV, (e) NS_NF3, (f) NF2, and (g) NF4.
(a)
(b)
(c)
(d)
(e)
(f)
(g)For thecar sequence, Figures 8 and 9 show the reconstructions with magnification factor P=2, using L=5 and L=15 LR frames, respectively. As expected, it can be observed that the reconstructions are better with 15 than 5 LR frames. From Figure 8, we can clearly see that the proposed method has better performance: the texts on the window of car are oversmoothed by the NS_NF3 method; by the red circle inserted in Figure 8(c), we can see the ringing artifact produced by the SAR method; near the vertical boundary of Figure 8(d), we can observe the unnatural and obvious artifacts produced by the TV method; also compared to TV and SAR method, our method has better noise removing performance with crisper edges. Although the visual performance difference in Figure 9 is not obvious as that in Figure 8, similar conclusions can be acquired by careful observation: especially, the NS_NF3 method oversmooths the texts on the window of the car; the artifacts near the vertical boundary in Figures 9(c) and 9(d) are obvious for SAR and TV method. Although the NS_NF3 method has the best noise removing performance, the edges are also oversmoothed by the l2 norm in the prior model, especially for the low number of frames when the prior has a strong effect on the reconstructed HR image; this can be observed by comparing Figures 8(e) and 9(e). To sum up, by the results in Figures 8 and 9, our proposed method has better trade-off between preserving the edges and removing the noise and artifacts.


## 4.1. Experimental Setup

This section presents experimental results of our proposed method on simulated and real images and compares it with(1) bicubic interpolation (BBC), that is, bicubic interpolation of the LR reference frame, (2) the variational SR method using simultaneous autoregressive model (SAR) prior [16] (denoted by SAR), (3) the variational SR method using TV prior model [16] (denoted by TV), and (4) the nonstationary image prior combination method [17] with the filter combination NF3, because it has the best performance in the combination selection (denoted by NS_NF3). It is easy to see that the l1 prior method in [20] is a special case of our proposed method. We have not compared our method with the one in [21], because how to determine the combination coefficients is not given in [21].For the sake of comparison, the filters in Section2.2 and the filter combinations are mainly designed as the ones in [17] except the Haar filters. The image filters used in this paper are the first-order difference (f.o.d.) filters f1, f2, f3, and f4, the second-order Laplacian filter f5, the second-order difference filters f6, f7, f8, and f9, the Prewitt filters f10 and f11, the Sobel filters f12 and f13, and the Haar filters f14, f15, and f16. All the corresponding convolution kernels are listed below:(36)f1=000-110000,f2=0-10010000,f3=00-1010000,f4=00001000-1,f5=0-10-14-10-10,f6=000-12-1000,f7=0-100200-10,f8=00-1020-100,f9=-10002000-1,f10=-1-1-1000111,f11=-101-101-101,f12=-1-2-1000121,f13=-101-202-101,f14=-110-110000,f15=-1-10110000,f16=1-10-110000.In our method, we use the combinations: NF2=f1,f2, NF3=f1,f2,f5, and NF4=f14,f15,f16, and our method with these combinations is denoted by NF2, NF3, and NF4, respectively. Our proposed method with NF2 is equivalent to the l1 prior method in [20].The initial values of the algorithm are chosen as follows: a noninformative hyperprior model assumption, that is,aθ=bθ=0,θ∈Θ. For the motion vectors, the initial values are estimated using the LK method [25], assuming that (Λ-kp)-1 is equal to zero matrix. The HR image estimate is initialized using bicubic interpolation of the reference frame. With the initialized HR image, the parameters wj, αj, and βk are initialized using the formula in (31) and (35) with the expectation operator removed.We usexn-xn-12/xn-12<10-5 as convergence criterion, where xn and xn-1 are the image estimates at iterations n and n-1, respectively.


## 4.2. Assessment of SR Algorithms

For experiments with simulated images, the objective methods have been used to assess the results of SR algorithms. With ground truth HR images, peak signal-to-noise ratio (PSNR) and structural similarity measure (SSIM) are used to assess the reconstructions of the SR algorithms. PSNR is defined as(37)PSNR=10log10NPx^-x2.SSIM is defined as(38)SSIM=2μxμx^+C12σxx^+C2μx2+μx^2+C1σx2+σx^2+C2,where x is the ground truth grayscale HR image, from which the sequence of LR images is generated, x^ is the reconstructed HR image using the SR algorithms, and C1,C2 are constant. In this paper, we set C1=C2=0.01. μx and μx^ are the mean of the component of vectors x and x^, respectively. σx and σx^ are the standard deviations of the associated images, and σxx^ is defined as(39)σxx^=1PN-1∑i=1PNxi-μxx^i-μx^.In this paper, high resolution images are all reconstructed iteratively for superresolution methods mentioned in this paper except the BBC. And for experiments with real images, the ground truth images are not available. So we choose the reconstruction images empirically such that they produce the visually most appealing results in the iteration process.


## 4.3. Experiments with Simulated Images

The simulated images are generated using the grayscale high resolution images in Figure2. For every image, synthetic sequences of 5 LR images are generated according to (1), that is, the following step: (1) warp, including transition and rotation, the transitions are(40)00,00.5,0.50,10,01pixels, respectively, and the rotation angles are(41)0∘,3∘,-3∘,5∘,-5∘;(2) blur using isotropic 3×3 Gaussian PSF with standard deviation 1; (3) downsample the row and column of the image by factor of P=2; (4) add independent identically distributed (i.i.d.) Gaussian noise of SNR (signal noise ratio): 10 dB, 15 dB, 20 dB, 25 dB, and 30 dB; SNR is defined as(42)SNR=10log10σx2σn2,where σx and σn are the standard derivation of the HR image and noise, respectively. We conduct the simulations with 10 different noise realizations at each SNR, and the average and variance of the performance are reported.Figure 2
Images used in the synthetic experiments.
(a)
(b)
(c)
(d)For our proposed method with different filter combinations, the performance comparisons in terms of PSNR and SSIM of the reconstructed image are shown in Figures3, 4, and 5. Here we take Figure 2(b) as an example. The filter combinations in Figure 3 are coupled by the distance between the positive and negative values of the filter, distance 1: (1) f1,f2, (2) f14,f15,f16; distance2: (3) f3,f4; distance 2: (4) f10,f11, and (5) f12,f13. From Figure 3, we can observe that the performance degrades as the distance increases. Figure 4 shows the reconstruction performance for the second-order filters combinations, distance 1: (1) f6,f7, distance 2: (2) f8,f9 and (3) f5. The performance degrades as the distance increases, but f5 has better performance than f8,f9. Figure 5 shows the performance comparison when combining the first- and second-order filters, the performance is mainly decided by the first-order filters, and the addition of the second-order filters slightly degrades the performance.Figure 3
(a) PSNR and (b) SSIM of the reconstructed images from the image in Figure2(b), for our proposed method, using first-order combinations: case 1 f1,f2, case 2 f14,f15,f16, case 3 f3,f4, case 4 f10,f11, and case 5 f12,f13.
(a)
(b)Figure 4
(a) PSNR and (b) SSIM of the reconstructed images from the image in Figure2(b), for our proposed method, using second-order combinations: case 1 f6,f7, case 2 f5, and case 3 f8,f9.
(a)
(b)Figure 5
(a) PSNR and (b) SSIM of the reconstructed images from the image in Figure2(b), for our proposed method, using first- and second-order combinations: case 1 f1,f2,f6,f7, case 2 f1,f2,f8,f9, case 3 f1,f2,f5, case 4 f14,f15,f16,f6,f7, case 5 f14,f15,f16,f8,f9, case 6 f14,f15,f16,f5, case 7 f3,f4,f6,f7, case 8 f3,f4,f8,f9, and case 9 f3,f4,f5.
(a)
(b)Figure6 shows the performance comparison in terms of PSNR and SSIM for the images presented in Figure 2 with different methods. It is observed that the proposed method has the highest PSNR and SSIM for all images with all SNR and this performance advantage is more obvious for low SNR. For our method, the best performance is achieved by filter combination of NF2 and NF4. In the following comparison, only NF2 and NF4 are used for our method.Figure 6
PSNR and SSIM for the reconstructed images for the images in Figure2 with different method: (a)-(b) for Figure 2(a), (c)-(d) for Figure 2(b), (e)-(f) for Figure 2(c), and (g)-(h) for Figure 2(d).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)To further compare the performance of different method, we set up another experiment. Here, we use the synthesized sequenceEIA from the MDSP benchmark dataset [27]. Detailed description about this synthesized sequence is available in [14]. We approximate the PSF using isotropic 3×3 Gaussian function with standard deviation 1, first 15 frames are used to reconstruct the HR images, and the first frame in the sequence is set as the reference frame. Figure 7 shows the HR reconstructions for different method. It is observed that our method has best overall performance. Ringing artifacts can be found for the SAR. The NS_NF3 has the best noise removing performance, but there are some unnatural artifacts around the edges, for example, in the red circle domain of the Figure 7(e). Compared to TV, our proposed method has better noise removing performance. So our proposed method has better trade-off between preserving the edges and removing noise and artifacts.Figure 7
Superresolution reconstructions on theEIA sequence (L = 15 frames, magnification factor 2). (a) Four low resolution frames (frames 1, 4, 7, and 11). Reconstruction using superresolution method: (b) BBC, (c) SAR, (d) TV, (e) NS_NF3, (f) NF2, and (g) NF4.
(a)
(b)
(c)
(d)
(e)
(f)
(g)


## 4.4. Experiments with Real Images

In this section, we demonstrate the performance comparison of our proposed method with real images. The color image sequencecar is downloaded from the website [28]. From Figures 8(a) and 9(a), we can observe that it is a challenging sequence because the LR car images are severely degraded by the blur and unknown noise model. We reconstruct each color channel separately using the superresolution algorithm. We approximate the PSF using isotropic 3×3 Gaussian function with standard deviation 1, and the first frame in the sequence is set as the reference frame.Figure 8
Superresolution reconstruction on thecar sequence (L = 5 frames, magnification factor 2). (a) The first four low resolution frames. Reconstruction using superresolution method: (b) BBC, (c) SAR, (d) TV, (e) NS_NF3, (f) NF2, and (g) NF4.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Figure 9
Superresolution reconstruction on thecar sequence (L = 15 frames, magnification factor 2). (a) Four low resolution frames (frames 1, 4, 7, and 11). Reconstruction using superresolution method: (b) BBC, (c) SAR, (d) TV, (e) NS_NF3, (f) NF2, and (g) NF4.
(a)
(b)
(c)
(d)
(e)
(f)
(g)For thecar sequence, Figures 8 and 9 show the reconstructions with magnification factor P=2, using L=5 and L=15 LR frames, respectively. As expected, it can be observed that the reconstructions are better with 15 than 5 LR frames. From Figure 8, we can clearly see that the proposed method has better performance: the texts on the window of car are oversmoothed by the NS_NF3 method; by the red circle inserted in Figure 8(c), we can see the ringing artifact produced by the SAR method; near the vertical boundary of Figure 8(d), we can observe the unnatural and obvious artifacts produced by the TV method; also compared to TV and SAR method, our method has better noise removing performance with crisper edges. Although the visual performance difference in Figure 9 is not obvious as that in Figure 8, similar conclusions can be acquired by careful observation: especially, the NS_NF3 method oversmooths the texts on the window of the car; the artifacts near the vertical boundary in Figures 9(c) and 9(d) are obvious for SAR and TV method. Although the NS_NF3 method has the best noise removing performance, the edges are also oversmoothed by the l2 norm in the prior model, especially for the low number of frames when the prior has a strong effect on the reconstructed HR image; this can be observed by comparing Figures 8(e) and 9(e). To sum up, by the results in Figures 8 and 9, our proposed method has better trade-off between preserving the edges and removing the noise and artifacts.


## 5. Conclusions

In this paper, a prior model combination method, or a class of prior model based on filter bank and thel1 norm, has been proposed for the variational Bayesian superresolution method. High resolution images and all parameters can be estimated automatically and cumbersome parameter-tuning by hand is no longer necessary. Example filters have been designed and used in the proposed method, and the proposed method outperforms the state of art methods in simulated and real data experiments. For the proposed method with different filter combination, filter combinations with small distance have better performance. For all image scenarios with different SNR, the proposed method has the highest PSNR and SSIM values and this performance advantage is more obvious for low SNR. By the visual effect in the simulated and real experiments, it is observed that the proposed method has better trade-off between preserving the edges and removing the noise and artifacts.



---

*Source: 2694638-2017-11-27.xml*