# A Robust Multiframe Image Super-Resolution Method in Variational Bayesian Framework

Authors: Lei Min; Xiangsuo Fan
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Multiframe image super-resolution (MISR) combines complementary information of a set of low-resolution (LR) images to reconstruct a high-resolution (HR) one. In this study, we propose a robust and fully data-driven MISR method in the variational Bayesian framework. Different from the existing variational super-resolution (SR) methods, we use the l1 norm-based observation model, which takes the acquisition noise, outliers, and impulse noise into account. Furthermore, we have evaluated three typical image prior models, and the most appropriate one is chosen for our proposed method. The proposed method has the following advantages: (1) the HR image and all parameters are automatically estimated in an optimal stochastic sense; (2) the algorithm is robust to impulse noise and outliers. Extensive experiments with synthetic and real images demonstrate the advantages of the proposed method.


---

## Body

## 1. Introduction

In most imaging applications, higher spatial resolution is usually desirable. By the Nyquist sampling theorem, the spatial resolution of the digital imaging system is determined by the bandwidth of the optics and sampling frequency of the image sensor. As a result, the most direct approaches to increase the spatial resolution are as follows: (1) increasing the bandwidth of the optics; (2) increasing the sampling frequency of the image sensor. For the former approach, a large aperture or complex lens combination is needed, which usually means big volume and high cost. For the latter solution, as the image sensor pixel patch decreases beyond a specific threshold (which has been reached by the current technologies), the amount of available light on the pixel cell also decreases, and the shot noise will severely degrade the signal-to-noise ratio (SNR). On the other hand, it is usually easy to get a set of LR images of the same scene by a single imaging system over a period of time or by multiple imaging sensors at the same time. The MISR is an algorithm-based approach to increase the image spatial resolution. It models the image degradation process and uses a signal processing algorithm to postprocess the set of LR images. Complementary information provided by the blur and subpixel motion between LR images is used to reconstruct the HR image. It trades off the computational cost and the hardware cost. SR methods have been widely utilized in several imaging fields, such as medical imaging [1], remote sensing [2, 3], surveillance [4], and optical imaging [5].Since Tsai and Huang proposed SR theory in the practice for the first time [2], the SR has been developed for more than three decades. Despite its long history, it is still an active research field and many problems remain unsolved. We refer the reader to [6–9] for a broad review in this area. The image SR methods can be divided into two categories: single-image super-resolution (SISR) and multiframe image super-resolution (MISR). SISR establishes a mapping between LR image and HR image and uses the mapping to map a single LR image to the desired HR image. Traditional machine learning algorithms have been developed to learn the mapping, such as methods based on self-example learning [10], neighborhood embedding [11], and sparse representation [12]. Currently, the mainstream of SISR research is the methods based on deep learning. It mainly includes SRCNN [13], SRGAN [14], STFGAN [15], LapSRN [16], SAN [17], RCAN [18], MAMSR [19–21], and so on, which has gradually achieved better subjective and objective performance.The MISR method combines complementary information of a set of LR images to reconstruct a HR one. The MISR method was first proposed by Tsai and Huang in the frequency domain [2]. However, these frequency approaches [22, 23] are sensitive to model errors and have difficulty in including spatial prior knowledge. But the prior knowledge has paramount importance in transferring the ill-posed SR problem to a well-posed one. Therefore, methods in the spatial domain have become more popular in recent years. These methods include projection onto convex sets (POCSs) [24, 25], iterative back projection (IBP) [26, 27], methods based on reaction-diffusion equation [28], and a group of stochastic methods, including maximum likelihood (ML) [29] and Bayesian methods [4, 30–32]. As pointed out in [6], robustness and flexibility in modeling noise characteristics and the prior knowledge about the solution are the major advantages of the stochastic SR approach. Hence, we use the Bayesian framework to solve the MISR problem.In the Bayesian paradigm, the performance of the SR algorithm is mainly determined by the HR image prior model and the observation model: the image prior model decides the trade-off between removing noise and preserving image details, while the observation model decides the robustness of the algorithm. By the Bayesian theorem, the image prior and the LR image observation are combined to get the posterior of the HR image and the parameters. However, the posterior is usually mathematically intractable. Babacan et al. [31] introduced the variational method to solve this difficulty. The image prior models have been extensively studied in the variational Bayesian framework [32–34]. However, all the mentioned works use a l2 norm-based observation model. As is well known, this observation model is sensitive to outliers and impulse noise. In many important applications, practical systems suffer from impulse noise and outliers [35–37], which are typically caused by motion outliers, transmission errors, camera sensor flaw, and occlusions of the scene. In these applications, l1 norm-based observation model is more appropriate. For image restoration, the l1 norm has been widely used and studied in the regularization framework. References [30, 36, 38–40] have used the l1 norm as data fidelity term, which is corresponding to the observation model in the Bayesian framework. However, in the regularization framework, the regularization parameters are usually determined by time-consuming manual selection. Moreover, an unexpected consequence of this modification (i.e., using the l1 norm-based observation model instead of the l2 one) is that it suggests a data-driven scale/regularization parameter selection technique [40]. In this study, we propose a robust and fully data-driven SR method utilizing the l1 norm-based observation model in the variational Bayesian framework, with the purpose of tackling the impulse noise and outliers and adaptively and efficiently selecting the optimal model parameters. Furthermore, by evaluating three typical prior models, the most appropriate prior model has been chosen for our proposed method.The rest of the study is organized as follows. Section2 introduces the imaging degradation model. In section 3, we model the unknowns in the hierarchical Bayesian framework. In section 4, we present the robust SR algorithm. Experimental results are given in section 5. Finally, section 6 concludes the study.


## 2. Imaging Degradation Model

The imaging degradation model describes the physical process of obtainingL LR images yk,k=1,...,L, from a HR image x. yk and x are arranged as N×1 and PN×1 vectors, where N=Nv×Nh (Nv and Nh are the pixel number of a LR image in the vertical and horizontal directions, respectively), and the integer P is the resolution enhancement factor. In this study, we adopt the matrix-vector notation: the image is represented in lexicographical order as a column vector, and the linear operators of the image, such as blurring and downsampling, are represented by corresponding matrices. The imaging degradation model consists of warping, blurring, downsampling, and adding noise, which is modeled as [41].(1)yk=DHkMkskx+nk=Bkx+nk,where nk, D, Hk, and Mk are, respectively, N×1 observation noise vector, N×PN downsampling matrix, PN×PN blurring matrix, and PN×PN warping matrix. Bk is a N×PN degradation matrix describing the whole effects of downsampling, blurring, and warping. Figure 1 illustrates the imaging degradation process.Figure 1
The block diagram of an imaging degradation process.In this study, we assumeHk is known, and Mksk is generated by the motion vector sk=sk1,sk2,sk3T=θk,ck,dkT, where θk is the rotation angle, and ck and dk are horizontal and vertical translations, respectively. The joint imaging degradation model for all LR images is given by(2)y=Bx+n,where y=y1T,...,yLTT, B=B1T,...,BLTT, and n=n1T,...,nLTT, and the superscript T denotes the transpose of a vector or matrix. Based on equation (2), we estimate the HR image x by all the information we have, including the LR observation y and the prior knowledge of B, x, and n. This estimation process is termed as the SR.


## 3. Hierarchical Bayesian Model

The degradation matrixB in equation (2) is determined by the unknown motion vector and is usually ill conditioned. So, the SR problem is an ill-posed inverse problem with unknown parameters. As a result, noise or small changes in LR images can result in large perturbations in the final estimation. The Bayesian method provides an efficiency approach to solve this inverse problem and the prior knowledge is introduced to transfer the ill-posed problem to a well one. In the Bayesian framework, all unknowns are treated as stochastic quantities and assigned corresponding probability distributions. The following subsections give the individual distribution models of the unknowns.


### 3.1. Proposed Observation Model

In this study, to enhance the robustness of the SR method and deal with impulse noise and outliers, we assume observation noisenk in equation (1) is zero-mean white Laplace noise with parameter βk, and the conditional distribution of yk is(3)pyk|x,βk,sk∝βkNexp−βk2yk−Bkskx1,where ||⋅||1 is the l1 norm of the vector, symbol ∝ means “proportional to,” and a nonessential normalization constant is ignored. For a I×1 vector a=a1,a2,...,aIT, the l1 norm and l2 norm are defined by(4)a1=∑i=1Iai,a2=∑i=1Iai2,, respectively.Assuming statistical independence of observation noise among LR frames, the conditional probability of LR imagesy given x, β=β1,...,βL and s=skk=1L is(5)py|x,β,s∝∏k=1Lpyk|x,βk,sk.This l1 norm-based observation model corresponds to a weighted median filter [17]. Hence, robust performance to the impulse noise and outliers can be expected for our proposed observation model. Another interpretation at the end of section 4 also sheds some light on the robustness of our method.


### 3.2. Image Prior Model

