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 [5–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 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–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.
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.
| ||||||||||||||||
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 [19–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 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