Two new preconditioners for mean curvature-based image deblurring problem

1.   Introduction

In the last two decades, the nonlinear variational methods have received a great deal of attention in the field of image deblurring. Two main difficulties arise while applying nonlinear variational techniques to large-scale noisy and blurrred images. One of them, of course, is the nonlinearty and other is the solution of the large system which arise from the discretization of these such problems. The main focus of this paper is to handle these two computational difficulties. The most well-known nonlinear variational image deblurring model is total variation (TV) model ^[1,10,14]. It is one of the most commonly used regularization model and it has nice properties like edge preserving. But the main drawback of the TV model is that the resulting images look blocky. Because this model converts smooth functions into piecewise constant functions which create staircase effects in resulting images. To reduce the staircase effects one remedy is to use the mean curvature (MC) ^{[4,12,16,18,19]} based regularization models.

The MC-based regularization models are widely used in all image processing problems. In image deblurring, the MC-based models are very effective. These models not only preserve edges but also remove staircase effect in the recovery of digital images. However, the discretization of the associated Euler-Lagrange equations produce a large nonlinear ill-conditioned system which affect the convergence of the numerical algorithms like Krylov subspace methods (GMRES etc.). Furthermore the Jacobian matrix of the MC-based nonlinear system is block banded matrix with large bandwidth. The MC-based regularization methods are effective, but due to the high nonlinearity and ill-conditioned system, robust and fast numerical solution is a crucial issue. To overcome these difficulties, in this paper we introduce two new symmetric positive definite (SPD) preconditioners. So instead of applying ordinary Generalized Minimal Residual (GMRES) method we use PGMRES (preconditioned Generalized Minimal Residual) method to solve the system. Fast convergence has shown in the numerical results by using proposed new preconditioners.

The contributions of the manuscript include the following:

(i) our work presents an efficient algorithm for mean curvature-based image deblurring problem which combines a fixed point iteration with new SPD preconditioned matrices to handle the nonlinearity and ill-conditioned nature of the large system;

(ii) presents a better treatment for the computationally expensive higher order and nonlinear mean curvature regularization functional.

The paper is organized in different sections. The first section includes introduction while the second section includes problem description of image deblurring model. In the third section, we present nonlinear system of first order equations, cell discretization and CCFD method. In fourth section, we introduce our proposed SPD preconditioners and our algorithmic technique. The properties and eigenvalue analysis of the proposed preconditioned matrices are discussed in the fourth section. The numerical experiments are also presented in the fourth section. The conclusions about the proposed PGMRES method is discussed in the last section of the paper.

2.   Problem description

The focus of the paper is on image debluring problem, so we start by presenting its concise description. Mathematically, the relationship between $u$ (original image) and $z$ (recorded image) is as follows;

where $\varepsilon$ is the noise function. The noise can be Gaussian noise, salt and pepper noise, Brownian noise etc. In this paper, we have considered Gaussian noise. The $\vec{K}$ is the blurring operator which is a Fredholm-integral operator of first kind;

having translation invariance property i.e., $k(x, y) = k(x - y)$. If $\vec{K} = I$ problem (2.1) is called image denoising problem. The problem (2.1) becomes ill-posed ^[1,13,14] because $\vec{K}$ is a compact operator. Let $\Omega$ be a square in $\text{R}^{2}$ and $u$ be an image intensity function defined on $\Omega$. The $\textbf{x} = (x, y)$ defines the position in $\Omega$. Let $\left|\textbf{x}\right| = \sqrt{x^2 + y^2}$ is an Euclidean norm and $\|. \|$ is $L_2(\Omega)$ norm. The problem (2.1) is an inverse problem. The recovering of $u$ from $z$ makes (2.1) an unstable problem ^[1,13,14]. To make it stable, one remedy is to use the MC regularization functional ^{[4,12,16,18,19]},

where $\kappa$ is called the mean curvature surface. Then the problem (2.1) takes the form, find $u$ that minimizes the following problem