Image prior model encapsulates our prior knowledge and indicates our belief in the HR image. Combining with the observation model, we compare three typical image prior models: the simultaneous autoregressive model (SAR) prior [42], which promotes smoothness, the total variation (TV) prior [43] for its good trade-off between preserving edges and removing noise and artifacts, and the nonstationary image prior combination method (denoted by NS) [32]. They are defined by(6)px|α∝αPN/2exp−α2Cx22,px|α∝αPN/2exp−α∑i=1PNΔihx2+Δivx2,px|α∝∏j=1d∏i=1PNαjiexp−12αjixji2,, respectively, where α is the prior model parameter. For SAR and TV prior model, α is a scalar; for NS prior model, α is a dPN×1 vector. C is Laplacian operator, Δih and Δiv are the horizontal and vertical first order difference at pixel i, respectively, xj=Fjx, xji is the i th component of the vector xj, and Fj is the PN×PN convolution filter matrix, which is decided by predefined image filter.


### 3.3. Hyperprior of the Hyperparameters

We assemble the parameters in the observation model and image prior model as hyperparameter setΘ=β1,...,βL,α. As in Reference [32], Θ will be modeled using(7)pΘ=∏θ∈Θpθ,where pθ is the gamma distribution(8)pθ=Γθ|aθ,bθ=bθaθΓaθθaθ−1exp−bθθ,where θ>0, aθ>0, bθ>0, and Γ⋅ is the gamma function and the expectation Eθ=aθ/bθ. By the variational Bayesian inference, the hyperparameters can be automatically estimated and do not require cumbersome parameters tuning by hand.


### 3.4. Prior of the Motion Vectors

As in Reference [32], we use the Gaussian distribution(9)psk=Nsk||s¯kp,Λ¯kp,and(10)ps=∏k=1Lpsk.To model the prior of the initial motion vectors,Λ¯kp is the covariance matrix. For example, we can use the algorithm in [44] to get s¯kp from the set of LR images.


## 3.1. Proposed Observation Model

In this study, to enhance the robustness of the SR method and deal with impulse noise and outliers, we assume observation noisenk in equation (1) is zero-mean white Laplace noise with parameter βk, and the conditional distribution of yk is(3)pyk|x,βk,sk∝βkNexp−βk2yk−Bkskx1,where ||⋅||1 is the l1 norm of the vector, symbol ∝ means “proportional to,” and a nonessential normalization constant is ignored. For a I×1 vector a=a1,a2,...,aIT, the l1 norm and l2 norm are defined by(4)a1=∑i=1Iai,a2=∑i=1Iai2,, respectively.Assuming statistical independence of observation noise among LR frames, the conditional probability of LR imagesy given x, β=β1,...,βL and s=skk=1L is(5)py|x,β,s∝∏k=1Lpyk|x,βk,sk.This l1 norm-based observation model corresponds to a weighted median filter [17]. Hence, robust performance to the impulse noise and outliers can be expected for our proposed observation model. Another interpretation at the end of section 4 also sheds some light on the robustness of our method.


## 3.2. Image Prior Model

Image prior model encapsulates our prior knowledge and indicates our belief in the HR image. Combining with the observation model, we compare three typical image prior models: the simultaneous autoregressive model (SAR) prior [42], which promotes smoothness, the total variation (TV) prior [43] for its good trade-off between preserving edges and removing noise and artifacts, and the nonstationary image prior combination method (denoted by NS) [32]. They are defined by(6)px|α∝αPN/2exp−α2Cx22,px|α∝αPN/2exp−α∑i=1PNΔihx2+Δivx2,px|α∝∏j=1d∏i=1PNαjiexp−12αjixji2,, respectively, where α is the prior model parameter. For SAR and TV prior model, α is a scalar; for NS prior model, α is a dPN×1 vector. C is Laplacian operator, Δih and Δiv are the horizontal and vertical first order difference at pixel i, respectively, xj=Fjx, xji is the i th component of the vector xj, and Fj is the PN×PN convolution filter matrix, which is decided by predefined image filter.


## 3.3. Hyperprior of the Hyperparameters

We assemble the parameters in the observation model and image prior model as hyperparameter setΘ=β1,...,βL,α. As in Reference [32], Θ will be modeled using(7)pΘ=∏θ∈Θpθ,where pθ is the gamma distribution(8)pθ=Γθ|aθ,bθ=bθaθΓaθθaθ−1exp−bθθ,where θ>0, aθ>0, bθ>0, and Γ⋅ is the gamma function and the expectation Eθ=aθ/bθ. By the variational Bayesian inference, the hyperparameters can be automatically estimated and do not require cumbersome parameters tuning by hand.


## 3.4. Prior of the Motion Vectors

As in Reference [32], we use the Gaussian distribution(9)psk=Nsk||s¯kp,Λ¯kp,and(10)ps=∏k=1Lpsk.To model the prior of the initial motion vectors,Λ¯kp is the covariance matrix. For example, we can use the algorithm in [44] to get s¯kp from the set of LR images.


## 4. Variational Bayesian SR Algorithm

We follow the method in [45] to infer the posterior distributions of x, Θ, and s in the Bayesian framework. Their Bayesian inference will be based on the posterior probability distribution(11)px,Θ,s|y=px,Θ,s,ypy.The joint probability distributionpx,Θ,s,y in equation (11) is expressed in terms of the observation model py|x,β,s, the image prior px|α, the hyperprior pΘ, and the prior of the motion vectors ps as(12)px,Θ,s,y=py|x,Θ,spx|ΘpΘps.Sincepx,Θ,s|y is intractable due to its complicated expression, we apply variational methods to approximate px,Θ,s|y by a simpler distribution qx,Θ,s by minimizing the Kullback-Leibler (KL) divergence, which is given by(13)CKLqx,Θ,s||px,Θ,s|y=∫qx,Θ,slogqx,Θ,spx,Θ,s|ydxdΘds=∫qx,Θ,slogqx,Θ,spx,Θ,s,ydxdΘds+logpy,where logpy is a fixed constant for a given observation y. The KL divergence is non-negative and equal zero only when qx,Θ,s=px,Θ,s|y. Using the mean-field approximation [45], we assume that qx,Θ,s can be factorized as(14)qx,Θ,s=qxqα∏k=1Lqβk∏k=1Lqsk.Due to the TV prior model and l1 norm in the observation model, the calculation of the KL divergence is still intractable. The majorization-minimization (MM) approach [46] is used to overcome this difficulty. Utilizing the inequality, which states that, for real numbers g≥0 and h>0,(15)g≤g2+h22h,we make the minimization of equation (13) tractable. For the observation model in equation (3), let g=yk−Bkskxi and h=vik, using inequality equation (16), then(16)pyk||x,βk,sk∝βkNexp−βkyk−Bkskx1=βkNexp−βk∑i=1Nyk−Bkskxi≥βkNexp−βk2∑i=1Nyk−Bkskxi2+vikvik=Tyk,x,βk,sk,vk,where vk is a N×1 vector with vik>0 as it is i th component. Similarly, for the TV prior model, let g=Δihx2+Δivx2 and h=wi, then(17)px|α∝αPN2exp−α∑i=1PNΔihx2+Δivx2≥αPN2exp−α2∑i=1PNΔihx2+Δivx2+wiwi=Mx,α,w,where w is a PN×1 vector with wi>0 as its i th component. For SAR and NS prior models, the above MM procedure is not needed because the l2 norm is tractable. However, in order to express our method in a unified paradigm, we let Mx,α,w=px|α. Using equations (14) and (15), a lower bound of the joint distribution px,Θ,s,y can be found(18)px,Θ,s,y≥pΘpsMx,α,w∏k=1LTyk,x,βk,sk,vk=Fx,Θ,s,y,w,v,where v=vkk=1L. px,Θ,s,y is substituted in equation (13) with Fx,Θ,s,y,w,v, and equation (18) is used, and we have(19)CKLqx,Θ,s||px,Θ,s|y≤∫qx,Θ,slogqx,Θ,sFx,Θ,s,y,w,vdxdΘds+logpy=CKLqx,Θ,s||Fx,Θ,s,y,w,v+logpy.Here, we can minimize the tight upper bound ofCKLqx,Θ,s||px,Θ,s|y by(20)argminqx,Θ,sCKLqx,Θ,s||Fx,Θ,s,y,w,v+logpy=argminqx,Θ,sCKLqx,Θ,s||Fx,Θ,s,y,w,v,to approximate(21)argminqx,Θ,sCKLqx,Θ,s||px,Θ,s|y.The minimization process of equation (20) has been described in [33], and it is summarized as below, for the n th iteration:(22)vn=minvCKLqn−1xqn−1α∏k=1Lqn−1βk∏k=1Lqn−1sk||Fx,Θ,s,y,wn−1,v,qnβk˜=minqβk˜CKLqβk˜qn−1xqn−1α∏k=1,k≠k˜Lqn−1βk∏k=1Lqn−1sk||Fx,Θ,s,y,wn−1,vn,qnx=minqxCKLqxqn−1α∏k=1Lqnβk∏k=1Lqn−1sk||Fx,Θ,s,y,wn−1,vn,wn=minwCKLqnxqn−1α∏k=1Lqnβk∏k=1Lqn−1sk||Fx,Θ,s,y,w,vn,qnα=minqαCKLqαqn−1x∏k=1Lqn−1βk∏k=1Lqn−1sk||Fx,Θ,s,y,wn,vn−1,qnsk˜=minqsk˜CKLqsk˜qnxqnα∏k=1Lqnβk∏k=1,k≠k˜Lqn−1sk||Fx,Θ,s,y,wn−1,vn.Hence, we can acquire the following iteration equation:(23)vikn=Eqn−1x,qn−1skyk−Bkskxi2,(24)qnβk∝βkaβk−1+N/2exp−βkbβk+∑i=1Nvikn.For the SAR prior model,(25)qnx∝exp−12Eqn−1ααCx22−12∑k=1LEqnβkβk∑i=1NEqn−1skyk−Bkskxi2vikn.(26)qnα∝αaα−1+PN/2exp−αbα+EqnxCx22.For the TV prior model,(27)qnx∝exp−12Eqn−1αα∑i=1PNΔihx2+Δivx2win−1−12∑k=1LEqnβkβk∑i=1NEqn−1skyk−Bkskxi2vikn.(28)win=EqnxΔihx2+Δivx2.(29)qnα∝αaα−1+PN/2exp−αbα+∑i=1PNwin.For the NS prior model,(30)qnx∝∏j=1d∏i=1PNEqn−1αjiαjiexp−12Eqn−1αjiαjixji2exp−12∑k=1LEqnβkβk∑i=1NEqn−1skyk−Bkskxi2vikn,qnαji∝αjiaαji−1exp−αjibαji+Eqnxxji2,qnsk∝exp−12Eqnβkβk∑i=1NEqnxyk−Bkskxi2vikn−12sk−s¯kpTΛ¯kp−1sk−s¯kp,where Ea⋅ is the expectation operator with respect to the random variable a. It should be noted that the difference made by different prior model is in the form of qnx. Once the posterior distribution of the HR image and the parameters are acquired at the end of the iteration, the optimal point estimation of the HR and parameters can be achieved by the maximum a posterior (MAP) principle.It is worth emphasizing here the profound influence made by our proposed observation model to the proposed SR algorithm. For the variational Bayesian SR algorithm in [33, 34], only one parameter βk is used to quantify the noise level of k th LR image. Different from those algorithms, besides the parameter βk, local noise property of the LR images, indicated by the auxiliary vector v=vkk=1L in equation (23), is naturally introduced into the estimation process of qnx and qnsk. From equation (23), it is clear that the auxiliary vector vk represents the local noise level of every LR image. From the estimation formula of qnx and qnsk, we can see that pixels with low noise level will have large contribution to the estimation, and pixels with higher noise level will have less impact to the estimation. So, we can expect that this local noise adaptivity makes the algorithm more robust to the observation noise. This will be demonstrated in the experiments.


