Efficient Algorithm for Isotropic and Anisotropic Total Variation Deblurring and Denoising

A new deblurring and denoising algorithm is proposed, for isotropic total variation-based image restoration.The algorithm consists of an efficient solver for the nonlinear system and an acceleration strategy for the outer iteration. For the nonlinear system, the split Bregmanmethod is used to convert it into linear system, and an algebraicmultigridmethod is applied to solve the linearized system. For the outer iteration, we have conducted formal convergence analysis to determine an auxiliary linear term that significantly stabilizes and accelerates the outer iteration. Numerical experiments demonstrate that our algorithm for deblurring and denoising problems is efficient.


Introduction
The purpose of image restoration is to recover original image  from an observed data  (noisy and blurred image) from the relation where  is a known linear blurring operator and  is a Gaussian white noise.The well-known minimization problem for image denoising and deblurring based on isotropic total variation (TV) is proposed by Rudin et al. [1]: Here  > 0 is the penalty parameter, and  > 0 is the diffusion regularization parameter and is typically small.The functional in (2) is strictly convex with a unique global minimizer.The paper [2] showed the well-posedness of problem (2) as  → 0 + and the existence of a unique global minimizer.The corresponding Euler-Lagrange equation for (2) is given by where  * is the adjoint operator of  with respect to standard  2 inner product.
In addition to the isotropic models, the anisotropic models for a qualitative improvement at corners are shown in [3] with that a transfer of the discretization from the anisotropic model to the isotropic setting results in an improvement of rotational invariance.The anisotropic TV defined in [3] is The existence and uniqueness of solutions of the anisotropic total variation flow are shown in [4].The explicit timemarching schemes are applied in [1,5].Zuo and Lin [6] proposed a generalized accelerated proximal gradient (GAPG) approach for solving TV-based image restoration problems by replacing the Lipschitz constant with an appropriate positive-definite matrix, resulting in faster convergence, and the TV regularization can be either isotropic or anisotropic.The convergence rate of ( 2 ) is maintained by GAPG, where  is the number of iterations.The introduction of an anisotropy to TV regularization [7] indeed leads to improved denoising: the staircasing effect is reduced while at the same time the creation of artifacts is suppressed.
There are different algorithms to solve (3).Earlier works include the time-marching scheme [1], the affine scaling algorithm [8], the Newton's method with a continuation procedure on  [9], a fixed point iteration with multigrid method [10][11][12], the combination of fixed point iteration and preconditioned conjugate gradient which method is proposed in [13], a fast algorithm based on the discrete cosine transform [14,15], a new modular solver [16], a proximity algorithm [17], and a nonlinear primal-dual method [18].Though the corresponding linear system of the nonlinear primal-dual method is twice as large as the primal system, it is shown to be better conditioned.In [19], the combination of the algebraic multigrid (AMG) methods [20,21], Krylov subspace acceleration, and extrapolation of initial data are successfully combined to give an efficient algorithm for the purely denoising problem.In [22], the authors proposed the split Bregman method to solve image denoising and compressed sensing problem rapidly.Chang et al. [23] added a linear matrix to speed up the computational process.
In this paper, we propose an efficient and rapid algorithm for minimization problem (2) rather than solving the Euler-Lagrange equation (3) directly which combines the split Bregman method, the algebraic multigrid method, and Krylov subspace acceleration.Our algorithm consists of two steps.One uses the split Bregman method to convert (2) into linear system in the outer iteration, and the other adopts an AMG method previously developed by one of the authors for other applications [23,24] to solve a linear system in the inner iteration.Moreover, the outer iteration is accelerated by Krylov subspace extrapolation.One -cycle per inner iteration is sufficient in all our simulations.For the outer iteration, we have developed a stabilizing technique adapted to the blur operator .This is done by adding a linear term on both sides of the equation [23].Motivated by that, a wider linear term different from that in [23] is considered here.This linear stabilizing term is obtained via formal convergence analysis and is expressed explicitly in terms of the mask of .The inclusion of the linear stabilizing term plays a crucial role in our scheme.The outer iteration indeed may diverge without a proper linear stabilizing term (Table 6, Section 5).
The rest of the paper is organized as follows.In Section 2, we discuss the purely image denoising problem by describing briefly the AMG method and the split Bregman method for the isotropic TV model and the anisotropic TV model.In Section 3, we discuss the image denoising and deblurring problem and present the formal convergence analysis and the framework of our scheme.In Section 4, we give out explicit formulae of the linear stabilizing term, together with other implementation details.Numerical results are also compared among several algorithms in Section 5. We end this paper with a brief summary in Section 6.

