# Performance of the Restarted Homotopy Perturbation Method and Split Bregman Method for Multiplicative Noise Removal

Authors: Yu Du Han; Jae Heon Yun
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

In this paper, we first propose restarted homotopy perturbation methods (RHPM) for multiplicative noise removal of the RLO and AA2 models. The main difficulty in applying the RHPM to the nonlinear denoising problem is settled by using binomial series techniques. We next propose the split Bregman methods for multiplicative noise removal of the RLO and AA2 models. The difficulty in applying the split Bregman method to the nonlinear denoising problem can be handled by transforming ill-conditioned linear systems into well-conditioned linear systems using splitting techniques of singular matrices. Lastly, numerical experiments for several test problems are provided to demonstrate the efficiency and reliability of the RHPM and split Bregman methods.


---

## Body

## 1. Introduction

Image denoising which is one of the fundamental problems in image processing is to recover an original image from a given noisy image. LetΩ be an open rectangular domain in R2, let z:Ω→R be an observed noisy image to be denoised, and let u be the original unknown image to restore. The mathematical models of noisy images can be classified into additive or multiplicative ones according to the relations of noises and clean images, i.e.,(1)z=u+η~orz=u·η,where η~ and η are unknown additive and multiplicative noise functions, respectively. In this paper, we focus on multiplicative noise removal for the observed images corrupted by multiplicative Gaussian white noise or Gamma noise with mean equal to one.There are a lot of investigations for additive model (1) using the following general nonlinear variational model [1, 2]:(2)minu⁡Eu=∫Ωφ∇udx+α2∫Ωu-z2dx,where the first term of the right-hand side is a regularizer in general form to ensure the smoothness of the desired image with edge preserving, the second term is a data fitting term which reflects the fidelity between the original image and observed noisy image, and α>0 is a regularization parameter.Many approaches have been proposed to remove the additive noise, like variational approaches [3, 4], the stochastic approach [5, 6], wavelet approaches [7, 8], framelet approaches [9, 10], and RHPM approaches [11], to name a few. It is well known that the model introduced by Rudin, Osher, and Fatemi [12] (called the ROF model) is a popular model for image denoising due to its property of preserving edges. If φ∇u=∇u in (2), then (2) reduces to the ROF model using Total Variation (TV): (3)minu⁡Eu=∫Ω∇udx+α2∫Ωu-z2dx.Since the images restored by the ROF model are in the space of bounded variation functions, the ROF model preserves sharp edges or object boundaries, and it is especially effective on restoring the piecewise constant images. However, the ROF model does not preserve the details and textures since it yields staircase effects.Multiplicative noise models occur in the study of several coherent imaging systems, such as synthetic aperture radar (SAR) and sonar (SAS), ultrasound or laser imaging, and magnetic resonance imaging (MRI). The researches on variational models of additive noise removal have been successful, but the investigations on variational models of multiplicative noise removal have not been done much. The difference between a variational model of additive and that of multiplicative noise removal is that the data fitting term depends on noise distributions, such as Gamma, Poisson, Rayleigh, and Gauss distributions due to different means of image acquisition. SAR images are Gamma distributions [13], medical ultrasound or resonance images fit Rayleigh distributions [14], and astronomical microscopy images and medical SPECT/PET images fit Poisson distributions [15, 16]. If the multiplicative noise is not too strong, the abovementioned noise can be considered approximately to the Gauss distributions [17].For the multiplicative Gaussian white noise with mean equal to one, Rudin, Lions, and Osher [17] proposed the following variational diffusion model (called the RLO model):(4)minu⁡Eu=∫Ω∇udx+∫Ωαzu+β2zu-12dx,where α and β are weighted parameters.For the multiplicative Gamma noise with mean equal to one, Jin and Yang [18] applied the exponential transformation u→eu introduced by Huang et al. [19] with the fitting term of the AA model [13] and proposed the following variational model (called the AA2 model in this paper):(5)minu⁡Eu=∫Ω∇udx+α∫Ωu+z·e-udx,where α is a weighted parameter.The investigations on the algorithm design of the classic TV model have attracted a lot of interest since it was invented to improve its computation efficiency and the quality of restored images, such as the artificial time marching method [12], fixed-point iterative method [20], primal-dual method [21], dual method [22], split method [23], Bregman iterative method [24], split Bregman method [25], and unbiased Box-Cox transformation method [26]. Among them, the split Bregman method combines the advantages of the split method in easy implementation and the Bregman iterative method in good quality of restored images. The variational models of multiplicative noise removal are more complex than the corresponding ones of additive noise removal, but its previous studies focused on models according to different noise distributions. Some researchers pay their attention to algorithm design, such as [12, 20, 27, 28]. The starting point of our investigation is the RLO model of (4) and AA2 model of (5) with the TV regularization term. The purpose of this paper is to propose restarted homotopy perturbation methods (RHPM) and split Bregman methods for multiplicative noise removal models (4) and (5).The paper is organized as follows. In Section2, we review the TM (time marching) method for multiplicative noise removal models. In Section 3, we briefly review a restarted homotopy perturbation method (RHPM). In Section 4, we propose RHPM algorithms for the multiplicative noise removal models. In Section 5, we propose split Bregman algorithms for the multiplicative noise removal models. In Section 6, we provide numerical experiments for several test problems in order to evaluate the performance of the RHPM and split Bregman methods. Lastly we provide concluding remarks.


## 2. Review of the TM Method for Multiplicative Noise Removal

