Abstract

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. The main 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 RHAM methods.

1. 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 [59].

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 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 [2326]. 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.

2. Review of TM Method for TV-Based Image Denoising Problem

In this section, we only review TM method for solving a TV-based 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 thatImage 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 TV-based minimization problem:where 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. The Euler-Lagrange equation corresponding to (2) leads to a nonlinear elliptic PDE with homogeneous Neumann boundary conditions as follows: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:Hence, one obtainsIn order to solve the Euler-Lagrange equation (5), Rudin et al. [9] proposed TM method for solving the following time-dependent nonlinear parabolic PDE

For numerical implementation, let us assume that the domain has been split into cells, where the grid points are located at , , , , where and , denote the time step and iteration number, respectively. We denote the values of at the grid points by and .

Without loss of generality, we can assume that . Then, the nonlinear PDE (6) can be approximated by the following finite difference formula:wherewithand boundary conditions areHence, we obtain TM method, called Algorithm 1, for solving (6), where is the tolerance for stopping criterion, is the maximum number of iterations, is the Noisy image, 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.

(1) for ,
(2) for   to maxit  do
(3)  Compute using (8) for ,
(4)  Compute for ,
(5)  if   or   then
(6)   Stop
(7)  end if
(8) end for

3. Description of Restarted Homotopy Perturbation Method (RHPM)

In this section, we first describe the homotopy perturbation method (HPM) for solving nonlinear Partial Differential Equationwith the boundary conditionswhere 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 asBy using homotopy technique, one can construct a homotopy which satisfiesor, equivalently,where is an embedding parameter and is the initial approximation of (11) which satisfies the boundary conditions. Clearly, we haveThe changing process of from zero to unity is just that of changing from to . This is called deformation and also and are called homotopic in topology. If the embedding parameter 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 , the solution of (11) can be expressed as the sum of an infinite seriesFurther description and convergence for the HPM can be found in [1922, 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 formwith the initial condition .

Method 1 (HPM). We first solve (19) using the traditional HPM. The operators , , and in (13) can be viewed asLet . From (15), we readily obtain a homotopy which satisfiesThen, the solution of (21) can be expressed as in (17). Substituting (17) into (21), one obtainsComparing coefficients of terms with identical powers of , we obtain thatHence, we have Therefore, we obtain the homotopy solution of (21) Letting , 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 as an initial approximation of (19). Let the operators , , and be defined as in (20). Comparing coefficients of terms and in (22), one obtainsThen, the first two term approximation of (18) is defined as follows: From (27), using as an initial approximation, one easily obtainsFrom (29), the new approximation is given by From (27) with an initial approximation , one obtainsHence, the new approximation is given by Continuing this process, it can be seen that converges to the exact solution of (19), which is the same as the result of Method 1, as .

Generally speaking, the computation of 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 . 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 be an initial approximation of a nonlinear equation. First, we compute the first th term approximation in the following way: Following the above computational process with as an initial approximation, a new approximation is computed as follows: By repeating the above process, a new approximation can be computed from the initial approximation , which was obtained at the th step, in the following way:Notice that are coefficients of , respectively. Then, the exact solution can be approximated by ; that is, under some appropriate conditions, (see Example 1). In RHPM, which computes only vectors , every step is called method. Clearly, RHPM introduced in Example 1 is RHPM.

4. Application of RHPM to TV-Based Image Denoising Problem

In this section, we describe an application of RHPM to a TV-based 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) (i.e., and ). In order to handle this difficulty, the following binomial series is used. For ,If in (36), thenLetting in (37), for ,If is large, then, by introducing a small parameter and using (38), can be approximated bySimilarly, if and in (36), then, for or ,If is large, then, by introducing a small parameter and using (40), can be approximated byNote that a small parameter is used to guarantee the convergence of the binomial series when 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 (i.e., ). Since we only use four terms of the binomial series for numerical computation of (45), division by zero does not occur and thus is used for RHPM (see (39) and (41)). Let the operators , , and be defined as follows:where From (15), one obtainswhere Since the solution of (44) is of the form (17), by simple computation From these equations, one easily obtains where The and in (44) can be expressed asSubstituting into (39), one obtainswhere Hence, the first term of the right hand side in (45) can be written asSubstituting into (41), one obtains where Hence, the second term of the right hand side in (45) can be written asFrom (45), (52), and (55), (44) can be expressed asComparing the coefficients of , , and in (56),Applying the inverse operator of to (58), one obtainswhere and . From (60), the partial derivatives of arewhere , , and . From (61), , , and can be expressed asSubstituting (60) and (62) into (59) and applying the inverse operator of to (59), one obtainsNow, update whose new value is set to . By repeating the aforementioned process with the new updated , we obtain RHPM method.

