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

Super-resolution is a fusion process for reconstructing a high-resolution image from a set of lowresolution 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.


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

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 L 1 N 1 × L 2 N 2 written in lexicographical notation as the vector z z 1 , z 2 , . . ., z N T , where N L 1 N 1 × L 2 N 2 .Let the parameters L 1 and L 2 be the subsampling factors in the horizontal and vertical directions, respectively; each observed LR image has the size N 1 × N 2 .Thus, the LR image can be represented as y k y k 1 , y k 2 , . . ., y k M T , for k 1, 2, . . ., K and M N 1 × N 2 .A popular model assumes that LR images {y k } K k 1 are generated from HR image z through a sequence of operations that includes i geometrical motions M k , ii a linear spaceinvariant blur B, iii a subsampling step represented by S, and finally iv an additive white Gaussian noise n k 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: where W k SBM k represents the imaging system.The recovery of z from {y k } K k 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 S, B, and M k are known, or the accuracy in estimating the degraded operators.Throughout this paper, we assume that S, B, and M k 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 indepth study of this simple case allows much insight to be gained about the problems inherent to SR image reconstruction.

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 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 {y k } K k 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 z is corrupted by additive white Gaussian noise, the MAP estimation will lead to the data fitting term where Ω L ⊂ Ω and Ω denote a bounded and open domain of continuous LR image and HR image in R 2 , 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: 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.

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][20][21][22][23][24][25][26][27][28][29][30][31][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 wellknown 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, doublebackward 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.

FBS and DRS
Let H be a real Hilbert space, and let A, B : H → 2 H be two set-valued operators.We assume that A and B are maximal monotone; that is, their resolvents J A : I A −1 and J B : 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 The idea of the forward-backward splitting algorithm is that, for any constant γ > 0, we have 3.2 This leads to the following result.
Theorem 3.1 FBS 43, Theorem 2.3.17 .Suppose that A : H → 2 H 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 x 0 and step size γ ∈ 0, 2η , the forward-backward splitting algorithm 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
This leads to the following result.Theorem 3.2 DBS 43, Theorem 2.3.21 .Let A, B : H → 2 H be maximal monotone, and assume that a solution of 3.1 exists.Then, for any elements y 0 and x 0 and any step size γ > 0, the following Douglas-Rachford splitting algorithm converges weakly to an element y: x n 1 J γB y n 1 .

3.5
Furthermore, it holds that x J γB y satisfies 0 ∈ A x B x .If H is finite-dimensional, then the sequence {x n } converges to x.

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 min Dz i,j is the discrete total variation of z.Here, • denotes Euclidean norm, and D is given by where D x and D y are forward difference operators with periodic boundary condition Consider the equivalent problem of 3.6 min E z φ w μΨ z , s.t.w Dz.

3.8
The Lagrangian for problem 3.8 is where the dual variable λ ∈ R 2L 1 N 1 ×L 2 N 2 can be thought of as a vector of Lagrange multipliers.Therefore, the dual problem of 3.8 is where φ * resp., Ψ * denotes conjugate function of φ resp., Ψ .

Mathematical Problems in Engineering 7
Define operators A and B by By Fermat's rule, solving the dual problem 3.10 is equivalent to finding λ such that 0 ∈ A λ B λ .

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

3.15
As pointed out by Setzer in 43, 49 , we note that this iteration scheme coincides with DRS algorithm 3.5 with y n γ αλ n Dz n , x n γαλ 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 {z n } 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 ∂φ w n 1 λ n 1 and ∂Ψ z n 1 −1/μ λ n 1 .

w-Subproblem
It is not difficult to show that the minimization of 3.13 with respect to w is equivalent to solving for which the unique minimizer is given by the following two-dimensional shrinkage formula: 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.

z-Subproblem
Subproblem 3.14 is quadratic in z, and the minimizer z n 1 is given by the normal equations where W T k is the transpose operator of W k .Note that D T D in 3.18 is BCCB matrix and can be diagonalized by FFT.Moreover, with periodic boundary condition, both B and M k are BCCB matrices.Exploiting the fact that the product order of two BCCB matrices can commute, we get that 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 14 that couples the variable z by the matrix W k makes the algorithm computationally expensive.There are some techniques to overcome these problems.In 50 , the authors introduced the three constrains: w 1 Bz, w 2 Dz, w 3 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 zsubproblem 3.14 , which can decouple the variable z and the constraint matrix W k and make the Hessian matrix of the energy functional have BCCB structure. Let Then, ∇f 2 z n B T RBz n − s with s : k M k T S T y k .Define operators A ∂f 1 and B ∇f 2 ; the FBS method 3.3 with γ as the step size applied to 3.14 leads to the iterative scheme

3.19
The minimizer z n 1 is given by the normal equations

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 B ∇f 2 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 z n is ensured by Theorem 3.1, while γ ∈ 0, 2/β .

Optimal Initial Guess
Firstly, borrowed the idea of 31 , we take initial data z 0 from maximum likelihood ML estimation; that is, z 0 is the minimizer of the optimization problem as follows:

3.21
The gradient descent algorithm suggests the following iterative equation for the solution of 3.21 :

3.22
where Δt denotes step size.Let us define the blurred super-resolution image by u n B z n .Multiplying both sides of 3.22 with B, we get where s k M T k S T y k and R k M T k S T SM k .According to 52 , the steady state solution of 3.23 is given by u ∞ R −1 s.

3.24
We note that M T k S T copies the values from the LR grid to the HR grid after proper shifting and zero filling and SM k copies a selected set of pixels in HR grid back on the LR grid Figure 1 illustrates the effect of upsampling and downsampling matrices S T 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 −1 s is none other than the pixel-wise average of the measurement.Therefore, the noise of u ∞ is reduced due to the averaging.Because u ∞ B z ∞ , Wiener filter was applied to u ∞ .Then, the restoration image z ∞ is taken as the initial data z 0 .
As stated previously, precondition design procedure of the initial data z 0 is summed up in Algorithm 1 as follows.

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

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 N 2 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 N 2 log N 2 O N 2 log N .

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  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.
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 10 log 10 where z and z * are the original and the SR reconstructed images, respectively, and L 1 N 1 × L 2 N 2 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 z 0 B −1 R −1 s.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: Mathematical Problems in Engineering In the first test, we apply Gaussian kernel with window size 3 × 3, standard deviation 0.5, and different noise level noise variance σ 2 5, 20 .One of LR frames is presented in Figures 2-5 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.In Figure 9, 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.

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 Mathematical Problems in Engineering 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.

Figure 1 :
Figure 1: Effect of S T and S on image resolution enhancement factor is three .

Figure 2 :Figure 3 :
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 .

Figure 4 :
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 .

Figure 5 :
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 .

Figure 6 :
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 .

Figure 7 :
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 .
, 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.

Figure 8 :
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 .

Figure 9 :
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.

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.