The Euler-Lagrange equations corresponding to RLO model (4) and AA2 model (5) lead to the following nonlinear elliptic PDEs with the homogeneous Neumann boundary conditions, respectively: (6)∇·∇ux,y∇ux,y+α-βzx,yux,y2+βzx,y2ux,y3=0,x,y∈Ω,∇ux,y·n→=0,x,y∈∂Ω,(7)∇·∇ux,y∇ux,y-α1-zx,ye-ux,y=0,x,y∈Ω,∇ux,y·n→=0,x,y∈∂Ω,where n→ is the unit normal vector exterior to the boundary ∂Ω.We briefly describe the TM (time marching) method [12] only for AA2 model (7). To avoid division by zero in numerical implementation, we replace the nondifferentiable term ∇u(x,y) in (7) with a smooth approximation term (8)∇ux,yε=∇ux,y2+εfor  a  smallε>0.Then (7) is transformed into(9)∇·∇ux,y∇ux,yε-α1-zx,ye-ux,y=0,x,y∈Ω,∇ux,y·n→=0,x,y∈∂Ω.In order to apply the TM method to (9), we consider the following time-dependent nonlinear parabolic PDE corresponding to (9): (10)∂ux,y;t∂t=∇·∇ux,y;t∇ux,y;tε-α1-zx,ye-ux,y;t,x,y∈Ω∇u·n→=0,x,y∈∂Ω,where u(x,y,0) is given.For numerical implementation, let us assume that the domainΩ has been split into m×n cells where the grid points are located at (xi=ihx,yj=jhy), 1≤i≤m,1≤j≤n,tk=k△t, where △t and k=1,2,… denote the time step and iteration number, respectively. We denote the values of u(x,y;t) at the grid points (xi,yj;tk) by ui,jk and ui,j0=z(xi,yj).Without loss of generality, we can assume thathx=hy=1. Then, nonlinear PDE (10) can be approximated by the following finite difference formula:(11)ui,jk+1=ui,jk+△tκui,jk-α1-zi,je-ui,jk,where for 1≤i≤m,1≤j≤n(12)κui,jk=∇·∇u∇uεi,jk=∂∂xux∇uε+∂∂yuy∇uεi,jk=△-x△+xui,jk△+xui,jk2+△+yui,jk2+ε+△-y△+yui,jk△+xuijk2+△+yui,jk2+εi,jwith(13)△+xui,jk=ui+1,jk-ui,jk,△-xui,jk=-ui-1,jk+ui,jk△+yui,jk=ui,j+1k-ui,jk,△-yui,jk=-ui,j-1k+ui,jk,u0,jk=u1,jk,um+1,jk=um,jk,ui,0k=ui,1k,ui,n+1k=ui,nk.In a similar way as was done for the AA2 model, we can obtain the following formula for the RLO model:(14)ui,jk+1=ui,jk+△tκui,jk+zi,jα-βui,jk+βzi,jui,jk3+ε,where κ(ui,jk) and the boundary conditions are defined the same as those in the AA2 model.Hence, we obtain TM Algorithms1 and 2 for the AA2 and RLO models. maxit denotes the maximum number of iterations, z denotes the noisy image, uk denotes the restored image at the kth iteration, and MAE(k) denotes themean absolute error at the kth iteration, i.e., (15)MAEk=uk-uk-11m·n,where ·1 stands for the l1-norm.Algorithm 1:TM algorithm for the AA2 model.
1: MAE(0)=255, ui,j0=zi,j for 1≤i≤m,1≤j≤n 2:for  k=0 to maxit  do 3:Compute κ(ui,jk) using Equation (12) for 1≤i≤m,1≤j≤n 4:Compute ui,jk+1 using Equation (11) for 1≤i≤m,1≤j≤n 5:if  MAE(k)<MAE(k+1)  then 6:Stop 7:end if 8:end forAlgorithm 2:TM algorithm for the RLO model.
1: MAE(0)=255, ui,j0=zi,j for 1≤i≤m,1≤j≤n 2:for  k=0 to maxit  do 3:Compute κ(ui,jk) using Equation (12) for 1≤i≤m,1≤j≤n 4:Compute ui,jk+1 using Equation (14) for 1≤i≤m,1≤j≤n 5:if  MAE(k)<MAE(k+1)  then 6:Stop 7:end if 8:end for


## 3. Review of the Restarted Homotopy Perturbation Method (RHPM)

We briefly describe the homotopy perturbation method (HPM) for solving nonlinear partial differential equation(16)Au-fr=0,r∈Ωwith the boundary conditions (17)Bu,∂u∂n=0,r∈Γ,where A is a general differential operator, B is a boundary operator, f(r) is a known analytic function, and Γ is the boundary of the domain Ω. The operator A can be divided into two parts L and N, where L is linear and N is nonlinear. Therefore, (16) can be written as(18)Lu+Nu-fr=0.By using the homotopy technique, one can construct a homotopy v(r,p):Ω×[0,1]→R which satisfies (19)Hv,p=1-pLv-Lu0+pAv-fr=0,or, equivalently,(20)Hv,p=Lv-Lu0+pLu0+pNv-fr=0,where p∈[0,1] is an embedding parameter and u0 is the initial approximation of (16) which satisfies the boundary conditions. Clearly, we have (21)Hv,0=Lv-Lu0=0,Hv,1=Av-fr=0.The changing process of p from zero to unity is just that of v(r,p) changing from u0(r) to u(r). This is called deformation and also L(v)-L(u0) and A(v)-f(r) are called homotopic in topology. If the embedding parameter 0≤p≤1 is considered as a "small parameter" applying the classical perturbation technique, then the solution of (20) can be expressed as a power series in p, that is,(22)v=vr,p=v0+pv1+p2v2+⋯.Setting p→1, the solution u of (16) can be expressed as the sum of an infinite series(23)u=limp→1⁡vr,p=v0+v1+v2+⋯.Generally speaking, the computation ofvk in the HPM becomes more complicated as k increases. To circumvent this problem, the restarted homotopy perturbation method (RHPM) proposed by Han and Yun [11] repeats the HPM process by computing only the first few terms instead of computing infinite terms of vk. This is a big advantage of the RHPM as compared with the HPM. The general procedure of the RHPM for solving a nonlinear equation is described below (see [11] for more details).Letu00 be an initial approximation of a nonlinear equation. First, we compute the first k-th term approximation u01 in the following way: (24)v00=u00Computev10,v20,…,vk0using20and22u01=v00+v10+v20+⋯+vk0.Following the above computational process with u01 as an initial approximation, a new approximation u02 is computed as follows: (25)v01=u01Computev11,v21,…,vk1using20and22u02=v01+v11+v21+⋯+vk1.By repeating the above process, a new approximation u0l+1 can be computed from the initial approximation u0l, which was obtained at the lth step, in the following way: (26)v0l=u0lComputev1l,v2l,…,vklusing20and22u0l+1=v0l+v1l+v2l+⋯+vkl.Notice that v0l,v1l,…,vkl are coefficients of p0,p1,…,pk, respectively. Then the exact solution u can be approximated by u0l+1, that is, under some appropriate conditions u=liml→∞⁡u0l. The RHPM which computes only k vectors v1l,v2l,…,vkl every step is called the RHPM(k) method.


## 4. RHPM Algorithms for Multiplicative Noise Removal

