Abstract

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

1. Introduction

Superresolution (SR) technique [1–5] is an important branch in image fusion technology which targets the reconstruction of a high-resolution (HR) image from a set of degraded low-resolution (LR) images. The LR images are usually affected by warping, blurring, downsampling, and noising. Due to advancement of the technology, superresolution has been widely used in a broad range of applications such as computer vision, medical imaging, public safety, and military reconnaissance.

SR is a typically ill-posed inverse problem that is hard to solve without the introduction of some prior image information [6–8]. Thus, a number of regularization-based SR approaches have been proposed [9–11] by incorporating the prior knowledge of the unknown high-resolution image in the regularization strategy [12]. The authors in [9] proposed a regularization-based SR approach based on the Tikhonov regularizer using norm. It can remove noises effectively; however, it also blurs the image edges. Farsiu et al. [10] proposed the bilateral total variation (BTV) regularizer, which penalized gradient magnitudes measured by norm, with the intention of preserving edges in the reconstruction. However, the BTV regularizer often produces artifacts in the smoothed regions due to the use of norm [11]. In [11], the authors improved BTV with the norm , where is the gradient value and is a positive scale parameter, and proposed the bilateral edge-preserving (BEP) regularizer. This new norm can adaptively use and norms, and the transition from to can be controlled by modifying the positive scale parameter . Unfortunately, the BEP regularizer employed a fixed scale parameter for the whole image and ignored the local image features. With this drawback, it cannot well control the level of punishment to local gradients.

It is well known that accurate registration which refers to the estimation of motion information is a very important factor in the success of the SR reconstruction method [13]. The above SR methods only registered the observations once before reconstructing the HR image. This is not a robust method. The error caused by an inaccurate registration would produce poor effects in the reconstruction process [13]. Other than the motion parameters, the hyperparameters related to the image prior model and the noise model could intimately affect the image reconstructing quality [12]. Recently, the variational Bayesian method has been used for SR [14–16]. This method approximates the joint posterior distribution of the unknowns, including the HR image, the motion parameters, and the hyperparameters, using a tractable distribution by minimizing the Kullback-Leibler (KL) divergence and by estimating the unknowns simultaneously in each iteration.

The Bayesian approach uses a prior model to estimate a priori knowledge of the unknown HR image. For edge-preservation, the variational Bayesian approaches based on TV prior model [14] and prior model [16] have been proposed. In these two prior models, the gradient magnitudes are measured by the norm, and consequently both of them can preserve edges. Compared with the TV prior, the prior contains two model parameters, which constrain certain preferred edge directions [15, 16]. However, the norm may not always distinguish the real image features from the effects of errors caused by inaccurate registration and additive Gaussian noise. Therefore, some artifacts will be produced in the smoothed image regions. Villena et al. made progress by combining the TV prior model and simultaneous autoregressive (SAR) (i.e., based on norm) prior model [15]. A similar approach is also applied to the prior model. The problem is that it is difficult to choose a proper parameter to balance the priors in the combination.

In this paper, a variational Bayesian estimate of the HR image, the motion parameters, and the hyperparameters are proposed and derived and then given a set of LR images. In order to overcome the difficulty of aforementioned SR methods, an adaptive image prior model is proposed, in which the degree of interaction between pixels is adjusted adaptively by using the norm in order to preserve edges and remove noises. It is important to note that we propose a method for automatically estimating the scale parameters (i.e., ) to the proposed prior model. Information needed to determine these scale parameters is updated iteratively based on the available estimated HR image. Experiment results show the effectiveness of the proposed SR method.

The rest of this paper is organized as follows: the observation model and a Bayesian framework for image SR are described in Section 2. The adaptive image prior model is sketched in Section 3. The proposed SR approach is presented in Section 4. Experimental results are demonstrated in Section 5. Finally, the conclusion is given in Section 6. In addition, Appendix is presented at the end of the paper showing the solving process in detail.

2. The Problem Formulation

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

2.1. The Observation Model

Formulating an observation model that relates the original HR image to the LR observations is the first step to fully analyze the SR reconstruction problem.