The Algebraic Multigrid
In this section, we will show how to directly use AMG in solving (3).We first show some denotations to discretize (3) (cf.[19]).We partition the domain Ω = (0, 1) × (0, 1) into  ×  uniform cells with mesh size ℎ = 1/.The cell centers are Following the ideas of [12,19], we discretize (3) by standard five-point finite difference scheme to get with homogeneous Neumann boundary condition: where with and so forth.
To simplify the notation, we abbreviate (6) as where which is fully nonlinear with wildly varying coefficient.Consider the general system (10) in an  2 ×  2 : with the matrix  = () in (10).In general,  varies wildly near areas of high contrast of the image and need not be diagonally dominant.Nevertheless,  is symmetric and positive definite.Moreover, the matrix  *  is widebanded, and the spectra of the matrices () and  *  are quite differently distributed.It is difficult to compute (10).In [13], Vogel and Oman adopted the combination of a fixed point iteration and a product PCG to handle the nonlinear term and the linear system, respectively.In [15], Chan et al.
proposed another preconditioner based on the fast cosine transform.Algebraic multigrid method is considered in [19,24,25] as a linear solver.In contrast to geometric multigrid method, the construction of the coarse grids Ω  is solely based on the algebraic information provided by the matrix  in an AMG setting.The details and effectiveness can be found, for example, in [23,24,26].For the restriction and coarse grid operators, a simple and popular approach is the Galerkin type construction.Namely,  +1  = (  +1 ) * and  +1 =  +1      +1 .How to construct the interpolation operators   +1 is the most important part in the method.The convergence proof for this improved AMG method was given in [25] when   is symmetric positive definite.The robustness of these interpolation formulae is supported by plenty of numerical evidence [24,25].The detailed algorithm using the AMG method to solve (10) is in Algorithm 1.
2.1.The Split Bregman Method for Isotropic TV Model.In [22], the split Bregman method is applied to the ROF model for image denoising problem.For completeness, we show the main progress about the split Bregman method applied to the ROF model.Firstly, consider discretizing ROF model into Setting   ≈ ∇   and   ≈ ∇  , the split Bregman formulation of the isotropic problem then becomes min where Consider it into the following three subproblems: The formula ( 16) is  2 -norm; it can be solved using the variation approach.The corresponding Euler-Lagrange equation for ( 16) is The difference scheme is When  = , the difference scheme of the previous formula is simplified as Using Gauss-Seidel method for (21), we have Problem ( 17) can be solved using the fast shrinkage operators: Initialize: The optimal value of  can be explicitly gotten using the above shrinkage operators.The formula (18) can be discretized directly.Thus the split Bregman algorithm of the isotropic TV denoising is as in Algorithm 2.

The Split Bregman Method Based on AMG for Isotropic TV
Model.It is easy to find that the accuracy of the split Bregman method is not high because the Gauss-Seidel iteration of  is executed only once per loop.At the same time, the AMG algorithm has fast convergence and high precision.It can be applied to large signal-to-noise ratio of image denoising, which has strong stability.So, we consider the combination of the two methods, and get a new denoising algorithm which is the split Bregman method based on AMG (here we call it as Algorithm 3).This algorithm contains four steps.
Moreover, we use Krylov technique to accelerate subspace extrapolation.

The Split Bregman Method for Anisotropic TV Model.
The anisotropic total variation of  [27] is represented by The anisotropic TV model for denoising can be formulated as the following minimization problem: where  is an appropriately chosen positive parameter.Define an operator [27] cut With the previous operator, we have Algorithm 4 to solve the anisotropic TV denoising problem.
It should be noticed that there is no need to solve any equations like the isotropic model.This algorithm need not solve any elliptic equation of  in [27], so the procedure is rapid, simple, and not costly.