In this section, we describe an application of the RHPM to multiplicative denoising problems of the RLO and AA2 models. Before doing this, we describe the binomial series which is used for an application of the RHPM. The difficulty in applying the RHPM to the TV-based image denoising problem comes from (36) and (62) (i.e., (vx2+vy2)-1/2, (vx2+vy2)-3/2, u-2, and u-3). In order to handle this difficulty, the following binomial series is used: for x<1,(27)1+xn=∑k=0∞nkxk=1+nx+nn-12!x2+nn-1n-23!x3+⋯.If n=-1/2 in (27), then(28)1+x-1/2=11+x=1-12x+38x2-516x3+⋯.Letting x=τ-1 in (28), for 0<τ<2  τ-1/2 can be approximated by(29)τ-1/2=1τ≈1-12τ-1+38τ-12=158-54τ+38τ2.If τ>0 is large, then, by introducing a small parameter ε>0 and using (29),(30)τ-1/2=1τ=εετ≈ε158-5ε4τ+3ε28τ2.Similarly, if n=-3/2 and x=τ-1 in (27), then, by introducing a small parameter ε>0,(31)τ-3/2=1τ3=ε3ετ3≈ε3358-21ε4τ+15ε28τ2.Note that a small parameter ε>0 in (30) and (31) is used to guarantee convergence of the binomial series when τ>0 is large.We first describe how to apply the RHPM to the following time-dependent nonlinear parabolic PDE corresponding to AA2 model (7): (32)∂ux,y;t∂t=∇·∇ux,y;t∇ux,y;t-α1-zx,ye-ux,y;t,x,y∈Ω∇u·n→=0,x,y∈∂Ω.For simplicity of exposition, we describe the computational process for RHPM(2) (i.e., k=2). We use three terms of binomial series (27) for numerical computation of (36) and (62), so division by zero does not occur. Let the operators L, N, and f be defined as (33)Lv=∂vx,y;t∂t,fx,y=0,Nv=-∇·∇vx,y;t∇vx,y;t+α1-zx,ye-vx,y;t,where (34)∇vx,y;t=∂v∂x,∂v∂yT,∇vx,y;t=vx2+vy2.From (20), one obtains(35)∂v∂t=∂u0∂t-p∂u0∂t+p∇·∇vx,y;t∇vx,y;t-pα1-zx,ye-vx,y;t,where(36)∇·∇v∇v=vxx+vyyvx2+vy2-vx2vxx+vy2vyyvx2+vy23.Using the Taylor series expansion of e-v with the second order and power series expansion of v of form (22), e-v can be approximated as(37)e-v≈1-v+12v2=1-v0+pv1+p2v2+⋯+12v0+pv1+p2v2+⋯2=1-v0+12v02-p1-v0v1+⋯.Since the solution v of (35) is of form (22), by simple computation (38)vx=v0x+pv1x+p2v2x+⋯,vxx=v0xx+pv1xx+p2v2xx+⋯,vy=v0y+pv1y+p2v2y+⋯,vyy=v0yy+pv1yy+p2v2yy+⋯,vx2=v0x2+2v0xv1xp+⋯,vy2=v0y2+2v0yv1yp+⋯.From these equations, one easily obtains (39)vx2+vy2=D1+2D2p+⋯,vx2+vy22=D12+4D1D2p+⋯,vxx+vyy=D3+D4p+⋯,vx2vxx+vy2vyy=D5+D6p+⋯,where(40)D1=←v0x2+v0y2,D2=←v0xv1x+v0yv1y,D3=←v0xx+v0yy,D4=←v1xx+v1yy,D5=←v0x2v0xx+v0y2v0yy,D6=←v0x2v1xx+v0y2v1yy+2v0xv0xxv1x+2v0yv0yyv1y.Substituting τ=(vx)2+(vy)2 into (30) and using (39), one obtains (41)1∇v=εετ=εεvx2+vy2≈ε158-5ε4vx2+vy2+3ε28vx2+vy22=ε158-5ε4D1+3ε28D12+ε3-52+3ε2D1D2p+⋯=A1+B1D2p+⋯,where (A1)=←ε{15/8-5ε/4(D1)+3ε2/8D1)2 and (B1)=←ε3-5/2+3ε/2(D1).Hence, the first term of the right-hand side in (36) can be written as(42)vxx+vyy∇v≈A1D3+A1D4p+B1D2D3p+⋯.Substituting τ=(vx)2+(vy)2 into (31) and using (39), one obtains (43)1∇v3=ε3ετ3=ε3ε3vx2+vy23≈ε3358-21ε4vx2+vy2+15ε28vx2+vy22=ε3358-21ε4D1+15ε28D12+ε5-212+15ε2D1D2p+⋯=A2+B2D2p+⋯,where (A2)=←ε3{35/8-21ε/4(D1)+15ε2/8(D1)2} and (B2)=←ε5-21/2+15ε/2(D1).Hence, the second term of the right-hand side in (36) can be written as(44)vx2vxx+vy2vyy∇v3≈A2D5+B2D2D5p+A2D6p+⋯.From (36), (37), (42), and (44), (35) can be expressed as(45)∂v0∂t+p∂v1∂t+p2∂v2∂t+⋯=∂u0∂t-p∂u0∂t+pvxx+vyy∇nv-vx2vxx+vx2vyy∇v3-p1-z1-v0+12v02α-p2z1-v0v1α≈∂u0∂t-p∂u0∂t-A1D3-A2D5+p2A1D4-A2D6+B1D3-B2D5D2-pα1-z1-v0+12v02-p2αz1-v0v1+⋯.Comparing coefficients of p0, p1, and p2 in (44), (46)p0:∂v0∂t=∂u0∂t⇒v0=u0=z:given an initial approximation(47)p1:∂v1∂t=-∂v0∂t+A1D3-A2D5-α1-z1-v0+12v02(48)p2:∂v2∂t=A1D4-A2D6+B1D3-B2D5D2-αz1-v0v1.Notice that the initial approximation u0 is set to z which is a noisy image. Applying the inverse operator of L to (47), one obtains(49)v1=z-v0+A1D3-A2D5-α1-z1-v0+12v02t=v10+v11t,where (50)v10=←z-v0,v11=←A1D3-A2D5-α1-z1-v0+12v02.From (49), the partial derivatives of v1 are(51)v1x=v10x+tv11x,v1y=v10y+tv11y,v1xx=v10xx+tv11xx,v1yy=v10yy+tv11yy.Using (51), (D2),(D4), and (D6) can be expressed as(52)D2=v0xv10x+v0yv10y+v0xv11x+v0yv11yt=D20+D21t,D4=v10xx+v10yy+v11xx+v11yyt=D40+D41t,D6=v0x2v10xx+v0y2v10yy+2v0xv0xxv10x+v0yv0yyv10y+v0x2v11xx+v0y2v11yy+2v0xv0xxv11x+v0yv0yyv11yt=D60+D61t,where (53)D20=←v0xv10x+v0yv10y,D21=←v0xv11x+v0yv11y,D40=←v10xx+v10yy,D41=←v11xx+v11yy,D60=←v0x2v10xx+v0y2v10yy+2v0xv0xxv10x+v0yv0yyv10y,D61=←v0x2v11xx+v0y2v11yy+2v0xv0xxv11x+v0yv0yyv11y.Substituting (49) and (52) into (48) and applying the inverse operator of L to (48), one obtains(54)v2=A1D40-A2D60+B1D3-B2D5D20t+12A1D41-A2D61+B1D3-B2D5D21t2-αz1-v0v10t-α2z1-v0v11t2=v21t+v22t2,where (55)v21=←A1D40-A2D60+B1D3-B2D5D20-αz1-v0v10,v22=←12A1D41-A2D61+B1D3-B2D5D21-α2z1-v0v11.Now update u0 whose new value is set to v0+v1+v2. By repeating the aforementioned process with the new updated u0, we obtain the RHPM(2) method for the AA2 model.For numerical implementation of the RHPM(2) method, let us assume that the domainΩ has been split into m×n cells where the grid points are located at (xi=ihx,yj=jhy), 1≤i≤m,1≤j≤n. Let tk=k△t, where △t and k=1,2,… refer to the time step and iteration number, respectively. We denote the values of u(x,y;t) at the grid points (xi,yj;tk) by (uk)i,j. Without loss of generality, we can assume that h=hx=hy=1. For simplicity of exposition, we define the following notation: (56)v0k=v0ki,j=v0xi,yj;tk∣1≤i≤m,1≤j≤n,k=1,2,…,v1k=v1ki,j=v1xi,yj;tk∣1≤i≤m,1≤j≤n,k=1,2,…,v2k=v2ki,j=v2xi,yj;tk∣1≤i≤m,1≤j≤n,k=1,2,….The boundary conditions are defined the same as those in the TM method. Assume that (u00)i,j=z(xi,yj)=zi,j and (u0k)i,j has been computed for all i and j. Then (v0k)i,j=(u0k)i,j and discretization of (49) is given by(57)v1ki,j=v10i,j+v11i,j·tk,where (58)v10i,j=zi,j-v0ki,j,v11i,j=A1i,jD3i,j-A2i,jD5i,j-α1-zi,j1-v0ki,j+12v0ki,j2.Also, discretization of (54) is given by(59)v2ki,j=v21i,j·tk+v22i,j·tk2,where (60)v21i,j=A1i,jD40i,j-A2i,jD60i,j+B1i,jD3i,j-B2i,jD5i,jD20i,j-αzi,j1-v0ki,jv10i,j,v22i,j=12A1i,jD41i,j-A2i,jD61i,j+B1i,jD3i,j-B2i,jD5i,jD21i,j-α2zi,j1-v0ki,jv11i,j.We now computeu0k+1=v0k+v1k+v2k. Repeating the above process with u0k+1 as an initial approximation, we can obtain the RHPM(2) method, called Algorithm 3, for multiplicative noise removal of the AA2 model.Algorithm 3:RHPM(2) algorithm for the AA2 model.
1: MAE(0)=1, (u00)i,j=zi,j for all i and j 2:for  k=0 to maxit  do 3:Set v0k=u0k 4:Compute v1k using (57) for all i and j 5:Compute v2k using (59) for all i and j 6:u0k+1=v0k+v1k+v2k 7:if  MAE(k)<MAE(k+1)  then 8:Stop 9:end if 10:end forWe next describe how to apply RHPM(2) to the following time-dependent nonlinear parabolic PDE corresponding to RLO model (6): (61)∂ux,y;t∂t=∇·∇ux,y;t∇ux,y;t+α-βzx,yux,y;t2+βzx,y2ux,y;t3,x,y∈Ω∇u·n→=0,x,y∈∂Ω,where ux,y,0 is given. Let the operators L, N, and f be defined as(62)Lv=∂vx,y;t∂t,fx,y=0,Nv=-∇·∇vx,y;t∇vx,y;t-α-βzx,yvx,y;t2-βzx,y2vx,y;t3.In order to handle the difficulty in applying the RHPM to the RLO model, we use three terms of binomial series (27) for v-2 and v-3 in (62). That is, for a small parameter ε>0, v-2 and v-3 can be approximated by (63)1v2=ε2εv2≈ε26-8εv+3ε2v2,1v3=ε3εv3≈ε310-15εv+6ε2v2.Using (62) and (63), the RHPM(2) method for the RLO model can be derived in a similar way as was done for the AA2 model as follows: (64)∂v0∂t+p∂v1∂t+p2∂v2∂t+⋯=∂u0∂t-p∂u0∂t+pvxx+vyy∇v-vx2vxx+vx2vyy∇v3+pα-βzv2+pβz2v3≈∂u0∂t-p∂u0∂t+pA1D3-A2D5+zα-βA3+βzA4+p2A1D4-A2D6+B1D3-B2D5D2+p2zα-βB3+βzB4+⋯,where (A1),(A2),(B1),(B2),(D2),(D3),(D4), and (D5) are defined the same as the AA2 model and (A3),(A4),(B3), and (B4) are defined as (65)A3=ε26-8εv0+3ε2v02,B3=-ε28-6εv0,A4=ε310-15εv0+6ε2v02,B4=-ε315-12εv0.Comparing coefficients of p0, p1, and p2 in (64),(66)p0:∂v0∂t=∂u0∂t⇒v0=u0=z:givenaninitialapproximation(67)p1:∂v1∂t=-∂v0∂t+A1D3-A2D5+zα-βA3+βzA4(68)p2:∂v2∂t=A1D4-A2D6+B1D3-B2D5D2+zα-βB3+βzB4.Applying the inverse operator of L to (67), one obtains(69)v1=z-v0+A1D3-A2D5+zα-βA3+βzA4t=v10+v11t,where (v10) is the constant term and (v11) is the coefficient of t.Substituting (69) and (52) into (68) and applying the inverse operator of L to (68), one obtains(70)v2=A1D40-A2D60+B1D3-B2D5D20t+12A1D41-A2D61+B1D3-B2D5D21t2+zα-βB3+βzB4t+z2α-βB3+βzB4t2=v21t+v22t2,where (v21) and (v22) are coefficients of t and t2, respectively.For numerical implementation of the RHPM(2) method, assume that(u00)i,j=z(xi,yj)=zi,j and (u0k)i,j has been computed for all i and j. Then (v0k)i,j=(u0k)i,j, and discretizations of (69) and (70) are given by (71)v1ki,j=v10i,j+v11i,j·tk,(72)v2ki,j=v21i,j·tk+v22i,j·tk2,where (73)v10i,j=zi,j-v0ki,jv11i,j=A1i,jD3i,j-A2i,jD5i,j+zi,jα-βA3i,j+βzi,jA4i,jv21i,j=A1i,jD40i,j-A2i,jD60i,j+B1i,jD3i,j-B2i,jD5i,jD20i,j+zi,jα-βB3i,j+βzB4i,jv22i,j=12A1i,jD41i,j-A2i,jD61i,j+B1i,jD3i,j-B2i,jD5ijD21i,j+12zi,jα-βB3i,j+βzB4i,j.We now computeu0k+1=v0k+v1k+v2k. Repeating the above process with u0k+1 as an initial approximation, we can obtain the RHPM(2) method, called Algorithm 4, for multiplicative noise removal of the RLO model.Algorithm 4:RHPM(2) algorithm for the RLO model.
1: MAE(0)=1, (u00)i,j=zi,j for all i and j 2:for  k=0 to maxit  do 3:Set v0k=u0k 4:Compute v1k using (71) for all i and j 5:Compute v2k using (72) for all i and j 6:u0k+1=v0k+v1k+v2k 7:if  MAE(k)<MAE(k+1)  then 8:Stop 9:end if 10:end forFrom RHPM(2) Algorithms3 and 4, we can easily obtain RHPM(1) Algorithms 5 and 6 whose computational costs per iteration are much cheaper than those of Algorithms 3 and 4. In Algorithms 3–6, maxit, z, and u0k denote the maximum number of iterations, the noisy image, and the restored image at the kth iteration, respectively. Notice that all images are assumed to have an intensity range of [0,1] to guarantee numerical stability of the RHPM methods.Algorithm 5:RHPM(1) algorithm for the AA2 model.
1: MAE(0)=1, (u00)i,j=zi,j for all i and j 2:for  k=0 to maxit  do 3:Set v0k=u0k 4:Compute v1k using (57) for all i and j 5:u0k+1=v0k+v1k 6:if  MAE(k)<MAE(k+1)  then 7:Stop 8:end if 9:end forAlgorithm 6:RHPM(1) algorithm for the RLO model.
1: MAE(0)=1, (u00)i,j=zi,j for all i and j 2:for  k=0 to maxit  do 3:Set v0k=u0k 4:Compute v1k using (71) for all i and j 5:u0k+1=v0k+v1k 6:if  MAE(k)<MAE(k+1)  then 7:Stop 8:end if 9:end for


