A New Extrapolation Economy Cascadic Multigrid Method for Image Restoration Problems

Zhaoteng Chu, Ziqi Yan, Chenliang Li

School of Mathematics and Computational Science, Guilin University of Electronics Technology, Guilin, China.

DOI: 10.4236/ajcm.2023.132016                

Abstract

In this paper, a new extrapolation economy cascadic multigrid method is proposed to solve the image restoration model. The new method combines the new extrapolation formula and quadratic interpolation to design a nonlinear prolongation operator, which provides more accurate initial values for the fine grid level. An edge preserving denoising operator is constructed to remove noise and preserve image edges. The local smoothing operator reduces the influence of staircase effect. The experiment results show that the new method not only improves the computational efficiency but also ensures good recovery quality.

Keywords

Extrapolation Economy Cascadic Multigrid Method, New Extrapolation Formula, Edge Preserving Denoising Operator, Local Smoothing Operator

1. Introduction

In the process of image formation, transmission, recording and processing, it is easy to be affected by equipment, environment, human factors and so on. It is very important to recover the noisy and blurred image. Therefore, it is very meaningful to construct a fast and effective image restoration algorithm, which has become one of the research hotspots in the field of image processing.

Consider the following two-dimensional image degradation model,

$g\left(x,y\right)=h\left(x,y\right)\ast u\left(x,y\right)+\eta \left(x,y\right),$ (1)

where $h\left(x,y\right)$ is an impulse response function, which is called the point spread function in the spatial domain, $\ast$ represents the two-dimensional convolution operation, $u\left(x,y\right)$ represents the input image to be processed, $\eta \left(x,y\right)$ represents the additive noise, and $g\left(x,y\right)$ represents the degraded image.

Ignoring the noise $\eta \left(x,y\right)$ , Equation (1) can be expressed as

$g\left(x,y\right)=h\left(x,y\right)\ast u\left(x,y\right)$ . (2)

In practical problems, it is usually necessary to restore the blur- and noise-contaminated image. Therefore, we usually consider the following problem,

$h\left(x,y\right)\ast u\left(x,y\right)=g^{\delta }\left(x,y\right),$ (3)

where $g^{\delta }\left(x,y\right)$ represents the blur- and noise-contaminated image.

Equation (3) is an ill-posed problem. In the case of direct solution, a small disturbance of the right-hand side may lead to great changes, so the recovered image is usually of low quality. The Tikhonov regularization method can eliminate this instability by solving the following model

$\underset{u}{\min }\frac{1}{2}{\displaystyle {\int }_{-\infty }^{+\infty }{\displaystyle {\int }_{-\infty }^{+\infty }{\left(h\ast u-g^{\delta }\right)}^2\text{d}x\text{d}y}}+\alpha \psi \left(u\right),$ (4)

where the first term is the fidelity term, $\psi \left(u\right)$ is the regularization term, also known as constraints, and α is the regularization parameter, which is used to balance the fidelity term and regularization term. The discrete regularization term is denoted by $D\left(u\right)$ , and the selection of the regularization operator $D\left(u\right)$ is crucial for the quality of image restoration. The most representative one is the nonlinear anisotropic diffusion equation proposed by Perona and Malik in [1] , and the regularization operator is expressed as follows

$D\left(u\right)=div\left[c\left(\left|\nabla u\right|\right)\nabla u\right],$ (5)

the diffusion coefficients are chosen as follows,

$c\left(\left|\nabla u\right|\right)=\frac{1}{1+{\left(\frac{\left|\nabla u\right|}{k}\right)}^2},$ (6)

where k is the threshold.

For nonlinear Equation (4), [2] and [3] propose the explicit and semi-implicit discretization schemes, which have good preserving effect on the edges of the image and can provide high quality restoration of images. However, whether it is the explicit or semi-implicit discretization schemes its computation cost is very large. In addition, the selection of the regularization parameters and the time intervals are not easy.

Equation (2) is discretized to obtain the following linear system

$g=Hu,H\in {ℝ}^{n\times n},u,g\in {ℝ}^n,$ (7)

where H represents the fuzzy operator, which is usually a real symmetric block Teoplitiz matrix, u is the original image, and g represents the blur and noise-free image.

The blur- and noise-contaminated image g is obtained by adding random normally distributed noise e with mean 0 and standard deviation 1, such as

$g^{\delta }=g+e.$ (8)

Thus, Equation (3) is discretized to obtain the following linear system

$Hu=g^{\delta },H\in {ℝ}^{n\times n},u,g^{\delta }\in {ℝ}^n.$ (9)

we set the solution is $u^{\delta }$ .

For the linear Equation (9), the matrix H is usually a large-scale singular matrix. The system of linear equations is usually an ill-posed problem. Although we can solve Equation (9) by some iterative method, due to the loss of high-frequency information, these iterative methods may produce ringing effect and cannot accurately recover edges, which seriously affects the quality of image restoration.

In order to solve the problems lots of effective methods have been proposed.

In [4] , a standard iterative method is proposed to compute several approximate solutions and then extract a new candidate solution from the linear subspace expanded by these approximate solutions. This method can be used for large-scale problems. In [5] , LSQR algorithm and RRGMRES algorithm are proposed to solve the large-scale linear ill-posed problems, and corresponding stopping criteria are proposed, which can obtain meaningful approximate solutions. In [6] , the regularization properties of the global GMRES method are analyzed, and a regularized global GMRES method is proposed for solving linear discrete ill-posed problems with multiple right-hand sides contaminated by errors. In [7] , NL-means algorithm is proposed to correct the ladder effect.

