Iterative Methods for Nonlinear Equations or Systems and Their ApplicationsView this Special Issue
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 Bregman method is used to convert it into linear system, and an algebraic multigrid method 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.
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. : Here is the penalty parameter, and is the diffusion regularization parameter and is typically small. The functional in (2) is strictly convex with a unique global minimizer. The paper  showed the well-posedness of problem (2) as 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 inner product.
In addition to the isotropic models, the anisotropic models for a qualitative improvement at corners are shown in  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  is The existence and uniqueness of solutions of the anisotropic total variation flow are shown in . The explicit time-marching schemes are applied in [1, 5]. Zuo and Lin  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 is maintained by GAPG, where is the number of iterations. The introduction of an anisotropy to TV regularization  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 , the affine scaling algorithm , the Newton’s method with a continuation procedure on , a fixed point iteration with multigrid method [10–12], the combination of fixed point iteration and preconditioned conjugate gradient which method is proposed in , a fast algorithm based on the discrete cosine transform [14, 15], a new modular solver , a proximity algorithm , and a nonlinear primal-dual method . 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 , 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 , the authors proposed the split Bregman method to solve image denoising and compressed sensing problem rapidly. Chang et al.  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 . Motivated by that, a wider linear term different from that in  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.
2. 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. ). We partition the domain into uniform cells with mesh size . 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 : 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 wide-banded, and the spectra of the matrices and are quite differently distributed. It is difficult to compute (10). In , 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 , 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, and . How to construct the interpolation operators is the most important part in the method. The convergence proof for this improved AMG method was given in  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 , 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 where Consider it into the following three subproblems:
The formula (16) is -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:
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.
2.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.
2.3. The Split Bregman Method for Anisotropic TV Model
The anisotropic total variation of  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  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 , so the procedure is rapid, simple, and not costly.
3. 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 with . 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).
Theorem 1. If the matrix satisfies the equation (30) is locally convergent.
Let be the solution of (29); that is, We then define , and the error equation is which is equivalent to with It is easy to show that the matrices and are symmetric. Taking the inner product on both sides of (34) with gives and we will now rearrange the formula (36) as Rewrite the previous formula as If and are symmetric and positive definite, it follows from (38) that where with positive definite matrix . The behavior of (39) leads to In other words, the iteration (30) is locally convergent if both are symmetric and positive definite which means 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  is no longer required when (31) is satisfied.
4. Numerical Experiments and Discussion
4.1. Images and Blurring Operators
The numerical examples given in this paper are based on the three images shown in Figure 1. The first one is a satellite image (Image I), the second one is made up of simple shapes (Image II) used in the literature , and the last one is Lena image (Image III). We have experiments on restoring the three images blurred by the following three different blurring operators .
Blur operator I:
Blur operator II: an out-of-focus blur with the kernel of convolution given by the kernel where .
Blur operator III: a truncated Gaussian blur given by where and .
4.2. The Linear Stabilizing Term
4.2.1. Diagonal Matrix
To proceed with the estimate of , 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 . See Section 5 for more details.
4.2.2. 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 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.
4.3. The Stopping Criterion and Acceleration Technique
The outer iteration is stopped upon a relative decrease of the (normalized) residual by a factor of for the blurring operators, namely, when the condition is reached. Here the normalization of the residual is given by where is an approximate solution of (28) obtained with one -cycle iteration (see also Section 5) and . Note that the larger entries of the normalizing factor correspond to regions where is less smooth. It amplifies the (unnormalized) residual on regions where 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 precise, we take two fixed integers , and for any integer , let The residual of can be approximated by
One can minimize with respect to to get and reset to before going to the outer iteration. The numerical results demonstrate that the Krylov acceleration method is an efficient way to accelerate the outer iteration. See Section 5.
5. Numerical Results
To generate an observed image , a Gaussian white noise with mean 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 and 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 and all images are contaminated with small noise and large noise . We choose according to (52) with 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 . The rate of convergence of the outer iteration can be significantly improved by the Krylov subspace technique.
5.1. 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 (SNR = 14.90). Table 3 shows the contrast results of Image III containing pixels. Figure 3 shows the denoising image for Image III () with noise contamination (SNR = 26.725). It is clear that our algorithm (Algorithm 3) is efficient and suitable for elliptic equation.
5.2. 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 and blurring operator III (). Figure 5 shows the denoising and deblurring images for Image II contaminated by noise with and blurring operator I ().
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 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 and 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 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.
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.
The authors thank the referees for very useful remarks and suggestions. The research is supported by NSFC (no. 10801049 and no. 11271126) and Foundation of North China Electric Power University.
L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D, vol. 60, no. 1–4, pp. 259–268, 1992.View at: Google Scholar
Y. Li and F. Santosa, “An affine scaling algorithm for minimizing total variation in image enhancement,” Tech. Rep. 12/94, Center for Theory and Simulation in Science and Engineering, Cornell University, 1994.View at: Google Scholar
R. Chan, T. Chan, and H. Zhou, “Advanced signal processing algorithms,” in Proceedings of the International Society of Photo-Optical Instrumentation Engineers, F. Luk, Ed., pp. 314–325, SPIE, 1995.View at: Google Scholar
M. E. Oman, “Fast multigrid techniques in total variation-based image reconstruction,” in Proceedings of the Copper Mountain Conference on Multigrid Methods, 1995.View at: Google Scholar
R. H. Chan, T. F. Chan, and C. K. Wong, “Cosine transform based preconditioned for total variation deblurring,” IEEE Transactions on Image Processing, vol. 8, no. 10, pp. 1472–1478, 1999.View at: Google Scholar
Q. Chang, L. Chien, W. Wang, and J. Xu, “A robust algorithm variation deblurring and denoising,” in Proceedings of the 2nd International Conference on Scale Space and Variational Methods in Computer Vision (SSVM '09), vol. 5567 of Lecture Notes in Computer Science, Springer, 2009.View at: Google Scholar