Multiframe super-resolution based on half-quadratic prior with artifacts suppress

J. Vis. Commun. Image R. 40 (2016) 656–670

Contents lists available at ScienceDirect

J. Vis. Commun. Image R. journal homepage: www.elsevier.com/locate/jvci

Multiframe super-resolution based on half-quadratic prior with artifacts suppress Renchao Jin a, Shengrong Zhao b,a, Xiangyang Xu a, Enmin Song a,⇑ a b

School of Computer Science and Technology, Huazhong University of Science and Technology, Wuhan 430074, China School of Information, Qilu University of Technology, Jinan 250353, China

a r t i c l e

i n f o

Article history: Received 6 January 2016 Revised 8 August 2016 Accepted 10 August 2016 Available online 11 August 2016 Keywords: Super-resolution Variational Bayesian inference Regularization Discontinuity-preserving Artifacts-suppressing

a b s t r a c t The multiframe super-resolution (SR) technique aims to obtain a high-resolution (HR) image by using a set of observed low-resolution (LR) images. In the reconstruction process, artifacts may be possibly produced due to the noise, especially in presence of stronger noise. In order to suppress artifacts while preserving discontinuities of images, in this paper a multiframe SR method is proposed by involving the reconstruction properties of the half-quadratic prior model together with the quadratic prior model using a convex combination. Moreover, by analyzing local features of the underlined HR image, these two prior models are combined by using an automatically calculated weight function, making both smooth and discontinuous pixels handled properly. A variational Bayesian inference (VBF) based algorithm is designed to efficiently and effectively seek the solution of the proposed method. With the VBF framework, motion parameters and hyper-parameters are all determined automatically, leading to an unsupervised SR method. The efficiency of the hybrid prior model is demonstrated theoretically and practically, which shows that our SR method can obtain better results from LR images even with stronger noise. Extensive experiments on several visual data have demonstrated the efficacy and superior performance of the proposed algorithm, which can not only preserve image details but also suppress artifacts. Ó 2016 Elsevier Inc. All rights reserved.

1. Introduction Super-resolution (SR) [1] is an important branch in image fusion technology, which aims to reconstruct a high-resolution (HR) image from a single or a set of low-resolution (LR) images. The SR concept was first proposed by Tsai and Huang in [2]. This technique has been widely used in many applications, such as computer vision, medical imaging, public safety, and military reconnaissance. In some applications, only one LR image is available. Inspired by this, the single-frame SR method attracts increasing research attention. The simplest methods are the interpolation-based approaches, such as the bicubic interpolation. However, these methods often lead to blurring the edges and introducing artifacts in the results. More advanced single-frame SR methods are proposed. In [3], an iterative multiscale semilocal interpolation method is proposed to generate the HR image by using a maximum a posteriori (MAP) estimation. In this estimation, the top texturerelevant LR pixels are used to construct a data fidelity term, and a bilateral total variation is used as the regularization term. ⇑ Corresponding author. E-mail address: [email protected] (E. Song). http://dx.doi.org/10.1016/j.jvcir.2016.08.006 1047-3203/Ó 2016 Elsevier Inc. All rights reserved.

Learning-based SR methods are popular in many applications recently. In [4], a statistical learning method for SR with both global and local constraints is proposed. The global parametric constraint guarantees the super-resolved global image to agree with the sparse property of natural images, and the local nonparametric constraint is used to infer the residues between the image derived from the global constraint and the HR image. A statistical prediction model based on sparse representations of LR and HR image patches for single image SR is proposed in [5], which goes beyond the standard assumption of sparse representation invariance over a low and high resolution dictionary pair. Inference with their model leads to a low-complexity scale-up scheme that has the useful interpretation of a feedforward neural network. In [6], a self-similarity based single-frame SR is proposed. The internal patch search space is expanded by allowing geometric variations, and the additional affine transformations is incorporated to accommodate local shape variations. In [7], a joint SR and smoothing framework is proposed with an L0 gradient smooth constraint. A robust coupled dictionary learning method with locality coordinate constraints is introduced to reconstruct the HR depth map. And an adaptively regularized shock filter is incorporated to simultaneously reduce the jagged noise and sharpen the edges. Deep learning techniques have been successfully applied in image SR

R. Jin et al. / J. Vis. Commun. Image R. 40 (2016) 656–670

657

Fig. 1. The proposed framework.

problems. In [8], the authors argue that domain expertise represented by the conventional sparse coding model is still valuable, and it can be combined with the key ingredients of deep learning to achieve further improved results. They show that a sparse coding model particularly designed for super-resolution can be incarnated as a neural network, and trained in a cascaded structure from end to end. The method in [9] learns an end-to-end mapping between low- and high-resolution images, with little extra pre/ post-processing beyond the optimization. This method also shows that conventional sparse-coding-based image SR methods can be reformulated into a deep convolutional neural network. The multi-frame SR method requires a set of LR image as the input. And the regularization approach has become more popular due to its perfect mathematical background, and it can introduce a prior knowledge concerning the solution into the reconstruction process. Such as, the Tikhonov prior model proposed in [10] is quadratic, which assumes global smoothness, hence it often blurs image edges. On the other hand, the total variation (TV) type [11–13] model could preserve edges due to their non-quadratic form. However, this form has low penalty for the noise gradients, thus the noise would not be suppressed effectively and lead to artifacts [14]. The widely used markov random field (MRF) model [15] includes the Gaussian MRF model [16] and the Huber MRF model [17]. As in the Tikhonov model, the Gaussian MRF model is also quadratic, which may blur image edges. There exists a threshold in the Huber model, which can separate the quadratic and linear regions [18]. However, it is a difficult matter to choose a suitable threshold in actual experiment condition. Thus, the concept of discontinuity adaptive MRF model or edge-preserving regularization term [15,19,20] is becoming more adopted due to the usage of half-quadratic potential functions. In [21], the authors proposed a deterministic SR approach based on the bilateral edge-preserving (BEP) model and using such a potential function (called Yang potential function here). In [22], the traditional TV model is combined with the Frobenius norm regularization. And this method was proposed to preserve edges while avoiding staircase effects in the homogeneous regions of images without considering the motion estimation. In [23], Villena et al. considered the simultaneous auto regressive (SAR) model with quadratic form involving the TV type based SR method, which has been proved effective for the SR problem. However, the weighting parameter used to balance these two prior models has to be adjusted manually, which is a difficult matter in actual experiment condition. In [24], a prior model based on the Yang potential function was proposed, which was called the anisotropic adaptive (AA) prior model here. The aforementioned prior models are all half-quadratic, hence, they may