Convergence Analysis
We consider the discretized system with general noise and blur: To simplify the notation, we abbreviate (27) as In general, the matrix  *  would be dense.There is one way to completely avoid the matrix inversion.A natural approach would then be It turns out that the iteration ( 29) is not robust and may diverge even for weak blur operator corresponding to the mask (1/64)(1, 1, 4, 1, 1)  (1, 1, 4, 1, 1) with  = 10.To overcome this instability, we propose to add a linear term  +  on both sides of (29) and get where  is arbitrary small positive number and  is a symmetric and positive definite matrix to be determined through the following formal analysis.Now we will analyze the convergence property of (29).
Step 2: To get  by solving (22) using the AMG method.
It is easy to show that the matrices  and  are symmetric.
Taking the inner product on both sides of (34) with  (+1) −  () gives and we will now rearrange the formula (36) as Rewrite the previous formula as Since −Δ and  *  are symmetric and positive definite, a sufficient condition for (42) is given by (31).
Remark.We can see that (31) indeed is a good guideline for devising the iteration of .We will elaborate on the choice of  better suited for numerical purposes in Section 4. Indeed, the iteration (30) may fail to converge when (31) is not satisfied.See Section 4.2 for more details.In addition, the preprocessing of initial data proposed in [19] is no longer required when (31) is satisfied.second one is made up of simple shapes (Image II) used in the literature [15], and the last one is Lena image (Image III).

Numerical Experiments and Discussion
We have experiments on restoring the three images blurred by the following three different blurring operators [23].
Blur operator I: Blur operator II: an out-of-focus blur with the kernel of convolution given by the kernel where  = 1.2188.Blur operator III: a truncated Gaussian blur given by where  = 0.176 and  = 0.38.

The Linear Stabilizing
Term.The condition (31) serves as a general guideline for the choice of the matrix.In this section, we will consider a quinta-diagonal matrix  to be compared with a diagonal matrix  [23].

Diagonal Matrix 𝐵.
In order to solve (30) efficiently, an obvious candidate is diagonal matrices of the form The condition (31) then reads To proceed with the estimate of  max (), we notice that our sample blurring operators I-III can be represented by a mask of the form with It is easy to see that, under the Neumann boundary conditions, the largest eigenvalue for such  is given by where the second term of (49) is independent and the corresponding eigenfunction is a constant.In summary, we have the following sufficient condition for (31): Similar conditions can be derived for variants of (50), such as  or In our examples, the diagonal entries of  *  are constant-valued except for the near-boundary pixels.We expect the performances of (50), (51), and (52) to be comparable to each other.In this paper, most of our numerical experiments are conducted using (52) with  = 1.See Section 5 for more details.

Quinta-Diagonal Matrix 𝐵.
In addition to diagonal matrices, we have also explored quinta-diagonal matrices  , given by the mask , has the same support as laplace operator Δ; therefore the computational cost in the AMG step is comparable to diagonal .In view of (33) and (31), it is clear that the optimal  (,) * is the one that solves min ,  ((−Δ +  (,) ) −1 ( (,) −  * )) where () is the spectral radius of .The actual minimizer of (54) depends on the true image  (∞) and there is no simple way of finding it.An accessible approximation is given by the following modified minimization problem: In our examples, the common eigenbasis of  , and  is given by It is straightforward to compute the corresponding eigenvalues for ( (,) −  * ) and the numerical solution of (55) by varying  and  over a suitable range.The optimal (, ) * , together with the critical values of , , and , is given in Table 1.It is not clear a priori why (55) should result in better performance.Nevertheless, we find from our experiment that  (,) * indeed performs no worse than diagonal  and can be significantly faster in some cases.See more details in Section 5.

The Stopping Criterion and Acceleration Technique.
The outer iteration is stopped upon a relative decrease of the (normalized) residual by a factor of 10 −4 for the blurring operators, namely, when the condition       ()      2      (1)     2 ≤ 10 −4 (57) is reached.Here the normalization of the residual is given by where V (+1) is an approximate solution of (28) obtained with one -cycle iteration (see also Section 5) and  (+1) = diag(−Δ + ).Note that the larger entries of the normalizing factor ( (+1) ) −1 correspond to regions where V (+1) is less smooth.It amplifies the (unnormalized) residual on regions where V (+1) either has a jump or requires further denoising, therefore doing a better job measuring the image quality.Numerical experiments confirmed that the normalized residual is indeed a good indicator for the quality of denoising and deblurring.
To accelerate and stabilize the outer iteration, we have incorporated the Krylov subspace method [19,28,29] in our implementation.The approximate solution is optimized every  steps on a subspace of dimension  ≤ .To be more

