Efficient and Effective Total Variation Image Super-Resolution: A Preconditioned Operator Splitting Approach

Abstract

Super-resolution is a fusion process for reconstructing a high-resolution image from a set of low-resolution images. This paper proposes a novel approach to image super-resolution based on total variation (TV) regularization. We applied the Douglas-Rachford splitting technique to the constrained TV-based variational SR model which is separated into three subproblems that are easy to solve. Then, we derive an efficient and effective iterative scheme, which includes a fast iterative shrinkage/thresholding algorithm for denoising problem, a very simple noniterative algorithm for fusion part, and linear equation systems for deblurring process. Moreover, to speed up convergence, we provide an accelerated scheme based on precondition design of initial guess and forward-backward splitting technique which yields linear systems of equations with a nice structure. The proposed algorithm shares a remarkable simplicity together with a proven global rate of convergence which is significantly better than currently known lagged diffusivity fixed point iteration algorithm and fast decoupling algorithm by exploiting the alternating minimizing approach. Experimental results are presented to illustrate the effectiveness of the proposed algorithm.

1. Introduction

Multiframe image super-resolution (SR) is one of the promising techniques in image processing community since it enables us to obtain an image with a resolution that exceeds the hardware limitation, for example, the number of pixels in a charge-coupled device (CCD). Super-resolution is the process of combining a sequence of low-resolution (LR) noisy blurred images to produce a high-resolution (HR) image of sequence. It overcomes the inherent resolution limitation by bringing together the additional information from each image.

Generally, SR techniques can be divided into two broad categories: frequency domain methods and spatial domain methods. Most of the earlier SR work was developed in frequency domain using discrete Fourier transform (DFT), such as the work of Tsai and Huang [1], Kim et al. [2, 3], where high-frequency information is extracted from low-frequency data in the given LR frames. Many other popular frequency domain methods were proposed in Discrete cosine transform (DCT) domain and wavelet domain [4–7]. Although the frequency domain methods are intuitively simple and computationally cheap, they are extremely sensitive to model error [8], limiting their use. In a spatial domain alternative, Patti et al. [9] and Stark and Oskoui [10] proposed projection onto convex sets (POCSs) algorithm. A related method, the iterative back projection, was developed in [11, 12]. Although projection-based algorithms are usually robust to noise and allow some modeling flexibility, they are also known for their low rate of convergence. The hybrid maximum likelihood (ML)-maximum a posteriori (MAP)-POCS method was proposed in [13]. Based on these basic reconstruction methods, researchers have produced many extended algorithms, such as nonlocal-means (NLMs) based approach [14], multidimensional kernel regression-based approach [15], the joint formulation of reconstruction and registration [16–19], algorithms for multi-spectral and color [20, 21], hyperspectral [22], and compressed [23, 24] sequence.

The spatial domain methods discussed so far are generally confronted with the problem of slow convergence and expensive computation. To apply the SR algorithm to practical situations, many novel and powerful algebraic techniques have been proposed to reduce the computation complexity. For example, authors of [25–29] proposed efficient preconditioners to accelerate convergence of a conjugate gradient minimization algorithm. Xiao et al. [30] proposed an efficiency SR reconstruction algorithm employing the Armijo rule to identify the step length instead of the exact line search and replaced numerical approximation of the gradient of the MAP object function by analytic approximation. For pure translational motion and common space invariant blurring model, Elad and Hel-Or [31] proposed a novel SR algorithm that separates fusion and deblurring. The fusion method is a very simple noniterative algorithm, while preserving its optimality in ML sense. Farsiu et al. [32] proposed an efficient two-stage method for minimizing a novel framework combining a robust 𝑙₁ norm fidelity term and a bilateral prior, leading to an initial Median Shift-And-Add operation on Bayer-filtered LR data followed by a deblurring and interpolating stage. Huang et al. [33] proposed a fast decoupling algorithm by exploiting the alternating minimization approach.

In this paper, we propose a general framework for multiple shifted and linear space-invariant noisy blurred LR image frames which subsume several well-known SR models. The proposed model combines the TV regularization to formulate the SR image reconstruction as an optimization problem.

The purpose of this paper is to study an efficient TV-based SR reconstruction algorithm. There are two major contributions in this paper. As the first contribution, we propose an efficient algorithm that takes full advantage of the problem structures; that is, geometrical motion matrices, blur matrix, and the first-order finite-difference matrix all have block-circulant-circulant-block (BCCB) structure under periodic boundary condition. As such, we propose to compute the minimizer of our SR model by applying Douglas-Rachford splitting (DBS) techniques, respectively, alternating direction methods of multipliers (ADMM), which separate the SR model into three subproblems that are easy to solve. As the second contribution, we provide an accelerated scheme based on precondition design of initial value and forward-backward splitting (FBS) to speed up convergence. Our method can separate the SR treatment into measurements fusion, denoising, and deblurring. The fusion part is shown to be a very simple noniterative algorithm. The denoising problem can be solved by a linear time shrinkage operation. Fast Fourier transform (FFT) is employed to solve the deblurring problem. Finally, experimental results are presented to illustrate the effectiveness of the proposed algorithm.