where $\alpha > 0$ is a regularization parameter. The well-posedness of the problem (2.2) for a particular case (synthetic image denoising problem) is explained in ^[18]. Then the Euler-Lagrange equations associated with (2.2) are,

where $\vec{K}^*$ denotes the adjoint operator of $\vec{K}$ and $\beta > 0$ is used to avoid non-differentiability at zero. The equation (2.3) is a nonlinear fourth order partial differential equation.

As it is known that the MC-based model does not only preserve edges but also removes staircase effect in the recovery of digital images. However, fourth order derivatives appear in the Euler-Lagrange equations, which create problems in developing an efficient numerical algorithm. One key problem in presenting the method is to give a proper approximation to the nonlinear mean curvature functional. We have treated this difficulty by reducing the nonlinear fourth order Euler-Lagrange equation into a system of first order equations.

3.   The first order nonlinear system

The Eq (2.3) can be expressed as first order nonlinear system,

where

We will take advantage of this special structure to derive our proposed algorithm.

3.1.   Cell discretization

For MC-based image deblurring problem, the domain $\Omega = (0, 1)\times (0, 1)$ is partitioned by $\delta_x \times \delta_y,$ ^[11], where

where $n_x$ represents the number of equispaced partitions in the $x$ or $y$ directions and $(x_i, y_j)$ denotes centers of the cells. The

where $h = \frac{1}{n_x}$. The $(x_{i \pm \frac{1}{2}}, y_j)$ and $(x_i, y_{j \pm \frac{1}{2}})$ are representing midpoints of cell edges,

The set

represents a cell with $(x_i, y_j)$ as a center. Let

And

Approximations of $u$ and $w$ are

and

respectively, where $u_{ij} = \overline {U}(x_i, y_j)$ and $w_{ij} = \overline{W}(x_i, y_j)$. The representation of the data $z$ is

where $z_{ij}$ can be calculated at cell averages. By applying midpoint quadrature approximation, we have

Denote $\overrightarrow{ v } = (v^x, v^y), \overrightarrow{ p } = (p^x, p^y)$ and $\overrightarrow{ t } = (t^x, t^y)$. The approximation of $x$ and $y$ components of $\overrightarrow{ v }$ are

respectively. $\overline{V} = [\overline{V}^x$ $\overline{V}^y]^t$ denotes the discretization of $\overrightarrow{ v }$. Similarly, approximation of the components of $\overrightarrow{ p }$ and $\overrightarrow{ t }$ are

and

respectively. The $\overline{P} = [\overline{P}^x$ $\overline{P}^y]^t$ and $\overline{T} = [\overline{T}^x$ $\overline{T}^y]^t$ denote the discretization of the vectors $\overrightarrow{ p }$ and $\overrightarrow{ t }$ respectively.

3.2.   The CCFD method

Here, we present the cell-centered finite difference (CCFD) method for MC-based image deblurring problem. By considering lexicographical ordering of the unknowns,

Now by applying CCFD method to (3.1)–(3.5) we obtain the following system,

where midpoint quadrature rule is used for the integral term. The matrices $K_h, A_h$ and $I_h$ (the identity matrix) each one is of size $n_x^2 \times n_x^2$. The matrix $B_h$ is of size $2n_x(n_x - 1) \times n_x^2$. The matrices $C_h$ and $D_h$ are of size $2n_x(n_x - 1) \times 2n_x(n_x - 1)$. So we have the following system

The matrix $K_h$ is block Toeplitz with Toeplitz blocks (BTTB) and $K^*_h K_h$ is SPD. The matrix $A_h$ is a diagonal matrix having following structure,

where both $A_1$ and $A_2$ are of size $n_x^2 \times n_x^2$.

where $\otimes$ is a tensor product. The size of identity matrix $\tilde{I}$ is $n_x \times n_x$. The matrix

is of size $n_x\times n_x$. The matrix $B_h$ has the following structure,