Let denote the original HR image of size , and each observed LR image is of size . and denote the downsampling factors in the horizontal and vertical directions, respectively. The HR image and LR images are written in lexicographical notations as the vectors of sizes and , respectively. In this work, the image acquisition process is modeled by geometric transformation, blurring, downsampling, and adding with white Gaussian noise. Thus, the following popular matrix notation [15] is adopted to describe the acquisition process:where is one of a set () of LR images, (the dimension of is ) is the downsampling matrix, (the dimension of is ) is the blur matrix, (the dimension of is ) is the warp matrix, and (the dimension of is ) is the additive white Gaussian noise.

The motion that occurs during the image acquisition is represented by warp matrix . In this paper, we assume that the motion includes global rotation and translation; that is, , where , , and are the rotation angle and the displacement in horizontal direction and vertical direction of the th HR image, respectively, relative to the reference frame. In this study, we assume that the blur matrices are determined by sensors. The sensor PSF is usually modeled as a spatial averaging operator, and the characteristics of this blur are usually assumed to be known. Under these assumptions, the matrices and have a block circulant with circulant block structure. The downsampling matrix is determined by the downsampling factors and , and it generates LR images from the warped and blurred HR image.

For convenience, we define and , where is a set.

2.2. The Hierarchical Bayesian Framework

We use a hierarchical Bayesian framework [17] to model the acquisition process, the HR image , the motion parameters , and the hyperparameters and . Parameters and are the model parameters of our proposed prior model and the noise model, respectively. Thus, we can obtain the following joint posterior distribution of , , , and given by using the Bayes rule:where represents the conditional distribution of the LR images and , , , and are prior distributions for the unknown image , the motion parameters , and the hyperparameters and , respectively.

Suppose that the noise in the LR image () is the white Gaussian noise with a zero-mean distribution , and we can writewhere denotes the norm of a given vector .

Due to , we can obtain the conditional distribution of the LR image ,Moreover, if we assume the statistical independence of noises among the LR images, we can obtain the conditional distribution of the LR images

By using the available registration algorithm, the initial values of motion parameters can be obtained. However, the obtained initial values are often inaccurate. In this paper, we refine the values of motion parameters by modeling them as stochastic variables following Gaussian distributions , which is similar to the model used in [15].

The prior distributions for the hyperparameters are defined to be Gamma distributions (, ). Gamma distributions for the hyperparameters were selected because they are conjugate for the variance of the Gaussian distribution; therefore, the posteriors will also have Gamma distributions in the Bayesian formulation [18].

Our proposed adaptive image prior model, denoted by , will be presented in the next section.

3. The Proposed Adaptive Image Prior Model

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

3.1. The Adaptive Image Prior Model

The adaptive norm is used to measure the horizontal and vertical gradients in which is the gradient value and is a positive scale parameter. These parameters are to determine the severity of penalizing image gradients. A new adaptive image prior model is proposed and defined as follows:where is the number of pixels in the HR image. Symbols and represent the horizontal and vertical gradient components, respectively, for the pixel , is the hyperparameter of this prior controlling the degree of regularization, and and are the scale parameters of the adaptive norms measuring the horizontal and vertical gradient component, respectively, of pixel .

If we take the logarithm, (6) becomeswhere is proportional to a linear combination of and , . Equation (7) is strictly a convex function because is strictly convex [11].

In addition, is adaptive, because the following approximations of can be obtained with the assumption that parameter is fixed:Thus, or , , in (7) can use either norm (i.e., (9)) or norm (i.e., (8)) adaptively based on the evaluation of the local image features.

Thus, the expression can be approximated as

For example, if there exists an edge along the horizontal direction in an LR image, of the pixel at the edge is relatively large while is small. Thus, approximates to and approximates to . Therefore, approximates to (10). In the SR process based on (10), when the LR pixel splits in both horizontal and vertical directions into a block of HR pixels, the large gradient magnitude in the vertical direction will be preserved due to the norm, and the small gradient magnitude in the horizontal direction will be smoothed due to the norm. Thus whenever there exists an edge along horizontal or vertical direction, can act as the norm in the direction perpendicular to the edge, preserving large gradients, and, at the same time, can act as the norm in the direction along the edge, obtaining relatively ideal smoothing effects. In (12), when and , the norm acts on both horizontal and vertical directions, to preserve the edges. In (13), when and , the noise will be effectively removed.