The outline of this paper is as follows. In Section 2, we present the image observation model of the SR problem, then propose a TV-based SR model. In Section 3, we present an efficient SR reconstruction algorithm. Experimental results are provided in Section 4. Finally, concluding remarks are given in Section 5.

2. Problem Formulation

2.1. Observation Model

The image observation model is employed to relate the desired referenced HR image to all the observed LR images. Consider the desired HR image of size 𝐿₁𝑁₁×𝐿₂𝑁₂ written in lexicographical notation as the vector 𝑧=[𝑧₁,𝑧₂,…,𝑧_𝑁]^𝑇, where 𝑁=𝐿₁𝑁₁×𝐿₂𝑁₂. Let the parameters 𝐿₁ and 𝐿₂ be the subsampling factors in the horizontal and vertical directions, respectively; each observed LR image has the size 𝑁₁×𝑁₂. Thus, the LR image can be represented as 𝑦_𝑘=[𝑦_𝑘1,𝑦_𝑘2,…,𝑦_𝑘𝑀]^𝑇, for 𝑘=1,2,…,𝐾 and 𝑀=𝑁₁×𝑁₂. A popular model assumes that LR images {𝑦_𝑘}^𝐾_𝑘=1 are generated from HR image 𝑧 through a sequence of operations that includes (i) geometrical motions 𝑀_𝑘, (ii) a linear space-invariant blur 𝐵, (iii) a subsampling step represented by 𝑆, and finally (iv) an additive white Gaussian noise 𝑛_𝑘 with zero mean that represents both measurements noise and model mismatch [32]. All these are linear operators, represented by a matrix multiplying the image they operate on. We assume hereafter that 𝐵 and 𝑆 are identical for all images in the sequence. This model leads to the following set of equations, where all images are ordered lexicographically:𝑦_𝑘=𝑆𝐵𝑀_𝑘𝑧+𝑛_𝑘∶=𝑊_𝑘𝑧+𝑛_𝑘for𝑘=1,2,…,𝐾,(2.1) where 𝑊_𝑘=𝑆𝐵𝑀_𝑘 represents the imaging system.

The recovery of 𝑧 from {𝑦_𝑘}^𝐾_𝑘=1 is thus an inverse problem, combining motion compensation, denoising, deblurring, scaling-up operation, and fusion of the different images, which all merged to one. The quality of the desired SR image depends on the assumption that 𝑆,𝐵, and 𝑀_𝑘 are known, or the accuracy in estimating the degraded operators. Throughout this paper, we assume that 𝑆,𝐵, and 𝑀_𝑘 are known. The decimation 𝑆 is dependent on the resolution scale factor that we aim to achieve, and, as such, it is easily fixed. In this work, we shall assume that this resolution factor is an integer 𝑠⩾1 on both axes. In most cases, the blur 𝐵 refers to the camera point spread function (PSF), and, therefore, it is also accessible. Even if this is not the case, the blurring is typically dependent on few parameters, and those, in the worst case, can be manually set. To be identical with the work of Elad and Farsiu [31, 32], we focus on the simplest of the motion models, namely, the translational model. Reference [34] detailed the several reasons for this. We believe that an in-depth study of this simple case allows much insight to be gained about the problems inherent to SR image reconstruction.

2.2. TV Regularization-Based SR Model

In general, SR is an ill-posed problem either because the information contained in the observed LR images is not sufficient or because it has great sensitivity to the noise. Procedures adopted to stabilize the inversion of ill-posed problem are called regularization. In such stabilization scheme, we reconstruct the original HR image by finding the minimizer of some appropriate functional𝐸(𝑧)∶=Φ(𝑧)+𝜇Ψ(𝑧),(2.2) where Φ(𝑧) is a regularization term which includes prior information about the original image, Ψ(𝑧) denotes the data fitting term depending on the given LR image {𝑦_𝑘}^𝐾_𝑘=1, and 𝜇 is a positive parameter which controls the tradeoff between the two terms for minimization.

In general, the data fitting term can be deduced from an MAP estimation. If 𝑧 is corrupted by additive white Gaussian noise, the MAP estimation will lead to the data fitting term Ψ(𝑧)=(1/2)∑^𝐾_𝑘=1‖𝑊_𝑘𝑧−𝑦_𝑘‖²_{𝐿2(Ω𝐿)}, where Ω_𝐿⊂Ω and Ω denote a bounded and open domain of continuous LR image and HR image in 𝑅², respectively. For regularization term, the popular choice is total variation seminorm ‖𝑧‖_TV=∫_Ω|∇𝑧|d𝑥, which was first proposed for image denoising [35], because TV norm can better preserve sharp edges or object boundaries that are usually the most important features to recover.