## 5. Experiments

### 5.1. Experimental Setup

In this section, we will show the experimental results on simulated and real images. The proposed method was compared with the state-of-the-art SR methods: BBC, bicubic interpolation of the LR reference frame; RSR [30], the robust SR method, which is based on l1-norm data fidelity term and bilateral TV prior; SAR [31], TV [31], and NS [32]: the variational SR methods, which are based on Gaussian observation model and SAR, TV, and NS prior, respectively; and RDESR [28], the super-resolution method based on reaction-diffusion equation. For the NS method, filter combination NF3 was used.The initial values of the iteration process are chosen as follows: the initial motion vectors are estimated using the LK method [44], and assuming Λ¯kp−1 equal to zero, a noninformative hyperprior model assumption, i.e., aθ=bθ=0,θ∈Θ, the HR image is initialized using bicubic interpolation of the LR reference frame; with the above initialization and Σsk=0 and Σx=0, the parameters w and vk and the distribution of βk are initialized using equation (28), and equations (23) and (24), respectively. Then, α is initialized using equation (26), equation (29), and equation (33) for SAR, TV, and NS prior models, respectively.We usedxn−xn−122/xn−122<10−7 as convergence criterion, where xn and xn−1 are acquired by using MAP to qnx and qn−1x, respectively.


### 5.2. Experiments with Simulated Images

The images in Figure2 were chosen to test the SR algorithms. The Lena and EIA image have a size of 120×120, and the size for the other two is 256×256. We normalized the observed LR images to the interval [0, 1]. Synthetic LR image sequences were generated from HR images according to equation (1). For the additive Gaussian noise nk, independent identically distributed (i.i.d.) zero-mean Gaussian noise of different SNR was generated. SNR is defined as(31)SNR=10log10σx2σn2,where σx and σn are the standard derivation of the HR image and noise, respectively. To test the robustness of the SR algorithm and simulate the outliers, we added the impulse noise defined as below to the above-degraded LR images. Let y˜k=Bkx+nk, for the impulse noise, then yk=Nry˜k, where Nr represents the corruption by impulse noise with a corruption rate r. There are mainly two types of impulse noise: (a) salt-and-pepper noise(32)y=Nry˜=y˜i,withprobability1−r,0,withprobabilityr2,1,withprobabilityr2,Figure 2
Original HR images for simulated experiments. (a) Lena, (b) EIA, (c) monarch, and (d) parrot.
(a)(b)(c)(d)(b) random-valued impulse noise(33)y=Nry˜=y˜i,witprobability1−r,nk,withprobabilityr,where nk is an i.i.d. uniformly distributed random variable with values in [0, 1].With ground-truth HR images, the quality of the estimated HR image was quantitatively measured by the peak signal-to-noise ratio (PSNR) [47] in decibels (dB) and structural similarity measure (SSIM) [48]. For the variational Bayesian methods, including SAR, TV, NS, and in our proposed method, the estimated HRs were chosen with the highest PNSR and SSIM in the iteration process. For the RSR and RDESR methods, we tuned the regularization parameter by hand to produce the estimated HR image with the highest PSNR and SSIM. At each noise level, we conducted simulations with 10 different noise realizations and presented the mean PSNR and SSIM.


#### 5.2.1. Effect of Image Prior Model

Before comparing with other SR methods, we study the effect of the image prior model for our proposed observation model. Here, we compared three typical prior models, including SAR using equations (25) and (26), TV using equations (27)-(29), and NS using equations (30)-(33), and denoted by La + SAR, La + TV, and La + NS, respectively. We evaluate the SR reconstruction on EIA and Lena for different noise models in cases where motion vector is available. According to equation (1), five LR images were generated by shifting the HR image with 0,0T, 0,0.5T, 0.5,0T, 1,0T, and 0,1T pixels, rotating it with 0∘, 3∘, −3∘, 5∘, and −5∘, blurring it with 3×3 Gaussian blur kernel with a standard deviation equal to 1, downsampling it by factor of 2, and corrupting it with Gaussian noise and salt-and-pepper noise. The PSNR and SSIM results for the different prior models are shown in Table 1 and 2, respectively. For the low SNR Gaussian noise, La + TV has a slight PSNR and SSIM decrease compared to La + NS. However, for other cases, La + TV has the highest PSNR and SSIM values and improves the PSNR and SSIM by at least 1.5 dB and 0.014, respectively. To see the visual effect of the different prior models, for example, HR restorations on Lena are shown in Figure 3 for the SNR = 35 dB Gaussian noise and r = 1% salt-and-pepper noise. From the red circle inserted in Figure 3, La + TV and La + NS have crisper edges compared with La + SAR; La + NS has better noise removing performance and favors a smooth solution, but some textures are also oversmoothed. The above PSNR, SSIM, and the virtual performance demonstrate that La + TV is more favorable, and we choose it as the image prior model for our proposed SR method.Table 1
PSNR (dB) results for different prior models.
Noise caseEIALenaLa + SARLa + TVLa + NSLa + SARLa + TVLa + NSGaussian 15 dB22.2123.1923.3128.3529.2529.55Gaussian 35 dB26.6731.0029.4733.0335.8133.68Salt-and-pepper 0.123.1825.6523.3330.8732.4030.48Salt-and-pepper 0.0126.0231.6428.8932.3936.0634.25Table 2
SSIM results for different prior models.
Noise caseEIALenaLa + SARLa + TVLa + NSLa + SARLa + TVLa + NSGaussian 15 dB0.96570.97270.97280.98420.98750.9880Gaussian 35 dB0.98810.99570.99390.99470.99720.9954Salt-and-pepper 0.10.97250.98470.97300.99120.99380.9902Salt-and-pepper 0.010.98610.99630.99300.99380.99740.9960Figure 3
Super-resolution on the Lena with PSNR and SSIM in the brackets: for 35 dB Gaussian noise, results of (a) La + SAR (33.02 dB, 0.9947), (b) La + TV (35.78 dB, 0.9972), and (c) La + NS (33.66 dB, 0.9954); for 1% salt-and-pepper noise, results of (d) La + SAR (32.34 dB, 0.9937), (e) La + TV (36.04 dB, 0.9973), and (f) La + NS (34.27 dB, 0.9960).
(a)(b)(c)(d)(e)(f)


#### 5.2.2. Robust to Noise Model

