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

Authors: Li-Li Huang; Liang Xiao; Zhi-Hui Wei
Journal: Mathematical Problems in Engineering
 (2011)

---

## 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.


---

## Body

## 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 l1 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 Section2, 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 sizeL1N1×L2N2 written in lexicographical notation as the vector z=[z1,z2,…,zN]T, where N=L1N1×L2N2. Let the parameters L1 and L2 be the subsampling factors in the horizontal and vertical directions, respectively; each observed LR image has the size N1×N2. Thus, the LR image can be represented as yk=[yk1,yk2,…,ykM]T, for k=1,2,…,K and M=N1×N2. A popular model assumes that LR images {yk}k=1K are generated from HR image z through a sequence of operations that includes (i) geometrical motions Mk, (ii) a linear space-invariant blur B, (iii) a subsampling step represented by S, and finally (iv) an additive white Gaussian noise nk 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 B and S are identical for all images in the sequence. This model leads to the following set of equations, where all images are ordered lexicographically:(2.1)yk=SBMkz+nk:=Wkz+nkfork=1,2,…,K,
where Wk=SBMk represents the imaging system.The recovery ofz from {yk}k=1K 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 S,B, and Mk are known, or the accuracy in estimating the degraded operators. Throughout this paper, we assume that S,B, and Mk are known. The decimation S 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 s⩾1 on both axes. In most cases, the blur B 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)E(z):=Φ(z)+μΨ(z),
where Φ(z) is a regularization term which includes prior information about the original image, Ψ(z) denotes the data fitting term depending on the given LR image {yk}k=1K, 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. Ifz is corrupted by additive white Gaussian noise, the MAP estimation will lead to the data fitting term Ψ(z)=(1/2)∑k=1K∥Wkz-yk∥L2(ΩL)2, where ΩL⊂Ω and Ω denote a bounded and open domain of continuous LR image and HR image in R2, respectively. For regularization term, the popular choice is total variation seminorm ∥z∥TV=∫Ω|∇z|dx, 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:(2.3)min⁡z∈BV(Ω){E(z)=‖z‖TV+μ2∑k=1K‖Wkz-yk‖L2(ΩL)2},
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 E(z) is globally strictly convex and possesses a unique minimizer. In terms of optimization, these are desirable properties.


## 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 sizeL1N1×L2N2 written in lexicographical notation as the vector z=[z1,z2,…,zN]T, where N=L1N1×L2N2. Let the parameters L1 and L2 be the subsampling factors in the horizontal and vertical directions, respectively; each observed LR image has the size N1×N2. Thus, the LR image can be represented as yk=[yk1,yk2,…,ykM]T, for k=1,2,…,K and M=N1×N2. A popular model assumes that LR images {yk}k=1K are generated from HR image z through a sequence of operations that includes (i) geometrical motions Mk, (ii) a linear space-invariant blur B, (iii) a subsampling step represented by S, and finally (iv) an additive white Gaussian noise nk 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 B and S are identical for all images in the sequence. This model leads to the following set of equations, where all images are ordered lexicographically:(2.1)yk=SBMkz+nk:=Wkz+nkfork=1,2,…,K,
where Wk=SBMk represents the imaging system.The recovery ofz from {yk}k=1K 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 S,B, and Mk are known, or the accuracy in estimating the degraded operators. Throughout this paper, we assume that S,B, and Mk are known. The decimation S 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 s⩾1 on both axes. In most cases, the blur B 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)E(z):=Φ(z)+μΨ(z),
where Φ(z) is a regularization term which includes prior information about the original image, Ψ(z) denotes the data fitting term depending on the given LR image {yk}k=1K, 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. Ifz is corrupted by additive white Gaussian noise, the MAP estimation will lead to the data fitting term Ψ(z)=(1/2)∑k=1K∥Wkz-yk∥L2(ΩL)2, where ΩL⊂Ω and Ω denote a bounded and open domain of continuous LR image and HR image in R2, respectively. For regularization term, the popular choice is total variation seminorm ∥z∥TV=∫Ω|∇z|dx, 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:(2.3)min⁡z∈BV(Ω){E(z)=‖z‖TV+μ2∑k=1K‖Wkz-yk‖L2(ΩL)2},
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 E(z) 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 B 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