not be possible to suppress noise, especially in presence of stronger noise. The regularization approach includes the deterministic and stochastic regularization approaches. The basic idea of the deterministic regularization approach is to use the regularization strategy to incorporate the prior knowledge of the unknown HR image. The stochastic regularization approach, typically a Bayesian approach, provides a flexible and convenient way to model a prior knowledge concerning the solution, such as the MAP approach. In [25], the HR image is reconstructed by minimizing a MAP-MRF energy function with the linear truncated prior. And the method in [26] employs the MAP estimator and the sparse image priors to keep useful salient structures while suppressing noise. The estimation of hyper-parameters related to the noise model and the prior image model can also affect the reconstruction quality. The joint posteriori probability distribution of the HR image, motion parameters and hyper-parameters has a complex expression and the MAP can not solve this issue. Recently, the variational Bayesian inference (VBF) estimator has been introduced into the SR reconstruction technique [22,24,23,27,28], which approximates this complex joint posteriori probability distribution analytically by a tractable distribution. In this paper, we extend the AA prior model using an extra quadratic function, and propose a variational Bayesian SR algorithm with a hybrid prior model using a convex combination scheme for HR image estimation. This approach efficiently interlaces the HR image estimation, the motion parameter estimation and the hyper-parameter estimation in a complete framework (Fig. 1). The numerical results show that our method achieves better reconstruction compared with the SR method under comparison. The contributions of this work are summarized as follows: (1) Involve the reconstruction properties of the AA prior model together with the Tikhonov prior model using a convex combination. (2) By analyzing local features of the underlined HR image, the two prior models are combined by using an automatically calculated weighting function, making both smooth and discontinuous pixels handled properly. (3) To demonstrate the efficiency of the combination strategy, the theoretical and practical analysis is presented in this paper. (4) The variational Bayesian inference (VBF) estimator has been used in this paper, and hence our approach is fully formulated in the framework of Bayesian statistics by the utilization of the Kullback-Leibler (KL) divergence. With VBF, the

658

R. Jin et al. / J. Vis. Commun. Image R. 40 (2016) 656–670

motion parameters and the hyper-parameters related to the composite prior and noise statistics are all determined automatically, leading to an unsupervised SR algorithm. The rest of this paper is organized as follows: the hierarchical Bayesian model is described briefly in Section 2. In Section 3, a hybrid prior model is proposed and a local-spatial-adaptive selection method is introduced to establish the new prior model. Also, the solving process is presented briefly in Section 3. In Section 4, the theory analysis of the proposed method is provided. The performance of the proposed SR method by several numerical experiments is demonstrated in Section 5. Finally, the conclusion is drawn in Section 6. And an appendix is presented at the end of the paper showing the solving process in detail. 2. The hierarchical Bayesian model In this paper, a hierarchical Bayesian framework consisting of two stages is adopted. The first stage is used to model the acquisition process, the unknown HR image and the motion parameters, and the second stage is used to model the hyper-parameters. 2.1. The observation model Our task is to estimate the HR image u 2 RM1 M2 from a set of

observed LR images v = fv k 2 RN1 N2 ; k ¼ 1; 2; . . . ; Lg. The total number of pixels in the LR image and HR image is N ¼ N 1  N 2 and M ¼ M1  M 2 , respectively. In this work, images u and v k are rearranged as M  1 and N  1 vectors, respectively. The HR image is geometrically warped, blurred, downsampled, and corrupted by noise to generate an LR image. Hence, the observation model can be defined mathematically as follows,

v k ¼ Dk Hk Cðsk Þu þ nk ; NM

where Dk 2 R

k ¼ 1; 2; . . . ; L ¼ Bðsk Þu þ nk ;

ð1Þ MM

is the downsampling matrix; Hk 2 R

is the

blurring matrix; Cðsk Þ 2 RMM is the warping matrix, where sk denotes the motion vector. The rotation and translation are consid-