Multigrid method is one of the most effective methods for solving elliptic boundary value problems. The cascadic multigrid method is proposed in [8] . The advantage of this method lies in its simplicity and efficiency. In [9] presents an economical cascadic multigrid method. Compared with the conventional cascadic multigrid method, the proposed method can get the same efficient solutions with fewer iterations. [10] [11] [12] have proposed some extrapolation cascadic multigrid methods. The methods use the new extrapolation with quadratic interpolation, which can provide a better initial value for the fine grid level.

Based on the advantages of these methods, many scholars have applied them to image restoration. [13] [14] have proposed multigrid method for image restoration. The step effect is reduced by introducing wavelet soft threshold denoising and smoothing. In [15] , a cascadic multilevel method is proposed to solve a linear discrete ill-posed problem with error contamination on the right side. Numerical experiments show that the algorithm is effective in image restoration. In [16] , a class of nonlinear edge-preserving continuation operators are proposed to be used as prolongation operators for cascadic multigrid methods. This method can reduce the ringing effect and restore good image edges. In [17] , the selection of fuzzy operator and the smooth method are further studied, and it is shown that the cascadic multilevel method has good efficiency for image restoration.

Combining the advantages of the new extrapolation method and the cascadic multigrid method, a new extrapolation economy cascadic multigrid method is proposed for image restoration in this paper. The experimental results show that this method not only improves computational efficiency but also ensures good recovery quality. The new method has the following advantages:

1) The new extrapolation formula and the quadratic interpolation are combined to provide more accurate values for the fine mesh level.

2) Fewer iterations are costed on each grid level without affecting the accuracy.

3) The edge preserving denoising operator can preserve the edge of image while denoising.

4) The local smoothing operator reduces the ladder effect and ensures good recovery quality.

The structure of this paper is as follows. In Section 2, introduces the basic framework of the new extrapolation economy cascadic multigrid. In Section 3, the linear prolongation operator, the new extrapolation formula prolongation operator combined with quadratic interpolation and the smoothing operator generated by the weighted least squares approximation is introduced. In Section 4, the edge preserving denoising operator is introduced. In Section 5, some numerical results are given to illustrate the effectiveness of the new method. In Section 6, a summary is presented.

2. New Extrapolation Economy Cascadic Multigrid Method

Selecting a Template

Let $V={\left[v^{\left(1\right)},v^{\left(2\right)},\cdots ,v^{\left(m\right)}\right]}^{\text{T}}$ , we define the following norm

$\Vert V\Vert ={\left(\frac{1}{m}{\displaystyle {\sum }_{i=1}^m{\left(v^{\left(i\right)}\right)}^2}\right)}^{1/2},$ (10)

For the Equation (9), we can solve it by using the conjugate gradient method (CG) or the minimal residual method (MR). Let the initial iterative value $u_0^l=0$ , here l represents the grid level. We use the following stopping rule to terminate the iteration.

Let $\delta =\Vert e\Vert$ , for a given $\delta \ge 0$ , and $c>1$ is a constant and independent of $\delta$ , when $\Vert g^{\delta }-H{u}_k^l\Vert \le c\delta$ , stop the iteration.

By using CG or MR to solve the linear model directly, although it cannot effectively restore the edge which seriously affects the quality of image restoration, it is simple and the time is short. In view of these two advantages, we can use CG or MR as the smoothing operator in the new extrapolation economic cascadic multigrid method.

For the sake of discussion, we set u as an $n^2\times 1$ column vector and $N=n^2$ , $G_1\subset \cdots \subset G_l\subset \cdots \subset G_L$ is a nested subspace. The subscript l indicates the number of the grid level, $G_1$ denotes the coarsest grid level, and $G_L$ indicates the finest grid level. Accordingly, we have $N_1<\cdots <N_l<\cdots <N_L=N$ . When N is odd, $N_{l-1}={\left(\sqrt{N_l}+1\right)}^2/4$ , when N is even, $N_{l-1}=N_l/4$ .

In the cascadic multigrid algorithm, the following linear models need to be solved iteratively at each grid level

$H_l{u}_l=g_l^{\delta },$ (11)

where $H_l$ is the representation of the fuzzy operator H in the $G_l$ .

Set ${\left(g_l^{\delta }\right)}^{\left(k\right)}$ be the k^th element of $g_l^{\delta }$ . Let

${\left(g_l^{\delta }\right)}^{\left(k\right)}={\left(g_{l+1}^{\delta }\right)}^{\left(2k-1\right)},1\le k\le N_l$ (12)

We define the restriction operator $R_l$ , which satisfies

$g_l^{\delta }=R_l{g}_{l+1}^{\delta }.$ (13)

When $l=L$ , $g_l^{\delta }=g^{\delta }$ .

$H_l$ is generated in the following way,

$H_l=R_l{H}_{l+1}{R}_l^{*}.$ (14)

where $R_l^{*}$ is the adjoint of $R_l$ , when $l=L$ , $H_l=H$ .