LetH be a real Hilbert space, and let A,B:H→2H be two set-valued operators. We assume that A and B are maximal monotone; that is, their resolvents JA:=(I+A)-1 and JB:=(I+B)-1 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 problem(3.1)0∈A(x̂)+B(x̂).The idea of the forward-backward splitting algorithm is that, for any constantγ>0, we have(3.2)0∈A(x̂)+B(x̂)⟺x̂-γB(x̂)∈x̂+γA(x̂)⟺x̂∈JγA(I-γB)(x̂).
This leads to the following result.Theorem 3.1 (FBS [43, Theorem 2.3.17]).
Suppose thatA:H→2H is maximal monotone and B:H→H is a monotone operator such that ηB is firmly nonexpansive for some η>0. Furthermore, assume that a solution of (3.1) exists. Then, for every start element x0 and step size γ∈(0,2η), the forward-backward splitting algorithm
(3.3)xn+1∈JγA(I-γB)(xn)
converges weakly to an element of the set of solutions (A+B)-1(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:(3.4)0∈A(x̂)+B(x̂)⟺x̂∈JγA(I-γB)(x̂)⟺x̂+γB(x̂)∈JγA(I-γB)(x̂)+γB(x̂)⟺x̂∈JγB(JγA(I-γB)(x̂)+γB(x̂)).
This leads to the following result.Theorem 3.2 (DBS [43, Theorem 2.3.21]).
LetA,B:H→2H be maximal monotone, and assume that a solution of (3.1) exists. Then, for any elements y0 and x0 and any step size γ>0, the following Douglas-Rachford splitting algorithm converges weakly to an element ŷ:
(3.5)yn+1=JγA(2xn-yn)+yn-xn,xn+1=JγB(yn+1).
Furthermore, it holds that x̂=JγB(ŷ) satisfies 0∈A(x̂)+B(x̂). If H is finite-dimensional, then the sequence {xn} converges to x̂.


### 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 form(3.6)min⁡z∈RL1N1×L2N2{E(z)=ϕ(Dz)︸Φ(z)+μΨ(z)}
with Ψ(z)=(1/2)∑k=1K∥Wkz-yk∥2, ϕ(Dz)=∑i,j∥(Dz)i,j∥ is the discrete total variation of z. Here, ∥·∥ denotes Euclidean norm, and D is given by (3.7)(Dz)i,j=((Dx+z)i,j;(Dy+z)i,j)=(zi,j+1-zi,j;zi+1,j-zi,j),for1≤i≤L1N1,1≤j⩽L2N2,
where Dx+ and Dy+ are forward difference operators with periodic boundary condition zi,L2N2+1=zi,1 and zL1N1+1,j=z1,j. Therefore, D is a BCCB matrix.Consider the equivalent problem of (3.6)(3.8)min⁡z∈RL1N1×L2N2w∈R2L1N1×L2N2{E(z)=ϕ(w)+μΨ(z)},s.t.w=Dz.
The Lagrangian for problem (3.8) is(3.9)L(w,z,λ)=ϕ(w)+μΨ(z)-〈λ,w-Dz〉,
where the dual variable λ∈R2L1N1×L2N2 can be thought of as a vector of Lagrange multipliers. Therefore, the dual problem of (3.8) is(3.10)max⁡λinf⁡w,zL(w,z,λ)=-min⁡λ{ϕ*(λ)+μΨ*(-μ-1D*λ)},
where ϕ*(resp., Ψ*) denotes conjugate function of ϕ (resp., Ψ).Define operatorsA and B by(3.11)A(λ)=∂ϕ*(λ),B(λ)=μ∂Ψ*(-μ-1D*λ).
By Fermat's rule, solving the dual problem (3.10) is equivalent to finding λ such that(3.12)0∈A(λ)+B(λ).By formally applying DRS (3.5) to (3.12) with α as the step size, respectively, the ADMM iterations [44–51] are given by  (3.13)wn+1=argmin⁡w{ϕ(w)-〈λn,w〉+α2‖w-Dzn‖2},(3.14)zn+1=argmin⁡z{μΨ(z)+〈λn,Dz〉+α2‖wn+1-Dz‖2},(3.15)λn+1=λn+α(Dzn+1-wn+1).
As pointed out by Setzer in [43, 49], we note that this iteration scheme coincides with DRS algorithm (3.5) with yn=γ(αλn+Dzn), xn=γαλn, and γ=α. Since operators ϕ* and Ψ* are proper, lower semicontinuous, and convex, the operators A and B are maximal monotone see [43]. According to Theorem 3.2, the sequence {zn} 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 ∂ϕ(wn+1)=λn+1 and ∂Ψ(zn+1)=(-1/μ)λn+1.


#### 3.2.1.w-Subproblem

It is not difficult to show that the minimization of (3.13) with respect to w is equivalent to solving L1N1×L2N2 two-dimensional problem of the form(3.16)min⁡wi,j{‖wi,j‖+α2‖wi,j-((Dzn)i,j+λi,jnα)‖2},
for which the unique minimizer is given by the following two-dimensional shrinkage formula:(3.17)wi,jn+1=max⁡{‖(Dzn)i,j+λi,jnα‖-1α,0}(Dzn)i,j+λi,jn/α‖(Dzn)i,j+λi,jn/α‖,
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.z-Subproblem

Subproblem (3.14) is quadratic in z, and the minimizer zn+1 is given by the normal equations(3.18)(DTD+μα∑kWkTWk)zn+1=DT(wn+1-λnα)+μα∑kWkTyk,
where WkT is the transpose operator of Wk.Note thatDTD in (3.18) is BCCB matrix and can be diagonalized by FFT. Moreover, with periodic boundary condition, both B and Mk are BCCB matrices. Exploiting the fact that the product order of two BCCB matrices can commute, we get that MkTBT=BTMkT, BMk=MkB. Then, ∑kWkTWk=BTRB (R:=∑kMkTSTSMk is a diagonal matrix [31]). However, ∑kWkTWk=BTRB 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 Ψ(z)=(1/2)∑k=1K∥Wkz-yk∥2 in (3.14) that couples the variable z by the matrix Wk makes the algorithm computationally expensive. There are some techniques to overcome these problems. In [50], the authors introduced the three constrains: w1=Bz, w2=Dz, w3=z 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 z-subproblem (3.14), which can decouple the variable z and the constraint matrix Wk and make the Hessian matrix of the energy functional have BCCB structure.Letf1(z)=(α/2μ)∥wn+1-Dz∥2+〈λn/μ,Dz〉 and f2(z)=(1/2)∑k=1K∥Wkz-yk∥2. Then, ∇f2(zn)=BT(RBzn-s) with s:=∑kMkTSTyk. Define operators A=∂f1 and B=∇f2; the FBS method (3.3) with γ as the step size applied to (3.14) leads to the iterative scheme(3.19)zn+1=argmin⁡z{α2‖wn+1-Dz‖2+〈λn,Dz〉+μ2γ‖z-zn+γ∇f2(zn)‖2}.
The minimizer zn+1 is given by the normal equations(3.20)(DTD+μγαI)zn+1=DT(wn+1-λnα)+μγα(zn-γ∇f2(zn)).
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 thatB=∇f2 is Lipschitz continuous with Lipschitz constant β=∥R∥∥B∥2, where ∥B∥ (resp., ∥R∥) is the matrix norm of B (resp., R). From [41,Theorem 2.3.19], (1/β)B is firmly nonexpansive. Then, the convergence of zn 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 z0 from maximum likelihood (ML) estimation; that is, z0 is the minimizer of the optimization problem as follows:(3.21)ẑ=argmin⁡z{∑k‖SBMkz-yk‖2}.
The gradient descent algorithm suggests the following iterative equation for the solution of (3.21):(3.22)ẑn+1=ẑn-Δt∑kMkTBTST(SBMkẑn-yk)=ẑn-ΔtBT∑kMkTST(SMkBẑn-yk),
where Δt denotes step size. Let us define the blurred super-resolution image by ûn=Bẑn. Multiplying both sides of (3.22) with B, we get(3.23)ûn+1=ûn+ΔtBBT[s-Rûn],
where s=∑kMkTSTyk and R=∑kMkTSTSMk. According to [52], the steady state solution of (3.23) is given by(3.24)û∞=R-1s.
We note that MkTST copies the values from the LR grid to the HR grid after proper shifting and zero filling and SMk copies a selected set of pixels in HR grid back on the LR grid (Figure 1 illustrates the effect of upsampling and downsampling matrices ST and S). Neither of these two operations changes the pixel value. It is easy to show that R is a diagonal matrix. Each diagonal entry in R 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 s is simply the addition of the measurements after proper zero-filling interpolation and motion compensation. Thus, û∞=R-1s is none other than the pixel-wise average of the measurement. Therefore, the noise of û∞ is reduced due to the averaging.Figure 1
Effect ofST and S on image (resolution enhancement factor is three).Becauseû∞=Bẑ∞, Wiener filter was applied to û∞. Then, the restoration image ẑ∞ is taken as the initial data z0.As stated previously, precondition design procedure of the initial dataz0 is summed up in Algorithm 1 as follows:Algorithm 1
Input: s=∑kMkTSTyk,R=∑kMkTSTSMk,B.Output: z0=B-1R-1s.


#### 3.2.4. Algorithm Description

To sum up the above arguments, the complete resulting algorithm is summarized in Algorithm2 as follows.Algorithm 2
Input: s=∑kMkTSTyk,R=∑kMkTSTSMk,B,toleranceζstop,γ∈]0,2/β[, and α,μ>0.Initialize: z0=B-1R-1s,λ0=0, and n=0.Do(1) Compute wn+1 according to (3.17) for fixed (zn,λn).(2) Compute zn+1 according to (3.20) for fixed (wn+1,λn).(3) Update λn+1 according to (3.15).(4)n=n+1.Until∥zn+1-zn∥/∥zn+1∥≤ζstop.


#### 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 O(N2) for an N×N image. Hence, the w-subproblem (3.16) and λ can be solved quickly. The solution of the z-subproblem (3.20) requires two FFTs (including one inverse FFT), which has a total complexity in the order of O(N2log(N2))=O(N2log(N)).


## 3.1. FBS and DRS

LetH be a real Hilbert space, and let A,B:H→2H be two set-valued operators. We assume that A and B are maximal monotone; that is, their resolvents JA:=(I+A)-1 and JB:=(I+B)-1 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 problem(3.1)0∈A(x̂)+B(x̂).The idea of the forward-backward splitting algorithm is that, for any constantγ>0, we have(3.2)0∈A(x̂)+B(x̂)⟺x̂-γB(x̂)∈x̂+γA(x̂)⟺x̂∈JγA(I-γB)(x̂).
This leads to the following result.Theorem 3.1 (FBS [43, Theorem 2.3.17]).
Suppose thatA:H→2H is maximal monotone and B:H→H is a monotone operator such that ηB is firmly nonexpansive for some η>0. Furthermore, assume that a solution of (3.1) exists. Then, for every start element x0 and step size γ∈(0,2η), the forward-backward splitting algorithm
(3.3)xn+1∈JγA(I-γB)(xn)
converges weakly to an element of the set of solutions (A+B)-1(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:(3.4)0∈A(x̂)+B(x̂)⟺x̂∈JγA(I-γB)(x̂)⟺x̂+γB(x̂)∈JγA(I-γB)(x̂)+γB(x̂)⟺x̂∈JγB(JγA(I-γB)(x̂)+γB(x̂)).
This leads to the following result.Theorem 3.2 (DBS [43, Theorem 2.3.21]).
LetA,B:H→2H be maximal monotone, and assume that a solution of (3.1) exists. Then, for any elements y0 and x0 and any step size γ>0, the following Douglas-Rachford splitting algorithm converges weakly to an element ŷ:
(3.5)yn+1=JγA(2xn-yn)+yn-xn,xn+1=JγB(yn+1).
Furthermore, it holds that x̂=JγB(ŷ) satisfies 0∈A(x̂)+B(x̂). If H is finite-dimensional, then the sequence {xn} converges to x̂.


## 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 form(3.6)min⁡z∈RL1N1×L2N2{E(z)=ϕ(Dz)︸Φ(z)+μΨ(z)}
with Ψ(z)=(1/2)∑k=1K∥Wkz-yk∥2, ϕ(Dz)=∑i,j∥(Dz)i,j∥ is the discrete total variation of z. Here, ∥·∥ denotes Euclidean norm, and D is given by (3.7)(Dz)i,j=((Dx+z)i,j;(Dy+z)i,j)=(zi,j+1-zi,j;zi+1,j-zi,j),for1≤i≤L1N1,1≤j⩽L2N2,
where Dx+ and Dy+ are forward difference operators with periodic boundary condition zi,L2N2+1=zi,1 and zL1N1+1,j=z1,j. Therefore, D is a BCCB matrix.Consider the equivalent problem of (3.6)(3.8)min⁡z∈RL1N1×L2N2w∈R2L1N1×L2N2{E(z)=ϕ(w)+μΨ(z)},s.t.w=Dz.
The Lagrangian for problem (3.8) is(3.9)L(w,z,λ)=ϕ(w)+μΨ(z)-〈λ,w-Dz〉,
where the dual variable λ∈R2L1N1×L2N2 can be thought of as a vector of Lagrange multipliers. Therefore, the dual problem of (3.8) is(3.10)max⁡λinf⁡w,zL(w,z,λ)=-min⁡λ{ϕ*(λ)+μΨ*(-μ-1D*λ)},
where ϕ*(resp., Ψ*) denotes conjugate function of ϕ (resp., Ψ).Define operatorsA and B by(3.11)A(λ)=∂ϕ*(λ),B(λ)=μ∂Ψ*(-μ-1D*λ).
By Fermat's rule, solving the dual problem (3.10) is equivalent to finding λ such that(3.12)0∈A(λ)+B(λ).By formally applying DRS (3.5) to (3.12) with α as the step size, respectively, the ADMM iterations [44–51] are given by  (3.13)wn+1=argmin⁡w{ϕ(w)-〈λn,w〉+α2‖w-Dzn‖2},(3.14)zn+1=argmin⁡z{μΨ(z)+〈λn,Dz〉+α2‖wn+1-Dz‖2},(3.15)λn+1=λn+α(Dzn+1-wn+1).
As pointed out by Setzer in [43, 49], we note that this iteration scheme coincides with DRS algorithm (3.5) with yn=γ(αλn+Dzn), xn=γαλn, and γ=α. Since operators ϕ* and Ψ* are proper, lower semicontinuous, and convex, the operators A and B are maximal monotone see [43]. According to Theorem 3.2, the sequence {zn} 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 ∂ϕ(wn+1)=λn+1 and ∂Ψ(zn+1)=(-1/μ)λn+1.


### 3.2.1.w-Subproblem

It is not difficult to show that the minimization of (3.13) with respect to w is equivalent to solving L1N1×L2N2 two-dimensional problem of the form(3.16)min⁡wi,j{‖wi,j‖+α2‖wi,j-((Dzn)i,j+λi,jnα)‖2},
for which the unique minimizer is given by the following two-dimensional shrinkage formula:(3.17)wi,jn+1=max⁡{‖(Dzn)i,j+λi,jnα‖-1α,0}(Dzn)i,j+λi,jn/α‖(Dzn)i,j+λi,jn/α‖,
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.z-Subproblem

Subproblem (3.14) is quadratic in z, and the minimizer zn+1 is given by the normal equations(3.18)(DTD+μα∑kWkTWk)zn+1=DT(wn+1-λnα)+μα∑kWkTyk,
where WkT is the transpose operator of Wk.Note thatDTD in (3.18) is BCCB matrix and can be diagonalized by FFT. Moreover, with periodic boundary condition, both B and Mk are BCCB matrices. Exploiting the fact that the product order of two BCCB matrices can commute, we get that MkTBT=BTMkT, BMk=MkB. Then, ∑kWkTWk=BTRB (R:=∑kMkTSTSMk is a diagonal matrix [31]). However, ∑kWkTWk=BTRB 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 Ψ(z)=(1/2)∑k=1K∥Wkz-yk∥2 in (3.14) that couples the variable z by the matrix Wk makes the algorithm computationally expensive. There are some techniques to overcome these problems. In [50], the authors introduced the three constrains: w1=Bz, w2=Dz, w3=z 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 z-subproblem (3.14), which can decouple the variable z and the constraint matrix Wk and make the Hessian matrix of the energy functional have BCCB structure.Letf1(z)=(α/2μ)∥wn+1-Dz∥2+〈λn/μ,Dz〉 and f2(z)=(1/2)∑k=1K∥Wkz-yk∥2. Then, ∇f2(zn)=BT(RBzn-s) with s:=∑kMkTSTyk. Define operators A=∂f1 and B=∇f2; the FBS method (3.3) with γ as the step size applied to (3.14) leads to the iterative scheme(3.19)zn+1=argmin⁡z{α2‖wn+1-Dz‖2+〈λn,Dz〉+μ2γ‖z-zn+γ∇f2(zn)‖2}.
The minimizer zn+1 is given by the normal equations(3.20)(DTD+μγαI)zn+1=DT(wn+1-λnα)+μγα(zn-γ∇f2(zn)).
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 thatB=∇f2 is Lipschitz continuous with Lipschitz constant β=∥R∥∥B∥2, where ∥B∥ (resp., ∥R∥) is the matrix norm of B (resp., R). From [41,Theorem 2.3.19], (1/β)B is firmly nonexpansive. Then, the convergence of zn 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 z0 from maximum likelihood (ML) estimation; that is, z0 is the minimizer of the optimization problem as follows:(3.21)ẑ=argmin⁡z{∑k‖SBMkz-yk‖2}.
The gradient descent algorithm suggests the following iterative equation for the solution of (3.21):(3.22)ẑn+1=ẑn-Δt∑kMkTBTST(SBMkẑn-yk)=ẑn-ΔtBT∑kMkTST(SMkBẑn-yk),
where Δt denotes step size. Let us define the blurred super-resolution image by ûn=Bẑn. Multiplying both sides of (3.22) with B, we get(3.23)ûn+1=ûn+ΔtBBT[s-Rûn],
where s=∑kMkTSTyk and R=∑kMkTSTSMk. According to [52], the steady state solution of (3.23) is given by(3.24)û∞=R-1s.
We note that MkTST copies the values from the LR grid to the HR grid after proper shifting and zero filling and SMk copies a selected set of pixels in HR grid back on the LR grid (Figure 1 illustrates the effect of upsampling and downsampling matrices ST and S). Neither of these two operations changes the pixel value. It is easy to show that R is a diagonal matrix. Each diagonal entry in R 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 s is simply the addition of the measurements after proper zero-filling interpolation and motion compensation. Thus, û∞=R-1s is none other than the pixel-wise average of the measurement. Therefore, the noise of û∞ is reduced due to the averaging.Figure 1
Effect ofST and S on image (resolution enhancement factor is three).Becauseû∞=Bẑ∞, Wiener filter was applied to û∞. Then, the restoration image ẑ∞ is taken as the initial data z0.As stated previously, precondition design procedure of the initial dataz0 is summed up in Algorithm 1 as follows:Algorithm 1
Input: s=∑kMkTSTyk,R=∑kMkTSTSMk,B.Output: z0=B-1R-1s.


### 3.2.4. Algorithm Description

To sum up the above arguments, the complete resulting algorithm is summarized in Algorithm2 as follows.Algorithm 2
Input: s=∑kMkTSTyk,R=∑kMkTSTSMk,B,toleranceζstop,γ∈]0,2/β[, and α,μ>0.Initialize: z0=B-1R-1s,λ0=0, and n=0.Do(1) Compute wn+1 according to (3.17) for fixed (zn,λn).(2) Compute zn+1 according to (3.20) for fixed (wn+1,λn).(3) Update λn+1 according to (3.15).(4)n=n+1.Until∥zn+1-zn∥/∥zn+1∥≤ζstop.


### 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 O(N2) for an N×N image. Hence, the w-subproblem (3.16) and λ can be solved quickly. The solution of the z-subproblem (3.20) requires two FFTs (including one inverse FFT), which has a total complexity in the order of O(N2log(N2))=O(N2log(N)).


## 3.2.1.w-Subproblem

It is not difficult to show that the minimization of (3.13) with respect to w is equivalent to solving L1N1×L2N2 two-dimensional problem of the form(3.16)min⁡wi,j{‖wi,j‖+α2‖wi,j-((Dzn)i,j+λi,jnα)‖2},
for which the unique minimizer is given by the following two-dimensional shrinkage formula:(3.17)wi,jn+1=max⁡{‖(Dzn)i,j+λi,jnα‖-1α,0}(Dzn)i,j+λi,jn/α‖(Dzn)i,j+λi,jn/α‖,
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.z-Subproblem

Subproblem (3.14) is quadratic in z, and the minimizer zn+1 is given by the normal equations(3.18)(DTD+μα∑kWkTWk)zn+1=DT(wn+1-λnα)+μα∑kWkTyk,
where WkT is the transpose operator of Wk.Note thatDTD in (3.18) is BCCB matrix and can be diagonalized by FFT. Moreover, with periodic boundary condition, both B and Mk are BCCB matrices. Exploiting the fact that the product order of two BCCB matrices can commute, we get that MkTBT=BTMkT, BMk=MkB. Then, ∑kWkTWk=BTRB (R:=∑kMkTSTSMk is a diagonal matrix [31]). However, ∑kWkTWk=BTRB 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 Ψ(z)=(1/2)∑k=1K∥Wkz-yk∥2 in (3.14) that couples the variable z by the matrix Wk makes the algorithm computationally expensive. There are some techniques to overcome these problems. In [50], the authors introduced the three constrains: w1=Bz, w2=Dz, w3=z 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 z-subproblem (3.14), which can decouple the variable z and the constraint matrix Wk and make the Hessian matrix of the energy functional have BCCB structure.Letf1(z)=(α/2μ)∥wn+1-Dz∥2+〈λn/μ,Dz〉 and f2(z)=(1/2)∑k=1K∥Wkz-yk∥2. Then, ∇f2(zn)=BT(RBzn-s) with s:=∑kMkTSTyk. Define operators A=∂f1 and B=∇f2; the FBS method (3.3) with γ as the step size applied to (3.14) leads to the iterative scheme(3.19)zn+1=argmin⁡z{α2‖wn+1-Dz‖2+〈λn,Dz〉+μ2γ‖z-zn+γ∇f2(zn)‖2}.
The minimizer zn+1 is given by the normal equations(3.20)(DTD+μγαI)zn+1=DT(wn+1-λnα)+μγα(zn-γ∇f2(zn)).
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 thatB=∇f2 is Lipschitz continuous with Lipschitz constant β=∥R∥∥B∥2, where ∥B∥ (resp., ∥R∥) is the matrix norm of B (resp., R). From [41,Theorem 2.3.19], (1/β)B is firmly nonexpansive. Then, the convergence of zn 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 z0 from maximum likelihood (ML) estimation; that is, z0 is the minimizer of the optimization problem as follows:(3.21)ẑ=argmin⁡z{∑k‖SBMkz-yk‖2}.
The gradient descent algorithm suggests the following iterative equation for the solution of (3.21):(3.22)ẑn+1=ẑn-Δt∑kMkTBTST(SBMkẑn-yk)=ẑn-ΔtBT∑kMkTST(SMkBẑn-yk),
where Δt denotes step size. Let us define the blurred super-resolution image by ûn=Bẑn. Multiplying both sides of (3.22) with B, we get(3.23)ûn+1=ûn+ΔtBBT[s-Rûn],
where s=∑kMkTSTyk and R=∑kMkTSTSMk. According to [52], the steady state solution of (3.23) is given by(3.24)û∞=R-1s.
We note that MkTST copies the values from the LR grid to the HR grid after proper shifting and zero filling and SMk copies a selected set of pixels in HR grid back on the LR grid (Figure 1 illustrates the effect of upsampling and downsampling matrices ST and S). Neither of these two operations changes the pixel value. It is easy to show that R is a diagonal matrix. Each diagonal entry in R 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 s is simply the addition of the measurements after proper zero-filling interpolation and motion compensation. Thus, û∞=R-1s is none other than the pixel-wise average of the measurement. Therefore, the noise of û∞ is reduced due to the averaging.Figure 1
Effect ofST and S on image (resolution enhancement factor is three).Becauseû∞=Bẑ∞, Wiener filter was applied to û∞. Then, the restoration image ẑ∞ is taken as the initial data z0.As stated previously, precondition design procedure of the initial dataz0 is summed up in Algorithm 1 as follows:Algorithm 1
Input: s=∑kMkTSTyk,R=∑kMkTSTSMk,B.Output: z0=B-1R-1s.


## 3.2.4. Algorithm Description

To sum up the above arguments, the complete resulting algorithm is summarized in Algorithm2 as follows.Algorithm 2
Input: s=∑kMkTSTyk,R=∑kMkTSTSMk,B,toleranceζstop,γ∈]0,2/β[, and α,μ>0.Initialize: z0=B-1R-1s,λ0=0, and n=0.Do(1) Compute wn+1 according to (3.17) for fixed (zn,λn).(2) Compute zn+1 according to (3.20) for fixed (wn+1,λn).(3) Update λn+1 according to (3.15).(4)n=n+1.Until∥zn+1-zn∥/∥zn+1∥≤ζstop.


## 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 O(N2) for an N×N image. Hence, the w-subproblem (3.16) and λ can be solved quickly. The solution of the z-subproblem (3.20) requires two FFTs (including one inverse FFT), which has a total complexity in the order of O(N2log(N2))=O(N2log(N)).


## 4. Experimental Results

In this section, we present some experimental examples to demonstrate the performance of our method. We use the three256×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.Results of different resolution enhancement methods applied to “Lena” image degraded by Gaussian blur with support size3×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)(e)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(4.1)PSNR=10log⁡10{L1N1×L2N2max⁡{z}2‖z*-z‖2},ReErr=‖z*-z‖2‖z‖2,
where z and z* are the original and the SR reconstructed images, respectively, and L1N1×L2N2 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 z0=B-1R-1s. 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: (4.2)ReDiff=‖zn+1-zn‖‖zn+1‖≤10-4.In the first test, we apply Gaussian kernel with window size3×3, standard deviation 0.5, and different noise level (noise variance σ2=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.
ExperimentsFigure2(a)Figure3(a)Figure4(a)Figure5(a)Figure6(a)Figure7(a)Figure8(a)LDFP methodPSNR (dB)27.413727.196825.771225.688124.898526.361821.9450ReErr0.00560.00590.00940.00960.01150.00710.0164Iterations46 (20)44 (20)49 (20)49 (20)43 (20)49 (20)44 (20)Time(s)25.515624.531327.312527.296930.125034.203125.1406AMM methodPSNR (dB)28.503728.003526.350826.124724.368026.394722.0199ReErr0.00440.00490.00820.00870.01300.00710.0161Iterations189 (7)176 (7)239 (7)228 (7)209 (7)184 (7)131 (7)Time(s)20.390618.953125.875024.625079.062569.593814.4531OSM methodPSNR28.882028.296826.564026.308624.878226.532222.0065ReErr0.00400.00460.00790.00830.01160.00690.0162Iterations909310010119418358Time(s)5.20315.34385.76565.843818.953117.82813.3750Results of different resolution enhancement methods applied to “Lena” image degraded by Gaussian blur with support size3×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)Results of different resolution enhancement methods applied to “Cameraman” image degraded by Gaussian blur with support size3×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)(e)Results of different resolution enhancement methods applied to “Cameraman” image degraded by Gaussian blur with support size3×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)Results of different resolution enhancement methods applied to “Lena” image degraded by Gaussian blur with support size15×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)Results of different resolution enhancement methods applied to “Cameraman” image degraded by Gaussian blur with support size15×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)Results of different resolution enhancement methods applied to “Fingerprint” image degraded by Gaussian blur with support size3×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)(e)In Figure9, we show that the convergence of the proposed method is faster than LDFP and AMM methods. The x-axis refers to the number of iterations. The y-axis in Figures 9(a), 9(b) refers to the PSNR and the y-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.Convergence performance of three methods in the Gaussian blur with support size3×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.
(a)(b)(c)(d)


## 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.



---

*Source: 380807-2011-02-21.xml*