Circulant preconditioners for mean curvature-based image deblurring problem

Abstract

The mean curvature-based image deblurring model is widely used to enhance the quality of the deblurred images. However, the discretization of the associated Euler–Lagrange equations produces a nonlinear ill-conditioned system which affects the convergence of the numerical algorithms such as Krylov subspace methods (generalized minimal residual etc.) To overcome this difficulty, in this paper, we present three new circulant preconditioners. An efficient algorithm is presented for the mean curvature-based image deblurring problem, which combines a fixed point iteration with new preconditioned matrices to handle the nonlinearity and ill-conditioned nature of the large system. The eigenvalues analysis is also presented in the paper. Fast convergence has shown in the numerical results by using the proposed new circulant preconditioners.

Introduction

In the last two decades, 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, blurred images. One of them, of course, is nonlinearity and the other is the solution of the large system which arises from the discretization after linearization. The main focus of this paper is to handle these two computational difficulties. The most well-known nonlinear variational image deblurring model is the total variation (TV) model.^1–5 It has nice properties such as 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)^6–18-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 the staircase effect in the recovery of digital images. However, the discretization of the associated Euler–Lagrange equations produces a large nonlinear ill-conditioned system which affects the convergence of the numerical algorithms such as Krylov subspace methods (generalized minimal residual (GMRES), etc.). Furthermore, the Jacobian matrix of the MC-based nonlinear system is a block banded matrix with a large bandwidth. The MC-based regularization methods are effective, but due to the high nonlinearity and ill-conditioned system, a robust and fast numerical solution is a crucial issue. To overcome these difficulties, in this paper, we introduce three new circulant symmetric positive definite (SPD) preconditioners. So instead of applying the ordinary GMRES method (without preconditioner), we use preconditioned GMRES (PGMRES) method (with new preconditioners) for the solution of the system. A fast convergence has been shown in the numerical results by using the proposed new preconditioners.

The contributions of the study include the following: (i) our work presents an efficient algorithm for an MC -based image deblurring problem, which combines a fixed point iteration (FPI) with new circulant 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 MC regularization functional. The paper is organized into different sections. The first section includes an introduction while the second section includes a problem description of an image deblurring model. In the third section, we present a nonlinear system of first-order equations and cell discritization. In the fourth section, we introduced our proposed circulant preconditioners for the PGMRES method. The properties and eigenvalues analysis of the proposed preconditioned matrices are discussed in the fourth section. The numerical experiments are also in the fourth section. The conclusions about the proposed PGMRES method are discussed in the last section of the paper.

Problem description

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

$$z=\overrightarrow{K}u+ϵ$$

(1)

where $\epsilon$ is the noise function. The noise can be Gaussian noise, salt and pepper noise, Brownian noise, etc. The $\overrightarrow{K}$ is the blurring operator, which is a Fredholm-integral operator of the first kind:

$$(\overrightarrow{K}u)(x)={\int }_{\mathrm{\Omega }}k(x,y)u(y)\mathrm{d}y,\mathbf{x}\in \mathrm{\Omega }$$

having translation invariance property, that is, $k(x,y)=k(x-y)$ . The problem (1) becomes ill-posed^1,19,4 because $\overrightarrow{K}$ is a compact operator. Let $\mathrm{\Omega }$ be a square in ${\text{R}}^2$ . $u\in \mathrm{\Omega }$ is an image intensity function. The $\mathbf{x}=(x,y)$ defines the position in $\mathrm{\Omega }$ . Let $|\mathbf{x}|=\sqrt{x^2+y^2}$ is an Euclidean norm and $\Vert .\Vert$ is $L_2(\mathrm{\Omega })$ norm. Equation (1) is an inverse problem. The recovering of $u$ from $z$ makes (1) an unstable problem.^1,19,4 To make it stable, one remedy is to use the MC regularization functional^6,14,16–18:

$$J(u)={\int }_{\mathrm{\Omega }}\kappa {(u)}^2\mathbf{d}\mathbf{x}={\int }_{\mathrm{\Omega }}{(▽.\frac{▽u}{|▽u|})}^2\mathbf{d}\mathbf{x}$$

Then problem (1) takes the form, find $u$ that minimize

$$T(u)=\frac{1}{2}\Vert \overrightarrow{K}u-z{\Vert }^2+\frac{\alpha }{2}J(u)$$

(2)

where $\alpha >0$ is a regularization parameter. The well-posedness of problem (2) for a particular case (synthetic image denoising problem) is explained in Zhu and Chan. ¹⁷ Then, the Euler–Lagrange equations of (2) are as follows:

$${\overrightarrow{K}}^{*}(\overrightarrow{K}u-z)+\alpha ▽.[\frac{▽\kappa }{\sqrt{{|▽u|}^2+{\beta }^2}}-\frac{▽\kappa .▽u}{(\sqrt{{|▽u|}^2+{\beta }^2}{)}^3}▽u]=0\text{ in }\mathrm{\Omega }$$

(3)

$$\frac{\mathrm{\partial }u}{\mathrm{\partial }n}=0\text{ in }\mathrm{\partial }\mathrm{\Omega }$$

(4)

$$\kappa (u)=0\text{ in }\mathrm{\partial }\mathrm{\Omega }$$

(5)

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