where both $B_1$ and $B_2$ are of size $n_x(n_x - 1) \times n_x^2$, and

is a matrix of size $(n_x - 1) \times n_x$. The matrix

is a diagonal matrix and its entries are obtained by the discretization of the expression $(\bigtriangledown w. \overrightarrow{ v })$. The matrix $C^x$ is of size $(n_x - 1) \times n_x,$ and the matrix $C^y$ is of size $n_x \times (n_x - 1)$. The matrix $D_h$ is also a diagonal matrix with positive diagonal entries and its entries are obtained by the discretization of the expression $\sqrt{\left|\bigtriangledown u\right|^2+\beta^2}$. The matrix $D_h$ has the following structure,

where $D^x$ is of size $(n_x - 1) \times n_x,$ and $D^y$ is of size $n_x \times (n_x - 1)$. Note that on horizontal and vertical edges of each cell $e_{ij}$, the values of the all unknowns are not available, so average operators are used to calculate their values.

Now if we eliminate $W, V, P$ and $T$ from (3.6)–(3.10), then we have the following nonlinear primal system,

where

The $L_h$ is block pentadiagonal matrix according to the lexicographical ordering of the unknowns. The main diagonal blocks are pentadiagonal while the off-diagonal blocks are tridiagonal matrices.

In the literature, one can find a number of numerical methods that have been investigated to mean curvature-based nonlinear image minimization problems ^{[2,4,8,12,14,15,16,18,19]}.

Since MC-based models produce a large and ill-conditioned nonlinear system, so almost all these standard methods get quite slow convergence. Furthermore the presence of higher order and nonlinear mean curvature regularization functional in the governing equation of the models makes these highly nonlinear systems harder for calculation. MC is very much computationally expensive, that is why, most of the existing methods performs quite poorly. To make a robust numerical method for MC-based nonlinear image deblurring problem, now we present a new preconditioned numerical method.

4.   Numerical implementation

Here we introduce the algorithms to solve the MC-based nonlinear system (3.11). First we apply discrete version of the FPI (fixed point iteration) to (3.11) to handle the nonlinearity of MC. The approach taken here is called "lagged diffusivity" scheme ^[14]. Its rate is just linear but in practice it has a quite rapid convergence. Furthermore, this scheme does not depends on initial guess to converge globally. This is why globalization is not an issue for this scheme. So by FPI we have a following linear system;

4.1.   Properties

Before proceeding further, we discuss the some important properties of our system (4.1).

(1) The Hessian matrix $K^*_h K_h +\alpha L_h$ is extremely large for practical application. When $\alpha$ is small, the Hessian matrix becomes quite ill-conditioned. This happens because the eigenvalues of the blurring operator $\vec{K}$ cluster to zero ^[14].

(2) The first term, $K^*_hK_h$, in the Hessian matrix is symmetric positive definite. Although, $K^*_hK_h$ is full but the blurring operator $\vec{K}$ has translation invariant property, which permit the use of FFT (Fast Fourier transformation) to evaluate $K^*_hK_h u$ in $O(n \log n)$ operations ^[14].

(3) The second term, $L_h$, in the Hessian matrix is sparse but not symmetric. The $L_h$ in (3.11) consists of three terms. The first and the last term in $L_h$ are symmetric positive semidefinite ^[14] but the middle term is not symmetric. Hence the system (4.1) is not symmetric positive definite.

4.2.   The Preconditioner

According to the properties of our system (4.1), mentioned above, GMRES (Generalized Minimal Residual) method, is appropriate for the solution of the system (4.1). Due to ill-conditioned system, GMRES can be quite slow to convergence. So we use preconditioned Generalized Minimal Residual (PGMRES) method. For effective solution, preconditioning matrix $P$, must be symmetric positive definite (SPD). Here we introduce following two SPD preconditioned matrices namely $P_1$ and $P_2$. For the first preconditioner $P_1$, we consider the following matrix