In this subsection, we utilized the same setup as the previous subsection to evaluate the robustness of the proposed method (La + TV in the last subsection) on simulated images, using the dataset EIA. Besides the Gaussian noise and salt-and-pepper noise, the mixed noise model is also used to test the robustness of the SR method. For the mixed noise, it is expressed as “SNR +r,” which means that the LR images are first corrupted by SNR dB Gaussian noise and then corrupted by the random-valued impulse noise with a corruption rate r. The PSNR and SSIM results for different SR methods are shown in Table 3 and Table 4, respectively. As expected, all SR algorithms result in better performance than BBC because it only uses the reference LR frame to reconstruct the HR image. For the low SNR Gaussian noise, our method has a slight PSNR and SSIM decrease compared to the TV and RDESR methods. However, for other cases, our proposed method produces the reconstructed HR images with the highest PSNR and SSIM values. In particular, an obvious performance advantage can be observed for the noise case with impulse noise. To see the visual effect of different methods, for example, HR restorations are shown in Figure 4 for the SNR = 35 dB Gaussian noise, in Figure 5 for r = 0.01 salt-and-pepper noise, and in Figure 6 for 35+0.001 mixed noise. In Figure 4, our method provides comparable visual performance with the TV, NS, and RDESR methods. In Figure 5, we can see that our method is efficient to remove salt-and-pepper noise, the RSR produces obvious ringing artifacts near the edges, and the SAR, TV, and RDESR methods are very sensitive to outliers because the l2 norm-based observation model is not robust against outliers. The NS and RDESR methods can moderately suppress the outliers, but the details are also oversmoothed. Our method provides the most visually appealing results without artifacts while still preserve edges well. In Figure 6, our method has better blur and ringing artifacts removing performance than RSR, only small fractional outliers can severely degrade the performance of SAR, TV, NS, and RDESR, and the specks are obvious in the reconstruction HR images. Our method is efficient, stable, and robust with various noise cases.Table 3
PSNR (dB) results with the accurate motion vector.
Noise casePSNRBBCRSRSARTVNSRDESRProposedGaussian 15 dB16.0823.5022.5223.5323.3423.4223.17Gaussian 35 dB16.5126.5327.0530.9729.4230.1030.99Salt-and-pepper 0.0115.8125.3922.2823.3522.4122.5031.64Salt-and-pepper 0.113.4023.9217.8717.9817.3618.1025.65Mixed 35 + 0.00116.4726.5025.6128.3026.6428.2830.94Mixed 30 + 0.0515.3925.4320.8221.4421.0721.6027.88Table 4
SSIM results with the accurate motion vector.
Noise casePSNRBBCRSRSARTVNSRDESRProposedGaussian 15 dB0.83330.97500.96830.97520.97430.97500.9726Gaussian 35 dB0.84910.98710.98920.99570.99380.99440.9957Salt-and-pepper 0.010.82520.98370.96640.97390.96820.97010.9963Salt-and-pepper 0.10.70750.97660.89840.89700.87590.90230.9847Mixed 35 + 0.0010.84780.98750.98470.99180.98810.99070.9956Mixed 30 + 0.050.80440.98350.95100.95740.95440.95890.9910Figure 4
Super-resolution on the EIA with PSNR and SSIM in the brackets for 35 dB Gaussian noise when the motion information is exact. Results of (a) BBC (16.52 dB, 0.8489), (b) RSR (26.81 dB, 0.9885), (c) SAR (27.04 dB, 0.9891), (d) TV (30.96 dB, 0.9956), (e) NS (29.52 dB, 0.9940), (f) RDESR (30.12 dB, 0.9953), and (g) proposed (30.98 dB, 0.9957).
(a)(b)(c)(d)(e)(f)(g)Figure 5
Super-resolution on the EIA with PSNR and SSIM in the brackets for 1% salt-and-pepper noise when the motion information is exact. Results of (a) BBC (16.18 dB, 0.8378), (b) RSR (25.55 dB, 0.9846), (c) SAR (22.17 dB, 0.9655), (d) TV (23.22 dB, 0.9730), (e) NS (22.21 dB, 0.9665), (f) RDESR (22.31 dB, 0.9689), and (g) proposed (31.55 dB, 0.9962).
(a)(b)(c)(d)(e)(f)(g)Figure 6
Super-resolution on the EIA with PSNR and SSIM in the brackets for 35 + 0.001 mixed noise when the motion information is exact. Results of (a) BBC (16.47 dB, 0.8475), (b) RSR (26.53 dB, 0.9878), (c) SAR (25.37 dB, 0.9839), (d) TV (27.66 dB, 0.9906), (e) NS (26.26 dB, 0.9872), (f) RDESR (26.33 dB, 0.9878), and (g) proposed (30.91 dB, 0.9956).
(a)(b)(c)(d)(e)(f)(g)


#### 5.2.3. Effect of Motion Vector Uncertainty

In this subsection, we conduct experiments to investigate the influence of inaccurate motion vector on the SR methods. We utilized the same setup and dataset as 5.2.2. To simulate the errors in motion estimation, we added a deliberate error in its translation parameters by i.i.d. zero-mean white Gaussian noise with standard deviation equal to 0.2 pixels, and the rotation parameters with noise uniformly distributed in−0..1∘,0.1∘. The mean PSNR and SSIM results for different SR methods are shown in Table 5 and 6, respectively. As expected, the performance of all SR algorithms decreases due to the motion error. The motion error has a slight effect on the variational Bayesian methods. However, the performance degradation for RSR and RDESR is severe (about 2–3 dB PSNR loss). Similarly, for the low SNR Gaussian noise, the TV method has a slight performance advantage (about 0.4 dB for PSNR). For the noise case with impulse noise, the proposed method improves the PSNR by at least 3–4 dB. Example HR reconstructions are shown in Figure 7 for 35+0.001 mixed noise. Compared with Figure 6, similar visual results are obtained for the variational Bayesian method, but the RSR has obvious ringing artifacts near the edge, and the RDESR method gets blurred edges and textures.Table 5
PSNR (dB) results with the inaccurate motion vector.
Noise casePSNRBBCRSRSARTVNSRDESRProposedGaussian 15 dB16.0522.4922.5123.2922.9222.5522.88Gaussian 35 dB16.5123.1627.0530.4627.3426.4130.49Salt-and-pepper 0.0116.1122.9522.2623.0521.9821.1731.03Salt-and-pepper 0.115.3921.7517.8617.9017.2417.0125.26Mixed 35 + 0.00116.4723.3325.6027.4825.7523.5430.57Mixed 30 + 0.0517.8623.2320.8121.2620.7119.3027.73Table 6
SSIM results with the inaccurate motion vector.
Noise casePSNRBBCRSRSARTVNSRDESRProposedGaussian 15 dB0.83330.96790.96830.97370.97170.97310.9709Gaussian 35 dB0.84910.97020.98910.99510.98980.98490.9952Salt-and-pepper 0.010.82520.97010.96640.97190.96490.96370.9957Salt-and-pepper 0.10.70750.95960.89810.89520.87260.86960.9832Mixed 35 + 0.0010.84780.97310.98470.99000.98540.97460.9952Mixed 30 + 0.050.80440.97300.95100.95550.95050.94830.9907Figure 7
Super-resolution on the EIA with PSNR and SSIM in the brackets for 35 + 0.001 mixed noise when the motion information is inexact. Results of (a) BBC (16.47 dB, 0.8475), (b) RSR (23.65 dB, 0.9753), (c) SAR (25.37 dB, 0.9839), (d) TV (27.63 dB, 0.9905), (e) NS (26.07 dB, 0.9866), (f) RDESR (23.58 dB, 0.9749), and (g) proposed (30.86 dB, 0.9955).
(a)(b)(c)(d)(e)(f)(g)


#### 5.2.4. Effect of the Number of Frames with Outliers

In this subsection, we study the effect of the number of frames with outliers on the performance of the SR methods. We generated 17 LR images from monarch by a downsampling factor of 4 and a3×3 Gaussian blurring kernel with a standard deviation equal to 1. All frames were corrupted by SNR dB Gaussian noise, and the first O2≤O≤17 frames were also corrupted by random-valued impulse noise with corruption rate r (denoted by SNR + r(O)). The SR reconstruction was performed with an exact motion vector. The mean PSNR and SSIM results for 25 + 0.05(O) mixed noise are shown in Figure 8. We observe that RSR are robust to outliers but with low PSNR and SSIM values, while SAR, TV, NS, and RDESR are not robust against outliers and the PSNR and SSIM values rapidly decrease when the number of frames with outliers increases. Our method has the highest PSNR and SSIM values and is robust to outliers. Example HR reconstructions are shown in Figure 9 for 25 + 0.05(9). Images in Figure 10 are the cropped regions corresponding to the rectangular area in Figure 9. It is observed that the proposed method has an advantage in achieving a balance between removing noise and preserving the image edges.Figure 8
PSNR and SSIM on the monarch dataset for 25 + 0.05(O).
(a)(b)Figure 9
Super-resolution on the monarch with PSNR and SSIM in the brackets for 25 + 0.05(9) mixed noise. (a) Original HR image, (b) one LR image, and results of (c) BBC (16.46 dB, 0.8147), (d) RSR (24.24 dB, 0.9714), (e) SAR (26.71 dB, 0.9841), (f) TV (27.62 dB, 0.9872), (g) NS (28.65 dB, 0.9900), (h) RDESR (28.72, 0.9908), and (i) proposed (29.73 dB, 0.9922).
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 10
The corresponding cropped images of Figure9.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


#### 5.2.5. Effect of the Sequence Length