In [16] , a cascadic multiresolution method is proposed. Numerical examples show that the method has good performance in image restoration. It is explained that for Equation (9), lots of iterations will affect the recovery effect. Therefore, the number of iterations is particularly critical, and the cascadic multigrid method proposed by [9] is of reference significance. We combine the advantages of the algorithms of [9] and [16] propose the following improved economic cascadic multigrid method.

Algorithm 1: Improved Economical Cascadic Multigrid Method for Image Restoration (IECMG)

Step 1: Set $u_1^0=0$ , $u_1^{\delta }=C_1^{m_1}{u}_1^0$ .

Step 2: When $l=2,3,\cdots ,L$ , do:

1) Prolongation: ${\tilde{u}}_l^0=L_l{u}_{l-1}^{\delta }$ or ${\tilde{u}}_l^0=P_l{u}_{l-1}^{\delta }$ .

2) Local smoothing: $u_l^S=S_l{\tilde{u}}_l^0$ .

3) Preserving edge and removing noise: $u_l^{DS}=D_l{S}_l{\tilde{u}}_l^0$ .

4) Smoothing: $u_l^0=u_l^{DS}$ , $u_l^{\delta }=C_l^{m_l}{u}_l^0$ .

Step 3: $u^{\delta }=u_L^{\delta }$ .

In [9] [10] [11] [12] , the economic cascadic multigrid method and the new extrapolation cascadic multigrid method are proposed respectively. The new extrapolation interpolation method can provide better initial values for the next level and has better convergence. Therefore, we propose the following new extrapolation economy cascadic multigrid method.

Algorithm 2: New Extrapolation Economy Cascadic Multigrid method for Image Restoration (EECMG)

Step 1: Set $u_l^0=0$ , $u_l^{\delta }=C_l^{m_l}{u}_l^0$ , $l=1,2$ .

Step 2: When $l=2,3,\cdots ,L$ , do:

1) New extrapolation: ${\tilde{u}}_l^0={\Pi }_l\left(u_{l-1}^{\delta },u_l^{\delta }\right)$ .

2) Quadratic interpolation prolongation: ${\tilde{u}}_{l+1}^0=P_l{\tilde{u}}_l^0$ .

3) Local smoothing: $u_{l+1}^S=S_{l+1}{\tilde{u}}_{l+1}^0$ .

4) Preserving edge and removing noise: $u_{l+1}^{DS}=D_{l+1}{S}_{l+1}{\tilde{u}}_{l+1}^0$ .

5) Smoothing: Let $u_{l+1}^0=u_{l+1}^{DS}$ , $u_{l+1}=C_{l+1}^{m_{l+1}}{u}_{l+1}^0$

Step 3: $u^{\delta }=u_L^{\delta }$

The number of iterations $m_l$ on each level is selected as follows [9] .

1) If $l>L_0$ , then $m_l=\left[m_0{\left(L-L_0\right)}^2{\beta }^{L-l}\right]$ .

2) If $l\le L_0$ , then $m_l=\left[m_{\ast }^{1/2}\left(L-\left(2-{\epsilon }_0\right)l\right)h_l^{-2}\right]$

Where $l=1,2,\cdots ,L$ , $L_0=\frac{1}{2}L$ , $m_0=1$ , $\beta =4$ , $m_{\text{*}}=1$ , ${\epsilon }_0=\frac{1}{2}$ , $h_l={\left(\frac{1}{2}\right)}^{l-1}$ , $\left[t\right]$ represents the smallest positive integer not less than t.

Note: In the above algorithms, $C_l^{m_l}$ denotes smoothing $m_l$ times by using CG or MR. In Sections 3 and Sections 4, we will introduce the linear prolongation operator $L_l$ , nonlinear prolongation operator $P_l$ , new extrapolation operator ${\Pi }_l$ , local smoothing operator $S_l$ and edge preserving denoising operator $D_l$ in Algorithm 1 and Algorithm 2 in detail.

3. Prolongation Operator and Local Smoothing Operator

3.1. Linear Prolongation Operator

Let $u_{\left(i,j\right)}^l$ denote the pixel value in $\left(i,j\right)$ on the level l. The linear operator $L_l$ can be defined as follows,

$\begin{array}{l}{u}_{\left(i,j\right)}^l=u_{\left(\left(i+1\right)/2,\left(j+1\right)/2\right)}^{l-1},\text{when}i,j\text{arebothodd};\\ u_{\left(i,j\right)}^l=\frac{1}{2}\left(u_{\left(\left(i+1\right)/2,j/2\right)}^{l-1}+u_{\left(\left(i+1\right)/2,j/2+1\right)}^{l-1}\right),\text{when}i\text{is}\text{odd},j\text{is}\text{even};\\ u_{\left(i,j\right)}^l=\frac{1}{2}\left(u_{\left(i/2,\left(j+1\right)/2\right)}^{l-1}+u_{\left(i/2+1,\left(j+1\right)/2\right)}^{l-1}\right),\text{when}i\text{iseven},j\text{is}\text{odd};\\ u_{\left(i,j\right)}^l=\frac{1}{2}\left(u_{\left(i/2,j/2\right)}^{l-1}+u_{\left(i/2+1,j/2+1\right)}^{l-1}\right),\text{when}i,j\text{arebotheven}.\end{array}$ (15)

Because the accuracy of the image on the fine grid level by piecewise linear interpolation is not high, a nonlinear prolongation operator, combined a new extrapolation formula with quadratic interpolation, is constructed to solve this problem.