3.2. Adaptive Calculation of Parameter

It is known that the gradients produced by edges should be preserved, while the gradients produced by errors, which suffer from inaccurate registration and additive Gaussian noise, need to be smoothed. Thus, the severity of penalizing image gradients should be determined according to the local image features. In , the scale parameter, , can control the severity of penalizing local gradients. For the gradients produced by errors, norm is a good choice. In this case, the scale parameter should be set as a large value. When the value of the scale parameter increases, accepts a larger range of errors and handles them like norm [11]. Otherwise, it should be set as a small value. Therefore, the property of scale parameter should be constrained with the following conditions: (1) its value is determined according to the local image features, (2) the value is inversely proportional to gradient magnitudes, and (3) the value is larger than zero.

We use the well-known monotonically decreasing function [19] to calculate , and in (6),where and represent the magnitudes of horizontal and vertical gradient components, respectively, for the pixel .

For convenience, we use the symbol to represent , and ; that is, .

4. Variational Bayesian Superresolution

The variational Bayesian inference [14] is utilized to estimate the unknowns. This variational Bayesian technique approximates the true posterior distribution analytically by a tractable distribution that minimizes the KL divergence, which can measure the difference between the two distributions and . That is to say, by minimizing this KL divergence, we can obtain the optimal approximation distribution. Thus, we can obtain the following expression: where and .

Thus, (11) can be rewritten asBecause (16) cannot be evaluated due to the form of (6), the inequality with , is employed in (6). This inequality was used in [14, 15]. Then, we can obtain where , the auxiliary variables and are related to the unknown HR image, and they need to be estimated (see (22)). Thus, we can obtain the upper bound of where .

As shown in [20], the minimization of can be replaced by the minimization of its upper bound . Thus, we can obtain the distributions , , , and and by minimizing ; that is,

The details of equations derivation are given in Appendix at the end of this paper. We obtain for the posterior distribution the multivariate Gaussian with mean valueand covariance where denotes the expected value of , denotes the transposed operator, is a 3 × 3 region in the covariance matrix in (24), , , where , , and denote partial derivatives of for , , and , respectively, and , is a diagonal matrix, in which the element on the diagonal is .

Then, the following expressions are obtained for:where denotes the expected value of using .

We also obtain for the posterior distribution the multivariate Gaussian with mean valueand covarianceIn (23) and (24), ,  ,   = with , and , for .

Finally, we obtain the mean values of the hyperparameter distributions, which are used to estimate hyperparameters,

At each iteration step, after obtaining the HR image, the scale parameters are updated using (10). We now conclude our algorithm (Algorithm 1).

5. Experimental Results

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

5.1. General Description

We compare the performance of our proposed approach with the bicubic approach, the BEP model [11] based approach, the model [16] based approach, the TV model [14] based approach, and the SAR model [15] based approach in terms of Peak Signal-to-Noise Ratio (PSNR), Structural Similarity (SSIM) values [21], and thorough visual inspection. The images presented in Figure 1 were chosen for the simulation test. We used MATLAB R2009a on a 3.0 GHz Pentium Dual core computer with 4.0 GB RAM.

(a)

(b)

(a)
(b)

Figure 1 

The 128 × 128 test images: (a) “zebra” and (b) “car.”