As stated previously, to invert the degradation process in (2.1), we can formulate a TV regularization model which requires solving the variational problem:min_{𝑧∈BV(Ω)}⎧⎨⎩𝐸(𝑧)=‖𝑧‖_TV+𝜇2^𝐾_𝑘=1‖‖𝑊_𝑘𝑧−𝑦_𝑘‖‖²_{𝐿2Ω𝐿}⎫⎬⎭,(2.3) where BV(Ω) is the space of functions of bounded variation. Note that the fitting term in (2.3) is strictly convex and coercive and the TV regularity term is also convex (though not strictly so) and lower semicontinuous. So the objective function 𝐸(𝑧) is globally strictly convex and possesses a unique minimizer. In terms of optimization, these are desirable properties.

3. An Operator Splitting Approach to TV-Based Super-Resolution

To solve the desired HR image of (2.3), commonly used method is the gradient descent method [19–32, 36]. Although this approach is simple, the nonlinearly and poor conditioning of the problem make this approach very slow. A more efficient class of solvers are those based on a linearized gradient method which solves the associated Euler-Lagrange equation via a lagged diffusivity fixed-point iteration [28, 29]. In each iteration of the linearized gradient method, a linear system needs to be solved, which becomes more and more difficult as 𝐵 becomes more ill-conditioned. Another group of algorithms is based on the well-known variable-splitting and penalty techniques in optimization. These ideas have gained wide application in image processing, such as works in [37–41]. Recently, Huang et al. [33] modified the SR model (2.3) by adding a quadratic term to get a simpler alternating minimization algorithm. The drawback of this method is the same as the lagged diffusivity fixed point method.

To efficiently solve the SR problem (2.3), in this section, we will show that the operator splitting method can be used to divide the problem (2.3) into subproblems that can be solved in sequence, and each of them permits a closed form solution. Among the current splitting methods, the most prominent splitting schemes are forward-backward splitting, double-backward splitting, Peaceman-Rachford splitting, and Douglas-Rachford splitting. In this paper, we will focus on the forward-backward splitting and Douglas-Rachford splitting. One may refer to [42, 43] and the references therein for more details.

3.1. FBS and DRS

Let 𝐻 be a real Hilbert space, and let 𝐴,𝐵∶𝐻→2^𝐻 be two set-valued operators. We assume that 𝐴 and 𝐵 are maximal monotone; that is, their resolvents 𝐽_𝐴∶=(𝐼+𝐴)⁻¹ and 𝐽_𝐵∶=(𝐼+𝐵)⁻¹ exist and are firmly nonexpansive. The problem which we will describe as a fundamental problem can be written in the form of a common zero inclusion problem0∈𝐴(𝑥̂)+𝐵(𝑥̂).(3.1)

The idea of the forward-backward splitting algorithm is that, for any constant 𝛾>0, we have0∈𝐴(𝑥̂)+𝐵(𝑥̂)⟺𝑥̂−𝛾𝐵(𝑥̂)∈𝑥̂+𝛾𝐴(𝑥̂)⟺𝑥̂∈𝐽_𝛾𝐴(𝐼−𝛾𝐵)(𝑥̂).(3.2) This leads to the following result.

Theorem 3.1 (FBS [43, Theorem 2.3.17]). Suppose that 𝐴∶𝐻→2^𝐻 is maximal monotone and 𝐵∶𝐻→𝐻 is a monotone operator such that 𝜂𝐵 is firmly nonexpansive for some 𝜂>0. Furthermore, assume that a solution of (3.1) exists. Then, for every start element 𝑥⁰ and step size 𝛾∈(0,2𝜂), the forward-backward splitting algorithm 𝑥^𝑛+1∈𝐽_𝛾𝐴(𝐼−𝛾𝐵)(𝑥^𝑛)(3.3) converges weakly to an element of the set of solutions (𝐴+𝐵)⁻¹(0).

We now describe the Douglas-Rachford splitting scheme, which does exhibit general convergence, at least when used with a constant step size in finite-dimensional spaces. To introduce it, we can rewrite the fixed point relation (3.2) as follows:0∈𝐴(𝑥̂)+𝐵(𝑥̂)⟺𝑥̂∈𝐽_𝛾𝐴(𝐼−𝛾𝐵)(𝑥̂)⟺𝑥̂+𝛾𝐵(𝑥̂)∈𝐽_𝛾𝐴(𝐼−𝛾𝐵)(𝑥̂)+𝛾𝐵(𝑥̂)⟺𝑥̂∈𝐽_𝛾𝐵𝐽_𝛾𝐴(𝐼−𝛾𝐵)(𝑥̂)+𝛾𝐵(𝑥̂).(3.4) This leads to the following result.

Theorem 3.2 (DBS [43, Theorem 2.3.21]). Let 𝐴,𝐵∶𝐻→2^𝐻 be maximal monotone, and assume that a solution of (3.1) exists. Then, for any elements 𝑦⁰ and 𝑥⁰ and any step size 𝛾>0, the following Douglas-Rachford splitting algorithm converges weakly to an element 𝑦̂: 𝑦^𝑛+1=𝐽_𝛾A(2𝑥^𝑛−𝑦^𝑛)+𝑦^𝑛−𝑥^𝑛,𝑥^𝑛+1=𝐽_𝛾𝐵𝑦^𝑛+1.(3.5) Furthermore, it holds that 𝑥̂=𝐽_𝛾𝐵(𝑦)̂ satisfies 0∈𝐴(𝑥̂)+𝐵(𝑥̂). If 𝐻 is finite-dimensional, then the sequence {𝑥^𝑛} converges to 𝑥̂.