The MC-based model not only preserves edges but also removes the 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 MC functional. We have treated this difficulty by reducing the nonlinear fourth-order Euler-Lagrange equation into a system of first-order equations.

The first-order nonlinear system

Equation (3) can be expressed as a first-order nonlinear system:

$${\overrightarrow{K}}^{*}\overrightarrow{K}u+\alpha ▽.\overrightarrow{p}-\alpha ▽.\overrightarrow{t}={\overrightarrow{K}}^{*}z$$

(6)

$$-w+▽.\overrightarrow{v}=0$$

(7)

$$\sqrt{{|▽u|}^2+{\beta }^2}\overrightarrow{v}-▽u=0$$

(8)

$$\sqrt{{|▽u|}^2+{\beta }^2}\overrightarrow{p}-▽w=0$$

(9)

$$\sqrt{{|▽u|}^2+{\beta }^2}\overrightarrow{t}-(▽w.\overrightarrow{v})\overrightarrow{v}=0$$

(10)

where

$$\overrightarrow{v}=\frac{▽u}{\sqrt{{|▽u|}^2+{\beta }^2}},w=▽.\overrightarrow{v},\overrightarrow{p}=\frac{▽w}{\sqrt{{|▽u|}^2+{\beta }^2}}$$

and

$$\overrightarrow{t}=\frac{(▽w.\overrightarrow{v})\overrightarrow{v}}{\sqrt{{|▽u|}^2+{\beta }^2}}$$

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

Cell discretization

For an MC-based image deblurring problem, the domain $\mathrm{\Omega }=(0,1)\times (0,1)$ is partitioned by ${\delta }_x\times {\delta }_y$ , ⁴ where

$$\begin{array}{rl}{\delta }_x:0& =x_{1/2}<x_{3/2}<x_{5/2}<\cdots <x_{n_x-1/2}<x_{n_x+1/2}=1\\ {\delta }_y:0& =y_{1/2}<y_{3/2}<y_{5/2}<\cdots <y_{n_x-1/2}<y_{n_x+1/2}=1\end{array}$$

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

$$\begin{array}{rl}{x}_i& =(i-\frac{1}{2})hi=1,2,3,\dots ,n_x\\ y_j& =(j-\frac{1}{2})hj=1,2,3,\dots ,n_x\end{array}$$

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

$$\begin{array}{rl}{x}_{i\pm \frac{1}{2}}& =x_i\pm \frac{h}{2}i=1,2,3,\dots ,n_x\\ y_{j\pm \frac{1}{2}}& =y_j\pm \frac{h}{2}j=1,2,3,\dots ,n_x\end{array}$$

The set

$$e_{ij}=\{(x,y):x\in [x_i-\frac{1}{2},x_i+\frac{1}{2}],y\in [y_j-\frac{1}{2},y_j+\frac{1}{2}]\}$$

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