3.2. New Extrapolation Quadratic Interpolation Prolongation Operator

Let us take the one-dimensional triplet grid $G_l\left(l=1,2,3\right)$ as an example, the pixel is denoted by $p_i^l$ , and the corresponding pixel value is denoted by $u_i^l$ . Let $G_1$ be the coarsest grid level and contains pixels $I_1=\left(p_i^1,p_{i+1}^1\right)$ . The pixels contained in $G_2$ and $G_3$ are represented as follows,

$I_2=\left(p_i^2,p_{\left(i+1\right)/2}^2,p_{i+1}^2\right);I_3=\left(p_i^3,p_{\left(i+1\right)/4}^3,p_{\left(i+1\right)/2}^3,p_{\left(i+3\right)/4}^3,p_{i+1}^3\right),$ (16)

the corresponding pixel value are expressed as,

$\begin{array}{l}{G}_1=\left(u_i^1,u_{i+1}^1\right);G_2=\left(u_i^2,u_{\left(i+1\right)/2}^2,u_{i+1}^2\right);\\ G_3=\left(u_i^3,u_{\left(i+1\right)/4}^3,u_{\left(i+1\right)/2}^3,u_{\left(i+3\right)/4}^3,u_{i+1}^3\right).\end{array}$ (17)

We combine the new extrapolation formula with quadratic interpolation, which is called the new extrapolation quadratic interpolation method. This can provide better initial values on the fine grid level $G_3$ . The details are as follows (see Figure 1).

1) New extrapolation $u_3={\Pi }_2\left(u_1,u_2\right)$ ,

$u_i^3=\frac{1}{4}\left(5u_i^2-u_i^1\right);$

$u_{i+1}^3=\frac{1}{4}\left(5u_{i+1}^2-u_{i+1}^1\right);$

$u_{i+\frac{1}{2}}^3=u_{i+\frac{1}{2}}^2+\frac{1}{8}\left[\left(u_i^2-u_i^1\right)+\left(u_{i+1}^2-u_{i+1}^1\right)\right].$

2) Quadratic interpolation $u_3=P_3\left(u_1,u_2,u_3\right)$ ,

$\begin{array}{c}{u}_{i+\frac{1}{4}}^3=\frac{3}{8}{u}_i^3+\frac{3}{4}{u}_{i+\frac{1}{2}}^3-\frac{1}{8}{u}_{i+1}^3\\ =\frac{1}{16}\left[\left(9u_i^2+12u_{i+\frac{1}{2}}^2-u_{i+1}^2\right)-\left(3u_i^1+u_{i+1}^1\right)\right];\end{array}$

$\begin{array}{c}{u}_{i+3/4}^3=-\frac{1}{8}{u}_i^3+\frac{3}{4}{u}_{i+1/2}^3+\frac{3}{8}{u}_{i+1}^3\\ =\frac{1}{16}\left[\left(9u_{i+1}^2+12u_{i+1/2}^2-u_i^1\right)-\left(3u_{i+1}^1+u_i^1\right)\right]\text{.}\end{array}$

Next, we apply this idea to the two-dimension problems. As shown in Figure 2, the hollow circles represent the pixels on the $G_{l-2}$ , the solid black dots represent the pixel values on the $G_{l-1}$ grid level, and the rectangular boxes represent the pixel values on the $G_l$ grid level. We use the new extrapolation formula to get the image value at the black node and the quadratic interpolation to obtian the image value at the rectangular box.

Figure 1. New extrapolation quadratic interpolation method for the 1-dimensional.

Figure 2. New extrapolation quadratic interpolation method for the 2-dimensional.

1) New extrapolation $u_l={\Pi }_l\left(u_{l-1},u_{l-2}\right)$ ,

$u_{\left(i,j\right)}^l=\frac{1}{4}\left(5u_{\left(i,j\right)}^{l-1}-u_{\left(i,j\right)}^{l-2}\right);$

$u_{\left(i,j+1\right)}^l=\frac{1}{4}\left(5u_{\left(i,j+1\right)}^{l-1}-u_{\left(i,j+1\right)}^{l-2}\right);$

$u_{\left(i,j+\frac{1}{2}\right)}^l=u_{\left(i,j+\frac{1}{2}\right)}^{l-1}+\frac{1}{8}\left[\left(u_{\left(i,j\right)}^{l-1}-u_{\left(i,j\right)}^{l-2}\right)+\left(u_{\left(i,j+1\right)}^{l-1}-u_{\left(i,j+1\right)}^{l-2}\right)\right];$

$u_{\left(i+1/2,j\right)}^l=u_{\left(i+1/2,j\right)}^{l-1}+\frac{1}{8}\left[\left(u_{\left(i,j\right)}^{l-1}-u_{\left(i,j\right)}^{l-2}\right)+\left(u_{\left(i+1,j\right)}^{l-1}-u_{\left(i+1,j\right)}^{l-2}\right)\right].$

2) Quadratic interpolation $u_l=P_l\left(u_l,u_{l-1},u_{l-2}\right)$ ,

$\begin{array}{r}\hfill u_{\left(i,j+\frac{1}{4}\right)}^l=\frac{1}{16}\left[\left(9u_{\left(i,j\right)}^{l-1}+12u_{\left(i,j+\frac{1}{2}\right)}^{l-1}-u_{\left(i,j+1\right)}^{l-1}\right)-\left(3u_{\left(i,j\right)}^{l-2}+u_{\left(i,j+1\right)}^{l-2}\right)\right];\end{array}$