The $L^{TV}_h$ is the SPD matrix lies in the first term of $L_h$ matrix (4.1). In fact, $L^{TV}_h$ arise by the discretization of total variation (TV) based image deblurring problem ^[14]. So we introduce preconditioned matrix $P_1$ as

where $I_h$ is an identity matrix and $diag(L^{TV})$ is a diagonal matrix and its entries are the diagonal entries of matrix $L^{TV}_h$. The $\gamma$ is the positive parameter. For the second preconditioner $P_2$, we consider the following SPD matrix

The $L^s_h$ is the SPD part (first and the third term) of $L_h$ matrix (4.1). So the second preconditioned matrix is

where $diag(L^s_h)$ is also a diagonal matrix whose entries are the diagonal entries of matrix $L^s_h$. While applying PGMRES to (4.1), the inversion of preconditioner matrix, will be required. Since the both terms in our preconditioning matrices ($P_1$ and $P_2$) are diagonal matrices, so inversion can be done easily. We have summarize the PGMRES method in the following algorithm.

——————Algorithm 1: The PGMRES Method——————

On mesh $\Omega_h$,

Step-I: Initial iteration $U^0,$

Step-II: for $m = 1; max$

Step-III: Use PGMRES to solve

where $P = P_1, P_2$.

end

————————————

4.3.   Eigenvalues

Now, let the eigenvalues of $K^*_h K_h$ and $L_h$ be $\lambda_i^K$ and $\lambda_i^L$ respectively such that $\lambda_i^K \downarrow 0$ and $\lambda_i^L \uparrow \infty$. So the eigenvalues of $P^{-1}_1\bar{A}$ and $P^{-1}_2\bar{A}$ are

respectively. Here $\lambda_i^{L^{TV}}$ are eigenvalues of $L^{TV}_h$ and $\lambda_i^{L^{ns}}$ are eigenvalues of $L^{ns}_h$. The $\bar{A} = K^*_h K_h +\alpha L_h$ is the Hessian matrix of the system (4.1). One can notice that $\lambda_i^{L^{TV}} \leq \lambda_i^{L^{ns}} \leq \lambda_i^L$. So for $\gamma > \alpha$,

Hence, for large values of $\gamma$, $P^{-1}\bar{A}$ has more favourable spectrum as compared to the Hessian matrix $\bar{A}$, so PGMRES converge more rapidly than ordinary GMRES. This can also be shown in the numerical examples that PGMRES is getting rapid convergence with large $\gamma > \alpha$.

4.4.   Numerical experiments

Now we present four numerical examples for MC-based image deblurring problem. We have used different values of $n_x$ so the resulting system has $n_x^2$ unknowns. For numerical computations, MATLAB software is used to obtain the numerical results on Intel(R) Core(TM) i7-4510U with CPU @ 2.00 GHz 2.60 GHz. To measure the quality of the deblurred images, we have used PSNR (Peak Signal to Noise Ratio) and SSIM (Structured Similarity Index Measure).

Initial guess. Since the method is globally convergent, the fixed point iteration (FPI) can be started with any initial guess ^[14]. In all experiments we have used $u^0 = z$ (the blurred data) as initial guess. For instance, taking the zero initial guess $u^0 = 0$ means that one extra FPI may be required to attain the accuracy obtained with $u^0 = z$.

Selection of the parameters. Although the automated selection of the parameter $\alpha$ is beyond the scope of the present paper, that information can be found in ^[3,18,19]. The purpose of $\alpha$ is to select the amount of blur / noise to be removed, so there is no point in choosing either a very large or very small $\alpha$. The optimum range of $\alpha$ for MC-based model lies between $1e - 3$ and $1e - 8$. To understand the performance of the parameter $\beta$ on our algorithm (GMRES), we have run the algorithm with fixed $\alpha = 1e - 8$ on the benchmark Cameraman image of size $256 \times 256$. In this experiment we have varied $\beta$ from $1e - 3$ to $1$ and found that the smaller value of $\beta$ does not improve the PSNR and SSIM. The results of this experiments summarized in Table 1.