$$\begin{array}{rl}{\chi }_i(x)& =\{\begin{array}{cc}1& x\in (x_i-\frac{1}{2},x_i+\frac{1}{2})\\ 0& \text{otherwise}\end{array}\\ {\chi }_j(y)& =\{\begin{array}{cc}1& y\in (y_i-\frac{1}{2},y_i+\frac{1}{2})\\ 0& \text{otherwise}\end{array}\end{array}$$

And

$$\begin{array}{rl}{\varphi }_i(x_l+\frac{1}{2})& ={\delta }_{il}\\ {\varphi }_k(y_j+\frac{1}{2})& ={\delta }_{jk}\end{array}$$

An approximation of $u$ and $w$ is as follows:

$$u(x,y)\cong \overline{U}(x,y)=\sum _{i=1}^{n_x}\sum _{j=1}^{n_x}{u}_{ij}{\chi }_i(x){\chi }_j(y)$$

and

$$w(x,y)\cong \overline{W}(x,y)=\sum _{i=1}^{n_x}\sum _{j=1}^{n_x}{w}_{ij}{\chi }_i(x){\chi }_j(y)$$

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

$$z(x,y)\cong \overline{Z}(x,y)=\sum _{i=1}^{n_x}\sum _{j=1}^{n_x}{z}_{ij}{\chi }_i(x){\chi }_j(y)$$

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

$$(Ku)(x_i,y_j)\cong [K_{h}U{]}_{(ij)}$$

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}$ is

$${\overline{V}}^x(x,y)=\sum _{i=1}^{n_x-1}\sum _{j=1}^{n_x}{v}_{ij}^x{\varphi }_i(x){\chi }_j(y)\text{and}{\overline{V}}^y(x,y)=\sum _{i=1}^{n_x-1}\sum _{j=1}^{n_x}{v}_{ij}^y{\varphi }_i(y){\chi }_j(x)$$

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

$${\overline{P}}^x(x,y)=\sum _{i=1}^{n_x-1}\sum _{j=1}^{n_x}{p}_{ij}^x{\varphi }_i(x){\chi }_j(y),{\overline{P}}^y(x,y)=\sum _{i=1}^{n_x-1}\sum _{j=1}^{n_x}{p}_{ij}^y{\varphi }_i(y){\chi }_j(x)$$

and

$${\overline{T}}^x(x,y)=\sum _{i=1}^{n_x-1}\sum _{j=1}^{n_x}{t}_{ij}^x{\varphi }_i(x){\chi }_j(y),{\overline{T}}^y(x,y)=\sum _{i=1}^{n_x-1}\sum _{j=1}^{n_x}{t}_{ij}^y{\varphi }_i(y){\chi }_j(x)$$

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.

The cell-centered finite difference method

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

$$\begin{array}{rll}U& =[{\overline{U}}_{11}\cdots {\overline{U}}_{n_x{n}_x}{]}^t,W=[{\overline{W}}_{11}\cdots {\overline{W}}_{n_x{n}_x}{]}^t,V& =[{\overline{V}}_{11}^x\cdots {\overline{V}}_{n_x-1n_x-1}^x{\overline{V}}_{11}^y\cdots {\overline{V}}_{n_x-1n_x-1}^y{]}^t\\ P& =[{\overline{P}}_{11}^x\cdots {\overline{P}}_{n_x-1n_x-1}^x{\overline{P}}_{11}^y\cdots {\overline{P}}_{n_x-1n_x-1}^y{]}^t\\ \text{and}T& =[{\overline{T}}_{11}^x\cdots {\overline{T}}_{n_x-1n_x-1}^x{\overline{T}}_{11}^y\cdots {\overline{T}}_{n_x-1n_x-1}^y{]}^t\end{array}$$

Now by applying the CCFD method to (6) to (10), we obtain the following system:

$$K^{*}{K}_{h}U-\alpha A_{h}W+\alpha B_h^{*}P-\alpha B_h^{*}T=K_h^{*}Z$$

(11)

$$-I_{h}W+B_h^{*}V=O$$

(12)

$$D_{h}V+B_{h}U=O$$

(13)

$$D_{h}P+B_{h}W=O$$

(14)

$$D_{h}T-C_{h}V=O$$

(15)

where the midpoint quadrature rule is used for the integral term. The matrices $K_h,A_h$ , and $I_h$ are 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:

$$\left[\begin{array}{ccccc}{K}_h^{*}{K}_h& -\alpha A_h& O& \alpha B_h^{*}& -\alpha B_h^{*}\\ O& -I_h& B_h^{*}& O& O\\ B_h& O& D_h& O& O\\ O& B_h& O& D_h& O\\ O& O& -C_h& O& D_h\end{array}\right]\left[\begin{array}{c}U\\ W\\ V\\ P\\ T\end{array}\right]=\left[\begin{array}{c}{K}_h^{*}Z\\ O\\ O\\ O\\ O\end{array}\right]$$

The matrix $K_h$ is block Toeplitz with Toeplitz blocks and $K_h^{*}{K}_h$ is SPD. The matrix $A_h$ is a diagonal matrix having the following structure:

$$A_h=\frac{2}{\beta h}(A_1+A_2)$$

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

$$A_1=\tilde{I}\otimes E\text{and}{A}_2=E\otimes \tilde{I}$$

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

$$E=\left[\begin{array}{cccccc}1& & & & & \\ & 0& & & & \\ & & \ddots & & & \\ & & & & 0& \\ & & & & & 1\end{array}\right]$$

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

$$B_h=\frac{1}{h}\left[\begin{array}{c}{B}_1\\ B_2\end{array}\right]$$

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

$$B_1=F\otimes \tilde{I}\text{and}{B}_2=\tilde{I}\otimes F.$$

$$F=\left[\begin{array}{cccccc}1& -1& & & & \\ & 1& -1& & & \\ & & \ddots & \ddots & & \\ & & & \ddots & -1& \\ & & & & 1& -1\end{array}\right]$$

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

$$C_h=\left[\begin{array}{cc}{C}^x& 0\\ 0& C^y\end{array}\right]$$

is a diagonal matrix and its entries are obtained by the discretization of the expression $(▽w\cdot \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{|▽u{|}^2+{\beta }^2}$ . The matrix $D_h$ has the following structure:

$$D_h=\left[\begin{array}{cc}{D}^x& 0\\ 0& D^y\end{array}\right]$$

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 (11) to (15), then we have the following primal system:

$$(K_h^{*}{K}_h+\alpha L_h(U))U=K_h^{*}Z$$

(16)

where

$$L_h=(B_h^{*}{D}_h^{-1}{B}_h{)}^2+A_h(B_h^{*}{D}_h^{-1}{B}_h)+B_h^{*}{D}_h^{-1}{C}_h{D}_h^{-1}{B}_h$$

(17)

The $L_h$ is a 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 MC-based nonlinear image minimization problems.^{20,6,21,14,4,22,16–18} 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 MC 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 perform quite poorly. To make a robust numerical method for the MC-based nonlinear image deblurring problem, now we present a new preconditioned numerical method.

Numerical implementation

Here we introduce the algorithms to solve the MC-based nonlinear system (16). First, we apply a discrete version of the FPI to (16) to handle the nonlinearity of MC. The approach taken here is called the “lagged diffusivity” scheme. ⁴ Its rate is just linear but in practice, it has a quite rapid convergence. Furthermore, this scheme does not depend on the initial guess to converge globally. This is why globalization is not an issue for this scheme. So by FPI, we have the following linear system:

$$(K_h^{*}{K}_h+\alpha L_h(U^m))U^{m+1}=K_h^{*}Z$$

(18)

Properties

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

The Hessian matrix $K_h^{*}{K}_h+\alpha L_h$ is extremely large for practical applications. When $\alpha$ is small, the Hessian matrix becomes quite ill-conditioned. This happens because the eigenvalues of the blurring operator $\overrightarrow{K}$ cluster as zero. ⁴

The first term $K_h^{*}{K}_h$ in the Hessian matrix is SPD. Although, $K_h^{*}{K}_h$ is full but the blurring operator $\overrightarrow{K}$ has a translation invariant property, which permit the use of fast Fourier transformation (FFT) to evaluate $K_h^{*}{K}_{h}u$ in $O(n\log n)$ operations. ⁴

The second term $L_h$ in the Hessian matrix is sparse but not symmetric. The $L_h$ (16) consists of three terms. The first and the last term in $L_h$ are symmetric positive semidefinite ⁴ but the middle term is not symmetric. Hence, system (18) is not SPD.

Preconditioned matrices

According to the properties of our system (18), mentioned above, the GMRES method is suitable for the solution of system (18). Due to an ill-conditioned system, GMRES can be quite slow to convergence. So we use the PGMRES method. For an effective solution, the preconditioning matrix must be SPD. Here we introduce the following three SPD circulant preconditioned matrices. The first preconditioning matrix $P_1$ is

$$P_1=(\gamma \widetilde{K_h^{*}}\widetilde{K_h}+\alpha I_h{)}^{1/2}(\gamma L_h+I_h)(\gamma \widetilde{K_h^{*}}\widetilde{K_h}+\alpha I_h{)}^{1/2}$$

(19)

where $\gamma$ is the positive parameter. The $\widetilde{K_h}$ is a circulant approximation of matrix $K_h$ and $I_h$ is an identity matrix. The second preconditioning matrix $P_2$

$$P_2=(\gamma \widetilde{K_h^{*}}\widetilde{K_h}+\alpha I_h{)}^{1/2}(\gamma L_h^s+I_h)(\gamma \widetilde{K_h^{*}}\widetilde{K_h}+\alpha I_h{)}^{1/2}$$

(20)

where the matrix

$$L_h^s=(B_h^{*}{D}_h^{-1}{B}_h{)}^2+B_h^{*}{D}_h^{-1}{C}_h{D}_h^{-1}{B}_h$$

(21)

is the SPD part (the first and the third term) of $L_h$ matrix (18). The third preconditioning matrix $P_3$ is

$$P_3=(\gamma \widetilde{K_h^{*}}\widetilde{K_h}+\alpha I_h{)}^{1/2}(\gamma L_h^{\mathrm{TV}}+I_h)(\gamma \widetilde{K_h^{*}}\widetilde{K_h}+\alpha I_h{)}^{1/2}$$

(22)

where the matrix

$$L_h^{\mathrm{TV}}=B_h^{*}{D}_h^{-1}{B}_h,$$

(23)

is the SPD matrix that lies in the first term of $L_h$ matrix (18). In fact, $L_h^{\mathrm{TV}}$ arise by the discretization of the TV-based image deblurring problem. ⁴

The preconditioned matrices $P_1,P_2$ , and $P_3$ are “product preconditioners” based on the analogy of “operator splitting” techniques.^23,4 While applying PGMRES to (18), the inversion of preconditioner matrices ( $P_1,P_2$ , and $P_3$ ) will be required. The inversion of the preconditioning matrices ( $P_1,P_2,$ and $P_3$ ) requires the inversion of the middle terms of (17), (20), and (20), respectively. The middle terms of $P_1,P_2$ , and $P_3$ have the following forms: $\gamma L_h+I_h,\gamma L_h^s+I_h$ , and $\gamma L_h^{\mathrm{TV}}+I_h$ , respectively. Since these middle terms are sparse matrices, so inversion can be done easily. The inversion of the preconditioning matrices ( $P_1,P_2,$ and $P_3$ ) also require the inversion of the matrix $(\gamma \widetilde{K_h^{*}}\widetilde{K_h}+\alpha I_h{)}^{1/2}$ . For this, we need only $O(n\log n)$ floating point operations using FFTs. The detail is given in Vogel and Oman. ⁴ The PGMRES method is summarized in Algorithm 4.1:

Algorithm 4.1.

The PGMRES method

1: On mesh ${\mathrm{\Omega }}_h$ , Initial iteration $U^0,$

2:  for $m=1;max$

$A^m=K_h^{*}{K}_h+\alpha L_h(U^m),$

$b^m=K_h^{*}Z,$

3: Use PGMRES to solve

$A^m{U}^{m+1}=b^m,$

with preconditioned matrix $P=P_1,P_2,P_3$ .

end

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 \mathrm{\infty }$ . So the eigenvalues of $P_1^{-1}\overline{A}$ are as follows:

$${\theta }_i^1=\frac{{\lambda }_i^K+\alpha {\lambda }_i^L}{(\gamma {\lambda }_i^K+\alpha )(\gamma {\lambda }_i^L+1)}$$

(24)

where $\overline{A}=K_h^{*}{K}_h+\alpha L_h$ is the Hessian matrix of system (18). Clearly, if one can choose small $\gamma$ such that $\alpha \le \gamma \le \frac{{\lambda }_i^L-1}{{\lambda }_i^L}$ , then ${\theta }_i^1\to 1$ as $i\to \mathrm{\infty }$ . Now consider the eigenvalues of $P_2^{-1}\overline{A}$ and $P_3^{-1}\overline{A}$ , these are

$${\theta }_i^2=\frac{{\lambda }_i^K+\alpha {\lambda }_i^L}{(\gamma {\lambda }_i^K+\alpha )(\gamma {\lambda }_i^{L^s}+1)}\text{and}{\theta }_i^3=\frac{{\lambda }_i^K+\alpha {\lambda }_i^L}{(\gamma {\lambda }_i^K+\alpha )(\gamma {\lambda }_i^{L^{\mathrm{TV}}}+1)}$$

(25)

respectively. Here, ${\lambda }_i^{L^{\mathrm{TV}}}$ and ${\lambda }_i^{L^{ns}}$ are eigenvalues of $L_h^{\mathrm{TV}}$ and $L_h^{ns}$ , respectively. One can notice that ${\lambda }_i^{L^{\mathrm{TV}}}\le {\lambda }_i^{L^s}\le {\lambda }_i^L$ . Then

$${\theta }_i^2\ge \frac{{\lambda }_i^K+\alpha {\lambda }_i^L}{(\gamma {\lambda }_i^K+\alpha )(\gamma {\lambda }_i^L+1)}\text{and}{\theta }_i^3\ge \frac{{\lambda }_i^K+\alpha {\lambda }_i^L}{(\gamma {\lambda }_i^K+\alpha )(\gamma {\lambda }_i^L+1)}$$

(26)

Clearly, for $\alpha \le \gamma \le \frac{{\lambda }_i^L-1}{{\lambda }_i^L}$

$${\theta }_i^2,{\theta }_i^3\to 1\text{as}i\to \mathrm{\infty }$$

Hence, for $\alpha \le \gamma \le [({\lambda }_i^L-1)/{\lambda }_i^L]\le 1$ , $P_1^{-1}\overline{A}$ , $P_2^{-1}\overline{A}$ and $P_3^{-1}\overline{A}$ have more favorable spectrum as compared to the Hessian matrix $\overline{A}$ , so PGMRES ( $P_1$ GMRES, $P_2$ GMRES and $P_3$ GMRES) converges more rapidly than ordinary GMRES. This can also be shown in the numerical examples.

Now we present four numerical examples for the MC-based image deblurring problem. We have used different values of $n_x$ so the resulting system has $n_x^2$ unknowns. In all examples, the value of the parameters $\alpha$ and $\beta$ is used according to the literature.^24,17,18 For numerical computations, MATLAB software is used to obtain the numerical results on Intel(R) Core(TM) i7-4510U with CPU@2.00 GHz and 2.60 GHz. To measure the quality of the deblurred images, we have used peak signal-to-noise ratio (PSNR), signal-to-noise ratio (SNR), and structured similarity index measure (SSIM).

Example 1

In this example, we have used Lena image. This is a complicated image, because it has a small-scale texture part (hat) and a large-scale cartoon part (face). The different aspects of the Lena image are presented in Figure 1. The size of each subfigure is $512\times 512$ . These are (a) exact image, (b) blurry image, (c) deblurred image by GMRES, (d) deblurred image by $P_1$ GMRES with $\gamma =1e-2$ (e) deblurred image by $P_2$ GMRES with $\gamma =1e-6$ , and (f) deblurred image by $P_3$ GMRES with $\gamma =1e-6$ . The calculation of relative residual at each iteration against different values of the parameter $\gamma$ is presented in Figure 2. The $k{e}_{-}\mathrm{gen}(N,300,4)$ kernel is used for numerical calculations. The parameters $\alpha =1e-8$ and $\beta =1$ . For the stopping criteria of a numerical method we have used tolerance $\mathrm{tol}=1e-7$ .

Figure 1.

Lena image: (a) exact image, (b) blurry image, (c) deblurred image by generalized minimal residual (GMRES), (d) deblurred image by $P_1$ GMRES with $\gamma =1e-2$ , (e) deblurred image by $P_2$ GMRES with $\gamma =1e-6$ , and (f) deblurred image by $P_3$ GMRES with $\gamma =1e-6$ .

Figure 2.

GMRES and PGMRES convergence at an FPI $m=1$ for Example 1. The norm of the residual at each iteration is present on the vertical axis. The blue line represents GMRES iterations. (a) $P_1$ GMRES iteration against different values of $\gamma$ (b) $P_2$ GMRES iteration against different values of $\gamma$ and (c) $P_3$ GMRES iteration against different values of $\gamma$ . GMRES: generalized minimal residual; PGMRES: preconditioned GMRES; FPI: fixed point iteration.

Remarks

One can notice that Figure 1(d) to (f) is having much better quality as compared with Figure 1(c). This means the PGMRES ( $P_1$ GMRES, $P_2$ GMRES and $P_3$ GMRES) method is generating better quality results by using the ordinary GMRES method without preconditioning.

From Figure 2, one can clearly observe the effectiveness of preconditioning. Here, the result was presented using an FPI count $m=1$ . The number of PGMRES iteration is much lesser than as compared to GMRES (without preconditioning) to reach the required accuracy $\mathrm{tol}=1e-7$ . The later FPIs are also having similar results.

From Figure 2, this can also be observed that PGMRES is getting rapid convergence when the value of $\gamma$ is close to the value of $\alpha =1e-8$ for all preconditioners ( $P_1$ , $P_2$ , and $P_3$ ).

From Table 1 and Figure 2, one can observe that the PSNR and SSIM by using the PGMRES method is much better as compared with the ordinary GMRES method. But the PGMRES method is producing better PSNR and SSIM in quite less iterations.

The $P_1$ GMRES method achieves its best PSNR and SSIM with $\gamma =1e-2$ . The $P_2$ GMRES method needs only 11 iterations with $\gamma =1e-6$ to get his best PSNR and SSIM. Similar is the case for $P_3$ GMRES. This also justify the selection of the best $\gamma$ , that is, $\alpha \le \gamma \le 1$ , which we have observed in the eigenvalue analysis (section “Eigvenvalues”). Although, the value of $\gamma$ smaller than $\alpha$ further reduces the number of iterations but it also decrease the deblurred quality (PSNR and SSIM).

It is also observed that although the $P_1$ GMRES method achieves little higher PSNR and SSIM as compared with both $P_2$ GMRES and $P_3$ GMRES but $P_1$ GMRES needs much more iterations than $P_2$ GMRES and $P_3$ GMRES. So, in this example, the performance of preconditioners $P_2$ and $P_3$ is more effective than the preconditioner $P_1$ .

Table 1.

Comparison of GMRES and PGMRES for Example 1.

Blurred	Blurred			Deblurred	Deblurred
PSNR	SSIM	Method	$\gamma$	PSNR	SSIM	Iterations
24.9488	0.7902	GMRES		32.1210	0.9795	${50}^{+}$
24.9488	0.7902	$P_1$ GMRES	$1e-2$	44.5108	0.9631	${50}^{+}$
			$1e-4$	41.4644	0.9825	50
			$1e-6$	42.0135	0.9841	11
			$1e-8$	35.5159	0.9819	5
			$1e-10$	32.1702	0.9796	4
24.9488	0.7902	$P_2$ GMRES	$1e-2$	41.6700	0.9832	${50}^{+}$
			$1e-4$	41.9035	0.9830	50
			$1e-6$	42.0184	0.9841	11
			$1e-8$	35.5160	0.9819	5
			$1e-10$	32.1702	0.9796	4
24.9488	0.7902	$P_3$ GMRES	$1e-2$	41.8196	0.9830	${50}^{+}$
			$1e-4$	41.9048	0.9830	50
			$1e-6$	42.0184	0.9841	11
			$1e-8$	35.5160	0.9819	5
			$1e-10$	32.1702	0.9796	4

GMRES: generalized minimal residual; PGMRES: preconditioned GMRES; PSNR: peak signal-to-noise ratio; SSIM: structured similarity index measure.

Example 2

Here we have used a nontexture Peppers image. The different aspects of the Peppers image are presented in Figure 3. The size of each subfigure is $256\times 256$ . These are (a) exact image, (b) blurry image, (c) deblurred image by GMRES, (d) deblurred image by $P_1$ GMRES with $\gamma =1e-4$ , (e) deblurred image by $P_2$ GMRES with $\gamma =1e-5$ , and (f) deblurred image $P_3$ GMRES with $\gamma =1e-5$ . For numerical calculations, we have used the $k{e}_{-}\mathrm{gen}(N/4,300,3)$ kernel. Here, $\alpha =1e-6$ and $\beta =0.1$ . For comparison, we have used three different values of $n_x$ . These are 64, 128, and 256. For the stopping criteria of numerical methods, we have used tolerance $\mathrm{tol}=1e-7$ .

Figure 3.

Peppers image: (a) exact image, (b) blurry image, (c) deblurred image by generalized minimal residual (GMRES), (d) deblurred image by $P_1$ GMRES with $\gamma =1e-4$ , (e) deblurred image by $P_2$ GMRES with $\gamma =1e-5$ , and (f) deblurred image $P_3$ GMRES with $\gamma =1e-5$ .

Remarks

Figure 3(c) to (f) is almost similar, so all methods are generating same quality results.

From Figure 4, one can clearly notice the effectiveness of preconditioning. For all values of $n_x$ , the number of PGMRES ( $P_1$ GMRES, $P_2$ GMRES, and $P_3$ GMRES) iteration is much lesser than as compared to GMRES to reach the required accuracy $\mathrm{tol}=1e-7$ . The later FPIs are also having similar results.

From Table 2, it is observed that the PSNR and SSIM by the PGMRES method are almost the same as compared with the ordinary GMRES method for all values of $n_x$ . But the PGMRES method is generating this PSNR and SSIM in quite less iterations. To achieve the PSNR and SSIM, the PGMRES ( $P_1$ GMRES, $P_2$ GMRES and $P_3$ GMRES) method needs only 10 iterations for $n_x=64$ . But the GMRES method needs 100 plus iterations to get the same PSNR and SSIM. The same is the case for other values of $n_x$ . Which means that the PGMRES method is faster than the GMRES for nontexture images.

From Table 2 and Figure 4, it is observed that the performance of all the preconditioners ( $P_1,P_2$ , and $P_3$ ) is the same.

Figure 4.

The GMRES and PGMRES convergence at FPI $m=1$ for Example 2. The norm of the residual at each iteration is presenting on vertical axis. The blue line represents GMRES iterations, the orange line represents $P_1$ GMRES iterations, the yellow line represents $P_2$ GMRES iteration, and the purple line represents $P_3$ GMRES iteration. (a) $n_x=64$ .(b) $n_x=128$ and (c) $n_x=256$ . GMRES: generalized minimal residual; PGMRES: preconditioned GMRES; FPI: fixed point iteration.

Table 2.

Comparison of CG, GMRES, and PGMRES for Example 2.

Mesh size	Blurred	Blurred	Method	Deblurred	Deblurred	Iterations
$h$	PSNR	SSIM		PSNR	SSIM
$1/64$	16.7522	0.6459	GMRES	35.1023	0.9905	$100+$
	16.7522	0.6459	$P_1$ GMRES	35.1023	0.9905	10
	16.7522	0.6459	$P_2$ GMRES	35.1024	0.9905	10
	16.7522	0.6459	$P_3$ GMRES	35.1024	0.9905	10
$1/128$	18.3185	0.6775	GMRES	36.1033	0.9840	$200+$
	18.3185	0.6775	$P_1$ GMRES	36.0281	0.9838	25
	18.3185	0.6775	$P_2$ GMRES	36.0292	0.9838	25
	18.3185	0.6775	$P_{3}2$ GMRES	36.0262	0.9838	25
$1/256$	19.5651	0.7013	GMRES	33.2606	0.9616	$200+$
	19.5651	0.7013	$P_1$ GMRES	33.2261	0.9613	40
	19.5651	0.7013	$P_2$ GMRES	33.0346	0.9587	40
	19.5651	0.7013	$P_3$ GMRES	33.0314	0.9587	40

CG: conjugate gradient; GMRES: generalized minimal residual; PGMRES: preconditioned GMRES; PSNR: peak signal-to-noise ratio; SSIM: structured similarity index measure.

Example 3

In this example, we have used a Goldhills image. This is a real and synthetic image. Here we have compared our MC-based algorithm with TV-based algorithm. Since the TV-based model generates an SPD matrix system, so for the solution we have used the conjugate gradient (CG) method. The different aspects of the Goldhills image are presented in Figure 5. The size of each subfigure is $512\times 512$ . These are (a) blurry image, (b) deblurred image by CG, (c) deblurred image by GMRES, (d) deblurred image by $P_1$ GMRES (e) deblurred image by $P_2$ GMRES, and (f) deblurred image $P_3$ GMRES. For numerical calculations, we have used the $k{e}_{-}\mathrm{gen}(N/2,300,2)$ kernel. For the TV-based method, we have used $\alpha =5e-4$ and $\beta =1$ according to. ⁴ For MC-based method we have used $\alpha =1e-4$ and $\beta =1$ and $\gamma =1e-3$ . For the stopping criteria of a numerical methods we have used tolerance $\mathrm{tol}=1e-2$ .

Figure 5.

Goldhills image: (a) blurry image, (b) deblurred image by CG, (c) deblurred image by GMRES, (d) deblurred image by $P_1$ GMRES, (e) deblurred image by $P_2$ GMRES, and (f) deblurred image $P_3$ GMRES. CG: conjugate gradient; GMRES: generalized minimal residual.

Remarks

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

From Figure 6, 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 $\mathrm{tol}=1e-2$ . The later FPIs are also having similar results.

From Table 3, it is observed that the PSNR and SSIM by MC-based PGMRES method are almost the same as compared with the ordinary MC-based GMRES method. But PGMRES method is generating the same PSNR and SSIM in quite less iterations. To achieve the same PSNR and SSIM, the $P_1$ GMRES method needs only 20 iterations, the $P_2$ GMRES method needs only 40 iterations and the $P_3$ GMRES method needs only 10 iterations. But the GMRES method needs 100 plus iterations to get the same PSNR and SSIM. The TV-based CG method is also getting 100 plus iterations to get its PSNR and SSIM. Which means that the 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_3$ is more effective than the preconditioners $P_1$ and $P_2$ .

Figure 6.

The CG, GMRES, and PGMRES convergence at an FPI $m=1$ for Example 3. The norm of the residual at each iteration is presented on the vertical axis. The blue line represents CG iterations, the orange line represents GMRES iterations, the yellow line represents $P_1$ GMRES iteration, the purple line represents $P_2$ GMRES iteration, and the green line represents $P_3$ GMRES iteration. CG: conjugate gradient; GMRES: generalized minimal residual; PGMRES: preconditioned GMRES; FPI: fixed point iteration.

Table 3.

Comparison of CG, GMRES, and PGMRES for Example 3.

	CG	GMRES	$P_1$ GMRES	$P_2$ GMRES	$P_{3}2$ GMRES
Blurred PSNR	23.8381	23.8381	23.8381	23.8381	23.8381
Blurred SSIM	0.3316	0.3316	0.3316	0.3316	0.3316
Deblurred PSNR	34.1331	34.5917	34.1467	34.0028	34.5874
Deblurred SSIM	0.8680	0.8758	0.8639	0.8524	0.8758
Iterations	$100+$	$100+$	20	40	10

CG: conjugate gradient; GMRES: generalized minimal residual; PGMRES: preconditioned GMRES; PSNR: peak signal-to-noise ratio; SSIM: structured similarity index measure.

Example 4

In this example we have applied our MC-based (GMRES and PGMRES) algorithms on a blurry image with a high level of Gaussian noise ( $\mathrm{mean}=0.1,\mathrm{variance}=0.04$ ). We have also presented a comparison of our algorithms with the TV-based (CG) algorithm. For this, we have used a complicated Cameraman image. The different aspects of the Cameraman image are presented in Figure 7. The size of each subfigure is $512\times 512$ . These are (a) noisy image, (b) denoised image by CG, (c) denoised image by GMRES, (d) denoised image by $P_1$ GMRES, (e) denoised image by $P_2$ GMRES, and (f) denoised image $P_3$ GMRES. For numerical calculations, we have used the $k{e}_{-}\mathrm{gen}(N,300,1)$ kernel. For the TV-based method, we have used $\alpha =1e-1$ and $\beta =1$ according to Vogal and Oman. ⁴ For the MC-based method, we have used $\alpha =1e-2$ , $\beta =1$ , and $\gamma =1e-1$ . For the stopping criteria of numerical methods, we have used tolerance $\mathrm{tol}=1e-3$ .

Figure 7.

Cameraman image: (a) noisy image, (b) denoised image by CG, (c) denoised image by GMRES, (d) denoised image by $P_1$ GMRES, (e) denoised image by $P_2$ GMRES, and (f) denoised image $P_3$ GMRES. CG: conjugate gradient; GMRES: generalized minimal residual.

Remarks

From Table 4, it is observed that the SNR and SSIM by MC-based (GMRES and PGMRES) methods are quite higher than the TV-based CG method. The same comparison can be observed from Figure 7(b) to (f). So MC-based (GMRES and PGMRES) methods are generating better-denoised images.

From Figure 8, one can clearly observe the effectiveness of preconditioning. The number of MC-based PGMRES iterations is lesser than as compared to both MC-based GMRES and TV-based CG methods to reach the required accuracy $\mathrm{tol}=1e-3$ . The later FPIs are also having similar results. It is also observed that the TV-based CG method is showing abrupt behavior.

From Table 4, it is observed that the SNR and SSIM obtained by using the MC-based PGMRES method are almost the same as compared with the ordinary MC-based GMRES method. But the PGMRES method is generating the same PSNR and SSIM in quite less iterations. To achieve the same PSNR and SSIM, the $P_1$ GMRES and $P_2$ GMRES methods need only 37 iterations, while the $P_3$ GMRES method needs only 65 iterations. But the GMRES method needs 70 plus iterations to get the same PSNR and SSIM. The TV-based CG method is also getting 70 plus iterations to get its litter lower SNR and SSIM. Which means that the MC-based PGMRES method is faster than the MC-based GMRES and TV-based CG methods for noisy images also. Here, the performance of the preconditioners $P_1$ and $P_3$ is more effective than the preconditioner $P_2$ .

Figure 8.

The CG, GMRES, and PGMRES convergence at FPI $m=1$ for Example 4. The norm of the residual at each iteration is presented on the vertical axis. The blue line represents CG iterations, the orange line represents GMRES iterations, the yellow line represents $P_1$ GMRES iteration, the purple line represents $P_2$ GMRES iteration, and the green line represents $P_3$ GMRES iteration. CG: conjugate gradient; GMRES: generalized minimal residual; PGMRES: preconditioned GMRES; FPI: fixed point iteration.

Table 4.

Comparison of CG, GMRES, and PGMRES for Example 4.

	CG	GMRES	$P_1$ GMRES	$P_2$ GMRES	$P_{3}2$ GMRES
Noisy SNR	6.5714	6.5714	6.5714	6.5714	6.5714
Noisy SSIM	0.5264	0.5264	0.5264	0.5264	0.5264
Denoised PSNR	15.4995	17.6802	17.6840	17.6985	17.6800
Denoised SSIM	0.7957	0.8331	0.8247	0.8118	0.8333
Iterations	$70+$	$70+$	37	65	37

CG: conjugate gradient; GMRES: generalized minimal residual; PGMRES: preconditioned GMRES; SNR: signal-to-noise ratio; PSNR: peak signal-to-noise ratio; SSIM: structured similarity index measure.

Conclusion

A numerical algorithm (PGMRES) is presented to solve the primal form of MC-based nonlinear image deblurring problem. Three new circulant preconditioned matrices ( $P_1,P_2,$ and $P_3$ ) are introduced. Four examples are tested by PGMRES using our new preconditioned matrices. Different kinds of images (complicated, real, synthetic, and nontexture) are tested by PGMRES using our new circulant preconditioned matrices. In Examples 3 and 4, we have compared our MC-based algorithm (PGMRES) with TV-based algorithm (CG) and observed the effectiveness of our algorithm on real and nontexture images. The convergence rates for a solution to the linear system and the norm of the residual at each iteration are presented in each example. Numerical experiments have shown the rapid convergence of the PGMRES method using new circulant preconditioners. Which means that the PGMRES method is robust, faster, and effective for an MC-based nonlinear image deblurring problem.