$u_{\left(i,j+\frac{3}{4}\right)}^l=\frac{1}{16}\left[\left(9u_{\left(i,j+1\right)}^{l-1}+12u_{\left(i,j+\frac{1}{2}\right)}^{l-1}-u_{\left(i,j\right)}^{l-1}\right)-\left(3u_{\left(i,j+1\right)}^{l-2}+u_{\left(i,j\right)}^{l-2}\right)\right];$

$u_{\left(i+\frac{1}{4},j\right)}^l=\frac{1}{16}\left[\left(9u_{\left(i,j\right)}^{l-1}+12u_{\left(i+\frac{1}{2},j\right)}^{l-1}-u_{\left(i+1,j\right)}^{l-1}\right)-\left(3u_{\left(i,j\right)}^{l-2}+u_{\left(i+1,j\right)}^{l-2}\right)\right];$

$u_{\left(i+\frac{3}{4},j\right)}^l=\frac{1}{16}\left[\left(9u_{\left(i+1,j\right)}^{l-1}+12u_{\left(i+\frac{1}{2},j\right)}^{l-1}-u_{\left(i,j\right)}^{l-1}\right)-\left(3u_{\left(i+1,j\right)}^{l-2}+u_{\left(i,j\right)}^{l-2}\right)\right];$

$u_{\left(i+1/4,j+1/4\right)}^l=\frac{3}{8}{u}_{\left(i+1/4,j\right)}^l+\frac{3}{4}{u}_{\left(i+1/4,j+1/2\right)}^l-\frac{1}{8}{u}_{\left(i+1/4,j+1\right)}^l;$

$u_{\left(i+1/4,j+3/4\right)}^l=-\frac{1}{8}{u}_{\left(i+1/4,j\right)}^l+\frac{3}{4}{u}_{\left(i+1/4,j+1/2\right)}^l+\frac{3}{8}{u}_{\left(i+1/4,j+1\right)}^l\text{.}$

3.3. Local Smoothing Operator

In order to get a smoother image, the local weighted least squares approximation is used to smooth the image in [2] . As shown in Figure 3, the pixel value $u_{\left(i,j\right)}^l$ is related to the image value of the surrounding eight points. The process is represented by the operator $S_l$ as follows.

First, we introduce the weighting function,

${\omega }_{\left(i,j\right)}^l\left(s,t\right)=\exp \left(-{\left(u_{\left(i,j\right)}^l-u_{\left(i+s,j+t\right)}^l\right)}^2\right),s,t\in \left\{0,\pm 1\right\},$ (18)

then we solve the locally weighted least squares problem,

$\underset{a_0,a_1,a_2}{\min }{\displaystyle {\sum }_{s,t\in \left\{0,\pm 1\right\}}{\left(u_{\left(i+s,j+t\right)}^l-\left(a_0+a_{1}s+a_{2}t\right)\right)}^2{\omega }_{\left(i,j\right)}^l}\left(s,t\right).$ (19)

Set ${\hat{a}}_0,{\hat{a}}_1,{\hat{a}}_2$ represent the solutions of the above problem, then the updated pixel value is $u_{\left(i,j\right)}^l={\hat{a}}_0$ .

4. Edge Preserving Denoising Operator

The isotropic diffusion equation diffuses equally along the tangential direction and the normal direction of the image edge, so it cannot protect the edge, nor can it preserve the original fine structure of the image well, so that the image becomes blurred. In order to deal with these shortcomings, an anisotropic diffusion equation is proposed. One is the nonlinear anisotropic diffusion equation proposed by Perona and Malik, abbreviated as P-M equation [1] , which is as follows,