For numerical implementation of RHPM method, let us assume that the domain has been split into cells, where the grid points are located at , , , . Let , where and , refer to the time step and iteration number, respectively. We denote the values of at the grid points by . Without loss of generality, we can assume that . For simplicity of exposition, we define the following notation: Boundary conditions are defined the same as those in TM method. Assume that and has been computed for all and . Then, and discretization of (60) is given byAlso, discretization of (63) is given by where We now compute . Repeating the above process with as an initial approximation, we can obtain RHPM method, called Algorithm 2, for solving (6).

(1) for all and
(2) for   to maxit do
(3)  Set
(4)  Compute using (65) for all and
(5)  Compute using (66) for all and
(6)  
(7)  if   or   then
(8)   Stop
(9)  end if
(10) end for

Since the computation of in RHPM method is so complicated, we also provide RHPM method, called Algorithm 3, whose computational cost is much cheaper than Algorithm 2. In RHPM and RHPM methods, , , , 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 may have large values which cause amplification of some incorrect results. Hence, we have used an intensity range for all test images used in RHPM methods to guarantee numerical stability.

(1) for all and
(2) for   to maxit do
(3)  Set
(4)  Compute using (65) for all and
(5)  
(6)  if   or   then
(7)   Stop
(8)  end if
(9) end for

5. Numerical Experiments

In this section, numerical experiments for several test images are carried out to compare the performance of RHPM and RHPM 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 . For example, adding of Gaussian white noise to the true image can be done as follows: 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 where and are the Original and restored images, respectively. Also stands for the value of Original image at the grid point and is the total number of pixels. It is clear that the larger PSNR is, the better the quality denoising image is.

We have tested the denoising algorithms for several images with an intensity range for TM and RHAM2 methods and for RHPM and RHPM methods, and Noisy images are generated by adding or of Gaussian white noise to true images. For all of the 4 algorithms, parameters and have been used. In TM and RHAM2 methods, parameters and are used. In RHPM and RHPM methods, parameters and are used. Note that the parameters are not optimal and are chosen as the best one by tries.

Tables 14 provide numerical results for TM, RHAM2, RHPM, and RHPM methods corresponding to , , , and , respectively. In Tables 14, column PSNR contains the PSNR values, column Iter contains the number of iterations, and column CPU contains the elapsed CPU time in seconds. Figures 1 to 7 show the images denoised by TM, RHAM2, RHPM, and RHPM methods with for Noisy images with of Gaussian white noise. Note that the denoised images for other values of and of Gaussian white noise are almost the same. The first row images are Original image and Noisy image with of Gaussian white noise and the image denoised by TM, and the second row images are the images denoised by RHAM2, RHPM, and RHPM methods. In Figures 1 to 4, Cameraman, Jet-plane, Lake, and Mandrill images with size are used. In Figure 5, Nodules image with size is used. In Figures 6 and 7, Rice and Caribou images with size are used.

As can be seen in Tables 14, yields the best performance of 4 different cases of for TM, RHPM, and RHPM methods, but yields the best performance for RHAM2 method. It can be also seen that RHPM methods perform much faster than TM and RHAM2 methods. On the other hand, PSNR values are close to one another. Overall, RHPM methods give a little bit higher PSNR values than TM and RHAM2 methods. In particular, RHPM gives about 1 higher PSNR value than TM and RHAM2 methods for Rice and Caribou images which have small sizes compared to the other images. According to numerical experiments, this phenomenon occurs because of the shape of image, not the size of image. That is, PSNR values are image-dependent, not size-dependent.

Observe that RHPM computing 2 terms , every iteration performs more efficiently than RHPM computing 3 terms of , , every iteration. The reason why RHPM performs worse than RHPM is as follows: the computational cost for computing is expensive as compared with , and its computation involves several approximation formulas, which use only a few terms of the infinite series, due to the difficulty in handling highly nonlinear terms. To this end, convergence rate of RHPM is slower than RHPM.

6. Conclusion

In this paper, we proposed a restarted homotopy perturbation method (RHPM) for solving a nonlinear PDE problem, and we presented an application of RHPM to TV-based image denoising problem. We also provided numerical methods RHPM and RHPM which use binomial series techniques to settle the main difficulty in handling nonlinear terms. Numerical experiments show that RHPM and RHPM are effective, feasible, and reliable for TV-based image denoising problem as compared with TM and RHAM2 methods. Notice that RHPM is very efficient compared to other methods since it has smaller iteration number and less CPU time than others. Also RHPM is much simpler to implement than RHPM. So, we recommend the use of RHPM with and over RHPM for TV-based image denoising problem. This work only studied the effectiveness, feasibility, and reliability of RHPM. So, future work will include comparison results between the proposed RHPM methods and other good methods in the literature.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors would like to thank the anonymous referees for their useful comments and for pointing out many related references, which greatly improved this paper. This work was supported by the research Grant of the Agency for Defense Development and the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by Ministry of Education (NRF-2013R1A1A2005722).