Example 1. In this example we have used Cameraman image. This is a complicated image, because it has a small scale texture part (shirt) and large scale cartoon part (face). The different aspects of Cameraman image are presented in Figure 1. The size of each subfigure is $256 \times 256$. These are (a) blurry image (b) deblurred image by GMRES (c) deblurred image by $P_1$GMRES with $\gamma = 1000$ and (d) deblurred image by $P_2$GMRES with $\gamma = 1000$. The calculation of relative residual at each iterations against different values of the parameter $\gamma$ is presented in Figure 2. The $fspecial('gaussian', [n_x \quad n_x], 1.5)$ kernel is used for numerical calculations. It is a $gaussian$ filter of size $n_x \times n_x$ with standard deviation $(\sigma = 1.5)$. The parameters $\alpha = 1e - 8$ and $\beta = 1$. For the stopping criteria of a numerical method we have used tolerance $tol = 1e-2$.

Remarks

(1) The Figures 1(b)–(d) are almost similar, this means both GMRES method (without preconditioning) and PGMRES ($P_1$GMRES and $P_2$GMRES) method are generating same quality results.

(2) From the Figure 2, one can clearly observe the effectiveness of preconditioning. Here, result were presented using fixed point iteration count $m = 1$. The number of PGMRES iteration is much lesser than as compared to GMRES (without preconditioning) to reach the required accuracy $tol = 1e-2$. The results for subsequent fixed point iterations are also similar.

(3) From the Figure 2, this can also be observed that PGMRES is getting rapid convergence for large values of $\gamma$ for both preconditioners ($P_1$ and $P_2$).

(4) From Table 2, one can observe that the PSNR and SSIM by PGMRES method is same as compared with the ordinary GMRES method. But PGMRES method is producing same PSNR and SSIM in quite less iterations. For example, to achieve the same PSNR and SSIM, the $P_1$GMRES method needs only 80 iterations with $\gamma = 10$ and the $P_2$GMRES method needs only 42 iterations with $\gamma = 10$. But the ordinary GMRES method needs more than 150 iterations to get the same PSNR and SSIM. The number of iterations further decreases with higher values of $\gamma$. Which means that PGMRES method is faster than the ordinary GMRES method.

(5) From Table 2, it is also observed that for both preconditioners ($P_1$ and $P_2$) the PSNR and SSIM do not get improvement with the increase in the value of $\gamma$. The higher value of $\gamma$ only decreases the number of iterations for both preconditioners ($P_1$ and $P_2$). But the decrease in the number of iterations is not much significant for the value of $\gamma$ between $500$ and $1000$. So we may suggest that suitable value of $\gamma$ lies between $500$ and $1000$.

Example 2. In this example we have used Goldhills image. This is real and synthetic image. Here we have compared our MC based algorithm with TV (total variation) based algorithm. Since TV-based model generates a SPD matrix system, so for the solution we have used CG (Conjugate Gradient) method. The different aspects of Goldhills image are presented in Figure 3. The size of each subfigure is $512 \times 512$. These are (a) exact image (b) blurry image (c) deblurred image by CG (d) deblurred image by GMRES (e) deblurred image by $P_1$GMRES and (f) deblurred image $P_2$GMRES. For numerical calculations, we have used the $ke_{-}gen (n_x, 300, 5)$ kernel ^[5,6,7,9]. It is a circular $gaussian$ filter of size $n_x \times n_x$ with radius $(r = 300)$ standard deviation $(\sigma = 5)$. The $ke_{-}gen (120, 40, 4)$ kernel is depicted in Figure 4. For TV-based method we have used $\alpha = 5e - 5$ and $\beta = 1$ according to ^[14]. For MC-based method we have used $\alpha = 1e - 8$ and $\beta = 1$ and $\gamma = 1000$. For the stopping criteria of a numerical methods we have used tolerance $tol = 1e-4$.