In this subsection, we evaluate the influence of the number of LR images on the quality of SR reconstruction. We generatedL5≤L≤19 LR images from the parrot by a downsampling factor of 4 and a 3×3 Gaussian blurring kernel with a standard deviation equal to 1. All frames were corrupted by SNR dB Gaussian noise and then by the random-valued impulse noise with corruption rate r. The SR restoration was performed with an exact motion vector. The PSNR and SSIM results for 30+0.01 mixed noise are shown in Figure 11. The performance of all algorithms increases as the number of LR images increases, as expected. The proposed method clearly outperforms the other methods independent of the number of LR images. Example HR reconstructions are shown in Figure 12 for 30+0.01 with L=9. Notice that this is an underdetermined case for the SR problem (L<16). Figure 13 presents the cropped regions corresponding to the rectangular area in Figure 12. Our proposed method produces the most visually appealing results and is very efficient in removing outliers and preserving detailed texture.Figure 11
Effect of number of LR images to the PSNR and SSIM on the parrot dataset.
(a)(b)Figure 12
Super-resolution on the parrot with PSNR and SSIM in the brackets for 30 + 0.01 mixed noise. (a) Original HR image, (b) one LR image, Results of (c) BBC (22.02 dB, 0.9423), (d) RSR (28.57 dB), (e) SAR (27.41 dB, 0.9835), (f) TV (27.99 dB, 0.9856), (g) NS (27.97 dB, 0.9856), and (h) proposed (31.51 dB, 0.9937).
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 13
The corresponding cropped images of Figure12.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


### 5.3. Experiments with Real Images

In this subsection, we present the performance comparison of SR methods on the challenging real image sequence called alpaca, which was downloaded from the MDSP benchmark dataset [49]. The sequence contains 55 LR frames (32×70), and 9 frames (frame 45 to 53) are used in our experiment. The 9 frames follow a complicated motion model and contain a high level of outliers. We approximated the PSF using an isotropic 3×3 Gaussian kernel with a standard deviation equal to 1.Since no ground truth image is available for the real sequence, we empirically chose the SR reconstructions such that they produce the visually most appealing results in the iteration process for the variational Bayesian method. For RSR and RDESR, the regularization parameters were manually tuned to produce the best performance. Figure14 shows the SR results with magnification factor P=3 using the 9 LR frames. Figure 15 shows the underdetermined case. It is observed that our proposed method produces the most visually appealing results. The SAR and TV methods suffer from visible ghosting artifacts caused by the alpaca’s movements. The NS and RDESR prefer a smooth solution, but important edges and texture are also oversmoothed, especially for the underdetermined case. Compared to the robust RSR method, our proposed method has better noise removing performance with crisper details, and most importantly, parameter tuning by hand is no longer necessary.Figure 14
Super-resolution on the 9 LR frames with magnification factor 3. Frame 49 is set as the reference frame. Results of (a) BBC, (b) RSR, (c) SAR, (d) TV, (e) NS, (f) RDESR, and (g) proposed.
(a)(b)(c)(d)(e)(f)(g)Figure 15
Super-resolution on the five LR frames (frame 46 to 50) with magnification factor 3. Frame 48 is set as the reference frame. Results of (a) BBC, (b) RSR, (c) SAR, (d) TV, (e) NS, (f) RDESR, and (g) proposed.
(a)(b)(c)(d)(e)(f)(g)


## 5.1. Experimental Setup

In this section, we will show the experimental results on simulated and real images. The proposed method was compared with the state-of-the-art SR methods: BBC, bicubic interpolation of the LR reference frame; RSR [30], the robust SR method, which is based on l1-norm data fidelity term and bilateral TV prior; SAR [31], TV [31], and NS [32]: the variational SR methods, which are based on Gaussian observation model and SAR, TV, and NS prior, respectively; and RDESR [28], the super-resolution method based on reaction-diffusion equation. For the NS method, filter combination NF3 was used.The initial values of the iteration process are chosen as follows: the initial motion vectors are estimated using the LK method [44], and assuming Λ¯kp−1 equal to zero, a noninformative hyperprior model assumption, i.e., aθ=bθ=0,θ∈Θ, the HR image is initialized using bicubic interpolation of the LR reference frame; with the above initialization and Σsk=0 and Σx=0, the parameters w and vk and the distribution of βk are initialized using equation (28), and equations (23) and (24), respectively. Then, α is initialized using equation (26), equation (29), and equation (33) for SAR, TV, and NS prior models, respectively.We usedxn−xn−122/xn−122<10−7 as convergence criterion, where xn and xn−1 are acquired by using MAP to qnx and qn−1x, respectively.


## 5.2. Experiments with Simulated Images

The images in Figure2 were chosen to test the SR algorithms. The Lena and EIA image have a size of 120×120, and the size for the other two is 256×256. We normalized the observed LR images to the interval [0, 1]. Synthetic LR image sequences were generated from HR images according to equation (1). For the additive Gaussian noise nk, independent identically distributed (i.i.d.) zero-mean Gaussian noise of different SNR was generated. SNR is defined as(31)SNR=10log10σx2σn2,where σx and σn are the standard derivation of the HR image and noise, respectively. To test the robustness of the SR algorithm and simulate the outliers, we added the impulse noise defined as below to the above-degraded LR images. Let y˜k=Bkx+nk, for the impulse noise, then yk=Nry˜k, where Nr represents the corruption by impulse noise with a corruption rate r. There are mainly two types of impulse noise: (a) salt-and-pepper noise(32)y=Nry˜=y˜i,withprobability1−r,0,withprobabilityr2,1,withprobabilityr2,Figure 2
Original HR images for simulated experiments. (a) Lena, (b) EIA, (c) monarch, and (d) parrot.
(a)(b)(c)(d)(b) random-valued impulse noise(33)y=Nry˜=y˜i,witprobability1−r,nk,withprobabilityr,where nk is an i.i.d. uniformly distributed random variable with values in [0, 1].With ground-truth HR images, the quality of the estimated HR image was quantitatively measured by the peak signal-to-noise ratio (PSNR) [47] in decibels (dB) and structural similarity measure (SSIM) [48]. For the variational Bayesian methods, including SAR, TV, NS, and in our proposed method, the estimated HRs were chosen with the highest PNSR and SSIM in the iteration process. For the RSR and RDESR methods, we tuned the regularization parameter by hand to produce the estimated HR image with the highest PSNR and SSIM. At each noise level, we conducted simulations with 10 different noise realizations and presented the mean PSNR and SSIM.


### 5.2.1. Effect of Image Prior Model

Before comparing with other SR methods, we study the effect of the image prior model for our proposed observation model. Here, we compared three typical prior models, including SAR using equations (25) and (26), TV using equations (27)-(29), and NS using equations (30)-(33), and denoted by La + SAR, La + TV, and La + NS, respectively. We evaluate the SR reconstruction on EIA and Lena for different noise models in cases where motion vector is available. According to equation (1), five LR images were generated by shifting the HR image with 0,0T, 0,0.5T, 0.5,0T, 1,0T, and 0,1T pixels, rotating it with 0∘, 3∘, −3∘, 5∘, and −5∘, blurring it with 3×3 Gaussian blur kernel with a standard deviation equal to 1, downsampling it by factor of 2, and corrupting it with Gaussian noise and salt-and-pepper noise. The PSNR and SSIM results for the different prior models are shown in Table 1 and 2, respectively. For the low SNR Gaussian noise, La + TV has a slight PSNR and SSIM decrease compared to La + NS. However, for other cases, La + TV has the highest PSNR and SSIM values and improves the PSNR and SSIM by at least 1.5 dB and 0.014, respectively. To see the visual effect of the different prior models, for example, HR restorations on Lena are shown in Figure 3 for the SNR = 35 dB Gaussian noise and r = 1% salt-and-pepper noise. From the red circle inserted in Figure 3, La + TV and La + NS have crisper edges compared with La + SAR; La + NS has better noise removing performance and favors a smooth solution, but some textures are also oversmoothed. The above PSNR, SSIM, and the virtual performance demonstrate that La + TV is more favorable, and we choose it as the image prior model for our proposed SR method.Table 1
PSNR (dB) results for different prior models.
Noise caseEIALenaLa + SARLa + TVLa + NSLa + SARLa + TVLa + NSGaussian 15 dB22.2123.1923.3128.3529.2529.55Gaussian 35 dB26.6731.0029.4733.0335.8133.68Salt-and-pepper 0.123.1825.6523.3330.8732.4030.48Salt-and-pepper 0.0126.0231.6428.8932.3936.0634.25Table 2
SSIM results for different prior models.
Noise caseEIALenaLa + SARLa + TVLa + NSLa + SARLa + TVLa + NSGaussian 15 dB0.96570.97270.97280.98420.98750.9880Gaussian 35 dB0.98810.99570.99390.99470.99720.9954Salt-and-pepper 0.10.97250.98470.97300.99120.99380.9902Salt-and-pepper 0.010.98610.99630.99300.99380.99740.9960Figure 3
Super-resolution on the Lena with PSNR and SSIM in the brackets: for 35 dB Gaussian noise, results of (a) La + SAR (33.02 dB, 0.9947), (b) La + TV (35.78 dB, 0.9972), and (c) La + NS (33.66 dB, 0.9954); for 1% salt-and-pepper noise, results of (d) La + SAR (32.34 dB, 0.9937), (e) La + TV (36.04 dB, 0.9973), and (f) La + NS (34.27 dB, 0.9960).
(a)(b)(c)(d)(e)(f)