## 5. Split Bregman Algorithms for Multiplicative Noise Removal

From (4) and (5), the general variational models based on TV regularization for multiplicative noise removal can be written as(74)minuEu=∫Ω∇udx+∫ΩHudx.To derive the alternating split Bregman algorithm for general multiplicative noise removal problem (74), we first introduce an auxiliary vector w=(wx,wy)T so that u and w minimize the following unconstrained problem:(75)minu,w⁡∫Ωwdx+θ2∫Ωw-∇u2dx+∫ΩHudx,where w=(wx)2+(wy)2 and θ>0 is a penalty parameter. Then unconstrained minimization problem (75) can be solved by the following split Bregman algorithm using an auxiliary vector b=(bx,by)T:(76)uk+1,wk+1=arg⁡minu,w⁡∫Ωw+θ2w-∇u-bk2+Hudx,(77)bk+1=bk+∇uk+1-wk+1,where b0=w0=0 and u0 is set to z (noisy image).Minimizing (76) alternatingly with respect to u and w, one obtains thealternating split Bregman algorithm introduced in [25]:(78)uk+1=arg⁡minu⁡∫ΩHu+θ2wk-∇u-bk2dx,(79)wk+1=arg⁡minw∫Ωw+θ2w-∇uk+1-bk2dx.The Euler-Lagrange equations corresponding to (78) and (79) are as follows:(80)∂∂uHuk+1-θ∇·∇uk+1+θ∇·wk-bk=0(81)1θ·wk+1wk+1+wk+1-∇uk+1-bk=0,where wk=(wxk,wyk)T, bk=(bxk,byk)T, and wk+1=(wxk+1)2+(wyk+1)2.First we consider how to solve (80) for uk+1. For the case of RLO model (4), H(u)=αz/u+β/2z/u-12 and thus (82)∂∂uHu=β-αzu2-βz2u3.To avoid division by zero or near-zero for z/u2 and z2/u3, u-2 and u-3 are approximated using four terms of binomial series (27). That is, (83)1u2≈10-20u+15u2-4u3,1u3≈20-45u+36u2-10u3,where u need to be scaled so that the value of u lies between 0 and 2 to guarantee convergence of the binomial series. Thus ∂/∂uH(u) can be approximated as (84)∂∂uHu≈β-α10-20u+15u2-4u3z-β20-45u+36u2-10u3z2.From this approximate equation, ∂/∂uHuk+1 in (80) can be approximated as (85)∂∂uHuk+1≈-10α-βz-20βz2+20α-βz+45βz2uk+1-15α-βz+36βz2uk2+4α-βz+10βz2uk3.Substituting this approximate equation into (80), one obtains the following approximate PDE for the RLO model:(86)20α-βz+45βz2uk+1-θ∇·∇uk+1=10α-βz+20βz2+15α-βz+36βz2uk2-4α-βz+10βz2uk3-θ∇·wk-bk.For the case of AA2 model (5), H(u)=αu+ze-u. Using four terms of the Taylor series for e-u, ∂/∂uH(u) can be approximated as (87)∂∂uHu=α1-ze-u≈α1-z1-u+12u2-16u3.From this approximate equation,∂/∂uHuk+1 in (80) can be approximated as (88)∂∂uHuk+1≈α1-z+zuk+1-12zuk2+16zuk3.Substituting this approximate equation into (80), we can obtain the following approximate PDE for the AA2 model:(89)αzuk+1-θ∇·∇uk+1=-α1-z-12zuk2+16zuk3-θ∇·wk-bk.Now we consider the numerical approach for solving PDE problems of RLO model (86) and AA2 model (89). Let ui,j, zi,j, ((bx)i,j,(by)i,j)T, and ((wx)i,j,(wy)i,j)T be the (i,j)-components of u, z, b, and w, respectively, where 1≤i≤m and 1≤j≤n. For an m×n array v, let vec(v) denote a long vector of size mn which is defined by vec(v)=(v∗1T,v∗2T,…,v∗nT)T and v∗i=(v1,i,v2,i,…,vm,i)T. Let ∇~=(∇~x,∇~y)T:Rm×n→R2m×n be a discrete gradient operator defined by (90)∇~ui,j=∇~xui,j,∇~yui,jT,∇~xui,j=0,i=1ui,j-ui-1,j,1<i≤m,∇~yui,j=0,j=1ui,j-ui,j-1,1<j≤n.Then the discrete gradient operator ∇~ can be represented by a matrix operator B of size 2mn×mn: (91)B=In⊗DmDn⊗Im,where ⊗ denotes the Kronecker product, Im denotes the identity matrices of order m, and Dm denotes the following matrix of order m: (92)Dm=00⋯⋯0-11⋯⋯0⋮⋱⋱⋱⋮0⋯-1100⋯0-11m×m.Using the matrix operator B, it is easy to show that ∇·∇u and ∇·w in (86) and (89) can be approximated as (93)∇·∇u≈-BTBvecu,∇·w≈BTvecwxvecwy.Hence, we obtain the finite difference approximate equations for RLO model (86) and AA2 model (89) which are given by the following linear systems, respectively:(94)diag20α-βz+45βz2+θBTBvecuk+1=vec10α-βz+20β2+15α-βz+36βz2uk2-4α-βz+10βz2uk3-θBTvecwxk-bxkvecwyk-byk,(95)diagαz+θBTBvecuk+1=vec-α1-z-12zuk2+16zuk3-θBTvecwxk-bxkvecwyk-byk,where z2, (uk)2, and (uk)3 denote elementwise exponentiations, uk denotes the restored image at the kth iteration, all multiplications of two m×n arrays denote elementwise multiplications, and diag(z) is an mn×mn diagonal matrix whose diagonal elements are the same as vec(z).Next, we consider how to solve (81) for wk+1. Equation (81) can be expressed as(96)1+1θ·1wk+1wxk+1=∇xuk+1+bxk1+1θ·1wk+1wyk+1=∇yuk+1+byk,where ∇uk+1=(∇xuk+1,∇yuk+1)T. By squaring and adding two equations in (96), one can obtain (97)1+1θ·1wk+12wk+12=∇uk+1+bk2,which implies(98)wk+1=∇uk+1+bk-1θ.Substituting (98) into (96), we have a generalized shrinkage formula(99)wxk+1=∇uk+1+bk-1θ∇xuk+1+bxk∇uk+1+bkwyk+1=∇uk+1+bk-1θ∇yuk+1+byk∇uk+1+bk.Generalized shrinkage formula (99) can be also found in [25, 29]. To make this paper self-contained, we provided a detailed proof for (99). Numerical implementation for formula (99) can be done as follows: for each (i,j)-component of wk+1=(wxk+1,wyk+1)T, (100)wxk+1i,j=max∇~uk+1+bki,j-1θ,0∇~xuk+1i,j+bxki,j∇~uk+1+bki,j,000=0wyk+1i,j=max∇~uk+1+bki,j-1θ,0∇~yuk+1i,j+byki,j∇~uk+1+bki,j,000=0.where (∇~uk+1+bk)i,j=((∇~xuk+1)i,j+(bxk)i,j,(∇~yuk+1)i,j+(byk)i,j)T.Lastly, numerical implementation for (77) can be done as follows: (101)vecbxk+1vecbyk+1=vecbxkvecbyk+B·vecuk+1-vecwxk+1vecwyk+1.Linear systems (94) and (95) are singular or ill-conditioned problems, so direct or iterative methods fail to get an approximate solution for these problems. For this reason, we propose a new technique for solving linear systems (94) and (95). Let (102)A1=diag20α-βz+45βz2+θBTB,A2=diagαz+θBTBr1k=right-hand  side  vector  of94,r2k=right-hand  side  vector  of95uk+1=vecuk+1,uk=vecukD=diagonal  matrix  whose  diagonal  elements  are  the  same  as  those  ofθBTB.Then linear systems (94) and (95) can be rewritten as(103)A1uk+1=r1k,A2uk+1=r2k.In order to solve ill-conditioned linear systems (103) for uk+1, we propose the following iterative method:Algorithm: SOLVER(i)Choosey1=ukforl=1 to maxlSolveAi+δDyl+1=δDyl+rik for yl+1end foruk+1=yl+1In iti=1 or 2, uk refers to the restored image computed at the previous kth step, and δ>0 is a parameter which should be chosen so that the coefficient matrix Ai+δD is well conditioned. To study semiconvergence of SOLVER(i), we first introduce some important properties for an iterative method.Theorem 1 (see [30, 31]).
LetS∈RN×N be a symmetric positive semidefinite matrix and b∈RN. Let S=M-N with M nonsingular and let G=M-1N. If MT+N is symmetric positive definite, then G is semiconvergent, and an iterative method yl+1=Gyl+M-1b, or equivalently Myl+1=Nyl+b, is semiconvergent. In other words, yl generated by the iterative method converges to some solution to Sy=b for each y1.Lemma 2.
Letα≥β≥0 and δ>0. Then Ai(i=1,2) are symmetric positive semidefinite, and D and Ai+δD(i=1,2) are symmetric positive definite.Proof.
It can be easily shown that all diagonal elements ofBTB are positive. Hence all properties of this lemma clearly hold from the assumptions.Theorem 3.
Letα≥β≥0 and δ>0. Then yl generated by SOLVER(i) converges to a solution uk+1 to Aiu=rik.Proof.
Note thatAi=(Ai+δD)-δD. Let M=Ai+δD and N=δD. From Lemma 2, Ai is symmetric positive semidefinite, and M and N are symmetric positive definite. Hence, M is nonsingular and MT+N is symmetric positive definite. From Theorem 1, the iterative method Myl+1=Nyl+rik, which is equivalent to (Ai+δD)yl+1=δDyl+rik, is semiconvergent. Hence the proof is complete.SinceAi+δD is symmetric positive definite for α≥β≥0 and δ>0, the linear systems in SOLVER(i) can be solved using the CG (Conjugate Gradient) method [32, 33]. Notice that convergence rate of the CG method depends on the condition number of coefficient matrix Ai+δD. So we next provide some conditional properties for Ai+δD in SOLVER(i). For a square matrix G, let κ(G) denote the condition number of G. For a symmetric matrix S∈RN×N, let λk(S) denote the kth largest eigenvalue of the matrix S for 1≤k≤N.Lemma 4 (see [34]).
LetS∈RN×N and S+E∈RN×N be symmetric matrices. Then for each 1≤k≤N(104)λkS+λNE≤λkS+E≤λkS+λ1E.Lemma 5.
Letα≥β≥0, δ>0, and N=mn. Then for each i=1,2(105)λ1Ai+δλNDλNAi+δλ1D≤κAi+δD≤λ1Ai+δλ1DλNAi+δλND.Proof.
SinceAi and Ai+δD are symmetric, from Lemma 4 one obtains(106)λ1Ai+δλND≤λ1Ai+δD≤λ1Ai+δλ1DλNAi+δλND≤λNAi+δD≤λNAi+δλ1D.From Lemma 2, Ai+δD is symmetric positive definite and thus(107)κAi+δD=λ1Ai+δDλNAi+δD.From (106) and (107), this lemma follows.Lemma 6.
Letα≥β≥0, δ>0, and N=mn. Then for each i=1,2(108)λ1Ai+2δλNAi+4δ≤κAi+δD≤λ1Ai+4δλNAi+2δ.Proof.
Note that all diagonal elements of the diagonal matrixD belong to a set {2,3,4}. Hence λ1(D)=4 and λN(D)=2. From Lemma 5, this lemma follows.Theorem 7.
Letα≥β≥0, δ>0, and N=mn. If 2λN(Ai)<λ1(Ai), then, for a given ϵ>0, there exists δ>0 such that 1/2<κ(Ai+δD)<2+ϵ.Proof.
Forδ>0, let (109)fδ=λ1Ai+4δλNAi+2δ,gδ=λ1Ai+2δλNAi+4δ.From Lemma 6, g(δ)<κ(Ai+δD)<f(δ). By simple calculation, one has(110)f′δ=4λNAi-2λ1AiλNAi+2δ2,g′δ=2λNAi-4λ1AiλNAi+4δ2.Since 2λN(Ai)<λ1(Ai), from (110) f′(δ)<0 and g′(δ)<0. Also note that limδ→∞f(δ)=2 and limδ→∞g(δ)=1/2. Thus f is a decreasing function which converges to 2, and g is a decreasing function which converges to 1/2. It follows that g(δ)>1/2 for any δ>0, and for a given ϵ>0 there exists a large δ>0 such that 2<f(δ)<2+ϵ. Since g(δ)<κ(Ai+δD)<f(δ), for a given ϵ>0 there exists a large δ>0 such that 1/2<κ(Ai+δD)<2+ϵ. Hence the proof is complete.Notice that the condition2λN(Ai)<λ1(Ai) in Theorem 7 usually holds when Ai is a singular or ill-conditioned matrix. Theorem 7 implies that the condition number of Ai+δD is small if a large number of δ>0 is chosen. In other words, for a suitably chosen large δ>0 the coefficient matrix Ai+δD in SOLVER(i) is well conditioned, so that the CG method for solving the linear systems in SOLVER(i) converges fast to the exact solution. The split Bregman method is an iterative method which monotonically decreases the residual norm between the original image and noisy image, so we do not have to get an accurate solution every iteration. For this reason, we have used maxl=1 and the CG method with tolerance set to 1×10-3 to get an approximate solution of the linear system in SOLVER(i). Notice that the CG method always achieves the desired tolerance within 15 iterations.Finally, one can obtain the split Bregman algorithms for the RLO and AA2 models. In Algorithms7 and 8,maxit denotes the maximum number of outer iterations, ·2 stands for the l2-norm, andtol denotes the tolerance of the stopping criterion which is set to 3×10-3 for the test problems used in this paper. Notice that all images are assumed to have an intensity range of [0,1] to guarantee convergence of the binomial series used in the split Bregman algorithm.Algorithm 7:Split Bregman algorithm for the RLO model.
1:b0=w0=0 and u0=z 2:for  k=0 to maxit  do 3:Compute uk+1 using SOLVER(1) 4:Compute wk+1 using Equation (100) 5:Compute bk+1 using Equation (101) 6:if  uk-uk+12uk+12<tol  then 7:Stop 8:end if 9:end forAlgorithm 8:Split Bregman algorithm for the AA2 model.
1:b0=w0=0 and u0=z 2:for  k=0 to maxit  do 3:Compute uk+1 using SOLVER(2) 4:Compute wk+1 using Equation (100) 5:Compute bk+1 using Equation (101) 6:if  uk-uk+12uk+12<tol  then 7:Stop 8:end if 9:end for