Remarks

(1) From Table 3, it is observed that the PSNR and SSIM by MC-based (GMRES and PGMRES) methods are little higher than TV-based CG method. Same comparison can be observed from Figures 3(c)–(f). So MC-based (GMRES and PGMRES) methods are generating better quality results.

(2) From the Figure 5, one can clearly observe the effectiveness of preconditioning. The number of MC-based PGMRES iteration is much lesser than as compared to both MC-based GMRES and TV-based CG methods to reach the required accuracy $tol = 1e-4$. The results for subsequent fixed point iterations are also similar.

(3) From Table 3, it is observed that the PSNR and SSIM by MC-based PGMRES method is almost same as compared with the ordinary MC-based GMRES method. But PGMRES method is generating same PSNR and SSIM in quite less iterations. To achieve the same PSNR and SSIM, the $P_1$GMRES method needs only 90 iterations and the $P_2$GMRES method needs only 40 iterations. But the GMRES method needs more than 200 iterations to get the same PSNR and SSIM. The TV-based CG method is also getting more than 200 iterations to get its PSNR and SSIM. Which means that MC-based PGMRES method is faster than the MC-based GMRES and TV-based CG method for real and synthetic images. Here, the performance of the preconditioner $P_2$ is again more effective than the preconditioner $P_1$.

Example 3. Here we have used nontexture Moon image. Here we have also compared our MC based algorithm with TV (total variation) based algorithm. For TV we have used CG (Conjugate Gradient) method. The different aspects of Moon image are presented in Figure 6. The size of each subfigure is $256 \times 256$. These are (a) exact image (b) blurry image (c) deblurred image by CG (d) deblurred image by GMRES (e) deblurred image by $P_1$GMRES and (f) deblurred image $P_2$GMRES. For numerical calculations, we have also used the $ke_{-}gen (n_x, 300, 5)$ kernel. For TV-based method we have used $\alpha = 1e - 4$ and $\beta = 1$ according to ^[14]. For MC-based method we have used $\alpha = 1e - 6, \beta = 1$ and $\gamma = 1000$. For comparison we have used three different values of $n_x$. These are 64,128 and 256. For the stopping criteria of a numerical methods we have used tolerance $tol = 1e-5$.

Remarks

(1) The Figures 6(b)–(d) are almost similar, so all methods are generating same quality results.

(2) From the Figure 7, one can clearly notice the effectiveness of preconditioning. For all values of $n_x$, the number of PGMRES ($P_1$GMRES and $P_2$GMRES) iteration is much lesser than as compared to MC-based GMRES and TV-based CG to reach to the required accuracy $tol = 1e-5$. The results for subsequent fixed point iterations are also similar.

(3) From Table 4, it is observed that the PSNR and SSIM by MC-based PGMRES method is almost same as compared with the ordinary MC-based GMRES method but much higher than TV-based CG method for all values of $n_x$. But PGMRES method is generating this better PSNR and SSIM in quite less iterations. For example, to achieve the better PSNR and SSIM, the $P_1$GMRES method needs only 40 iterations and the $P_2$GMRES method needs only 24 iterations for $n_x = 64$. But the GMRES method needs more than 150 iterations to get the same PSNR and SSIM. The TV-based CG method is also getting more than 150 iterations to get its lower PSNR and SSIM. Same is the case for other values of $n_x$. Which means that MC-based PGMRES method is faster than the MC-based GMRES and TV-based CG method for nontexture images. Here, the performance of the preconditioner $P_2$ is again more effective than the preconditioner $P_1$