The AA model: pðujaÞ / ða1 a2 Þ 4 ( M

 exp 

h r i N k pðv k ju; sk ; rk Þ / rk2 exp  jjv k  Bðsk Þujj22 : 2

ð2Þ

Since noise among the LR images is assumed to be mutually independent, thus,

pðv ju; fsk g; frk gÞ ¼

L Y

pðv k ju; sk ; rk Þ:

ð3Þ

k¼1

2.2. The image prior distribution In the Bayesian method, prior knowledge about the underlined HR image is represented by the prior image distribution. The estimation distinguishes between possible solutions by using different prior image models [29], and the following prior models are used in this paper. The Tikhonov model:

n c o M pðujcÞ / c 2 exp  jjrujj22 ; 2

ð4Þ

where c is the hyper-parameter of this model, r ¼ ðr x ; r y Þ; r x and r y represent the horizontal and vertical gradient operators, respectively.



qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  



qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi

a1 a1i a21i þ rix u 2  a21i þa2 a2i a22i þ riy u 2  a22i

 ;

i

ð5Þ

where for pixel i; r and r represent the horizontal and vertical gradient operators respectively, a1i and a2i are scale parameters for pixel i on horizontal and vertical directions respectively, and a ¼ fa1 ; a2 g; a1 and a2 are the hyper-parameters of this prior distribution. The scale parameter can control the severity of penalizing gradients. With the value of scale parameter increasing, a larger range of errors (sufferring from inaccurate registration and additive Gaussian noise) is accepted and handled using quadratic functions [15]. In theory, gradients produced by edges or corners should be preserved; while gradients produced by errors need to be smoothed. Thus, for the gradients produced by errors, the scale parameter should be set as a large value, otherwise, it should be set as a small value [24]. Based on the analysis, we design three conditions to constrain the scale parameter: (1) its value is determined according to the local image features, (2) the value is inversely proportional to gradient magnitudes, and (3) the value is larger than zero. Thus, a1i and a2i are estimated as follows, x i

y i

a1i ¼

1  2 ; 1 þ rix u

ð6Þ

a2i ¼

1  2 : 1 þ riy u

ð7Þ

For simplicity, denote R ¼ fa1i ; a2i ; i ¼ 1; 2; . . . ; Mg. 2.3. The motion prior distribution For the motion parameters, their corresponding prior distributions are defined as Gaussian distributions,

T

ered in this paper, thus sk ¼ ðhk ; dxk ; dyk Þ . nk 2 RN1 is additive white Gaussian noise. Assume that nk is white Gaussian noise with zero-mean and variance r1 k , and we can obtain,

X

pðfsk gÞ ¼

L

 Y N sk js0k ; R0sk ;

ð8Þ

k¼1

where s0k is the a prior mean vector and R0sk is the a prior covariance matrix. These two parameters can incorporate prior knowledge about the motion parameters into the estimation process. In this work, s0k is obtained by using the Lucas-Kanade method [30], and

R0sk is set as zero matrices. The Lucas-Kanade method is an optical flow algorithm, which owns many advantages: (1) it could more precisely calculate the moving distance in sub-pixel; (2) it is efficient to calculate the complex motions; (3) it is more robust; (4) it is easy to use. For these reasons, in this paper, we use LucasKanade method to calculate the motions as the initials. As the results of the Lucas-Kanade method is very precise, we suppose that the exact value is beside the results of the Lucas-Kanade method. That is the prior knowledge for the range of the motion parameter, which is satisfied the Gaussian distribution. 2.4. The hyperprior distribution Hyper-parameter estimation is also a critical issue, which is approached from a Bayesian perspective in this paper. Generally, it is difficult to directly obtain the prior information about hyperparameters. Therefore, it is usually expressed using the conjugate prior distribution which is calculated conveniently. Moreover, the corresponding posterior distribution has the same functional form with the prior distribution and hence the analytic solution can be

659

R. Jin et al. / J. Vis. Commun. Image R. 40 (2016) 656–670

obtained [27]. It is well known that the inverse Gamma distributions are the conjugate priors for the variance of the Gaussian distribution whose mean value is known. Thus, we assume that the hyper-parameters obey Gamma distributions, i.e.,

 0 pðtÞ ¼ Gamma tja0t ; bt ; t 2 ffrk g; a; cg

ð9Þ 0 bt .

a0t

with the shape parameter and rate parameter The related hyper-parameters can incorporate prior knowledge about the variances of the HR image and noise among the observed LR images into the estimation process. In the following SR method, 0

a0t is set equal to 1 and bt is set equal to 0, which corresponds to utilizing flat prior distributions for the hyper-parameters, in this case, only the observed LR images are responsible for the estimation process. Fig. 2 presents the hierarchical Bayesian model based on the AA prior model and the Tikhonov model, respectively. 3. The variational Bayesian super-resolution

Fig. 3. The trend of ki .

3.1. The objective function In this paper, the VBF estimator is used to obtain the posteriori distribution of the HR image, motion parameters and hyperparameters, given the set of observed LR images, i.e.,

pðHjv Þ ¼

pðv jHÞpðHÞ ; pðv Þ

ð10Þ

where H includes the HR image, motion parameters and hyperparameters. In the VBF estimator, this posteriori distribution is approximated by a tractable distribution qðHÞ, which is obtained by minimizing the Kullback-Leibler (KL) divergence defined as follows,

C KL ðqðHÞjjpðHjv ÞÞ ¼

Z

qðHÞ log

qðHÞ dH: pðHjv Þ

ð11Þ

In order to combine the AA and Tikhonov prior models, the following convex combination is established,

h qðHÞ ¼ arg min kC AA KL ðqðHAA ÞjjpðHAA jv ÞÞ qðHÞ

þ ð1 

i HTikhonov ÞjjpðHTikhonov jv ÞÞ

v ðqð kÞC Tikhono KL

¼ arg minEðqðHÞ; v Þ; qðHÞ

ð12Þ

where HAA ¼ fu; fsk g; frk g; ag; HTikhonov ¼ fu; fsk g; frk g; cg; H ¼ fu; fsk g; frk g; a; cg, and k ¼ ðk1 ; . . . ; kM ÞT is a weight function. 3.2. Calculation of k The weight k is used to control the balance between the Tikhonov prior model and the AA prior model. For pixel i, the AA prior model should be used when jri uj is small or large enough, so that artifacts can be suppressed in the homogeneous region while discontinuities being preserved. Meanwhile, the Tikhonov prior model should be applied when jri uj is moderate, thus, the stronger noise would be smoothed. In this paper, the following function is used to estimate ki ,

(

ki ¼

if jri uj P c;   0:5cos 2cp jri uj þ 0:5; if jri uj < c;

1;

ð13Þ

where c is a given threshold, and the curve of ki is given in Fig. 3. From Fig. 3, we can see that, the Tikhonov prior model is used mainly in the middle domain of jruj such that the stronger noise can be suppressed in the homogeneous regions, avoiding producing artifacts. Moreover, by using the AA prior model simultane-

Fig. 2. The hierarchical Bayesian model.

660

R. Jin et al. / J. Vis. Commun. Image R. 40 (2016) 656–670

ously, the small structure information, such as weak edges and textures, could be avoided being penalized severely by the Tikhonov prior model. The AA prior model is applied when jruj is large enough or equal to zero, and it encourages piecewise smooth reconstructions. In this paper, the pixel values are normalized to lie in ½0; 1, and the threshold parameter c ¼ 0:5 is selected.

Algorithm 1. The iterative SR algorithm.

Input: The observations fv k g; the initial set R0 and vector k0 ; The initial distribution parameters u0 ; R0u ; s0k ; R0sk ; 0

The half-quadratic form of the AA model makes the corresponding objective function difficult to be evaluated. To solve this problem, the majorization-minimization (MM) method was used, and the following lower bound of pðujaÞ in Eq. (5) was found,

8 2 0 1  x 2 > < 2 1 X M r u þ a þ gi 6 B C  a21i A Pðu; a; g 1 ; g 2 Þ ¼ ða1 a2 Þ 4 exp  4a1 @a1i i qffiffiffiffiffi1i > 1 : i¼1 2 gi 0 139  y 2 > = r u þ a22i þ g 2i B C7 qffiffiffiffiffi  a22i A5 ; þa2 @a2i i ð14Þ > ; 2 g2 i

where g 1 and g 2 are M-dimensional vectors with components g 1i and g 2i , respectively. g 1i and g 2i are auxiliary variables which need to be calculated. We will see in Eq. (17), g 1i and g 2i have large values in pixel locations with large gradient, in this case, the functional has low smoothing effect. Consequently, an upper bound to the integral in Eq. (12) is obtained as follows,

Z qðHAA Þ qðHÞ ¼ arg min k qðHAA Þ log dHAA qðHÞ KðHAA ; v ; g 1 ; g 2 Þ  Z qðHTikhonov Þ þð1  kÞ qðHTikhonov Þ log dHTikhonov ; pðHTikhonov ; v Þ

0

0

ð15Þ

The initial auxiliary variables g 10 and g 20 . Output: The mean of the image posterior distribution; 1: repeat 2: Estimate the image posterior distribution using Eq. (16); 3: Estimate the auxiliary variables using Eq. (17); 4: Estimate the scale parameters using Eqs. (6) and (7); 5: Estimate the weighting parameters using Eq. (13); 6: Estimate the motion posterior distribution using Eq. (18); 7: Estimate the hyper-posterior distributions using Eqs. (19), (20), (21), and (22), respectively. 8: until the convergence criterion jjunþ1  un jj22 =jjun jj22 < 105 is met

4. Efficiency analysis In this section, we will first prove that for the function Pðu; a; g 1 ; g 2 Þ, the ability of preserving edges still exists. Then, we will demonstrate the proposed SR method could preserve edges while suppressing artifacts. 4.1. The ability of preserving edges

¼ arg minEðqðHÞ; g 1 ; g 2 ; v Þ;

We

qðHÞ

Q where KðH1 ; v ; g 1 ; g 2 Þ ¼ Pðu; a; g 1 ; g 2 Þpða1 Þpða2 ÞpðcÞ Lk¼1 pðsk Þpðrk Þ. QL Since qðHÞ ¼ qðuÞqða1 Þqða2 ÞqðcÞ k¼1 qðsk Þqðrk Þ, from Eq. (15) we can obtain

qðuÞ ¼ N ðujhui; Ru Þ; 2

ð17Þ

2

and

gi

2

gi

function, and the potential function should satisfy the following three conditions for edge-preserving:

ð16Þ

g 1i ¼ hðrix uÞ þ a21i iqðuÞ ; g 2i ¼ hðriy uÞ þ a22i iqðuÞ ;

2   ðr x uÞ þa2 þg1 / rix u ¼ a1i i pffiffiffi11iffi i  a21i

denote

2   ðr y uÞ þa2 þg2 / riy u ¼ a2i i pffiffiffi22iffi i  a22i . /ðÞ can be viewed as a potential

(i) 2

0

a0rk ; brk ; a0c ; bc ; a0a1 ; ba1 ; a0a2 ; ba2 ;

3.3. Optimization

/0 ðf Þ 2f

continues and strictly decreasing for f P 0; 0

(ii) limf !þ1 /2fðf Þ ¼ 0; 0

qðsk Þ ¼ N ðujhsk i; Rsk Þ;

ð18Þ

qðrk Þ ¼ Gammaðrk jark ; brk Þ;

ð19Þ

qða1 Þ ¼ Gammaða1 jaa1 ; ba1 Þ;

ð20Þ

qða2 Þ ¼ Gammaða2 jaa2 ; ba2 Þ;

ð21Þ

qðcÞ ¼ Gammaðcjac ; bc Þ:

ð22Þ

(iii) limf !0 /2fðf Þ ¼ constant; 0 < constant < 1. For convenience, we denote F 1;i ¼

/0 ðrix uÞ 2rix u

and F 2;i ¼

/0 ðri uÞ y

y

2ri u

. Then,

we can obtain

F 1;i

Updating these probability distributions is via updating the distribution parameters, and these parameters are keys to obtain the analytic expressions of the variables. The specific distributions and analytic expressions of the parameters in these distributions are described in the Appendix. We summarize below the iterative SR algorithm, Algorithm 1. When the convergence criterion is met, estimates for the HR image, the motion parameters, and hyperparameters, i.e., H ¼ fu; fsk g; frk g; a; cg, can be finally given by the corresponding posterior mean estimates: u ¼ EqðuÞ ½u; a1 ¼ Eqða1 Þ ½a1 ; a2 ¼ Eqða2 Þ ½a2 ; c ¼ EqðcÞ ½c; sk ¼ Eqðsk Þ ½sk , and rk ¼ Eqðrk Þ ½rk ; k ¼ 1; 2; . . . ; L.

  /0 rix u a1i ¼ qffiffiffiffiffi ; ¼ 2rix u 2 g1

ð23Þ

  /0 riy u a2i ¼ qffiffiffiffiffi : ¼ 2riy u 2 g2

ð24Þ

i

F 2;i

i

Thus, we can obtain,

lim F 1;i ¼ 0; lim F 1;i ¼ c1; y

0 < c1 < 1;

ð25Þ

lim F 2;i ¼ 0; lim F 2;i ¼ c2; y

0 < c2 < 1;

ð26Þ

rix u!1

rix u!1

ri u!0

ri u!0

where c1 and c2 are constants. The above formulas indicate that Pðu; a; g 1 ; g 2 Þ still has the ability of preserving edges.

661

R. Jin et al. / J. Vis. Commun. Image R. 40 (2016) 656–670

(a) panda

(b) EIAcen

(c) Barbara

(d) pepper

(e) cameraman

(f) tower

(g) house

(h) text

Fig. 4. The images used in simulation experiments.

4.2. Potentials for the AA prior model with artifacts suppress In Algorithm 1, the HR image is estimated by solving the following equation in an iterative fashion

" XXX X hrk iBðhsk iÞT Bðhsk iÞ þ hnkij Oki ðhsk iÞT Okj ðhsk iÞ k

k

i

j

i þ k½ha1 iðr Þ W 1 ðr Þ þ ha2 iðr Þ W 2 ðr y Þ þ ð1  kÞhcirT r u " # X T ð27Þ hrk iBðhsk iÞ v k ; ¼ x T

y T

x

lim Z 1;i ¼ 0;

ð32Þ

rix u!1

k

where W 1 and W 2 are M  M diagonal matrices, in which, the elea1i a2i and , respecments on the diagonal are 2 2 y hðrix uÞ þa21i iqðuÞ hðri uÞ þa22i iqðuÞ tively. For ease of description, we first denote

a1i

Z 1;i ¼  x 2 ; h ri u þ a21i iqðuÞ

ð28Þ

ð29Þ

T

gradient operator; ðr x Þ ðr x Þ is a matrix that represents a second order gradient operator on the horizontal direction; and T

ðr y Þ ðr y Þ is a matrix that represents a second order gradient operT

ator on the vertical direction. Thus, ðrÞT ðrÞ; ðr x Þ W 1 ðr x Þ, and y T

ðr Þ W 2 ðr y Þ are equaling to filters, moreover, the last two filters are weighted by W 1 and W 2 respectively. In Eqs. (28) and (29), Z 1;i and Z 2;i are the diagonal elements of W 1 and W 2 respectively. If there exists weak noise, for the pixel i in homogeneous regions, the difference between the pixel i and the adjacent pixels is small, and ki is almost equal to 1, thus

lim Z 1;i ¼ cost1; lim Z 2;i ¼ const2;

0 < const1 < 1; 0 < const2 < 1;

0 < const2 < 1;

ð33Þ

For the corner, the difference between the adjacent pixels is large, i.e., rix u and riy u are very large, thus,

lim Z 1;i ¼ 0;

ð34Þ

lim Z 2;i ¼ 0:

ð35Þ

riy u!1

In Eq. (27), ðrÞT ðrÞ is a matrix that represents a second order

rix u!0

lim Z 2;i ¼ const2;

riy u!0

rix u!1

a2i : Z 2;i ¼  y 2 h ri u þ a22i iqðuÞ

riy u!0

where const1 and const2 are constants. That is, the pixels in the homogeneous regions are imposed on a strong smooth constraint, which can suppress the noise. For the discontinuities, their gradient modules are very large, hence ki is almost equal to 1. For the edge, on the direction perpendicular to it, the difference between the adjacent pixels is larger, while along the edge direction, the difference is smaller. Taking any pixel i on the vertical edge as an example, rix u is large, and riy u is very small, and we can obtain

ð30Þ ð31Þ

For pixels on edges, a very weak smooth constraint is imposed on the direction perpendicular to the edge, and a strong smooth constraint is imposed on the edge direction. Therefore, the edge details can be preserved. The corner is with large gradient module on both horizontal and vertical directions, from Eqs. (34) and (35) we know that a very weak smooth constraint is imposed on the corner, thus, it can avoid blurring. Sometimes, there exists stronger noise in images. The noised pixel has relatively large gradient, leading to very small weighting T

functions, in this case, the smoothing effect of ðr x Þ W 1 ðr x Þ and y T

ðr Þ W 2 ðr y Þ is not very good. However, the relatively large gradient could make ki 2 ð0; 1Þ, thus, the filter ðrÞT ðrÞ can work and suppress noise. In essence, the following expressions are obtained by adding a quadratic term as the additional term,

662

R. Jin et al. / J. Vis. Commun. Image R. 40 (2016) 656–670

Table 1 PSNR and SSIM comparisons with UF = 2, SNR = 5 dB. EIAcen

Bicubic SSR MSR BEP AA l1-SAR TV-SAR Ours

Text

House

Cameraman

PSNR

SSIM

PSNR

SSIM

PSNR

SSIM

PSNR

SSIM

21.46 21.34 8.97 23.03 23.04 23.45 24.22 25.64

0.5463 0.5463 0.1574 0.6788 0.6998 0.7005 0.7065 0.7367

16.33 17.01 12.93 21.27 22.61 21.67 21.60 23.97

0.3707 0.3650 0.2162 0.4380 0.4579 0.4470 0.4464 0.4775

27.97 28.59 17.27 31.44 33.18 32.58 32.56 34.57

0.5528 0.5500 0.1559 0.5748 0.6384 0.6321 0.6319 0.6552

23.94 24.22 16.31 26.65 27.62 27.17 27.16 28.89

0.4874 0.4492 0.1386 0.4932 0.5619 0.5419 0.5420 0.6519

Table 2 PSNR and SSIM comparisons with UF = 2, SNR = 10 dB. EIAcen

Bicubic SSR MSR BEP AA l1-SAR TV-SAR Ours

Text

House

Cameraman

PSNR

SSIM

PSNR

SSIM

PSNR

SSIM

PSNR

SSIM

22.15 22.26 9.26 25.12 24.78 28.41 28.37 30.10

0.612 0.612 0.2011 0.7493 0.7526 0.7596 0.569 0.7867

16.35 17.03 12.91 22.62 26.10 23.62 25.39 24.37

0.5486 0.4600 0.5315 0.4942 0.4477 0.5209 0.5001 0.5114

28.07 28.95 16.76 33.52 35.77 34.78 34.79 37.40

0.6054 0.6429 0.1460 0.6488 0.7451 0.7464 0.7461 0.7876

24.10 24.74 16.50 28.66 29.78 28.95 28.93 31.67

0.6436 0.6438 0.1970 0.5580 0.7134 0.7139 0.7140 0.7588

Table 3 PSNR and SSIM comparisons with UF = 4, SNR = 5 dB. Pepper

Bicubic SSR MSR BEP AA l1-SAR TV-SAR Ours

Tower

Panda

Barbara

PSNR

SSIM

PSNR

SSIM

PSNR

SSIM

PSNR

SSIM

23.01 23.18 15.86 28.36 29.43 28.67 28.64 30.42

0.6278 0.6399 0.2469 0.7472 0.7886 0.7914 0.7915 0.8082

26.89 26.95 21.47 32.24 32.67 32.42 32.49 33.17

0.5474 0.5579 0.2548 0.7083 0.7748 0.7610 0.7634 0.7764

21.46 21.34 15.89 25.11 25.15 24.39 24.39 26.02

0.3749 0.3550 0.1625 0.4836 0.4730 0.4493 0.4499 0.5462

21.02 21.14 15.37 27.67 28.04 27.76 27.74 28.16

0.5704 0.5999 0.2341 0.8188 0.8473 0.8417 0.8413 0.8512

Table 4 PSNR and SSIM comparisons with UF = 4, SNR = 10 dB. Pepper

Bicubic SSR MSR BEP AA l1-SAR TV-SAR Ours

G1;i ¼

G2;i ¼

Tower

Panda

Barbara

PSNR

SSIM

PSNR

SSIM

PSNR

SSIM

PSNR

SSIM

23.07 23.30 15.71 30.62 32.25 30.27 30.23 34.22

0.6564 0.6835 0.2420 0.8340 0.8896 0.8818 0.8817 0.8978

26.94 27.03 21.55 33.66 33.87 33.54 33.56 34.36

0.5677 0.5951 0.2565 0.8267 0.8555 0.8538 0.8551 0.8631

22.15 22.26 15.51 26.90 27.86 27.26 27.26 28.05

0.4940 0.4911 0.1814 0.5877 0.6415 0.6250 0.6251 0.6157

21.06 21.21 15.33 29.78 29.54 29.13 29.11 29.78

0.5809 0.6151 0.2240 0.8890 0.9037 0.9011 0.9008 0.9120

h    2 i0 ki /0 rix u þ ð1  ki Þ rix u 2rix u h    2 i0 ki /0 riy u þ ð1  ki Þ riy u 2riy u

ki a1i ¼ qffiffiffiffiffi þ 1  ki ; 2 g 1i

ð36Þ

ki a2i ¼ qffiffiffiffiffi þ 1  ki : 2 g 2i

ð37Þ

Obviously, when ki ¼ 1; G1;i and G2;i are equal to F 1;i and F 2;i , respectively; When ki is close to zero, the new prior model is quadratic and strongly penalizes the high-frequency information; When ki 2 ð0; 1Þ, for the pixel with large gradients, G1;i and G2;i are not equal to zero, and the stronger noise may be possibly suppressed. The combination strategy can improve condition (ii)

 0 0 limf !þ1 /2fðf Þ ¼ 0 , i.e., limf !þ1 /2fðf Þ ¼ 0 for the discontinuities

663

R. Jin et al. / J. Vis. Commun. Image R. 40 (2016) 656–670

(a) Bicubic

(b) SSR

(c) MSR

(d) BEP

(e) Bicubic error

(f) SSR error

(g) MSR error

(h) BEP error

(i) AA

(j) TV-SAR

(k) l1-SAR

(l) Ours

(m) AA error

(n) TV-SAR error

(o) l1-SAR error

(p) Ours error

Fig. 5. The reconstructed images of ‘‘cameraman” and the corresponding error images with UF = 2 and SNR = 5 dB.

0

while limf !þ1 /2fðf Þ ¼ constant for the stronger noise. Therefore, the proposed combination strategy is effective. 5. Experimental results and discussion In this section, we will show our experimental results in detail. The experiment environment and parameter settings are presented in detail in subSection 5.1. And then the evaluation standards and the simulated experimental results obtained by different SR methods are shown in subSection 5.2. And some of the real observations are shown in subSection 5.3. 5.1. General description In this paper, we measure the performance of our proposed method in terms of both objective measurements and subjective measurements, and compare it with the bicubic interpolation

method, the singleframe SR (SSR) method [5], the multiframe SR (MSR) method [26], the BEP model based method [21], the AA model based method [24], the TV-SAR model based method [23], and the l1-SAR model based method [23]. The eight images presented in Fig. 4 are used for simulation experiments. In order to obtain the LR images, the following five steps are used: sub-pixel motion, rotation, blurring, downsampling and adding noise. In the experiments, we let hk 2 ð3 ; 3 Þ, xk 2 ð1; 1Þ, yk 2 ð1; 1Þ, (k = 1, 2, 3, 4, 5) and the transformation of the LR images be variant from each other. A 3  3 Gaussian filter with 0 mean, 1 variance is used for Gaussian blur. The unsample factors (UF) are 2 and 4, respectively. Finally, white Gaussian noise is added to the LR images with Signal-toNoise Ratio (SNR) levels of 5 dB and 10 dB, respectively. Note that this is a much stronger noise than usually did. Thus, under each level, 5 LR images are obtained from each original image with UF = 2, and 17 LR images are obtained from each original image with UF = 4.

664

R. Jin et al. / J. Vis. Commun. Image R. 40 (2016) 656–670

(a) Bicubic

(b) SSR

(c) MSR

(d) BEP

(e) Bicubic error

(f) SSR error

(g) MSR error

(h) BEP error

(i) AA

(j) TV-SAR

(k) l1-SAR

(l) Ours

(m) AA error

(n) TV-SAR error

(o) l1-SAR error

(p) Ours error

Fig. 6. The reconstructed images of ‘‘text” and the corresponding error images with UF = 2 and SNR = 10 dB.

For real experiments, the image sequences are downloaded from Milanfar’s web site (https://users.soe.ucsc.edu/milanfar/software/sr-datasets.html). In the following, the initial HR image u0 , the initial auxiliary 0 0 and g 20 variables g 10 i i , and the initial hyper-parameters a1 ; a2 and 0 0 c are given. And the following initial values are used: u is a bicubic interpolation of the first observed LR image; pffiffiffiffiffi g 10 ¼ jrix u0 j2 þ a21i ; g 20 ¼ jriy u0 j2 þ a22i ; a01 ¼ P M=2 ; a02 ¼ i i i

P i

pffiffiffiffiffi

M=2 ða1i

ða1i

g 10 a21i Þ i

M , and c0 ¼ jjv Bðs . 0 2 k k Þu jj2

g 20 a22i Þ i

For the TV-SAR and l1-SAR methods [23], they are implemented by running the MATLAB codes provided by the authors. To be fair, the BEP model [21] is applied to the variational Bayesian framework in this paper, and the other parameters used in it are manually adjusted for the best SR reconstruction performance. Typically, we use MATLAB R2009a on a 3.0 GHz Pentium Dual core computer with 2.0 GB RAM.

5.2. Simulation experiments In this subsection, we will show the comparison experiments between our proposed method and some other SR methods, such as the bicubic interpolation method, the SSR method, the MSR method, the BEP method, the AA method, the TV-SAR method, and the l1-SAR method. To show the comparison, the PSNR values, and the SSIM values are used. It should be noted that in our work SNR = 5 dB and SNR = 10 dB are used, which means that our experiments are with strong noise. The PSNR values and the SSIM values achieved by the above eight methods are shown in Tables 1–4, where Tables 1 and Table 2 are with UF = 2, SNR = 5 dB and SNR = 10 dB respectively, and Tables 3 and Table 4 are with UF = 4, SNR = 5 dB and SNR = 10 dB, respectively. From them it is clearly observed that, the MSR method has achieved the lowest PSNR and SSIM values in each reconstruction problem; in most cases, the bicubic interpolation method and the SSR method have almost the same PSNR and SSIM

665

R. Jin et al. / J. Vis. Commun. Image R. 40 (2016) 656–670

(a) Bicubic

(b) SSR

(c) MSR

(d) BEP

(e) Bicubic error

(f) SSR error

(g) MSR error

(h) BEP error

(i) AA

(j) TV-SAR

(k) l1-SAR

(l) Ours

(m) AA error

(n) TV-SAR error

(o) l1-SAR error

(p) Ours error

Fig. 7. The reconstructed images of ‘‘panda” and the corresponding error images with UF = 4 and SNR = 5 dB.

values, and they perform better than the MSR method; the BEP method perform much better than the MSR method, the bicubic interpolation method and the SSR method; the TV-SAR method, the l1-SAR method and the AA method perform much better than BEP method. It is observed that our method has achieved the highest PSNR and SSIM values in nearly all the reconstruction problems, demonstrating the importance and necessity of the proposed method. For visual comparison, the HR images reconstructed by different methods are provided in Figs. 5–8. Fig. 5 shows the reconstructed images of the ‘‘cameraman” image and the corresponding error images with UF = 2 and SNR = 5 dB, Fig. 6 shows the reconstructed images of the ‘‘text” image and the corresponding error images with UF = 2 and SNR = 10 dB, Fig. 7 shows the reconstructed images of the ‘‘panda” image and the corresponding error images with UF = 4 and SNR = 5 dB, and Fig. 8 shows the reconstructed images of the ‘‘pepper” image and the corresponding error images with

UF = 4 and SNR = 10 dB. The error images, i.e., the difference between the estimated HR image and the original image. It is clearly seen from these figures that, the bicubic interpolation method and the SSR method produce blur images because they do not exploit the underlying structures in natural images fully; the MSR method, the BEP method and the AA method are uneffective in edge-preserving and artifact-suppressing; the TV-SAR method and the l1-SAR method have achieved similar reconstruction performance, and obtain better reconstructions than the bicubic interpolation method, the BEP method and the AA method. As expected, our method achieves the best visual quality among all the methods, thus its advantage to other methods is obvious in the term of visual perception. In the error images, brighter pixels represent a large error. From these error images, the difference between different SR approaches is clearly observed. In general, our proposed approach can obtain the reconstructed images with the best visual quality.

666

R. Jin et al. / J. Vis. Commun. Image R. 40 (2016) 656–670

(a) Bicubic

(b) SSR

(c) MSR

(d) BEP

(e) Bicubic error

(f) SSR error

(g) MSR error

(h) BEP error

(i) AA

(j) TV-SAR

(k) l1-SAR

(l) Ours

(m) AA error

(n) TV-SAR error

(o) l1-SAR error

(p) Ours error

Fig. 8. The reconstructed images of ‘‘pepper” and the corresponding error images with UF = 4 and SNR = 10 dB.

5.3. Experiments on real data In this subsection, the performance of our proposed method has been discussed with real observations. We will show the comparison experiments between our proposed method and some other SR methods, such as the bicubic interpolation method, the BEP method, the AA method, the TV-SAR method, and the l1-SAR method. The data set used in this subsection is the popularly used video sequences, downloaded from Milanfar’s web site (https://users.soe.ucsc.edu/milanfar/software/sr-datasets.html). We used the first ten LR images to reconstruct the HR image. Figs. 9 and 10 present the estimated HR images for two real image sequences downloaded from internet using different SR methods. In Figs. 9 and 10, the reconstructed images obtained by the bicubic interpolation method are blurred. From Figs. 9(b), (c) and 10(b), (c), it is noted that there are some artifacts in the reconstructed HR images by using the BEP method and the AA method.

By incorporating a quadratic function, such as the TV-SAR method (i.e., Figs. 9(d) and 10(d)), and the l1-SAR method (i.e., Figs. 9(e) and 10(e)), such artifacts have been eliminated to a certain extent. The images reconstructed by our method (i.e., Figs. 9 and 10(f)) own clearer edges. These indicate that our method is superior to the SR methods under comparison. 6. Conclusion This paper proposes a multiframe SR method by involving the reconstruction properties of the anisotropic adaptive prior model together with the Tikhonov prior model. This method is based on a convex combination of the two prior models, and an iterative way to determine the combination function according to the local features of the underlined HR image. Moreover, the effectiveness of the proposed hybrid prior model has been demonstrated theoretically and practically. The hybrid prior model can benefit of both the anisotropic adaptive prior model and the Tikhonov prior

R. Jin et al. / J. Vis. Commun. Image R. 40 (2016) 656–670

(a) Bicubic

(b) BEP

(c) AA

(d) TV-SAR

(e) l1-SAR

(f) Ours

667

Fig. 9. The reconstructed images for ‘‘Alpaca” sequence obtained by using different methods.

model, surpassing each individually in image reconstruction. In this paper, a variational Bayesian inference estimate of the superresolved image is obtained, with the motion parameter and hyper-parameters are determined automatically. Numerous experimental results demonstrate the effectiveness of the proposed approach and its superiority to several other SR methods. The proposed image prior model could suppress artifacts while preserving edges, however, a lot of information or details cannot be well restored until now. In our future work, we will focus on estimating the missing information or the reconstruction error, and establish a more efficient reconstruction model.

Acknowledgement This research was supported by NSFC Grant No. 61370179, 71271125, 61502260, the Fundamental Research Funds for the Central Universities, HUST No. 2015YGYL012 and No. 2016YXMS086. Appendix A. The solving process in detail Differentiating the integral on the right hand side in Eq. (12) with respect to qðuÞ; qðsk Þ; qðrk Þ; qða1 Þ; qða2 Þ, and qðcÞ can result in the following expressions.

668

R. Jin et al. / J. Vis. Commun. Image R. 40 (2016) 656–670

(a) Bicubic

(b) BEP

(c) AA

(d) TV-SAR

(e) l1-SAR

(f) Ours

Fig. 10. The reconstructed images for ‘‘Emily” sequence obtained by using different methods.

8 2 0 1  x 2 > < 2 1 X r u þ a þ g 6 B C i i qffiffiffiffiffi1i  a21i A qðuÞ / exp k4ha1 i @a1i > : i 2 g 1i 0 13  y 2 2 2 XB ri u þ a2i þ g i C7 qffiffiffiffiffi  a22i A5  ð1  kÞha3 ikruk22 þha2 i @a2i 2 i 2 gi ) X 2 ðA:1Þ  hrk ihkv k  Bðsk Þuk2 iqðsk Þ ; k

where hiqðÞ denotes the expected value using qðÞ, i.e., hiqðÞ ¼ EqðÞ ðÞ, and hi is used for simplicity in the rest of this paper. g 1i and g 2i are obtained from Eq. (12) as follows, 2

2

g 1i ¼ hðrix uÞ þ a21i iqðuÞ ; g 2i ¼ hðriy uÞ þ a22i iqðuÞ ;   T 0 1   1 qðsk Þ / exp  hrk ihkv k  Bðsk Þuk22 iqðuÞ þ sk  s0k Rk sk  s0k ; 2 ðA:2Þ

669

R. Jin et al. / J. Vis. Commun. Image R. 40 (2016) 656–670 N1þa0

qðrk Þ / rk2

(

 exp rk brk þ 0

M1þa0 a1 4

qða1 Þ / a1

hkv k  Bðsk Þuk22 iqðuÞ;qðsk Þ

!) ;

2

( exp a1 ba1 0

ðA:3Þ

!) X qffiffiffiffiffi 2 1 ; þ a1i g i  a1i

ðA:4Þ

i

M1þa0

qða2 Þ / a24

a2

( exp a2 ba2 þ 0

!) X qffiffiffiffiffi a2i g 2i  a22i ;

( qðcÞ / c 2 1þac exp c bc þ 0

0

hkruk22 iqðuÞ 2

Rk ¼



R0k

1

þ hrk iðCk þ Wk Þ

1 ;

ðA:13Þ

and

hsk i ¼ Rk



R0k

1

!) :

ðA:6Þ

In Eqs. (A.1), (A.2), (A.3) and (A.6), hkv k  Bðsk Þuk22 iqðsk Þ and

hkv k  Bðsk Þuk22 iqðuÞ are not easy to calculate since Cðsk Þ is nonlinear with respect to sk . Therefore, Cðsk Þ is expanded using its first-order Taylor series around the mean value hsk i of qðsk Þ, similar to [23], resulting in

"

#

@C

@C

@C

ðsk  hsk iÞ: Cðsk Þ Cðhsk iÞ þ ; ; @hk ðhsk iÞ @dxk ðhsk iÞ @dyk ðhsk iÞ

ðA:7Þ

 s0k þ hrk i½Ck hsk i þ Wk hsk i þ Q k  Uk  ;

ðA:14Þ

where Ckij ¼ !Tki !kj , for i; j ¼ 1; 2; 3, and Q ki ¼ ðv k  Bðhsk iÞT Oki ðhsk iÞhui, for i ¼ 1; 2; 3. Finally, we obtain the means of the distributions in Eqs. (A.3)– (A.6),

hrk i ¼

N=2 þ a0rk hkv k Bðsk Þuk22 iqðuÞ;qðs 2

kÞ

0

þ brk

;

M=4 þ a0a1 ; ha1 i ¼ X qffiffiffiffiffi  0 a1i g 1i  a21i þ ba1

ðA:15Þ

ðA:16Þ

i

Thus, Bðsk Þ can be approximated by

Bðsk Þ Bðhsk iÞ þ ½Ok1 ðhsk iÞ; Ok2 ðhsk iÞ; Ok3 ðhsk iÞðsk  hsk iÞ;

ðA:8Þ

@C @C Ok1 ðhsk iÞ ¼ Dk Hk @h jðhsk iÞ ; Ok2 ðhsk iÞ ¼ Dk Hk @dx jðhsk iÞ ,

where

ðhsk iÞhui; Ok3 ðhsk iÞhui; Uki ¼

trace½Bðhsk iÞT Oki ðhsk iÞRu , and Wkij ¼ trace½Oki ðhsk iÞT Okj ðhsk iÞRu , for i; j ¼ 1; 2; 3. Thus, the covariance and mean value of qðsk Þ can be calculated as

ðA:5Þ

i

M

!k ¼ ½Ok1 ðhsk iÞhui; Ok2

where

rk

k

and

k

M=4 þ a0a2 ; ha2 i ¼ X qffiffiffiffiffi  0 a2i g 2i  a22i þ ba2

ðA:17Þ

i

@C jðhsk iÞ . Ok3 ðhsk iÞ ¼ Dk Hk @dy k

Then the following approximation of hkv k  obtained,

Bðsk Þuk22 iqðsk Þ

can be

hkv k  Bðsk Þuk22 iqðsk Þ kv k  Bðhsk iÞuk22 þ

3 X 3 X hrk inkij uT Oki ðhsk iÞT Okj ðhsk iÞu:

ðA:9Þ

hci ¼

M=2 þ a0c hjjrujj22 iqðuÞ 2

0

þ bc

:

ðA:18Þ

References

i¼1 j¼1

Substituting Eq. (A.9) into Eq. (A.1), the inverse covariance and mean value of qðuÞ can be calculated as

"

Ru ¼

X XXX hrk iBðhsk iÞT Bðhsk iÞ þ hnkij Oki ðhsk iÞT Okj ðhsk iÞ k

k

i

j

þ kðha1 iðr Þ W 1 ðr Þ x T

x

i1 T þ ha2 iðr y Þ W 2 ðr y ÞÞ þ ð1  kÞha3 irT r ;

ðA:10Þ

and

" # X T hui ¼ Ru hrk iBðhsk iÞ v k ;

ðA:11Þ

k

where nkij is a 3  3 region in the covariance matrix Rk of qðsk Þ; W 1 and W 2 are M  M diagonal matrices, in which, the elements on the diagonal are ag1i1 and ag2i2 , respectively. i

i

The calculation of the covariance matrix Ru is a nontrivial job. Therefore, the conjugate gradient method is used to solve Eq. (A.11). And the approximation of hkv k  Bðsk Þuk22 iqðuÞ can be obtained by using Eq. (A.8),

hkv k  Bðsk Þuk22 iqðuÞ kv k  Bðhsk iÞhui  !k ðsk  hsk iÞk22 þ trace½Bðhsk iÞT Bðhsk iÞRu  þ 2UTk ðsk  hsk iÞ þ ðsk  hsk iÞT Wk ðsk  hsk iÞ;

ðA:12Þ

[1] J. Tian, K. Ma, A survey on super-resolution imaging, Signal Image Video Process. 5 (3) (2011) 329–342. [2] R. Tsai, T. Huang, Multiframe image restoration and registration, Adv. Comput. Vision Image Process. 1 (1984) 317–339. [3] K. Guo, X. Yang, H. Zha, W. Lin, S. Yu, Multiscale semilocal interpolation with antialiasing, IEEE Trans. Image Process. 21 (2) (2012) 615–625. [4] K. Guo, X. Yang, W. Lin, R. Zhang, S. Yu, Learning based super-resolution method with a combining of both global and local constraints, IET Image Process. 6 (4) (2012) 337–344. [5] T. Peleg, M. Elad, A statistical prediction model based on sparse representations for single image super-resolution, IEEE Trans. Image Process. 23 (6) (2014) 2569–2582. [6] J.-B. Huang, A. Singh, N. Ahuja, Single image super-resolution from transformed self-exemplars, in: 2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), Boston, Massachusetts, USA, 2015, pp. 5197– 5206. [7] J. Xie, R. Feris, S. Yu, M. Sun, Joint super resolution and denoising from a single depth image, IEEE Trans. Multimedia 17 (9) (2015) 1525–1537. [8] Z. Wang, D. Liu, J. Yang, W. Han, T. Huang, Deep networks for image superresolution with sparse prior, in: 2015 IEEE International Conference on Computer Vision (ICCV), Santiago, Chile, 2015, pp. 370–378. [9] C. Dong, C. Loy, K. He, X. Tang, Image super-resolution using deep convolutional networks, IEEE Trans. Pattern Anal. Machine Intell. 38 (2) (2016) 295–307. [10] X. Zhang, E. Lam, E. Wu, K. Wong, Application of Tikhonov regularization to super-resolution reconstruction of brain mri images, in: Medical Imaging and Informatics - 2nd International Conference, 2008, pp. 51–56. [11] M. Ng, H. Shen, E. Lam, L. Zhang, A total variation regularization based superresolution reconstruction algorithm for digital video, Eurasip J. Adv. Signal Process. 2007 (2007) 298–311. [12] S. Farsiu, M. Robinson, M. Elad, P. Milanfar, Fast and robust multiframe super resolution, IEEE Trans. Image Process. 13 (10) (2004) 1327–1344. [13] Q. Yuan, L. Zhang, H. Shen, Multiframe super-resolution employing a spatially weighted total variation model, IEEE Trans. Circ. Syst. Video Technol. 22 (3) (2012) 379–392. [14] F. Li, C. Shen, J. Fan, C. Shen, Image restoration combining a total variational filter and a fourth-order filter, J. Visual Commun. Image Represent. 18 (4) (2007) 322–330.

670

R. Jin et al. / J. Vis. Commun. Image R. 40 (2016) 656–670

[15] P. Charbonnier, L. Blanc-Feraud, G. Aubert, M. Barlaud, Deterministic edgepreserving regularization in computed imaging, IEEE Trans. Image Process. 6 (2) (1997) 298–311. [16] H. Shekarforoush, M. Berthod, J. Zerubia, 3D super-resolution using generalized sampling expansion, IEEE International Conference on Image Processing, Washington, DC, USA, vol. 2, 1996, pp. 300–303. [17] R. Schultz, R. Stevenson, Extraction of high-resolution frames from video sequences, IEEE Trans. Image Process. 5 (6) (1996) 996–1011. [18] K. Suresh, G. Kumar, A. Rajagopalan, Superresolution of license plates in real traffic videos, IEEE Trans. Intell. Transport. Syst. 8 (2) (2007) 321–331. [19] K. Suresh, A. Rajagopalan, Robust and computationally efficient superresolution algorithm, J. Opt. Soc. Am. A: Opt. Image Sci. Vision 24 (4) (2007) 984–992. [20] M. Ceccarelli, A finite markov random field approach to fast edge-preserving image recovery, Image Vision Comput. 25 (6) (2007) 792–804. [21] X. Zeng, L. Yang, A robust multiframe super-resolution algorithm based on half-quadratic estimation with modified btv regularization, Digital Signal Process. 23 (1) (2013) 98–109. [22] W. Shao, H.S. Deng, Z. Wei, Kullback-leibler divergence based composite prior modeling for bayesian super-resolution, J. Sci. Comput. 60 (1) (2014) 60–78.

[23] S. Villena, M. Vega, S. Babacan, R. Molina, A. Katsaggelos, Bayesian combination of sparse and non sparse priors in image super resolution, Digital Signal Process. 23 (2) (2013) 530–541. [24] S. Zhao, R. Jin, X. Xu, E. Song, C. Hung, A variational bayesian superresolution approach using adaptive image prior model, Math. Problems Eng. (2015). [25] D. Zhang, P.-M. Jodoin, C. Li, Y. Wu, G. Cai, Novel graph cuts method for multiframe super-resolution, IEEE Signal Process. Lett. 22 (12) (2015) 2279–2283. [26] Z. Wang, D. Liu, J. Yang, W. Han, T. Huang, Handling motion blur in multiframe super-resolutions, in: 2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), Boston, Massachusetts, USA, 2015, pp. 5224– 5232. [27] S. Babacan, R. Molina, A. Katsaggelos, Variational bayesian super resolution, IEEE Trans. Image Process. 20 (4) (2011) 984–999. [28] S. Villena, M. Vega, R. Molina, A. Katsaggelos, Bayesian super-resolution image reconstruction using an l1 prior, in: 2009 Proceedings of 6th International Symposium on Image and Signal Processing and Analysis, Piscataway, NJ, USA, 2009, pp. 152–157. [29] C. Sung, K. Min, G. Moon, Super-resolution image reconstruction: a technical overview, IEEE Signal Process. Mag. 20 (3) (2003) 21–36. [30] L. Zhang, H. Zhang, H. Shen, P. Li, A super-resolution reconstruction algorithm for surveillance images, Signal Process. 90 (3) (2010) 848–859.