3.2. Proposed Algorithm

In this subsection, we will apply the DRS to dual problem of the minimization functional (2.3). We first rewrite the energy functional (2.3) in a discrete formmin_{𝑧∈𝑅𝐿1𝑁1×𝐿2𝑁2}⎧⎪⎨⎪⎩𝐸(𝑧)=𝜙(𝐷𝑧)_{Φ(𝑧)}+𝜇Ψ(𝑧)⎫⎪⎬⎪⎭(3.6) with Ψ(𝑧)=(1/2)∑^𝐾_𝑘=1‖𝑊_𝑘𝑧−𝑦_𝑘‖², 𝜙(𝐷𝑧)=∑_𝑖,𝑗‖(𝐷𝑧)_𝑖,𝑗‖ is the discrete total variation of 𝑧. Here, ‖⋅‖ denotes Euclidean norm, and 𝐷 is given by (𝐷𝑧)_𝑖,𝑗=𝐷⁺_𝑥𝑧_𝑖,𝑗;𝐷⁺_𝑦𝑧_𝑖,𝑗=𝑧_𝑖,𝑗+1−𝑧_𝑖,𝑗;𝑧_𝑖+1,𝑗−𝑧_𝑖,𝑗,for1≤𝑖≤𝐿₁𝑁₁,1≤𝑗⩽𝐿₂𝑁₂,(3.7) where 𝐷⁺_𝑥 and 𝐷⁺_𝑦 are forward difference operators with periodic boundary condition 𝑧_{𝑖,𝐿2𝑁2+1}=𝑧_𝑖,1 and 𝑧_{𝐿1𝑁1+1,𝑗}=𝑧_1,𝑗. Therefore, 𝐷 is a BCCB matrix.

Consider the equivalent problem of (3.6)min_{𝑧∈𝑅𝐿1𝑁1×𝐿2𝑁2𝑤∈𝑅2𝐿1𝑁1×𝐿2𝑁2}{𝐸(𝑧)=𝜙(𝑤)+𝜇Ψ(𝑧)},s.t.𝑤=𝐷𝑧.(3.8) The Lagrangian for problem (3.8) is𝐿(𝑤,𝑧,𝜆)=𝜙(𝑤)+𝜇Ψ(𝑧)−⟨𝜆,𝑤−𝐷𝑧⟩,(3.9) where the dual variable 𝜆∈𝑅^{2𝐿1𝑁1×𝐿2𝑁2} can be thought of as a vector of Lagrange multipliers. Therefore, the dual problem of (3.8) ismax_𝜆inf_𝑤,𝑧𝐿(𝑤,𝑧,𝜆)=−min_𝜆𝜙^∗(𝜆)+𝜇Ψ^∗−𝜇⁻¹𝐷^∗𝜆,(3.10) where 𝜙^∗(resp., Ψ^∗) denotes conjugate function of 𝜙 (resp., Ψ).

Define operators 𝐴 and 𝐵 by𝐴(𝜆)=𝜕𝜙^∗(𝜆),𝐵(𝜆)=𝜇𝜕Ψ^∗−𝜇⁻¹𝐷^∗𝜆.(3.11) By Fermat's rule, solving the dual problem (3.10) is equivalent to finding 𝜆 such that0∈𝐴(𝜆)+𝐵(𝜆).(3.12)

By formally applying DRS (3.5) to (3.12) with 𝛼 as the step size, respectively, the ADMM iterations [44–51] are given by 𝑤^𝑛+1=argmin_𝑤𝜙(𝑤)−⟨𝜆^𝑛,𝑤⟩+𝛼2‖𝑤−𝐷𝑧^𝑛‖²,(3.13)𝑧^𝑛+1=argmin_𝑧𝜇Ψ(𝑧)+⟨𝜆^𝑛,𝐷𝑧⟩+𝛼2‖‖𝑤^𝑛+1−𝐷𝑧‖‖²,(3.14)𝜆^𝑛+1=𝜆^𝑛+𝛼𝐷𝑧^𝑛+1−𝑤^𝑛+1.(3.15) As pointed out by Setzer in [43, 49], we note that this iteration scheme coincides with DRS algorithm (3.5) with 𝑦^𝑛=𝛾(𝛼𝜆^𝑛+𝐷𝑧^𝑛), 𝑥^𝑛=𝛾𝛼𝜆^𝑛, and 𝛾=𝛼. Since operators 𝜙^∗ and Ψ^∗ are proper, lower semicontinuous, and convex, the operators 𝐴 and 𝐵 are maximal monotone see [43]. According to Theorem 3.2, the sequence {𝑧^𝑛} converges to solution of (3.6). We also note that the above ADMM algorithm coincides with the alternating split Bregman algorithm proposed by Goldstein and Osher [40] with 𝜕𝜙(𝑤^𝑛+1)=𝜆^𝑛+1 and 𝜕Ψ(𝑧^𝑛+1)=(−1/𝜇)𝜆^𝑛+1.

