Performance of Restarted Homotopy Perturbation Method for TV-Based Image Denoising Problem

We first propose a restarted homotopy perturbation method (RHPM) for solving a nonlinear PDE problem which repeats HPM process by computing only the first few terms instead of computing infinite terms, and then we present an application of RHPM to TV(Total Variation-) based image denoising problem.Themain difficulty in applying RHPM to the nonlinear denoising problem is settled by using binomial series techniques. We also provide finite difference schemes for numerical implementation of RHPM. Lastly, numerical experiments for several test images are carried out to demonstrate the feasibility, efficiency, and reliability of RHPM by comparing the performance of RHPM with that of existing TM and recently proposed RHAMmethods.


Introduction
Digital images generated by various digital devices or used in practical applications usually contain some noises, so denoising is necessary before analyzing the images [1].Over the last decade, a number of image denoising techniques preserving important image information such as edges have been proposed to obtain the true image from the Noisy (known) image.These methods include filtering technique [2], wavelet techniques [3,4], and variational methods [5][6][7][8][9].
Denoising based on linear filters normally does not perform satisfactorily since both noise and edges contain high frequencies, while denoising based on the nonlinear models has been found to be successful.To preserve image edges, Rudin and Osher [8] and Rudin et al. [9] introduced a TV-(Total Variation-) based denoising model, which is based on a variational problem using the TV norm as a nonlinear nondifferentiable functional.Also Rudin et al. [9] used a nonlinear parabolic PDE (Partial Differential Equation) as a solution procedure for solving the Euler-Lagrange equation arising from the TV-based denoising model and proposed TM (Time Marching) method for solving the nonlinear parabolic PDE.This model is found to be a successful tool for image denoising and edge enhancement.High-order denoising models [10,11] improving the TV denoising model have been proposed, but Euler-Lagrange equation associated with the TV functional for those models is a nonlinear PDE [12].Due to the difficulties in obtaining analytic solutions of the nonlinear PDE, there has been many researches in developing implicit or explicit methods for the nonlinear PDE.Also, many good methods for image denoising problems such as high-order TV method [13,14], total generalized variation (TGV) method, block-matching and 3D filtering (BM3D) method [15,16], four-directional TV-based method [17], and an adaptive method based on Tikhonov and TV regularization [18] have been recently proposed.
For the last 15 years, HPM (homotopy perturbation method), first introduced by He [19], has been used by many mathematicians and engineers to solve various functional equations.In this method, the solution is considered as the sum of an infinite series which converges rapidly to the exact solutions.Using homotopy technique in topology, a homotopy is constructed with an embedding parameter  ∈ [0, 1] which is considered as a "small parameter."Many researches have been conducted in applying the HPM to a class of linear and nonlinear equations.In addition, HPM was further applied to nonlinear oscillators with discontinuities [20], nonlinear wave equations [21], boundary value problem [22], and many other subjects [23][24][25][26].Recently, Ghanbari et al. [27] proposed RHAM (restarted homotopy analysis method), based on HAM [28], for a nonlinear denoising model problem, which is a main motivation for this paper.The purpose of this paper is to propose a restarted homotopy perturbation method (RHPM), based on the HPM, for solving a nonlinear PDE problem and then to apply RHPM to the TV-based image denoising problem.
This paper is organized as follows.In Section 2, we review TM method for TV-based image denoising problem.In Section 3, we propose a restarted homotopy perturbation method (RHPM) for solving a nonlinear PDE problem.In Section 4, we present an application of RHPM to TV-based image denoising problem.The main difficulty in applying RHPM is settled by using binomial series techniques.In Section 5, numerical experiments for several test images are carried out to demonstrate the feasibility, efficiency, and reliability of RHPM by comparing the performance of RHPM with that of existing TM and recently proposed RHAM methods.Lastly, some conclusions are drawn.