We need to derive a set of LR images from original images in the experiments. These LR images should have been generated with subpixel motion, rotation, blurring, downsampling, and noise addition. In the experiments, we let ∈ (−3°, 3°), ∈ (−1, 1), ∈ (−1, 1), = 1, 2, 3, 4, 5, and the transformation of LR images be variant from each other. 3 × 3 average and Gaussian blurring operators with deviation equal to 1 were used in the simulation. The downsampling factors were = 2 and = 2. Finally, additive zero-mean white Gaussian noises were added to the LR images with Signal-to-Noise Ratio (SNR) levels of 1 dB, 5 dB, 10 dB, and 30 dB, respectively. Thus, in each level, five LR images were obtained from each of the original images. The initial parameters used in the experiments were set as follows: our proposed method: , , , , , and ; the method: , , , and ; the TV method: , ; the SAR method: with Laplacian operator ; the BEP method: , where is the scale parameter, and shift by and pixels in the horizontal and vertical directions, respectively, and these parameters are selected to obtain the best reconstructions; consider . For all the methods, a bicubic interpolation of is used as initial value , , the initial motion parameters were estimated using the method proposed in [22], , and , . In our work, the convergence of the algorithms was set with the same criterion .

5.2. Simulation Experiments

We used LR images derived from artificial HR images to test our proposed approach. The PSNR and SSIM are used to evaluate the reconstruction quality of different approaches. The first set of experiments was with the high level of noise (i.e., SNR = 1 dB and SNR = 5 dB). For the second set of experiments, white Gaussian noises with SNR = 10 dB and SNR = 30 dB were used and the motion parameters were estimated in each iteration.

In the first set of experiments, we will show the improvement of our proposed approach in heavy noise (i.e., SNR = 1 dB and SNR = 5 dB). Tables 1 and 2 show the PSNR values of all the estimated HR images in heavy noise with average blur and Gaussian blur. Based on the PSNR values, our proposed method shows a good performance. Moreover, Tables 1 and 2 give the corresponding SSIM values that also demonstrate the efficiency of our proposed method.

For the second set of experiments, the motion parameters were estimated. In order to further illustrate the effectiveness of our approach, we will present the reconstruction results for “zebra” and “car” images in noise with SNR = 10 dB and SNR = 30 dB, respectively. Tables 3 and 4 show the comparisons of PSNR values of bicubic approach, BEP approach, approach, SAR approach, TV approach, and our approach. Tables 3 and 4 also give the SSIM results of these approaches with average blur and Gaussian blur, respectively. From Tables 3 and 4, we can see the improvement of our proposed approach over other tested methods. Among all the approaches tested, the proposed method achieves the highest PSNR and SSIM values on both test images.

For visual quality comparison, the error images, that is, the difference between the estimated HR image and the original image, are shown in Figures 2 and 3. Our proposed approach can preserve edge details well. Figures 2(g)–2(l) and 3(g)–3(l) show the corresponding error images to the reconstructed images obtained with different approaches. Brighter pixels represent a large error. From these error images, the difference between different SR approaches is clearly observed. In general, our proposed approach can obtain the highest PSNR and SSIM values, with the best visual quality.

(a) Bicubic

(b) SAR

(c)

(d) TV

(e) BEP

(f) Ours

(g) Error image of bicubic interpolation

(h) Error image of SAR

(i) Error image of

(j) Error image of TV

(k) Error image of BEP

(l) Error image of ours

(a) Bicubic
(b) SAR
(c)
(d) TV
(e) BEP
(f) Ours
(g) Error image of bicubic interpolation
(h) Error image of SAR
(i) Error image of
(j) Error image of TV
(k) Error image of BEP
(l) Error image of ours

Figure 2 

Results obtained by applying different approaches to LR zebra images corrupted with Gaussian blur and white Gaussian noise with SNR = 10 dB and their corresponding error images. Brighter pixels represent a large error in error image.

(a) Bicubic

(b) SAR

(c)

(d) TV

(e) BEP

(f) Ours

(g) Error image of bicubic interpolation

(h) Error image of SAR

(i) Error image of

(j) Error image of TV

(k) Error image of BEP

(l) Error image of ours

(a) Bicubic
(b) SAR
(c)
(d) TV
(e) BEP
(f) Ours
(g) Error image of bicubic interpolation
(h) Error image of SAR
(i) Error image of
(j) Error image of TV
(k) Error image of BEP
(l) Error image of ours

Figure 3 

Results obtained by applying different approaches to LR car images corrupted with average blur and white Gaussian noise with SNR = 30 dB and their corresponding error images. Brighter pixels represent a large error in error image.