$\{\begin{array}{l}\frac{\partial u}{\partial t}=\text{div}\left[g\left(\left|\nabla u\right|\right)\nabla u\right],\hfill \\ {u\left(x,y,t\right)|}_{t=0}=u_0\left(x,y\right),\hfill \end{array}$ (20)

where the $u\left(x,y,t\right)$ is the image of the t moment, $g\left(\left|\nabla u\right|\right)\in \left[0,1\right]$ is the diffusion coefficient. The more classical diffusion coefficient is as follows,

$g\left(\left|\nabla u\right|\right)=\frac{1}{1+{\left(\frac{\left|\nabla u\right|}{k}\right)}^2}$ (21)

Figure 3. Local smoothing method.

Set $c\left(x,y\right)=g\left(\left|\nabla u\right|\right)$ , discretization of (20) is as follows,

$\begin{array}{c}\frac{\partial u}{\partial t}=\text{div}\left[g\left(\left|\nabla u\right|\right)\nabla u\right]=\nabla \left[c\left(x,y\right)\nabla u\right]={\left(c\left(x,y\right)u_y\right)}_x+{\left(c\left(x,y\right)u_y\right)}_y\\ =\frac{\partial }{\partial x}\left[c\left(x,y\right)\frac{u\left(x+\frac{\Delta x}{2},y\right)-u\left(x-\frac{\Delta x}{2},y\right)}{\Delta x}\right]\\ +\frac{\partial }{\partial y}\left[c\left(x,y\right)\frac{u\left(x,\frac{\Delta x}{2}+y\right)-u\left(x,y+\frac{\Delta x}{2}\right)}{\Delta y}\right]\end{array}$

$\begin{array}{l}=\frac{1}{\Delta x}[c\left(x+\frac{\Delta x}{2},y\right)\frac{1}{\Delta x}\left(u\left(x+\Delta x,y\right)-u\left(x,y\right)\right)\\ -c\left(x-\frac{\Delta x}{2},y\right)\frac{1}{\Delta x}\left(u\left(x,y\right)-u\left(x-\Delta x,y\right)\right)]\\ +\frac{1}{\Delta y}[c\left(x,\frac{\Delta x}{2}+y\right)\frac{1}{\Delta y}\left(u\left(x,y+\Delta y\right)-u\left(x,y\right)\right)\\ -c\left(x,y-\frac{\Delta x}{2}\right)\frac{1}{\Delta y}\left(u\left(x,y\right)-u\left(x,y-\Delta y\right)\right)],\end{array}$ (22)

Taking $\Delta x=1,\Delta y=1$ , we have,

$\begin{array}{c}\frac{u_{\left(i,j\right)}^{n+1}-u_{\left(i,j\right)}^n}{\tau }=\frac{c_{\left(i+1,j\right)}^n+c_{\left(i,j\right)}^n}{2}\left(u_{\left(i+1,j\right)}^n-u_{\left(i,j\right)}^n\right)+\frac{c_{\left(i-1,j\right)}^n+c_{\left(i,j\right)}^n}{2}\left(u_{\left(i-1,j\right)}^n-u_{\left(i,j\right)}^n\right)\\ +\frac{c_{\left(i,j+1\right)}^n+c_{\left(i,j\right)}^n}{2}\left(u_{\left(i,j+1\right)}^n-u_{\left(i,j\right)}^n\right)+\frac{c_{\left(i,j-1\right)}^n+c_{\left(i,j\right)}^n}{2}\left(u_{\left(i,j-1\right)}^n-u_{\left(i,j\right)}^n\right)\\ =\underset{\left(k,l\right)\in Q\left(i,j\right)}{{\displaystyle \sum }}\frac{c_{\left(k,l\right)}^n+c_{\left(i,j\right)}^n}{2}\left(u_{\left(k,l\right)}^n-u_{\left(i,j\right)}^n\right),\end{array}$ (23)

where $Q\left(i,j\right)$ denotes the set of four neighboring points centered on $\left(i,j\right)$ . So according to the above equation we have,

$\begin{array}{c}{u}_{\left(i,j\right)}^{n+1}=u_{\left(i,j\right)}^n+\tau [\frac{c_{\left(i+1,j\right)}^n+c_{\left(i,j\right)}^n}{2}{u}_{\left(i+1,j\right)}^n+\frac{c_{\left(i-1,j\right)}^n+c_{\left(i,j\right)}^n}{2}{u}_{\left(i-1,j\right)}^n+\frac{c_{\left(i,j+1\right)}^n+c_{\left(i,j\right)}^n}{2}{u}_{\left(i,j+1\right)}^n\\ +\frac{c_{\left(i,j-1\right)}^n+c_{\left(i,j\right)}^n}{2}{u}_{\left(i,j-1\right)}^n-\frac{c_{\left(i+1,j\right)}^n+c_{\left(i-1,j\right)}^n+c_{\left(i,j+1\right)}^n+c_{\left(i,j-1\right)}^n+4c_{\left(i,j\right)}^n}{2}{u}_{\left(i,j\right)}^n]\text{.}\end{array}$ (24)

Equation (24) can be expressed in matrix form as,

$u^{n+1}=u^n+\tau A\left(u^n\right)u^n,$ (25)

where the elements of $A\left(u^n\right)$ are as follows,

$a_{\left(i,j\right),\left(s,t\right)}=\{\begin{array}{l}c\frac{c_{\left(i,j\right)}^n+c_{\left(s,t\right)}^n}{2},\left(s,t\right)\in Q\left(i,j\right)\\ -{\displaystyle {\sum }_{\left(m,n\right)\in Q\left(i,j\right)}\frac{c_{\left(i,j\right)}^n+c_{\left(m,n\right)}^n}{2}},\left(s,t\right)=\left(i,j\right)\text{.}\\ 0,\text{other}\end{array}$ (26)

5. Numerical Examples

We use numerical examples to illustrate the effectiveness of the new extrapolation economy cascadic multigrid method for restoring blur- and noise-contaminated image. The elements in the fuzzy operator H are defined by the following function.

$h_{\left(i,j\right)}=\{\begin{array}{ll}\frac{1}{\sigma \sqrt{2\pi }}\exp \left(-\frac{{\left(i-j\right)}^2}{2{\sigma }^2}\right),\hfill & \text{if}\left|j-k\right|\le \text{band}\hfill \\ 0,\hfill & \text{other}\hfill \end{array}$

$h={\left[h_{\left(i,j\right)}\right]}_{i,j=1,\cdots ,n},$

$H=h\otimes h\text{.}$

where $\otimes$ denotes the Kronecker product. σ represents the variance of the Gaussian point diffusion function. The “band” represents the half bandwidth of h. The degree of blurring is denoted by $b_i$ , and the following cases are considered in this paper,

$b_1:\{\begin{array}{l}\sigma =1\hfill \\ \text{band}=7\hfill \end{array},b_2:\{\begin{array}{l}\sigma =2\hfill \\ \text{band}=9\hfill \end{array},b_3:\{\begin{array}{l}\sigma =3\hfill \\ \text{band}=11\hfill \end{array},b_4:\{\begin{array}{l}\sigma =4\hfill \\ \text{band}=13\hfill \end{array}\text{.}$

Defining the noise levels as follows,

$\frac{\Vert u-u^{\delta }\Vert }{\Vert u\Vert }=\frac{\Vert e\Vert }{\Vert u\Vert }\approx v,$

the noise level v is set to the following values,

$v_1=5\times {10}^{-2},v_2=1\times {10}^{-1},v_3=5\times {10}^{-1},v_4=8\times {10}^{-1}\text{.}$

In this experiment, the blur and noise level of the image is denoted as $b_i{v}_i$ ( $i=1,2,3,4$ ). For example, $b_1{v}_1$ denotes that $\sigma =1$ , $\text{band}=7$ , $v_1=5\times {10}^{-2}$ .

The image restoration effect is measured by the following peak signal ratio,

$\text{PSNR}\left(u,\hat{u}\right)=20{\log }_{10}\frac{255}{\Vert u-\hat{u}\Vert }\text{dB}\text{.}$

where u represents blur- and noise-free image, $\hat{u}$ represents the recovered image.

Figure 4(a) is the original image, Figures 4(b)-(e) are the images with four different levels of blur and noise.

Firstly, CG algorithm and MR algorithm are used for image restoration. When the stopping rule is met, the iteration times of both algorithms are 4. The restoration effect of four different levels of blur- and noise-contaminated images are shown in Figure 5.

Table 1 and Figure 5 show that the MR algorithm has a faster computational speed than the CG algorithm in the same case. For b₁v₁ and b₂v₂, the recovery effect of CG and MR is almost the same. For b₃v₃ and b₄v₄, the MR algorithm has better recovery effect.

(a) (b) (c) (d) (e)

Figure 4. Original images and contaminated images with different levels of blur and noise. (a) Original image; (b) Blur and noise image b₁v₁; (c) Blur and noise image b₂v₂; (d) Blur and noise image b₃v₃; (e) Blur and noise image b₄v₄.

Table 1. Numerical results of CG and MR.

In the following experiments, with different smoothers, Algorithm 1 and Algorithm 2 are abbreviated as shown in the Table 2. Algorithm 1 and Algorithm 2 with different smoother. Next, we compare the numerical results of restoration Algorithm 1 and Algorithm 2.

From Table 3 and Table 4, Figure 6 and Figure 7, it can be found that,

1) The quadratic interpolation can obtain better accuracy initial value. So IECMG-P is better than IECMG-L.

2) New extrapolation with quadratic interpolation can improve the accuracy of the initial values. So EECMG is better than IECMG.

3) With MR as smoother, EECMG (MR) and IECMG (MR) have the best restoration quality and spend the least recovery time.

Moreover, Figure 5, Figure 8 and Figure 9 show that EECMG (MR) is best, whether PSNR or CPU time.

Finally, we show the edge preserving ability of the new extrapolation cascadic multigrid method.

Figure 10(a) is the original image, Figure 10(b) is the third-level polluted image. Figure 10(c) is the restored image by using EECMG (MR), which PSNR = 27.9617. The good effect of edge preservation is shown in Figure 10(d).

(a) (b) (c) (d) (e) (f) (g) (h)

Figure 5. Numerical results of CG and MR. (a) CG, b₁v₁; (b) CG, b₂v₂; (c) CG, b₃v₃; (d) CG, b₄v₄; (e) MR, b₁v₁; (f) MR, b₂v₂; (g) MR, b₂v₂; (h) MR, b₄v₄.

Table 2. Algorithm 1 and Algorithm 2 with different smoother.

(a)(b)

Figure 6. PSNR and Time of IECMG (CG) and EECMG (CG). (a) PSNR; (b) Time.

Table 3. Numerical results of IECMG (CG) and EECMG (CG).

(a)(b)

Figure 7. PSNR and TIME of IECMG (MR) and EECMG (MR). (a) PSNR; (b) Time.

Table 4. Numerical results of IECMG (MR) and EECMG (MR).

(a) (b) (c) (d)

Figure 8. Recovery image of EECMG (MR) for four different blur and noise levels. (a) EECMG, b₁v₁, PSNR = 63.0973; (b) EECMG, b₂v₂, PSNR = 41.7086; (c) EECMG, b₃v₃, PSNR = 35.2055; (d) EECMG, b₄v₄, PSNR = 33.0595.

(a)(b)

Figure 9. Numerical results of CG, MR, IECMG and EECMG. (a) PSNR; (b) Time.

(a) (b) (c) (d)

Figure 10. Edge preserving effect of EECMG (MR). (a) Original image; (b) Blur and noise image with b₃v₃; (c) Recovery image; (d) Recovery of image edge.

6. Conclusions

Combining the advantages of the new extrapolation method and the economical cascadic multigrid method, this paper proposes a new extrapolation economical cascadic multigrid method for image restoration, which not only improves the computational efficiency but also ensures the good quality of restoration. Numerical experiments show that,

1) The edge preserving denoising operator can achieve the effect of preserving image edges and removing noise;

2) The weighted least squares approximation smooths the image locally and reduces the staircase effect;

3) The new extrapolation formula and quadratic interpolation construct a nonlinear continuation operator, which can provide more accurate pixel values for the fine grid level and ensure good recovery quality;

4) The iteration number selection method controls the iteration number of each grid level, which greatly reduces the computational work.