3.2.1. 𝑤-Subproblem

It is not difficult to show that the minimization of (3.13) with respect to 𝑤 is equivalent to solving 𝐿₁𝑁₁×𝐿₂𝑁₂ two-dimensional problem of the formmin_{𝑤𝑖,𝑗}⎧⎨⎩‖‖𝑤_𝑖,𝑗‖‖+𝛼2‖‖‖𝑤_𝑖,𝑗−⎛⎜⎝(𝐷𝑧^𝑛)_𝑖,𝑗+𝜆^𝑛_𝑖,𝑗𝛼⎞⎟⎠‖‖‖²⎫⎬⎭,(3.16) for which the unique minimizer is given by the following two-dimensional shrinkage formula:𝑤^𝑛+1_𝑖,𝑗=max⎧⎨⎩‖‖‖(𝐷𝑧^𝑛)_𝑖,𝑗+𝜆^𝑛_𝑖,𝑗𝛼‖‖‖−1𝛼,0⎫⎬⎭(𝐷𝑧^𝑛)_𝑖,𝑗+𝜆^𝑛_𝑖,𝑗/𝛼‖‖(𝐷𝑧^𝑛)_𝑖,𝑗+𝜆^𝑛_𝑖,𝑗/𝛼‖‖,(3.17) where the convention 0⋅(0/0)=0 is followed. Here, shrinkage formula (3.17) which serves as nonlinear low-pass filter to restored HR image is tantamount to denoising treatment.

3.2.2. 𝑧-Subproblem

Subproblem (3.14) is quadratic in 𝑧, and the minimizer 𝑧^𝑛+1 is given by the normal equations⎛⎜⎝𝐷^𝑇𝐷+𝜇𝛼_𝑘𝑊^𝑇_𝑘𝑊_𝑘⎞⎟⎠𝑧^𝑛+1=𝐷^𝑇𝑤^𝑛+1−𝜆^𝑛𝛼+𝜇𝛼_𝑘𝑊^𝑇_𝑘𝑦_𝑘,(3.18) where 𝑊^𝑇_𝑘 is the transpose operator of 𝑊_𝑘.

Note that 𝐷^𝑇𝐷 in (3.18) is BCCB matrix and can be diagonalized by FFT. Moreover, with periodic boundary condition, both 𝐵 and 𝑀_𝑘 are BCCB matrices. Exploiting the fact that the product order of two BCCB matrices can commute, we get that 𝑀^𝑇_𝑘𝐵^𝑇=𝐵^𝑇𝑀^𝑇_𝑘, 𝐵𝑀_𝑘=𝑀_𝑘𝐵. Then, ∑_𝑘𝑊^𝑇_𝑘𝑊_𝑘=𝐵^𝑇𝑅𝐵 (𝑅∶=∑_𝑘𝑀_𝑘^𝑇𝑆^𝑇𝑆𝑀_𝑘 is a diagonal matrix [31]). However, ∑_𝑘𝑊^𝑇_𝑘𝑊_𝑘=𝐵^𝑇𝑅𝐵 does not have BCCB structure. Therefore, it does not allow us to apply FFT implementation to (3.18) directly as done in [37, 38]. The quadratic term Ψ(𝑧)=(1/2)∑^𝐾_𝑘=1‖𝑊_𝑘𝑧−𝑦_𝑘‖² in (3.14) that couples the variable 𝑧 by the matrix 𝑊_𝑘 makes the algorithm computationally expensive. There are some techniques to overcome these problems. In [50], the authors introduced the three constrains: 𝑤₁=𝐵𝑧, 𝑤₂=𝐷𝑧, 𝑤₃=𝑧 and used the alternating split Bregman technique to maximally decouple the variables. In [51], the linear system was solved noniteratively by using Sherman-Morrison-Woodbury (SMW) matrix inversion formula and FFT to diagonalize the Hessian matrix of the energy functional. In this paper, we use the forward-backward splitting method (3.3) to efficiently solve the 𝑧-subproblem (3.14), which can decouple the variable 𝑧 and the constraint matrix 𝑊_𝑘 and make the Hessian matrix of the energy functional have BCCB structure.