Review of TM Method for TV-Based Image Denoising Problem
In this section, we only review TM method for solving a TVbased image denoising problem (see [27] for details of the recently proposed RHAM method).TM method has been proven to be a successful tool for image denoising and has been used extensively in many practical application.
Let Ω be a bounded domain, let an intensity function (, ) be the pixel values of an observed image, and let (, ) be an unknown noise function.Then, we would like to construct the desired clean image (, ) from the observed image (, ) such that  (, ) =  (, ) +  (, ) . (1) Image denoising problem is to reconstruct (, ) from (, ) which contains the additive noise (, ).It is well known that image denoising problem is an ill-conditioned problem, so regularization techniques should be used to approximate (, ) from (, ).Rudin et al. [9] proposed the use of the following TVbased minimization problem: where  > 0 is a regularization parameter.The first term is the TV of (, ), a regularizer, while the second term is a fidelity term ensuring that the denoised image (, ) will be close to the given image (, ).The regularization parameter  is important for balancing denoising and smoothing.That is, it can be observed that the parameter  may be small in the presence of high noise and  may be large for little noise.
where ⃗  is the unit normal vector exterior to the boundary Ω.
To avoid division by zero in numerical implementation, we replace the nondifferentiable term |∇(, )| by a smooth approximation as follows: for some small  > 0.
In order to solve the Euler-Lagrange equation ( 5), Rudin et al. [9] proposed TM method for solving the following timedependent nonlinear parabolic PDE For numerical implementation, let us assume that the domain Ω has been split into × cells, where the grid points are located at (  = ℎ  ,   = ℎ  ), 1 ≤  ≤ , 1 ≤  ≤ ,   = Δ, where Δ and  = 1, 2, . .., denote the time step and iteration number, respectively.We denote the values of (, ; ) at the grid points (  ,   ;   ) by    and  0  = (  ,   ).Without loss of generality, we can assume that ℎ  = ℎ  = 1.Then, the nonlinear PDE (6) can be approximated by the following finite difference formula: where with and boundary conditions are Hence, we obtain TM method, called Algorithm 1, for solving (6), where tol is the tolerance for stopping criterion, maxit is the maximum number of iterations,  is the Noisy image, V is the true image, and PSNR is the peak signal-to-noise ratio between the true and restored images which is defined in Section 5.

Description of Restarted Homotopy Perturbation Method (RHPM)
In this section, we first describe the homotopy perturbation method (HPM) for solving nonlinear Partial Differential Equation  () −  () = 0,  ∈ Ω (11) with the boundary conditions where  is a general differential operator,  is a boundary operator, () is a known analytic functions, and Γ is the boundary of the domain Ω.The operator  can be divided in two parts  and , where  is linear and  is nonlinear.Therefore, (11) can be written as By using homotopy technique, one can construct a homotopy V(, ) : Ω × [0, 1] → R which satisfies or, equivalently, where  ∈ [0, 1] is an embedding parameter and  0 is the initial approximation of (11) which satisfies the boundary conditions.Clearly, we have The changing process of  from zero to unity is just that of V(, ) changing from  0 () to ().This is called deformation and also (V) − ( 0 ) and (V) − () are called homotopic in topology.If the embedding parameter 0 ≤  ≤ 1 is considered as a "small parameter" applying the classical perturbation technique, then the solution of ( 15) can be expressed as a power series in ; that is, Setting  → 1, the solution  of ( 11) can be expressed as the sum of an infinite series  = lim Further description and convergence for the HPM can be found in [19][20][21][22]26].Before giving a general description of the restarted homotopy perturbation method (RHPM), we first provide a simple example of how RHPM can be applied to a PDE problem.

Example 1. Consider two-dimensional heat equation with variable coefficients of the form
with the initial condition (, , 0) =  2 .
Let  0 = (, , 0).From ( 15), we readily obtain a homotopy V(, ) which satisfies Then, the solution V of ( 21) can be expressed as in (17).Substituting ( 17) into ( 21), one obtains Comparing coefficients of terms with identical powers of , we obtain that Hence, we have Therefore, we obtain the homotopy solution of ( 21) Letting  → 1, we obtain the exact solution of (19) which is given by Method 2 (RHPM).We propose a new approach for solving (19), called RHPM, which is based on HPM.We choose  0 = (, , 0) =  2 as an initial approximation of (19).Let the operators , , and  be defined as in (20).Comparing coefficients of terms  0 and  1 in (22), one obtains Mathematical Problems in Engineering 5 Then, the first two term approximation of ( 18) is defined as follows: From ( 27), using  0 =  1 0 as an initial approximation, one easily obtains From ( 29), the new approximation  2 0 is given by From ( 27) with an initial approximation  0 =  2 0 , one obtains Hence, the new approximation  3 0 is given by Continuing this process, it can be seen that   0 converges to the exact solution  of (19), which is the same as the result of Method 1, as  → ∞.
Generally speaking, the computation of V  in the HPM becomes more complicated as  increases.To circumvent this problem, RHPM repeats the HPM process by computing only the first few terms instead of computing infinite terms of V  .This is a big advantage of RHPM as compared with HPM.From the idea of RHPM introduced in Example 1, the general procedure of RHPM for solving a nonlinear equation is described below.
Let  0 0 be an initial approximation of a nonlinear equation.First, we compute the first th term approximation  1 0 in the following way: 15) and ( 17) Following the above computational process with  1 0 as an initial approximation, a new approximation  2 0 is computed as follows: 15) and ( 17) By repeating the above process, a new approximation  ℓ+1 0 can be computed from the initial approximation  ℓ 0 , which was obtained at the ℓth step, in the following way: . ., V ℓ  using ( 15) and ( 17) Notice that V ℓ 0 , V ℓ 1 , . . ., V ℓ  are coefficients of  0 ,  1 , . . .,   , respectively.Then, the exact solution  can be approximated by  ℓ+1 0 ; that is, under some appropriate conditions,  = lim ℓ → ∞  ℓ 0 (see Example 1).In RHPM, which computes only