Example 4. In this example, we have compared our MC-based algorithm with MC-based augmented lagrangian methods of (Zhu, Tai and Chan) ^[19,20] and (Zhang, Deng, Shi, Wang and Yonggui) ^[17]. We have done the comparison for image denoising case (when blurring operator $\vec{K} = I$ identity operator) because in the given literature ^[17,19,20] augmented lagrangian method proposed only for image denoising problem. For simplicity, we have named augmented lagrangian methods in ^[19,20] and ^[17] by FFTALM, GSALM and LALM, respectively. For the comparison we have used Lena and Cameraman images. The different aspects of both images are presented in Figure 8 and 9. The size of each subfigure is $256 \times 256$. The MC-based augmented lagrangian methods (FFTALM, GSALM and LALM) use certain sets of parameters which we have used according to ^[17,19,20]. For FFTALM, the parameters are $\alpha = 8e - 1, r_1 = 3e-4, r_2 = 50, r_3 = 5e+1, r_4 = 5e+4$ and $\beta = 1$. For GSALM, the parameters are $\alpha = 8e - 1, r_1 = 4e-3, r_2 = 50, r_3 = 1e+3, r_4 = 5e+4$ and $\beta = 1$. For LALM, the parameters are $\alpha = 8e - 1, r_1 = 4e-3, r_2 = 50, r_3 = 1e+3, r_4 = 5e+4, \beta = 1, \delta_1 = 0.0038$ and $\delta_2 = 0.00012$. For our MC-based algorithm (GMRES and PGMRES) we have used $\alpha = 1e - 6, \beta = 1$ and $\gamma = 1000$. In both noisy images we have used random noise with mean $(\mu = 0)$ and standard deviation $(\sigma = 10)$. For the stopping criteria of a numerical methods we have used tolerance $tol = 9e-4$ for Lena image and $tol = 8e-4$ for Cameraman image.

Remarks

(1) From the Figures 8 and 9 one can notice that all methods are generating same quality results.

(2) From the Table 5, it is observed that the PSNR value of all methods is almost same for both images. But our methods (GMRES, $P_1$GMRES and $P_2$GMRES) are generating same PSNR in quite less CPU-Time as compared to all augmented lagrangian methods (FFTALM, GSALM and LALM). For example, the $P_1$GMRES method is generating its PSNR for Lena image in just 6.528 seconds. But FFTALM, GSALM and LALM are generating their PSNR in 19.515, 11.406 and 9.078 seconds, respectively. Which means $P_1$GMRES method is saving more than 60$\%$ of CPU-Time as compared to FFTALM, more than 50$\%$ of CPU-Time as compared to GSALM and more than 30$\%$ of CPU-Time as compared to LALM. Same is the case for $P_2$GMRES method. So we can say, our proposed MC-based methods (GMRES, $P_1$GMRES and $P_2$GMRES) are faster than the MC-based augmented lagrangian methods (FFTALM, GSALM and LALM).

5.   Conclusions

A numerical algorithm (PGMRES) is presented to solve the primal form of mean curvature-based nonlinear image deblurring problem. Two new SPD preconditioner matrices ($P_1$ and $P_2$) are introduced. Four examples are tested by PGMRES using our new preconditioner matrices. Different kinds of images (Complicated, real, synthetic and nontexture)are tested by PGMRES using our new preconditioner matrices. The convergence rates for solution to the linear system and the norm of the residual at each iteration are presented in each example. The comparison between our proposed MC-based methods (GMRES, $P_1$GMRES and $P_2$GMRES) and MC-based augmented lagrangian methods (FFTALM, GSALM and LALM) is also presented. Numerical experiments show the rapid convergence of PGMRES method using new preconditioners. Which means that PGMRES method is a robust and faster numerical method. It is also observed that the performance of the preconditioner $P_2$ is more effective than the preconditioner $P_1$.

Acknowledgments

The authors would like to thank King Fahd University of Petroleum and Minerals (KFUPM) for its continuous supports. The authors also thank the referees for their very careful reading and valuable comments. This work is funded by KFUPM under Project #SB201012.

Conflict of interest

The authors declare that there is no conflict of interest.