Let 𝑓₁(𝑧)=(𝛼/2𝜇)‖𝑤^𝑛+1−𝐷𝑧‖²+⟨𝜆^𝑛/𝜇,𝐷𝑧⟩ and 𝑓₂(𝑧)=(1/2)∑^𝐾_𝑘=1‖𝑊_𝑘𝑧−𝑦_𝑘‖². Then, ∇𝑓₂(𝑧^𝑛)=𝐵^𝑇(𝑅𝐵𝑧^𝑛−𝑠) with 𝑠∶=∑_𝑘𝑀_𝑘^𝑇𝑆^𝑇𝑦_𝑘. Define operators 𝐴=𝜕𝑓₁ and 𝐵=∇𝑓₂; the FBS method (3.3) with 𝛾 as the step size applied to (3.14) leads to the iterative scheme𝑧^𝑛+1=argmin_𝑧𝛼2‖‖𝑤^𝑛+1−𝐷𝑧‖‖²+⟨𝜆^𝑛,𝐷𝑧⟩+𝜇2𝛾‖‖𝑧−𝑧^𝑛+𝛾∇𝑓₂(𝑧^𝑛)‖‖².(3.19) The minimizer 𝑧^𝑛+1 is given by the normal equations𝐷^𝑇𝐷+𝜇𝛾𝛼𝐼𝑧^𝑛+1=𝐷^𝑇𝑤^𝑛+1−𝜆^𝑛𝛼+𝜇𝛾𝛼𝑧^𝑛−𝛾∇𝑓₂(𝑧^𝑛).(3.20) Under the periodic boundary condition, (3.20) can be solved by two FFTs, which simultaneously performs the LR measurements fusion and deblurring treatment. In the next subsection, the fusion part will be shown to be a simple noniterative.