Numerical Results
To generate an observed image , a Gaussian white noise with mean 0 and variance  is added to the blurred image.The level of noise can be measured by the signal-to-noise ratio (SNR): The signal-to-blur ratio (SBR) measures the strength of the blurring operator .The two kinds of error 2 = ‖ +1 −   ‖  2 and 0 = ‖ +1 −   ‖  ∞ are used to judge the iterative error as the iteration ends.The cost of the iterative process is denoted by "time." These indexes are captioned in each example for readers' convenience.
In all our examples, the iteration begins with the observed image  (0) =  and all images are contaminated with small noise  = 10 and large noise  = 40.We choose  according to (52) with  = 1.0 except in Table 6 where we demonstrate the necessity of the linear stabilizing term and the overall efficiency of our scheme by varying .It is enough to apply one single -cycle in the AMG step with the Gauss-Seidel iteration as the smoother.We apply the Krylov acceleration procedure every 4 steps with the Krylov dimension  = 2.The rate of convergence of the outer iteration can be significantly improved by the Krylov subspace technique.

Image Denoising. Table 2 shows the contrast results of
Image II with the split Bregman method for isotropic TV model (Algorithm 2), the split Bregman method based on AMG (Algorithm 3), and the split Bregman method for anisotropic TV model (Algorithm 4). Figure 2 shows the denoising image for Image II with noise contamination of  = 10 (SNR = 14.90%).Table 3 shows the contrast results of Image III containing 512 * 512 pixels.Figure 3 shows the denoising image for Image III (512 * 512) with noise contamination  = 40 (SNR = 26.725%).It is clear that our algorithm (Algorithm 3) is efficient and suitable for elliptic equation.

Image Denoising and Deblurring. The compared results
are given in Tables 4-5. Figure 4 shows the denoising and deblurring images for Image I contaminated by noise with  = 40 and blurring operator III (SNR = 51.34%,SBR = 45.31%).Figure 5 shows the denoising and deblurring images for Image II contaminated by noise with  = 10 and blurring operator I (SNR = 5.96%, SBR = 48.79%).We continue our numerical experiments by varying the parameter  and compare their performance.The number of iterations of a typical example is shown in Table 6.The combination of  ≈  * and  = 2 gives the optimal result among possible choices of  and  in general.The number of iterations is insensitive to the variation of  near  * and gradually grows as  increases.For simplicity of presentation, we take  = 1 and  = 2 in all other numerical examples.We have also included the results for  =  (,) * in Table 6 for comparison.In general,  =  (,) * performs no worse than (52) with  =  * .In some cases,  =  (,) * can be up to 25% faster.In addition, if the blur operator is mild, our method can be used to solve  =  directly.
According to the previous figures and tables, we get some conclusions here.First of all, in the second part of (42), we have neglected the −Δ term to obtain a sufficient condition (31), from which we further derived the critical value  * in (52).We therefore expect  * to be slightly overestimated.In other words, the outer iteration should converge for  ≥  * , but not necessarily for  <  * .This is in accordance with our numerical result.See the row with  = 0 in Table 6.Secondly, the Krylov subspace technique helps to stabilize the outer iteration for  <  * , but eventually fails if  becomes too small.This result demonstrates the necessity of the linear stabilizing term.

Conclusion
In this paper, we propose a new algorithm by adding a linear term for the total variation-based image denoising and deblurring which combine the split Bregman method, the algebraic multigrid method, and Krylov subspace acceleration.Through formal convergence analysis, we derived an explicit formula for a linear stabilizing term.Numerical experiments demonstrate that our algorithm is efficient and robust for deblurring and denoising problems.

Figure 1 :
Figure 1: Original images.Image I, Image II, and Image III.

Table 3 :
Results of Figure3.The numerical results demonstrate that the Krylov acceleration method is an efficient way to accelerate the outer iteration.See Section 5.