### 5.2.2. Robust to Noise Model

In this subsection, we utilized the same setup as the previous subsection to evaluate the robustness of the proposed method (La + TV in the last subsection) on simulated images, using the dataset EIA. Besides the Gaussian noise and salt-and-pepper noise, the mixed noise model is also used to test the robustness of the SR method. For the mixed noise, it is expressed as “SNR +r,” which means that the LR images are first corrupted by SNR dB Gaussian noise and then corrupted by the random-valued impulse noise with a corruption rate r. The PSNR and SSIM results for different SR methods are shown in Table 3 and Table 4, respectively. As expected, all SR algorithms result in better performance than BBC because it only uses the reference LR frame to reconstruct the HR image. For the low SNR Gaussian noise, our method has a slight PSNR and SSIM decrease compared to the TV and RDESR methods. However, for other cases, our proposed method produces the reconstructed HR images with the highest PSNR and SSIM values. In particular, an obvious performance advantage can be observed for the noise case with impulse noise. To see the visual effect of different methods, for example, HR restorations are shown in Figure 4 for the SNR = 35 dB Gaussian noise, in Figure 5 for r = 0.01 salt-and-pepper noise, and in Figure 6 for 35+0.001 mixed noise. In Figure 4, our method provides comparable visual performance with the TV, NS, and RDESR methods. In Figure 5, we can see that our method is efficient to remove salt-and-pepper noise, the RSR produces obvious ringing artifacts near the edges, and the SAR, TV, and RDESR methods are very sensitive to outliers because the l2 norm-based observation model is not robust against outliers. The NS and RDESR methods can moderately suppress the outliers, but the details are also oversmoothed. Our method provides the most visually appealing results without artifacts while still preserve edges well. In Figure 6, our method has better blur and ringing artifacts removing performance than RSR, only small fractional outliers can severely degrade the performance of SAR, TV, NS, and RDESR, and the specks are obvious in the reconstruction HR images. Our method is efficient, stable, and robust with various noise cases.Table 3
PSNR (dB) results with the accurate motion vector.
Noise casePSNRBBCRSRSARTVNSRDESRProposedGaussian 15 dB16.0823.5022.5223.5323.3423.4223.17Gaussian 35 dB16.5126.5327.0530.9729.4230.1030.99Salt-and-pepper 0.0115.8125.3922.2823.3522.4122.5031.64Salt-and-pepper 0.113.4023.9217.8717.9817.3618.1025.65Mixed 35 + 0.00116.4726.5025.6128.3026.6428.2830.94Mixed 30 + 0.0515.3925.4320.8221.4421.0721.6027.88Table 4
SSIM results with the accurate motion vector.
Noise casePSNRBBCRSRSARTVNSRDESRProposedGaussian 15 dB0.83330.97500.96830.97520.97430.97500.9726Gaussian 35 dB0.84910.98710.98920.99570.99380.99440.9957Salt-and-pepper 0.010.82520.98370.96640.97390.96820.97010.9963Salt-and-pepper 0.10.70750.97660.89840.89700.87590.90230.9847Mixed 35 + 0.0010.84780.98750.98470.99180.98810.99070.9956Mixed 30 + 0.050.80440.98350.95100.95740.95440.95890.9910Figure 4
Super-resolution on the EIA with PSNR and SSIM in the brackets for 35 dB Gaussian noise when the motion information is exact. Results of (a) BBC (16.52 dB, 0.8489), (b) RSR (26.81 dB, 0.9885), (c) SAR (27.04 dB, 0.9891), (d) TV (30.96 dB, 0.9956), (e) NS (29.52 dB, 0.9940), (f) RDESR (30.12 dB, 0.9953), and (g) proposed (30.98 dB, 0.9957).
(a)(b)(c)(d)(e)(f)(g)Figure 5
Super-resolution on the EIA with PSNR and SSIM in the brackets for 1% salt-and-pepper noise when the motion information is exact. Results of (a) BBC (16.18 dB, 0.8378), (b) RSR (25.55 dB, 0.9846), (c) SAR (22.17 dB, 0.9655), (d) TV (23.22 dB, 0.9730), (e) NS (22.21 dB, 0.9665), (f) RDESR (22.31 dB, 0.9689), and (g) proposed (31.55 dB, 0.9962).
(a)(b)(c)(d)(e)(f)(g)Figure 6
Super-resolution on the EIA with PSNR and SSIM in the brackets for 35 + 0.001 mixed noise when the motion information is exact. Results of (a) BBC (16.47 dB, 0.8475), (b) RSR (26.53 dB, 0.9878), (c) SAR (25.37 dB, 0.9839), (d) TV (27.66 dB, 0.9906), (e) NS (26.26 dB, 0.9872), (f) RDESR (26.33 dB, 0.9878), and (g) proposed (30.91 dB, 0.9956).
(a)(b)(c)(d)(e)(f)(g)


### 5.2.3. Effect of Motion Vector Uncertainty

In this subsection, we conduct experiments to investigate the influence of inaccurate motion vector on the SR methods. We utilized the same setup and dataset as 5.2.2. To simulate the errors in motion estimation, we added a deliberate error in its translation parameters by i.i.d. zero-mean white Gaussian noise with standard deviation equal to 0.2 pixels, and the rotation parameters with noise uniformly distributed in−0..1∘,0.1∘. The mean PSNR and SSIM results for different SR methods are shown in Table 5 and 6, respectively. As expected, the performance of all SR algorithms decreases due to the motion error. The motion error has a slight effect on the variational Bayesian methods. However, the performance degradation for RSR and RDESR is severe (about 2–3 dB PSNR loss). Similarly, for the low SNR Gaussian noise, the TV method has a slight performance advantage (about 0.4 dB for PSNR). For the noise case with impulse noise, the proposed method improves the PSNR by at least 3–4 dB. Example HR reconstructions are shown in Figure 7 for 35+0.001 mixed noise. Compared with Figure 6, similar visual results are obtained for the variational Bayesian method, but the RSR has obvious ringing artifacts near the edge, and the RDESR method gets blurred edges and textures.Table 5
PSNR (dB) results with the inaccurate motion vector.
Noise casePSNRBBCRSRSARTVNSRDESRProposedGaussian 15 dB16.0522.4922.5123.2922.9222.5522.88Gaussian 35 dB16.5123.1627.0530.4627.3426.4130.49Salt-and-pepper 0.0116.1122.9522.2623.0521.9821.1731.03Salt-and-pepper 0.115.3921.7517.8617.9017.2417.0125.26Mixed 35 + 0.00116.4723.3325.6027.4825.7523.5430.57Mixed 30 + 0.0517.8623.2320.8121.2620.7119.3027.73Table 6
SSIM results with the inaccurate motion vector.
Noise casePSNRBBCRSRSARTVNSRDESRProposedGaussian 15 dB0.83330.96790.96830.97370.97170.97310.9709Gaussian 35 dB0.84910.97020.98910.99510.98980.98490.9952Salt-and-pepper 0.010.82520.97010.96640.97190.96490.96370.9957Salt-and-pepper 0.10.70750.95960.89810.89520.87260.86960.9832Mixed 35 + 0.0010.84780.97310.98470.99000.98540.97460.9952Mixed 30 + 0.050.80440.97300.95100.95550.95050.94830.9907Figure 7
Super-resolution on the EIA with PSNR and SSIM in the brackets for 35 + 0.001 mixed noise when the motion information is inexact. Results of (a) BBC (16.47 dB, 0.8475), (b) RSR (23.65 dB, 0.9753), (c) SAR (25.37 dB, 0.9839), (d) TV (27.63 dB, 0.9905), (e) NS (26.07 dB, 0.9866), (f) RDESR (23.58 dB, 0.9749), and (g) proposed (30.86 dB, 0.9955).
(a)(b)(c)(d)(e)(f)(g)


### 5.2.4. Effect of the Number of Frames with Outliers

In this subsection, we study the effect of the number of frames with outliers on the performance of the SR methods. We generated 17 LR images from monarch by a downsampling factor of 4 and a3×3 Gaussian blurring kernel with a standard deviation equal to 1. All frames were corrupted by SNR dB Gaussian noise, and the first O2≤O≤17 frames were also corrupted by random-valued impulse noise with corruption rate r (denoted by SNR + r(O)). The SR reconstruction was performed with an exact motion vector. The mean PSNR and SSIM results for 25 + 0.05(O) mixed noise are shown in Figure 8. We observe that RSR are robust to outliers but with low PSNR and SSIM values, while SAR, TV, NS, and RDESR are not robust against outliers and the PSNR and SSIM values rapidly decrease when the number of frames with outliers increases. Our method has the highest PSNR and SSIM values and is robust to outliers. Example HR reconstructions are shown in Figure 9 for 25 + 0.05(9). Images in Figure 10 are the cropped regions corresponding to the rectangular area in Figure 9. It is observed that the proposed method has an advantage in achieving a balance between removing noise and preserving the image edges.Figure 8
PSNR and SSIM on the monarch dataset for 25 + 0.05(O).
(a)(b)Figure 9
Super-resolution on the monarch with PSNR and SSIM in the brackets for 25 + 0.05(9) mixed noise. (a) Original HR image, (b) one LR image, and results of (c) BBC (16.46 dB, 0.8147), (d) RSR (24.24 dB, 0.9714), (e) SAR (26.71 dB, 0.9841), (f) TV (27.62 dB, 0.9872), (g) NS (28.65 dB, 0.9900), (h) RDESR (28.72, 0.9908), and (i) proposed (29.73 dB, 0.9922).
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 10
The corresponding cropped images of Figure9.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