Notice that 𝐵=∇𝑓₂ is Lipschitz continuous with Lipschitz constant 𝛽=‖𝑅‖‖𝐵‖², where ‖𝐵‖ (resp., ‖𝑅‖) is the matrix norm of 𝐵 (resp., 𝑅). From [41,Theorem 2.3.19], (1/𝛽)𝐵 is firmly nonexpansive. Then, the convergence of 𝑧^𝑛 is ensured by Theorem 3.1, while 𝛾∈]0,2/𝛽[.

3.2.3. Optimal Initial Guess

Firstly, borrowed the idea of [31], we take initial data 𝑧⁰ from maximum likelihood (ML) estimation; that is, 𝑧⁰ is the minimizer of the optimization problem as follows:𝑧̂=argmin_𝑧⎧⎨⎩_𝑘‖‖𝑆𝐵𝑀_𝑘𝑧−𝑦_𝑘‖‖²⎫⎬⎭.(3.21) The gradient descent algorithm suggests the following iterative equation for the solution of (3.21):𝑧̂^𝑛+1=𝑧^𝑛̂−Δ𝑡_𝑘𝑀^𝑇_𝑘𝐵^𝑇𝑆^𝑇𝑆𝐵𝑀_𝑘𝑧^𝑛̂−𝑦_𝑘=𝑧^𝑛̂−Δ𝑡𝐵^𝑇_𝑘𝑀^𝑇_𝑘𝑆^𝑇𝑆𝑀_𝑘𝐵𝑧^𝑛̂−𝑦_𝑘,(3.22) where Δ𝑡 denotes step size. Let us define the blurred super-resolution image by 𝑢^𝑛̂=𝐵𝑧^𝑛̂. Multiplying both sides of (3.22) with 𝐵, we get𝑢̂^𝑛+1=𝑢^𝑛̂+Δ𝑡𝐵𝐵^𝑇𝑠−𝑅𝑢^𝑛̂,(3.23) where 𝑠=∑_𝑘𝑀^𝑇_𝑘𝑆^𝑇𝑦_𝑘 and 𝑅=∑_𝑘𝑀^𝑇_𝑘𝑆^𝑇𝑆𝑀_𝑘. According to [52], the steady state solution of (3.23) is given by𝑢_∞̂=𝑅⁻¹𝑠.(3.24) We note that 𝑀^𝑇_𝑘𝑆^𝑇 copies the values from the LR grid to the HR grid after proper shifting and zero filling and 𝑆𝑀_𝑘 copies a selected set of pixels in HR grid back on the LR grid (Figure 1 illustrates the effect of upsampling and downsampling matrices 𝑆^𝑇 and 𝑆). Neither of these two operations changes the pixel value. It is easy to show that 𝑅 is a diagonal matrix. Each diagonal entry in 𝑅 corresponds to one pixel in the super-resolution image. Its value is a nonnegative integer, counting the number of measurements contributing to it. The fusion image 𝑠 is simply the addition of the measurements after proper zero-filling interpolation and motion compensation. Thus, 𝑢_∞̂=𝑅⁻¹𝑠 is none other than the pixel-wise average of the measurement. Therefore, the noise of 𝑢_∞̂ is reduced due to the averaging.

Figure 1 

Effect of 𝑆𝑇 and 𝑆 on image (resolution enhancement factor is three).

Because 𝑢_∞̂=𝐵𝑧_∞̂, Wiener filter was applied to 𝑢_∞̂. Then, the restoration image 𝑧_∞̂ is taken as the initial data 𝑧⁰.

As stated previously, precondition design procedure of the initial data 𝑧⁰ is summed up in Algorithm 1 as follows:

Input: 𝑠=∑𝑘𝑀𝑇𝑘𝑆𝑇𝑦𝑘,𝑅=∑𝑘𝑀𝑇𝑘𝑆𝑇𝑆𝑀𝑘,𝐵. Output: 𝑧0=𝐵−1𝑅−1𝑠.

Algorithm 1 

3.2.4. Algorithm Description

To sum up the above arguments, the complete resulting algorithm is summarized in Algorithm 2 as follows.

Input: 𝑠=∑𝑘𝑀𝑇𝑘𝑆𝑇𝑦𝑘,𝑅=∑𝑘𝑀𝑇𝑘𝑆𝑇𝑆𝑀𝑘,𝐵,tolerance𝜁stop,𝛾∈]0,2/𝛽[, and 𝛼,𝜇>0. Initialize: 𝑧0=𝐵−1𝑅−1𝑠,𝜆0=0, and 𝑛=0. Do  (1) Compute 𝑤𝑛+1 according to (3.17) for fixed (𝑧𝑛,𝜆𝑛).  (2) Compute 𝑧𝑛+1 according to (3.20) for fixed (𝑤𝑛+1,𝜆𝑛).  (3) Update 𝜆𝑛+1 according to (3.15).  (4)𝑛=𝑛+1. Until ‖𝑧𝑛+1−𝑧𝑛‖/‖𝑧𝑛+1‖≤𝜁stop.

Algorithm 2 

3.2.5. Some Complexity Notes

It is clear that the complexity of the proposed algorithm mainly includes three parts. The calculation in (3.17) and (3.15) have linear-time complexity of order 𝑂(𝑁²) for an 𝑁×𝑁 image. Hence, the 𝑤-subproblem (3.16) and 𝜆 can be solved quickly. The solution of the 𝑧-subproblem (3.20) requires two FFTs (including one inverse FFT), which has a total complexity in the order of 𝑂(𝑁²log(𝑁²))=𝑂(𝑁²log(𝑁)).

4. Experimental Results

In this section, we present some experimental examples to demonstrate the performance of our method. We use the three 256×256 test images (“Lena” Figure 2(b), “Cameraman” Figure 4(b), and “Fingerprint” Figure 8(b)) for the synthetic test. A sequence of LR frames of 64×64 pixels from the original image is generated as follows. First, the original image was shifted by one pixel in the vertical direction. Then, to simulate the effect of camera PSF, this shifted image was convolved with a Gaussian blur kernel. The resulting image was subsampled by the factor 4 in each direction. The same approach with different motion vectors in the vertical and horizontal directions was used to produce eight LR images from the original scene. The resulting LR frames were corrupted with white Gaussian noise. All experiments were performed under Windows XP and MATLAB v7.1 running on a desktop with an Intel Core Dual Processor 3.00 GHz and 4.00 GB of memory.

(a)

(b)

(c)

(d)

(e)

Figure 2 

Results of different resolution enhancement methods applied to “Lena” image degraded by Gaussian blur with support size 3×3, standard deviation 0.5, and Gaussian noise with variance 5. (a) One LR image; (b) original image; reconstructed image (c) using LDFP method, (d) using AMM method, and (e) using the proposed method (𝜇=130).

For the objective comparison between the original HR and SR reconstructed images, we measure the peak signal-to-noise ratio (PSNR) and the relative error (ReErr) defined as PSNR=10log₁₀𝐿₁𝑁₁×𝐿₂𝑁₂max{𝑧}²‖𝑧^∗−𝑧‖²,ReErr=‖𝑧^∗−𝑧‖²‖𝑧‖²,(4.1) where 𝑧 and 𝑧^∗ are the original and the SR reconstructed images, respectively, and 𝐿₁𝑁₁×𝐿₂𝑁₂ represents the image size.

We compare the proposed method (operator splitting method, OSM) with the lagged diffusivity fixed point iteration (LDFP) [28, 29] and alternating minimization method (AMM) [33]. In all the tests, we set the initial guess 𝑧⁰=𝐵⁻¹𝑅⁻¹𝑠. The choice of parameters in three methods all base on the tradeoff between reconstruction effect and computing time. In the proposed method, the value of 𝛾 and 𝛼 are fixed to be 1.5 and 4, respectively. The stopping criterion of all the methods is that the relative difference (ReDiff) between the successive iterative of the SR reconstructed image should satisfy the following inequality: ReDiff=‖‖𝑧^𝑛+1−𝑧^𝑛‖‖‖‖𝑧^𝑛+1‖‖≤10⁻⁴.(4.2)

In the first test, we apply Gaussian kernel with window size 3×3, standard deviation 0.5, and different noise level (noise variance 𝜎²=5,20). One of LR frames is presented in Figures 2–5, 8(a), respectively. The corresponding reconstructed images of the three methods are shown in Figures 2, 4, 8(c)–8(e) and Figures 3, 5(b)–5(d), respectively. Our second test uses Gaussian kernel with support size 15×15, standard deviation 1.7, and noise variance 5. One of LR images is presented in Figures 6, 7(a), respectively. The corresponding reconstructed images by the three methods are shown in Figures 6, 7(b)–7(d). We can see theses SR reconstructed images by different methods are very similar in real visual perception. In Table 1, we compare their reconstruction performances in PSNRs and ReErrs. On one hand, we see from the table that both PSNRs and ReErrs of the reconstructed images by the proposed method are better than those by the LDFP and AMM method. On the other hand, it is clear from Table 1 that the proposed method is more efficient (in iterations and computation times) than the other two methods.

Table 1 

The PSNR, ReErr, number of iterations, and computational times of the reconstructed images using three methods. The numbers in the bracket in the “Iterations” row refer to the total number of inner conjugate gradient iterations.

(a)

(b)

(c)

(d)

Figure 3 

Results of different resolution enhancement methods applied to “Lena” image degraded by Gaussian blur with support size 3×3, standard deviation 0.5, and Gaussian noise with variance 20. (a) One LR image; reconstructed image (b) using LDFP method, (c) using AMM method, and (d) using the proposed method (𝜇=130).

(a)

(b)

(c)

(d)

(e)

Figure 4 

Results of different resolution enhancement methods applied to “Cameraman” image degraded by Gaussian blur with support size 3×3, standard deviation 0.5, and Gaussian noise with variance 5. (a) One LR image; (b) original image; reconstructed image (c) using LDFP method, (d) using AMM method, and (e) using the proposed method (𝜇=130).

(a)

(b)

(c)

(d)

Figure 5 

Results of different resolution enhancement methods applied to “Cameraman” image degraded by Gaussian blur with support size 3×3, standard deviation 0.5, and Gaussian noise with variance 20. (a) One LR image; reconstructed image (b) using LDFP method, (c) using AMM method, and (d) using the proposed method (𝜇=130).

(a)

(b)

(c)

(d)

Figure 6 

Results of different resolution enhancement methods applied to “Lena” image degraded by Gaussian blur with support size 15×15, standard deviation 1.7, and Gaussian noise with variance 5. (a) One LR image; reconstructed image (b) using LDFP method, (c) using AMM method, and (d) using the proposed method (𝜇=500).

(a)

(b)

(c)

(d)

Figure 7 

Results of different resolution enhancement methods applied to “Cameraman” image degraded by Gaussian blur with support size 15×15, standard deviation 1.7, and Gaussian noise with variance 5. (a) one LR image; reconstructed image (b) using LDFP method, (c) using AMM method, and (d) using the proposed method (𝜇=500).

(a)

(b)

(c)

(d)

(e)

Figure 8 

Results of different resolution enhancement methods applied to “Fingerprint” image degraded by Gaussian blur with support size 3×3, standard deviation 0.5, and Gaussian noise with variance 5. (a) One LR image, (b) original image; reconstructed image (c) using LDFP method, (d) using AMM method, and (e) using the proposed method (𝜇=130).

In Figure 9, we show that the convergence of the proposed method is faster than LDFP and AMM methods. The 𝑥-axis refers to the number of iterations. The 𝑦-axis in Figures 9(a), 9(b) refers to the PSNR and the 𝑦-axis in Figures 9(c), 9(d) refers to the relative difference between the successive iteration of the reconstructed image. These figures show that the proposed method can provide good quality of reconstructed images in an efficient manner.

(a)

(b)

(c)

(d)

Figure 9 

Convergence performance of three methods in the Gaussian blur with support size 3×3, standard deviation 0.5, and Gaussian noise with variance 5 case: measured by PSNR value of the SR reconstructed (a) “Lena” image, (b) “Cameraman” image, measured by ReDiff value of the SR reconstructed (c) “Lena” image, (d) “Cameraman” image.

5. Conclusion

This paper proposes a general framework for multiple shifted and linear space-invariant blurred LR image frames which subsume several well-known SR models. The proposed model combines total variation (TV) regularization to formulate the SR image reconstruction as an optimization problem. Then, we propose an efficient algorithm that takes full advantage of the problem structures. As such, we propose to compute the minimizer of our SR model by applying DRS techniques (resp., ADMM) which separated the SR model into three subproblems that can be easily solved. Moreover, to speed up convergence, we provide an accelerated scheme based on precondition design of initial value and FBS. The proposed algorithm reduces the computational complexity. The good performance of the proposed explicit algorithm has been tested for synthetic data sets of several images degraded with Gaussian blur and contaminated with Gaussian white noise. Numerical results indicate that the algorithm recovers well edges and small features not appearing in the original degraded images. The experimental results indicate that the proposed algorithm has considerable effectiveness in terms of both objective measurements and visual evaluation.

Acknowledgments

This work was supported in part by the Natural Science Foundation of China under Grant no. 60802039, by the Doctoral Foundation of Ministry of Education of China under Grant no. 20070288050, NUST Research Funding under Grant no. 2010ZDJH07, and also sponsored by Natural Science Foundation of Jiangsu under Grant no. BK2010488 and “Qing Lan Project” of Jiangsu Province. The authors would like to express their gratitude to the anonymous referees for making helpful and constructive suggestions.