## 6. Numerical Experiments

In this section, we provide numerical performance results of the TM, split Bregman, RHPM(2), and RHPM(1) methods for multiplicative noise removal problems. Numerical tests are performed using MATLAB R2016a on a personal computer with 3.6GHz CPU and 8GB RAM. The multiplicative noisy images are generated by multiplying Gaussian white noise or Gamma noise to the clean image using the built-in MATLAB functionrandn orrandg. More specifically, the multiplicative noisy images z are generated by (111)η=1+0.01×randnm,n;z=u·η,η=0.01×randg10.01,m,n;z=u·η,where u is the clean image, η=1+0.01×randn(m,n) generates the Gaussian white noise with mean 1 and variance 0.01, η=0.01×randg1/0.01,(m,n) generates the Gamma noise with mean 1 and variance 0.01, and u·η is componentwise multiplication of u and η.In order to illustrate the efficiency and reliability of the proposed denoising algorithms (i.e., Algorithms1–8), we provide numerical results for 4 test images such as the Caribou, Boat, Jet Plane, and Lake. The pixel size of the Caribou image is 256×256, and the pixel size of the other 3 test images is 512×512. To evaluate the quality of the restored images, we use the peak signal-to-noise ratio (PSNR) between the restored image and original image which is defined by (112)PSNR=10log10⁡m·n·maxi,j⁡ui,j2u-u~F2,where ·F refers to the Frobenius norm and u and u~ are the original and restored images, respectively. Also ui,j stands for the value of the original image u at the pixel point (i,j) and m·n is the total number of pixels. It is generally true that the larger PSNR stands for the better quality of the denoised image.For numerical experiments, the TM method uses the test images with an intensity range of[0,255], the split Bregman and RHPM methods use the test images with an intensity range of [0,1], and the noise η is Gaussian white noise or Gamma noise with mean 1 and variance 0.01 or 0.03. For all test problems, an initial image was set to the observed noisy image z, and we have used maxit=200, α=0.2, and β=0.1 for the RLO model and α=10-4 for the AA2 model. For the TM and RHPM methods (i.e., Algorithms 1–6), we have used △t=0.35, and the iteration was terminated when MAE(k), defined in Section 2, started to increase. For the split Bregman method (i.e., Algorithms 7 and 8), we have used θ=255, and the iteration was terminated when the following stopping criterion was satisfied: (113)uk-uk+12uk+12<tolwhere tol=3×10-3 for all test problems.Tables1 and 2 provide numerical results for the TM, split Bregman, RHPM(1), and RHPM(2) methods corresponding to the AA2 model with multiplicative Gamma noise and the RLO model with multiplicative Gaussian white noise, respectively. In Tables 1 and 2, P0 represents the PSNR values for the noisy images,Psnr represents the PSNR values for the restored images, parameters ϵ and δ are chosen as the best one by numerical tries,Iter denotes the number of iterations, andTime denotes the elapsed CPU time in seconds.Table 1
Numerical results of multiplicative noise removal for the AA2 model.
Variance=0.01 Variance=0.03 Method Image Caribou Boat Jet Plane Lake Caribou Boat Jet Plane Lake P 0 24.89 25.34 22.12 24.65 20.12 20.58 17.34 19.90 TM ε 0.012 0.011 0.011 0.033 0.05 0.049 0.072 0.2 Iter 41 39 53 48 63 67 99 91 Psnr 28.25 29.33 27.77 27.25 25.57 26.12 23.95 23.60 Time 1.27 4.81 6.55 6.08 1.75 7.32 10.71 11.75  SplitBregman δ 24 30 16 28 17 21 12 22 Iter 16 16 14 16 20 21 18 22 Psnr 30.18 29.82 28.94 29.38 27.82 27.48 26.50 27.03 Time 0.45 1.91 1.65 1.82 0.55 2.52 2.05 2.51  RHPM(1) ε 0.00062 0.0008 0.0008 0.00082 0.00014 0.00014 0.00014 0.00014 Iter 23 21 21 21 53 53 53 53 Psnr 29.67 29.58 27.82 29.06 26.18 26.34 23.76 25.76 Time 0.03 1.08 0.93 1.01 1.47 2.81 2.77 2.73  RHPM(2) ε 0.052 0.039 0.067 0.044 0.071 0.07 0.078 0.073 Iter 3 3 3 3 3 3 3 3 Psnr 30.38 29.92 29.05 29.49 27.53 27.36 25.58 26.94 Time 0.12 0.49 0.49 0.49 0.14 0.5 0.5 0.5Table 2
Numerical results of multiplicative noise removal for the RLO model.
Variance=0.01 Variance=0.03 Method Image Caribou Boat Jet Plane Lake Caribou Boat Jet Plane Lake P 0 24.83 25.35 22.10 24.67 20.05 20.58 17.33 19.90 TM ε 0.0004 0.00021 0.00019 0.00019 0.0004 0.00042 0.0004 0.038 Iter 40 39 53 47 63 66 97 84 Psnr 28.36 29.38 27.76 27.33 25.60 26.24 24.10 23.97 Time 1.57 6.50 8.69 8.00 2.70 11.24 16.53 14.33  SplitBregman δ 24 30 16 28 18 22 12 22 Iter 16 16 14 16 21 22 18 22 Psnr 30.20 29.84 28.92 29.39 27.84 27.53 26.52 27.03 Time 0.58 2.43 2.41 2.39 0.59 2.40 2.06 2.47  RHPM(1) ε 0.0005 0.00045 0.00045 0.00045 0.00014 0.00014 0.00014 0.00014 Iter 27 29 28 29 53 53 53 53 Psnr 29.71 29.60 27.75 29.06 26.14 26.38 23.78 25.76 Time 0.30 1.38 1.38 1.39 0.64 2.48 2.24 2.45  RHPM(2) ε 0.053 0.039 0.067 0.044 0.071 0.07 0.077 0.073 Iter 3 3 3 3 3 3 3 3 Psnr 30.40 29.95 29.03 29.50 27.51 27.40 25.59 26.93 Time 0.14 0.41 0.42 0.40 0.59 0.36 0.46 0.42Figures1–4 show the denoised images by the TM, split Bregman, RHPM(1), and RHPM(2) methods for the AA2 model with multiplicative Gamma noise of variance 0.03. Figures 5–8 show the denoised images by the TM, split Bregman, RHPM(1), and RHPM(2) methods for the RLO model with multiplicative Gaussian white noise of variance 0.03. The first-row images are the original image, noisy image, and the denoised image by the TM method. The second-row images are the images denoised by the split Bregman, RHPM(1), and RHPM(2) methods.Figure 1
Denoising of multiplicative Gamma noise for the Caribou image with a256×256 size.
(a)
True image (b)
Noisy image (c)
TM denoising (d)
Split Bregman denoising (e)
RHPM(1) denoising (f)
RHPM(2) denoisingFigure 2
Denoising of multiplicative Gamma noise for the Boat image with a512×512 size.
(a)
True image (b)
Noisy image (c)
TM denoising (d)
Split Bregman denoising (e)
RHPM(1) denoising (f)
RHPM(2) denoisingFigure 3
Denoising of multiplicative Gamma noise for the Jet Plane image with a512×512 size.
(a)
True image (b)
Noisy image (c)
TM denoising (d)
Split Bregman denoising (e)
RHPM(1) denoising (f)
RHPM(2) denoisingFigure 4
Denoising of multiplicative Gamma noise for the Lake image with a512×512 size.
(a)
True image (b)
Noisy image (c)
TM denoising (d)
Split Bregman denoising (e)
RHPM(1) denoising (f)
RHPM(2) denoisingFigure 5
Denoising of multiplicative Gaussian white noise for the Caribou image with a256×256 size.
(a)
True image (b)
Noisy image (c)
TM denoising (d)
Split Bregman denoising (e)
RHPM(1) denoising (f)
RHPM(2) denoisingFigure 6
Denoising of multiplicative Gaussian white noise for the Boat image with a512×512 size.
(a)
True image (b)
Noisy image (c)
TM denoising (d)
Split Bregman denoising (e)
RHPM(1) denoising (f)
RHPM(2) denoisingFigure 7
Denoising of multiplicative Gaussian white noise for the Jet Plane image with a512×512 size.
(a)
True image (b)
Noisy image (c)
TM denoising (d)
Split Bregman denoising (e)
RHPM(1) denoising (f)
RHPM(2) denoisingFigure 8
Denoising of multiplicative Gaussian white noise for the Lake image with a512×512 size.
(a)
True image (b)
Noisy image (c)
TM denoising (d)
Split Bregman denoising (e)
RHPM(1) denoising (f)
RHPM(2) denoisingAs can be seen in Tables1 and 2, the split Bregman method denoises best for variance 0.03, while RHPM(2) denoises best for variance 0.01. That is, the split Bregman method has the highest PSNR values for variance 0.03, while RHPM(2) has the highest PSNR values for variance 0.01. Observe that RHPM(1) computes 2 vectors v0,v1 every iteration, while RHPM(2) computes 3 vectors v0,v1,v2 every iteration. It means that the computational cost of RHPM(2) per iteration is more expensive than RHPM(1). However, the convergence rate of RHPM(2) is much faster than RHPM(1), so that total CPU time of RHPM(2) is much less than RHPM(1). In addition, RHPM(2) denoises much better than RHPM(1), and RHPM(2) has the smallest CPU time of all the methods considered in this paper. Overall, the split Bregman and RHPM(2) methods perform better than the RHPM(1) and TM methods. The split Bregman method denoises almost as well as RHPM(2) for most cases, while it denoises significantly better than RHPM(2) for the Jet Plane image with variance 0.03. In this regard, the split Bregman method performs more stably than RHPM(2).


## 7. Conclusion

In this paper, we have proposed restarted homotopy perturbation methods (RHPM) and split Bregman methods for multiplicative noise removal of the RLO and AA2 models. For the RHPM method, we have used binomial series techniques to settle the main difficulty in handling nonlinear terms. For the split Bregman method, we have used splitting techniques of singular matrices to handle the difficulty in solving ill-conditioned linear systems. Numerical experiments have shown that these techniques work well for all test problems. More specifically, the split Bregman and RHPM methods perform better than the TM method.Numerical experiments also showed that the split Bregman and RHPM(2) methods perform better than RHPM(1). In addition, the split Bregman method denoises almost as well as RHPM(2) for most cases, but it performs more stably than RHPM(2) for noisy images with large variance. Notice that RHPM(2) has the smallest execution time because of the fastest convergence rate. Based on these facts, the RHPM(2) and split Bregman methods are preferred over RHPM(1), and the split Bregman method is preferred over RHPM(2) for noisy images with large variance. The proposed RHPM and split Bregman methods are only compared with the TM method, so future work will include comparison studies between the proposed methods and other existing methods for multiplicative noise removal models.



---

*Source: 7696798-2018-12-06.xml*