### 5.2.5. Effect of the Sequence Length

In this subsection, we evaluate the influence of the number of LR images on the quality of SR reconstruction. We generatedL5≤L≤19 LR images from the parrot by a downsampling factor of 4 and a 3×3 Gaussian blurring kernel with a standard deviation equal to 1. All frames were corrupted by SNR dB Gaussian noise and then by the random-valued impulse noise with corruption rate r. The SR restoration was performed with an exact motion vector. The PSNR and SSIM results for 30+0.01 mixed noise are shown in Figure 11. The performance of all algorithms increases as the number of LR images increases, as expected. The proposed method clearly outperforms the other methods independent of the number of LR images. Example HR reconstructions are shown in Figure 12 for 30+0.01 with L=9. Notice that this is an underdetermined case for the SR problem (L<16). Figure 13 presents the cropped regions corresponding to the rectangular area in Figure 12. Our proposed method produces the most visually appealing results and is very efficient in removing outliers and preserving detailed texture.Figure 11
Effect of number of LR images to the PSNR and SSIM on the parrot dataset.
(a)(b)Figure 12
Super-resolution on the parrot with PSNR and SSIM in the brackets for 30 + 0.01 mixed noise. (a) Original HR image, (b) one LR image, Results of (c) BBC (22.02 dB, 0.9423), (d) RSR (28.57 dB), (e) SAR (27.41 dB, 0.9835), (f) TV (27.99 dB, 0.9856), (g) NS (27.97 dB, 0.9856), and (h) proposed (31.51 dB, 0.9937).
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 13
The corresponding cropped images of Figure12.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 5.2.1. Effect of Image Prior Model

Before comparing with other SR methods, we study the effect of the image prior model for our proposed observation model. Here, we compared three typical prior models, including SAR using equations (25) and (26), TV using equations (27)-(29), and NS using equations (30)-(33), and denoted by La + SAR, La + TV, and La + NS, respectively. We evaluate the SR reconstruction on EIA and Lena for different noise models in cases where motion vector is available. According to equation (1), five LR images were generated by shifting the HR image with 0,0T, 0,0.5T, 0.5,0T, 1,0T, and 0,1T pixels, rotating it with 0∘, 3∘, −3∘, 5∘, and −5∘, blurring it with 3×3 Gaussian blur kernel with a standard deviation equal to 1, downsampling it by factor of 2, and corrupting it with Gaussian noise and salt-and-pepper noise. The PSNR and SSIM results for the different prior models are shown in Table 1 and 2, respectively. For the low SNR Gaussian noise, La + TV has a slight PSNR and SSIM decrease compared to La + NS. However, for other cases, La + TV has the highest PSNR and SSIM values and improves the PSNR and SSIM by at least 1.5 dB and 0.014, respectively. To see the visual effect of the different prior models, for example, HR restorations on Lena are shown in Figure 3 for the SNR = 35 dB Gaussian noise and r = 1% salt-and-pepper noise. From the red circle inserted in Figure 3, La + TV and La + NS have crisper edges compared with La + SAR; La + NS has better noise removing performance and favors a smooth solution, but some textures are also oversmoothed. The above PSNR, SSIM, and the virtual performance demonstrate that La + TV is more favorable, and we choose it as the image prior model for our proposed SR method.Table 1
PSNR (dB) results for different prior models.
Noise caseEIALenaLa + SARLa + TVLa + NSLa + SARLa + TVLa + NSGaussian 15 dB22.2123.1923.3128.3529.2529.55Gaussian 35 dB26.6731.0029.4733.0335.8133.68Salt-and-pepper 0.123.1825.6523.3330.8732.4030.48Salt-and-pepper 0.0126.0231.6428.8932.3936.0634.25Table 2
SSIM results for different prior models.
Noise caseEIALenaLa + SARLa + TVLa + NSLa + SARLa + TVLa + NSGaussian 15 dB0.96570.97270.97280.98420.98750.9880Gaussian 35 dB0.98810.99570.99390.99470.99720.9954Salt-and-pepper 0.10.97250.98470.97300.99120.99380.9902Salt-and-pepper 0.010.98610.99630.99300.99380.99740.9960Figure 3
Super-resolution on the Lena with PSNR and SSIM in the brackets: for 35 dB Gaussian noise, results of (a) La + SAR (33.02 dB, 0.9947), (b) La + TV (35.78 dB, 0.9972), and (c) La + NS (33.66 dB, 0.9954); for 1% salt-and-pepper noise, results of (d) La + SAR (32.34 dB, 0.9937), (e) La + TV (36.04 dB, 0.9973), and (f) La + NS (34.27 dB, 0.9960).
(a)(b)(c)(d)(e)(f)


## 5.2.2. Robust to Noise Model

In this subsection, we utilized the same setup as the previous subsection to evaluate the robustness of the proposed method (La + TV in the last subsection) on simulated images, using the dataset EIA. Besides the Gaussian noise and salt-and-pepper noise, the mixed noise model is also used to test the robustness of the SR method. For the mixed noise, it is expressed as “SNR +r,” which means that the LR images are first corrupted by SNR dB Gaussian noise and then corrupted by the random-valued impulse noise with a corruption rate r. The PSNR and SSIM results for different SR methods are shown in Table 3 and Table 4, respectively. As expected, all SR algorithms result in better performance than BBC because it only uses the reference LR frame to reconstruct the HR image. For the low SNR Gaussian noise, our method has a slight PSNR and SSIM decrease compared to the TV and RDESR methods. However, for other cases, our proposed method produces the reconstructed HR images with the highest PSNR and SSIM values. In particular, an obvious performance advantage can be observed for the noise case with impulse noise. To see the visual effect of different methods, for example, HR restorations are shown in Figure 4 for the SNR = 35 dB Gaussian noise, in Figure 5 for r = 0.01 salt-and-pepper noise, and in Figure 6 for 35+0.001 mixed noise. In Figure 4, our method provides comparable visual performance with the TV, NS, and RDESR methods. In Figure 5, we can see that our method is efficient to remove salt-and-pepper noise, the RSR produces obvious ringing artifacts near the edges, and the SAR, TV, and RDESR methods are very sensitive to outliers because the l2 norm-based observation model is not robust against outliers. The NS and RDESR methods can moderately suppress the outliers, but the details are also oversmoothed. Our method provides the most visually appealing results without artifacts while still preserve edges well. In Figure 6, our method has better blur and ringing artifacts removing performance than RSR, only small fractional outliers can severely degrade the performance of SAR, TV, NS, and RDESR, and the specks are obvious in the reconstruction HR images. Our method is efficient, stable, and robust with various noise cases.Table 3
PSNR (dB) results with the accurate motion vector.
Noise casePSNRBBCRSRSARTVNSRDESRProposedGaussian 15 dB16.0823.5022.5223.5323.3423.4223.17Gaussian 35 dB16.5126.5327.0530.9729.4230.1030.99Salt-and-pepper 0.0115.8125.3922.2823.3522.4122.5031.64Salt-and-pepper 0.113.4023.9217.8717.9817.3618.1025.65Mixed 35 + 0.00116.4726.5025.6128.3026.6428.2830.94Mixed 30 + 0.0515.3925.4320.8221.4421.0721.6027.88Table 4
SSIM results with the accurate motion vector.
Noise casePSNRBBCRSRSARTVNSRDESRProposedGaussian 15 dB0.83330.97500.96830.97520.97430.97500.9726Gaussian 35 dB0.84910.98710.98920.99570.99380.99440.9957Salt-and-pepper 0.010.82520.98370.96640.97390.96820.97010.9963Salt-and-pepper 0.10.70750.97660.89840.89700.87590.90230.9847Mixed 35 + 0.0010.84780.98750.98470.99180.98810.99070.9956Mixed 30 + 0.050.80440.98350.95100.95740.95440.95890.9910Figure 4
Super-resolution on the EIA with PSNR and SSIM in the brackets for 35 dB Gaussian noise when the motion information is exact. Results of (a) BBC (16.52 dB, 0.8489), (b) RSR (26.81 dB, 0.9885), (c) SAR (27.04 dB, 0.9891), (d) TV (30.96 dB, 0.9956), (e) NS (29.52 dB, 0.9940), (f) RDESR (30.12 dB, 0.9953), and (g) proposed (30.98 dB, 0.9957).
(a)(b)(c)(d)(e)(f)(g)Figure 5
Super-resolution on the EIA with PSNR and SSIM in the brackets for 1% salt-and-pepper noise when the motion information is exact. Results of (a) BBC (16.18 dB, 0.8378), (b) RSR (25.55 dB, 0.9846), (c) SAR (22.17 dB, 0.9655), (d) TV (23.22 dB, 0.9730), (e) NS (22.21 dB, 0.9665), (f) RDESR (22.31 dB, 0.9689), and (g) proposed (31.55 dB, 0.9962).
(a)(b)(c)(d)(e)(f)(g)Figure 6
Super-resolution on the EIA with PSNR and SSIM in the brackets for 35 + 0.001 mixed noise when the motion information is exact. Results of (a) BBC (16.47 dB, 0.8475), (b) RSR (26.53 dB, 0.9878), (c) SAR (25.37 dB, 0.9839), (d) TV (27.66 dB, 0.9906), (e) NS (26.26 dB, 0.9872), (f) RDESR (26.33 dB, 0.9878), and (g) proposed (30.91 dB, 0.9956).
(a)(b)(c)(d)(e)(f)(g)