Although this paper puts forward several kinds of algorithms for linear model of image restoration problem and has made some research achievements, there are still many works worth further research, such as:

1) Whether the algorithm proposed in this paper can also have good restoration effect on other image restoration models?

2) Whether the algorithm in this paper can restore large-size images?

3) Can the algorithm in this paper be applied to the restoration of three-dimensional images?

Acknowledgements

This work was supported by Natural Science Foundation of China (11661027), the Guangxi Natural Science Foundation (2020GXNSFAA159143), and National Natural Science Foundation of China (12161027).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Perona, P. and Malik, J. (1990) Scale-Space and Edge Detection Using Anisotropic Diffusion. IEEE Transactions on Pattern Analysis and Machine Intelligence, 12, 629-639.
[2]	Handlovičová, A., Mikula, K. and Sgallari, F. (2003) Semi-Implicit Complementary Volume Scheme for Solving Level Set Like Equations in Image Processing and Curve Evolution. Numerische Mathematik, 93, 675-695. https://doi.org/10.1007/s002110100374
[3]	Weickert, J., Romeny, B.M.T.H. and Viergever, M.A. (1998) Efficient and Reliable Schemes for Nonlinear Diffusion Filtering. IEEE Transactions on Image Processing, 7, 398-410. https://doi.org/10.1109/83.661190
[4]	Hochstenbach, M.E. and Reichel, L. (2012) Combining Approximate Solutions for Linear Discrete Ill-Posed Problems. Journal of Computational and Applied Mathematics, 236, 2179-2185. https://doi.org/10.1016/j.cam.2011.09.040
[5]	Morigi, S., Reichel, L., Sgallari, F. and Zama, F. (2006) Iterative Methods for Ill-Posed Problems and Semiconvergent Sequences. Journal of Computational and Applied Mathematics, 193, 157-167. https://doi.org/10.1016/j.cam.2005.05.028
[6]	Zhang, H. and Dai, H. (2021) The Regularizing Properties of Global GMRES for Solving Large-Scale Linear Discrete Ill-Posed Problems with Several Right-Hand Sides. Applied Numerical Mathematics, 164, 57-71. https://doi.org/10.1016/j.apnum.2020.08.021
[7]	Buades, A., Coll, B. and Morel, J.M. (2006) The Staircasing Effect in Neighborhood Filters and Its Solution. IEEE Transactions on Image Processing, 15, 1499-1505. https://doi.org/10.1109/TIP.2006.871137
[8]	Bornemann, F.A. and Deuflhard, P. (1996) The Cascadic Multigrid Method for Elliptic Problems. Numerische Mathematik, 75, 135-152. https://doi.org/10.1007/s002110050234
[9]	Shi, Z.-C., Xu, X.-J., Huang, Y. and Huang, Q. (2007) Economical Cascadic Multigrid Method. Science in China Series A: Mathematics, 50, 1083-1098. https://doi.org/10.1007/s11425-007-0127-z
[10]	Liu, Y. (2010) A New Extrapolation Cascadic Multi-Grid Method on Quasiuniform Meshes. International Journal of Applied and Computational Mathematics, 5, 573-583. https://mathscinet.ams.org/mathscinet-getitem?mr=2738414
[11]	Chen, C.M., Hu, H.L., Xie, Z.Q. and Li, C.L. (2008) Analysis of Extrapolation Cascadic Multigrid Method. Science in China Series A: Mathematics, 51, 1349-1360. https://doi.org/10.1007/s11425-008-0119-7
[12]	Chen, C.M., Shi, Z.-C. and Hu, H.L. (2011) On Extrapolation Cascadic Multigrid Method. Journal of Computational Mathematics, 29, 684-697. https://doi.org/10.4208/jcm.1110-m11si05
[13]	Donatelli, M. (2012) An Iterative Multigrid Regularization Method for Toeplitz Discrete Ill-Posed Problems. Numerical Mathematics: Theory, Methods and Applications, 5, 43-61. https://doi.org/10.4208/nmtma.2011.m12si03
[14]	Donatelli, M. (2010) A Multigrid for Image Deblurring with Tikhonov Regularization. Numerical Linear Algebra with Applications, 12, 715-729. https://doi.org/10.1002/nla.446
[15]	Reichel, L. and Shyshkov, A. (2010) Cascadic Multilevel Methods for Ill-Posed Problems. Journal of Computational and Applied Mathematics, 233, 1314-1325. https://doi.org/10.1016/j.cam.2009.03.019
[16]	Morigi, S., Reichel, L., Sgallari, F. and Shyshkov, A. (2008) Cascadic Multiresolution Methods for Image Deblurring. SIAM Journal on Imaging Sciences, 1, 51-74. https://doi.org/10.1137/070694065
[17]	Morigi, S., Reichel, L. and Sgallari, F. (2010) Cascadic Multilevel Methods for Fast Nonsymmetric Blur-and Noise-Removal. Applied Numerical Mathematics, 60, 378-396. https://doi.org/10.1016/j.apnum.2009.07.007