Application of RHPM to TV-Based Image Denoising Problem
In this section, we describe an application of RHPM to a TVbased image denoising problem.Before doing this, we first describe the binomial series which is used for an application of RHPM.The difficulty in applying RHPM to the TV-based image denoising problem comes from (45 ).In order to handle this difficulty, the following binomial series is used.For || < 1, If  = −1/2 in (36), then Letting  =  − 1 in (37), for 0 <  < 2, If  > 0 is large, then, by introducing a small parameter  > 0 and using (40), 1/ √  3 can be approximated by Note that a small parameter  > 0 is used to guarantee the convergence of the binomial series when  > 0 is large.Now, we describe how to apply RHPM to the TV-based image denoising problem (6).For simplicity of exposition, we describe computational process for RHPM(2) (i.e.,  = 2).Since we only use four terms of the binomial series for numerical computation of (45), division by zero does not occur and thus  = 0 is used for RHPM (see (39) and (41)).Let the operators , , and  be defined as follows: where Since the solution V of ( 44) is of the form (17), by simple computation From these equations, one easily obtains where The (V) and (V − ) in (44) can be expressed as Substituting  = (V  ) 2 + (V  ) 2 into (39), one obtains where Hence, the first term of the right hand side in (45) can be written as Substituting  = (V  ) 2 + (V  ) 2 into (41), one obtains where Hence, the second term of the right hand side in (45) can be written as From (45), (52), and ( 55), (44) can be expressed as Comparing the coefficients of  0 ,  1 , and  2 in (56), = : given an initial approximation (57) Applying the inverse operator of  to (58), one obtains where  = V 0 −  and  = ( 1 )( 3 ) − ( 2 )( 5 ).From (60), the partial derivatives of V 1 are where ( 1 ) =   +   , ( 2 ) =   +   , ( 3 ) =   +   , and ( 4 ) =   +   .From (61), ( 2 ), ( 4 ), and ( 6 ) can be expressed as Substituting (60) and ( 62) into (59) and applying the inverse operator of  to (59), one obtains Now, update  0 whose new value is set to By repeating the aforementioned process with the new updated  0 , we obtain RHPM(2) method.
Boundary conditions are defined the same as those in TM method.Assume that ( 0 0 )  = (  ,   ) =   and (  0 )  has been computed for all  and .Then, (V  0 )  = (  0 )  and discretization of (60) is given by Also, discretization of (63) is given by where ( We now compute Repeating the above process with  +1 0 as an initial approximation, we can obtain RHPM(2) method, called Algorithm 2, for solving (6).Since the computation of V  2 in RHPM(2) method is so complicated, we also provide RHPM(1) method, called Algorithm 3, whose computational cost is much cheaper than Algorithm 2. In RHPM(2) and RHPM(1) methods, tol, maxit, , and  are the tolerance for stopping criterion, the maximum number of iterations, the Noisy image, and the true image, respectively.In (68), the coefficients of ( 1 ) may have large values which cause amplification of some incorrect results.Hence, we have used an intensity range [0, 1] for all test images used in RHPM methods to guarantee numerical stability.

Numerical Experiments
In this section, numerical experiments for several test images are carried out to compare the performance of RHPM(2) and RHPM(1) with those of existing TM and recently proposed RHAM2 in [27].Numerical tests are performed using MATLAB R2014a on a personal computer with a processer at 3.6 GHz and 8 GB RAM.The Noisy images are generated by adding Gaussian white noise to the clean images using the built-in MATLAB function randn.For example, adding 15% where ‖ ⋅ ‖  refers to the Frobenius norm.
In order to illustrate efficiency and reliability of the proposed denoising algorithms (i.e., Algorithms 2 and 3), we provide experimental results for several images such as artificial, synthetic, and real images.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 PSNR = 10log ) , where  and   are the Original and restored images, respectively.Also   stands for the value of Original image