## 5.2.3. Effect of Motion Vector Uncertainty

In this subsection, we conduct experiments to investigate the influence of inaccurate motion vector on the SR methods. We utilized the same setup and dataset as 5.2.2. To simulate the errors in motion estimation, we added a deliberate error in its translation parameters by i.i.d. zero-mean white Gaussian noise with standard deviation equal to 0.2 pixels, and the rotation parameters with noise uniformly distributed in−0..1∘,0.1∘. The mean PSNR and SSIM results for different SR methods are shown in Table 5 and 6, respectively. As expected, the performance of all SR algorithms decreases due to the motion error. The motion error has a slight effect on the variational Bayesian methods. However, the performance degradation for RSR and RDESR is severe (about 2–3 dB PSNR loss). Similarly, for the low SNR Gaussian noise, the TV method has a slight performance advantage (about 0.4 dB for PSNR). For the noise case with impulse noise, the proposed method improves the PSNR by at least 3–4 dB. Example HR reconstructions are shown in Figure 7 for 35+0.001 mixed noise. Compared with Figure 6, similar visual results are obtained for the variational Bayesian method, but the RSR has obvious ringing artifacts near the edge, and the RDESR method gets blurred edges and textures.Table 5
PSNR (dB) results with the inaccurate motion vector.
Noise casePSNRBBCRSRSARTVNSRDESRProposedGaussian 15 dB16.0522.4922.5123.2922.9222.5522.88Gaussian 35 dB16.5123.1627.0530.4627.3426.4130.49Salt-and-pepper 0.0116.1122.9522.2623.0521.9821.1731.03Salt-and-pepper 0.115.3921.7517.8617.9017.2417.0125.26Mixed 35 + 0.00116.4723.3325.6027.4825.7523.5430.57Mixed 30 + 0.0517.8623.2320.8121.2620.7119.3027.73Table 6
SSIM results with the inaccurate motion vector.
Noise casePSNRBBCRSRSARTVNSRDESRProposedGaussian 15 dB0.83330.96790.96830.97370.97170.97310.9709Gaussian 35 dB0.84910.97020.98910.99510.98980.98490.9952Salt-and-pepper 0.010.82520.97010.96640.97190.96490.96370.9957Salt-and-pepper 0.10.70750.95960.89810.89520.87260.86960.9832Mixed 35 + 0.0010.84780.97310.98470.99000.98540.97460.9952Mixed 30 + 0.050.80440.97300.95100.95550.95050.94830.9907Figure 7
Super-resolution on the EIA with PSNR and SSIM in the brackets for 35 + 0.001 mixed noise when the motion information is inexact. Results of (a) BBC (16.47 dB, 0.8475), (b) RSR (23.65 dB, 0.9753), (c) SAR (25.37 dB, 0.9839), (d) TV (27.63 dB, 0.9905), (e) NS (26.07 dB, 0.9866), (f) RDESR (23.58 dB, 0.9749), and (g) proposed (30.86 dB, 0.9955).
(a)(b)(c)(d)(e)(f)(g)


## 5.2.4. Effect of the Number of Frames with Outliers

In this subsection, we study the effect of the number of frames with outliers on the performance of the SR methods. We generated 17 LR images from monarch by a downsampling factor of 4 and a3×3 Gaussian blurring kernel with a standard deviation equal to 1. All frames were corrupted by SNR dB Gaussian noise, and the first O2≤O≤17 frames were also corrupted by random-valued impulse noise with corruption rate r (denoted by SNR + r(O)). The SR reconstruction was performed with an exact motion vector. The mean PSNR and SSIM results for 25 + 0.05(O) mixed noise are shown in Figure 8. We observe that RSR are robust to outliers but with low PSNR and SSIM values, while SAR, TV, NS, and RDESR are not robust against outliers and the PSNR and SSIM values rapidly decrease when the number of frames with outliers increases. Our method has the highest PSNR and SSIM values and is robust to outliers. Example HR reconstructions are shown in Figure 9 for 25 + 0.05(9). Images in Figure 10 are the cropped regions corresponding to the rectangular area in Figure 9. It is observed that the proposed method has an advantage in achieving a balance between removing noise and preserving the image edges.Figure 8
PSNR and SSIM on the monarch dataset for 25 + 0.05(O).
(a)(b)Figure 9
Super-resolution on the monarch with PSNR and SSIM in the brackets for 25 + 0.05(9) mixed noise. (a) Original HR image, (b) one LR image, and results of (c) BBC (16.46 dB, 0.8147), (d) RSR (24.24 dB, 0.9714), (e) SAR (26.71 dB, 0.9841), (f) TV (27.62 dB, 0.9872), (g) NS (28.65 dB, 0.9900), (h) RDESR (28.72, 0.9908), and (i) proposed (29.73 dB, 0.9922).
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 10
The corresponding cropped images of Figure9.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 5.2.5. Effect of the Sequence Length

In this subsection, we evaluate the influence of the number of LR images on the quality of SR reconstruction. We generatedL5≤L≤19 LR images from the parrot by a downsampling factor of 4 and a 3×3 Gaussian blurring kernel with a standard deviation equal to 1. All frames were corrupted by SNR dB Gaussian noise and then by the random-valued impulse noise with corruption rate r. The SR restoration was performed with an exact motion vector. The PSNR and SSIM results for 30+0.01 mixed noise are shown in Figure 11. The performance of all algorithms increases as the number of LR images increases, as expected. The proposed method clearly outperforms the other methods independent of the number of LR images. Example HR reconstructions are shown in Figure 12 for 30+0.01 with L=9. Notice that this is an underdetermined case for the SR problem (L<16). Figure 13 presents the cropped regions corresponding to the rectangular area in Figure 12. Our proposed method produces the most visually appealing results and is very efficient in removing outliers and preserving detailed texture.Figure 11
Effect of number of LR images to the PSNR and SSIM on the parrot dataset.
(a)(b)Figure 12
Super-resolution on the parrot with PSNR and SSIM in the brackets for 30 + 0.01 mixed noise. (a) Original HR image, (b) one LR image, Results of (c) BBC (22.02 dB, 0.9423), (d) RSR (28.57 dB), (e) SAR (27.41 dB, 0.9835), (f) TV (27.99 dB, 0.9856), (g) NS (27.97 dB, 0.9856), and (h) proposed (31.51 dB, 0.9937).
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 13
The corresponding cropped images of Figure12.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 5.3. Experiments with Real Images

In this subsection, we present the performance comparison of SR methods on the challenging real image sequence called alpaca, which was downloaded from the MDSP benchmark dataset [49]. The sequence contains 55 LR frames (32×70), and 9 frames (frame 45 to 53) are used in our experiment. The 9 frames follow a complicated motion model and contain a high level of outliers. We approximated the PSF using an isotropic 3×3 Gaussian kernel with a standard deviation equal to 1.Since no ground truth image is available for the real sequence, we empirically chose the SR reconstructions such that they produce the visually most appealing results in the iteration process for the variational Bayesian method. For RSR and RDESR, the regularization parameters were manually tuned to produce the best performance. Figure14 shows the SR results with magnification factor P=3 using the 9 LR frames. Figure 15 shows the underdetermined case. It is observed that our proposed method produces the most visually appealing results. The SAR and TV methods suffer from visible ghosting artifacts caused by the alpaca’s movements. The NS and RDESR prefer a smooth solution, but important edges and texture are also oversmoothed, especially for the underdetermined case. Compared to the robust RSR method, our proposed method has better noise removing performance with crisper details, and most importantly, parameter tuning by hand is no longer necessary.Figure 14
Super-resolution on the 9 LR frames with magnification factor 3. Frame 49 is set as the reference frame. Results of (a) BBC, (b) RSR, (c) SAR, (d) TV, (e) NS, (f) RDESR, and (g) proposed.
(a)(b)(c)(d)(e)(f)(g)Figure 15
Super-resolution on the five LR frames (frame 46 to 50) with magnification factor 3. Frame 48 is set as the reference frame. Results of (a) BBC, (b) RSR, (c) SAR, (d) TV, (e) NS, (f) RDESR, and (g) proposed.
(a)(b)(c)(d)(e)(f)(g)


## 6. Conclusions

In this study, a robust and fully data-driven SR reconstruction method has been proposed in the hierarchical Bayesian framework, using a l1 norm-based observation model. After evaluating three typical image prior models, we chose the TV prior model for our proposed method. The local noise adaptivity of our method makes it more robust to impulse noise and outliers. Compared with other competitive methods, experiments on simulated show that our method improves the PSNR by at least 2–4 dB for the noise case with impulse noise, and experiments on simulated and real images demonstrate that our method has a better trade-off between preserving detailed image information and removing noise and outliers.



---

*Source: 1497107-2022-03-23.xml*