HPM-Based Dynamic Wavelet Transform and Its Application in Image Denoising

Abstract

Wavelet-based multiscale interpolation operator is often employed to construct the adaptive numerical method for PDEs, in which the computational complexity of the wavelet transform is one of the main factors affecting the algorithm efficiency. As the wavelet transform just acts as the detector of the characteristic points in the interpolation operator, the multiscale wavelet interpolation operator can be viewed as a nonlinear problem. Based on this assumption, we construct an approximate dynamic interpolation operator with the homotopy perturbation method (HPM), which decreases the computational complexity of the wavelet transform appearing in the wavelet interpolation operator from 𝑂((1/3)4^2𝐽−1) to 𝑂(4^𝐽), where 𝐽 is the amount of the wavelet scales. Then an adaptive algorithm solving the Perona-Malik model on image denoising is constructed with the HPM-based interpolation operator. Last, the quasi-Shannon wavelet is employed to design the experiments on the medical image and some artificial images denoising. The experiment results show that the simplified wavelet interpolation operator based on HPM possesses the adaptability and nonsensitivity to the time step, which is helpful to improve the algorithm efficiency. This illustrates that the HPM-based wavelet interpolation operator is an effective tool to solve the problems in image processing.

1. Introduction

In recent years, the partial differential equation (PDE) method has been widely used in image denoising, especially in medical images and remote sensing images. Various variational models have been constructed to satisfy different requirements on image processing such as denoising and segmentation. One of the excellent models is the Perona-Malik model [1] which is widely used in image denoising. Perona-Malik model is expressed as a nonlinear 2-dimension partial differential equation, which overcomes the drawback of the scale-space technique introduced by Witkin that involves generating coarser resolution images by convolving the original image with a Gaussian kernel. In this approach, a new definition of scale space was suggested, and a class of algorithms was introduced; then, the accurate locations of the “semantically meaningful” edges at coarse scales using a diffusion process can be obtained; that is, a high quality edge detector which successfully exploits global information was obtained with this new method.

There is no doubt that Perona-Malik model with anisotropic diffusion is a very useful tool in many tasks of computer vision [2, 3]. As a nonlinear 2D partial differential equation, the bottleneck of the model is the computational efficiency. It is difficult to find an exact analytical solution of the nonlinear PDEs. Up to now, the infinite difference method is the main numerical algorithm for Perona-Malik model [4], which can bring artifact into the images due to the nonsmoothness of the basis function of the infinite difference method [5, 6]. Multilevel wavelet numerical method for the nonlinear PDEs has been proposed over ten years, which can take full advantage of the adaptability of the wavelet analysis. The artifacts in image can be eliminated with the wavelet numerical algorithm instead of the infinite difference method, as wavelet basis function possesses many excellent properties such as smoothness and compact support. But the support range of wavelet function is much wider than the basis function in the infinite difference method. This leads to a lower computational efficiency of wavelet transform in solving 2D nonlinear PDEs.

The multilevel wavelet transform usually appears in the wavelet interpolation operator [7, 8]. In regard to the evolvement of PDEs, the wavelet interpolation operator should be dynamic with the evolvement. That means each iteration step of the evolvement corresponds to one change of the interpolation operator. Therefore, the efficiency of the algorithm depends mainly on the wavelet transform. For larger images such as the remote sensing images and medical images, the computational complexity is higher. As a matter of fact, we could view the dynamic process of the wavelet interpolation operator as a nonlinear problem, and the effective approximate method can be employed to decrease the calculation amount. Among numerous linearization methods, the homotopy analysis method (HAM) [9–11] proposed by Shijun Liao is a simple and efficient method, which has been used widely in various applications [12–14]. The HAM enjoys great freedom in choosing the auxiliary linear operator and initial guess. In particular, it provides a convenient way to guarantee the convergence of solution series. In most cases, the solution series given by the HAM converge quickly. Based on HAM, Jihuan He proposed a more simple method with perturbation series instead of Tayor series, which is called homotopy perturbation method (HPM). The better improvement is adding an auxiliary parameter into the homotopy equation, which is helpful to eliminate the secular term in the perturbation solution. This can improve the rate of convergence greatly. Unlike analytical perturbation methods, HPM does not depend on a small parameter which is difficult to find [15–18] and has been widely used in solving various nonlinear problems [19–27]. The purpose of this research is to construct a dynamic wavelet interpolation operator with HPM and apply it to solve the Perona-Malik model which is a classical image denoising model.

2. Anisotropic Diffusion Model and Its Discretization Scheme

2.1. Perona-Malik Model

The anisotropic diffusion image denoising model was proposed by Perona and Malik, which has been widely used in various image processing fields. So, it is usually called Perona-Malik model, which can be expressed as the nonlinear partial differential equations: 𝜕𝑢𝑥,𝑦,𝑡(𝑐),𝑢,𝜕𝑡=div(|∇𝑢|)∇𝑢𝑥,𝑦,0=𝑓𝑥,𝑦 (1) where (𝑥,𝑦) denotes pixel position, 𝑡 is the time parameter, 𝑓(𝑥,𝑦) is the 2D image being processed, 𝑢(𝑥,𝑦,𝑡) is the image after diffusion processing, and 𝑢(𝑥,𝑦,0) is the initial value. div denotes the divergence operator, and the ∇𝑢 denotes the gradient operator; 𝑐(|∇𝑢|) denotes the diffusion coefficient, which is a nonnegative decreasing function of the image gradient modulus. It is usually taken as 𝑐=1(|∇𝑢|)1+(∇𝑢/𝑘)² (2) or 𝑐(|∇𝑢|)=exp−∇𝑢𝑘², (3) where 𝑘 is a constant.

2.2. Wavelet Numerical Discretization Scheme of Perona-Malik Model

Let the definition domain of the image be (𝑥_min,𝑥_max)×(𝑦_min,𝑦_max); the discretization points can be defined as (𝑥^𝑗_𝑘1,𝑦^𝑗_𝑘2), where 𝑗 is a scale parameter and 𝑘₁ and 𝑘₂ are position parameters. So, 𝑥^𝑗_𝑘1=𝑥_min+𝑘₁𝑥_max−𝑥_min2^𝑗,𝑦^𝑗_𝑘1=𝑦_min+𝑘₂𝑦_max−𝑦_min2^𝑗,𝑗,𝑘₁,𝑘₂∈ℤ. (4)

In addition, 𝑤_𝑘^{𝑗(𝑚,𝑛)}_1,𝑘2(𝑥,𝑦) denotes the multiscale wavelet function and the corresponding 𝑚th and 𝑛th derivatives with respect to 𝑥 and 𝑦, respectively. The level set function  𝜙(𝑥,𝑦,𝑡) and the corresponding derivative function can be descretized as follows: 𝑢^{𝐽(𝑚,𝑛)}=𝑥,𝑦,𝑡_{1𝑘011=0𝑘02=0}𝑢𝑥⁰_𝑘01,𝑦⁰_𝑘02𝑤_𝑘^0(𝑚,𝑛)_01,𝑘02+𝑥,𝑦_{𝐽−12𝑗=0𝑗−1𝑘112=0𝑗−1𝑘12=0}𝛼¹_{𝑗,𝑘11,𝑘12}(𝑡)𝑤^{𝑗+1(𝑚,𝑛)}_{2𝑘11+1,2𝑘12}𝑥,𝑦+𝛼²_{𝑗,𝑘11,𝑘12}(𝑡)𝑤^{𝑗+1(𝑚,𝑛)}_{2𝑘11,2𝑘12+1}𝑥,𝑦+𝛼³_{𝑗,𝑘11,𝑘12}(𝑡)𝑤^{𝑗+1(𝑚,𝑛)}_{2𝑘11+1,2𝑘12+1}𝑥,𝑦, (5) where 𝑗 and 𝐽 are constants, which denote the wavelet scale number and the maximum of the scale number, respectively. 𝛼¹_{𝑗,𝑘11,𝑘12}, 𝛼²_{𝑗,𝑘11,𝑘12}, and 𝛼³_{𝑗,𝑘11,𝑘12}are the wavelet coefficients at the point (𝑥^𝑗_𝑘1,𝑦^𝑗_𝑘2). According to the interpolation wavelet transform theory, the wavelet coefficients can be written as 𝛼¹_{𝑗,𝑘1,𝑘2}=𝑢𝑥_{𝑗+1,2𝑘1+1},𝑦_{𝑗+1,2𝑘2}−𝐼_𝑗𝑢𝑥_{𝑗+1,2𝑘1+1},𝑦_{𝑗+1,2𝑘2},𝛼²_{𝑗,𝑘1,𝑘2}=𝑢𝑥_{𝑗+1,2𝑘1},𝑦_{𝑗+1,2𝑘2+1}−𝐼_𝑗𝑢𝑥_{𝑗+1,2𝑘1},𝑦_{𝑗+1,2𝑘2+1},𝛼³_{𝑗,𝑘1,𝑘2}=𝑢𝑥_{𝑗+1,2𝑘1+1},𝑦_{𝑗+1,2𝑘2+1}−𝐼_𝑗𝑢𝑥_{𝑗+1,2𝑘1+1},𝑦_{𝑗+1,2𝑘2+1}, (6) where 𝐼_𝑗  denotes the multilevel interpolation operator. In order to obtain the multilevel interpolation operator, it is necessary to express the wavelet coefficients  𝛼¹_{𝑗,𝑘1,𝑘2}, 𝛼²_{𝑗,𝑘1,𝑘2}, and  𝛼³_{𝑗,𝑘1,𝑘2} as a weighted sum of 𝑢 in all of the collocation points in the 𝐽-level. Therefore, we should give the definition of the restriction operator as follows: 𝑅_𝑘^{𝑙,𝑙,𝑗,𝑗}_{1,𝑘2,𝑚1,𝑚2}=1,𝑥^𝑙_𝑘1=𝑥^𝑗_𝑚1,𝑦^𝑙_𝑘2=𝑦^𝑗_𝑚20,otherwise. (7) Using the restriction operator, 𝑢(𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2),𝑢(𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1),  and   𝑢(𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1) can be rewritten as 𝑢𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2=^2𝐽_𝑛1²₌₀^𝐽_𝑛2=0𝑅^{𝑗+1,𝑗+1,𝐽,𝐽}_{2𝑘1+1,2𝑘2,𝑛1,𝑛2}𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2,𝑢𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1=^2𝐽_𝑛1²₌₀^𝐽_𝑛2=0𝑅^{𝑗+1,𝑗+1,𝐽,𝐽}_{2𝑘1,2𝑘2+1,𝑛1,𝑛2}𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2,𝑢𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1=^2𝐽_𝑛1²₌₀^𝐽_𝑛2=0𝑅^{𝑗+1,𝑗+1,𝐽,𝐽}_{2𝑘1+1,2𝑘2+1,𝑛1,𝑛2}𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2. (8) Introducing the extension operators 𝐶1, 𝐶2, and 𝐶3 and substituting (8) into (6), the wavelet coefficients can be rewritten as 𝛼¹_{𝑗,𝑘1,𝑘2}=^2𝐽_𝑛1²₌₀^𝐽_𝑛2=0𝑅^{𝑗+1,𝑗+1,𝐽,𝐽}_{2𝑘1+1,2𝑘2,𝑛1,𝑛2}𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2−^2𝐽_𝑛1²₌₀^𝐽_𝑛2²₌₀^𝑗0_𝑘01²₌₀^𝑗0_𝑘02=0𝑅^𝑗0,𝑗0_𝑘^,𝐽,𝐽_{01,𝑘02,𝑛1,𝑛2}𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2×𝜙^𝑗0_{𝑘01,𝑘02}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+^𝑗−1_𝑗1=𝑗0^2𝐽_𝑛1²₌₀^𝐽_𝑛2²₌₀^𝑗1_𝑘11²₌₀^𝑗1_𝑘12=0𝐶1^𝑗1,𝑗1_𝑘^,𝐽,𝐽_{11,𝑘12,𝑛1,𝑛2}𝜙^𝑗1+1_{2𝑘11+1,2𝑘12}×𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2+𝐶2^𝑗1,𝑗1_𝑘^,𝐽,𝐽_{11,𝑘12,𝑛1,𝑛2}𝜙^𝑗1+1_{2𝑘11,2𝑘12+1}×𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2+𝐶3^𝑗1,𝑗1_𝑘^,𝐽,𝐽_{11,𝑘12,𝑛1,𝑛2}×𝜙^𝑗1+1_{2𝑘11+1,2𝑘12+1}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2×𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2=^2𝐽_𝑛1²₌₀^𝐽_𝑛2=0𝐶1_𝑘^{𝑗,𝑗,𝐽,𝐽}_{1,𝑘2,𝑛1,𝑛2}𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2. (9)𝛼²_{𝑗,𝑘1,𝑘2}and 𝛼³_{𝑗,𝑘1,𝑘2} are similar to 𝛼¹_{𝑗,𝑘1,𝑘2}. From the above equation, the extension operator can be obtained as 𝐶1_𝑘^{𝑗,𝑗,𝐽,𝐽}_{1,𝑘2,𝑛1,𝑛2}=𝑅^{𝑗+1,𝑗+1,𝐽,𝐽}_{2𝑘1+1,2𝑘2,𝑛1,𝑛2}−^2𝑗0_𝑘01²₌₀^𝑗0_𝑘02=0𝑅^𝑗0,𝑗0_𝑘^,𝐽,𝐽_{01,𝑘02,𝑛1,𝑛2}𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2×𝜙^𝑗0_{𝑘01,𝑘02}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+^𝑗−1_𝑗1=𝑗0^2𝐽_𝑛2²₌₀^𝑗1_𝑘11=0𝐶1^𝑗1,𝑗1_𝑘^,𝐽,𝐽_{11,𝑘12,𝑛1,𝑛2}𝜙^𝑗1+1_{2𝑘11+1,2𝑘12}×𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2+𝐶2^𝑗1,𝑗1_𝑘^,𝐽,𝐽_{11,𝑘12,𝑛1,𝑛2}𝜙^𝑗1+1_{2𝑘11,2𝑘12+1}×𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2+𝐶3^𝑗1,𝑗1_𝑘^,𝐽,𝐽_{11,𝑘12,𝑛1,𝑛2}𝜙^𝑗1+1_{2𝑘11+1,2𝑘12+1}×𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2×𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2. (10)𝐶2 and 𝐶3 can be obtained with the same method. Therefore, the calculation time complexity of the wavelet transform coefficient 𝛼¹_{𝑗,𝑘11,𝑘12}, 𝛼²_{𝑗,𝑘11,𝑘12}, and 𝛼³_{𝑗,𝑘11,𝑘22}is 𝑂((1/3)4^2𝐽−1).

Substituting 𝛼¹_{𝑗,𝑘11,𝑘12}, 𝛼²_{𝑗,𝑘11,𝑘12}, and 𝛼³_{𝑗,𝑘11,𝑘22} and 𝐶1, 𝐶2 and 𝐶3 into (5), the multilevel wavelet interpolation operator can be obtained as 𝐼_𝑛1,𝑛2=𝑥,𝑦_{2𝑗0𝑘012=0𝑗0𝑘02=0}𝑅^𝑗0,𝑗0_𝑘^,𝐽,𝐽_{01,𝑘02,𝑛1,𝑛2}𝜙^𝑗0_{𝑘01,𝑘02}+𝑥,𝑦_{𝐽−1𝑗=𝑗02𝑗𝑘12=0𝑗𝑘2=0}𝐶1_𝑘^{𝑗,𝑗,𝐽,𝐽}_{1,𝑘2,𝑛1,𝑛2}𝜙^𝑗+1_{2𝑘1+1,2𝑘2}×𝑥,𝑦𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2+𝐶2_𝑘^{𝑗,𝑗,𝐽,𝐽}_{1,𝑘2,𝑛1,𝑛2}𝜙^𝑗+1_{2𝑘1,2𝑘2+1}×𝑥,𝑦𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2+𝐶3_𝑘^{𝑗,𝑗,𝐽,𝐽}_{1,𝑘2,𝑛1,𝑛2}𝜙^𝑗+1_{2𝑘1+1,2𝑘2+1}×𝑥,𝑦𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2. (11) Then, (5) can be rewritten as 𝑢^{𝐽(𝑚,𝑛)}𝑥,𝑦,𝑡=^2𝐽_𝑛1^2𝐽_𝑛2𝐼_𝑛1,𝑛2𝑥,𝑦𝑢𝑥^𝐽_𝑛1,𝑦^𝐽_𝑛2. (12) Substituting (12) into (1), the multilevel wavelet discretization scheme of Perona-Malik model can be obtained.

3. HPM-Based Wavelet Interpolation Scheme

The purpose of constructing the multilevel wavelet collocation method is to decrease the amount of the collocation points and then improve the efficiency of the algorithm. But the efficiency will be eliminated if the computation complexity of the multilevel wavelet interpolation operator is too high. It is easy to understand that the interpolation wavelet coefficient is the error between the interpolation result and the exact result at the same collocation point. And so the wavelet coefficient must be the function of the parameter 𝑡. In other words, the wavelet coefficient should vary with the time parameter 𝑡. Then, the interpolation operator can be viewed as a nonlinear problem. HPM is an efficient and effective tool to solve nonlinear problem. Aiming to improve the efficiency of the multilevel wavelet interpolation operator, HPM would be employed to construct a novel interpolation operator in this section.

For convenience, 𝑢 and its derivative in (1) should be rewritten as 𝜕𝑢𝜕𝑡=𝐹𝑡,𝑥,𝑦,𝑢,𝜕𝑢,𝜕𝑥𝜕𝑢,𝜕𝜕𝑦²𝑢𝜕𝑥²,𝜕²𝑢,𝜕𝜕𝑥𝜕𝑦²𝑢𝜕𝑦²,(),𝑢𝑡>0𝑥,𝑦,0=𝑢₀𝑥,𝑦,𝑑𝑢^𝐽𝑥,𝑦,𝑡𝑑𝑡=𝐹𝑡,𝑥,𝑦,𝑢^𝐽𝑥,𝑦,𝑡,𝑢^𝐽(1,0),𝑢𝑥,𝑦,𝑡_𝐽(0,1)𝑥,𝑦,𝑡,𝑢^𝐽(2,0),𝑢𝑥,𝑦,𝑡_𝐽(1,1)𝑥,𝑦,𝑡,𝑢^𝐽(0,2)𝑥,𝑦,𝑡, (13) respectively, where 𝑢^𝐽=𝑥,𝑦,𝑡_{1𝑘011=0𝑘02=0}𝑢𝑥⁰_𝑘01,𝑦⁰_𝑘02∅⁰_{𝑘01,𝑘02}+𝑥,𝑦_{𝐽−12𝑗=0𝑗−1𝑘112=0𝑗−1𝑘12=0}𝛼¹_{𝑗,𝑘11,𝑘12}∅^𝑗+1_{2𝑘11+1,2𝑘12}𝑥,𝑦+𝛼²_{𝑗,𝑘11,𝑘12}∅^𝑗+1_{2𝑘11,2𝑘12+1}𝑥,𝑦+𝛼³_{𝑗,𝑘11,𝑘12}∅^𝑗+1_{2𝑘11+1,2𝑘12+1}𝑥,𝑦. (14) The value of 𝑢^𝐽(𝑥,𝑦,𝑡_𝑛) at 𝑡_𝑛 is denoted as 𝑢_𝑛, and 𝐹𝑡_𝑛,𝑥,𝑦,𝑢^𝐽𝑥,𝑦,𝑡_𝑛,𝑢^𝐽(1,0)𝑥,𝑦,𝑡_𝑛,𝑢^𝐽(0,1)𝑥,𝑦,𝑡_𝑛,𝑢^𝐽(2,0)𝑥,𝑦,𝑡_𝑛,𝑢^𝐽(1,1)𝑥,𝑦,𝑡_𝑛,𝑢^𝐽(0,2)𝑥,𝑦,𝑡_𝑛 (15) is denoted as 𝐹_𝑛. And then a linear homotopy function can be constructed as 𝑢^𝐽𝑥,𝑦,𝑡=(1−𝜀)𝐹_𝑛+𝜀𝐹_𝑛+1. (16) It is easy to identify the homotopy parameter as 𝜀(𝑡)=𝑡−𝑡_𝑛𝑡_𝑛+1−𝑡_𝑛,𝑡∈𝑡_𝑛,𝑡_𝑛+1∴𝜀∈[0,1]. (17) According to the perturbation theory, the solution of (16) can be expressed as the power series expansion of 𝜀𝑢^𝐽=𝑢^𝐽₀+𝜀𝑢^𝐽₁+𝜀²𝑢^𝐽₂+⋯. (18) Substituting (18) into (16) and rearranging based on powers of  𝜀-terms, we have 𝜀⁰∶𝑢^𝐽₀=𝐹_𝑛𝜀¹∶𝑢^𝐽₁=𝐹_𝑛+1−𝐹_𝑛⋮ (19) According to HPM, we obtain the wavelet coefficient 𝛼¹_{𝑗,𝑘1,𝑘2}(𝑡_𝑛+1),    𝛼²_{𝑗,𝑘1,𝑘2}(𝑡_𝑛+1),   and  𝛼³_{𝑗,𝑘1,𝑘2}(𝑡_𝑛+1) at 𝑡_𝑛  as follows: 𝛼¹_{𝑗,𝑘1,𝑘2}=𝑢𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2−𝐼_𝑗𝑢𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2=𝑢𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2−¹_𝑘01¹_=0𝑘02=0𝑢𝑥⁰_𝑘01,𝑦⁰_𝑘02∅⁰_{𝑘01,𝑘02}𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+^𝑗−1_𝑗1²₌₀^𝑗1_𝑘11²₌₀^𝑗1_𝑘12=0𝛼¹_{𝑗1,𝑘11,𝑘12}∅^𝑗1+1_{2𝑘11+1,2𝑘12}×𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+𝛼²_{𝑗1,𝑘11,𝑘12}∅^𝑗1+1_{2𝑘11,2𝑘12+1}×𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+𝛼³_{𝑗1,𝑘11,𝑘12}∅^𝑗1+1_{2𝑘11+1,2𝑘12+1}×𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2,𝛼²_{𝑗,𝑘1,𝑘2}=𝑢𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1−𝐼_𝑗𝑢𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1=𝑢𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1−¹_𝑘01¹_=0𝑘02=0𝑢𝑥⁰_𝑘01,𝑦⁰_𝑘02∅⁰_{𝑘01,𝑘02}×𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1+^𝑗−1_𝑗1²₌₀^𝑗1_𝑘11²₌₀^𝑗1_𝑘12=0𝛼¹_{𝑗1,𝑘11,𝑘12}∅^𝑗1+1_{2𝑘11+1,2𝑘12}×𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1+𝛼²_{𝑗1,𝑘11,𝑘12}∅^𝑗1+1_{2𝑘11,2𝑘12+1}×𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1+𝛼³_{𝑗1,𝑘11,𝑘12}∅^𝑗1+1_{2𝑘11+1,2𝑘12+1}×𝑥^𝑗+1_2𝑘1,𝑦^𝑗+1_2𝑘2+1,𝛼³_{𝑗,𝑘1,𝑘2}=𝑢𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1−𝐼_𝑗𝑢𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1=𝑢𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1−¹_𝑘01¹_=0𝑘02=0𝑢𝑥⁰_𝑘01,𝑦⁰_𝑘02∅⁰_{𝑘01,𝑘02}×𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1+^𝑗−1_𝑗1²₌₀^𝑗1_𝑘11²₌₀^𝑗1_𝑘12=0𝛼¹_{𝑗1,𝑘11,𝑘12}∅^𝑗1+1_{2𝑘11+1,2𝑘12}×𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1+𝛼²_{𝑗1,𝑘11,𝑘12}∅^𝑗1+1_{2𝑘11,2𝑘12+1}×𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1+𝛼³_{𝑗1,𝑘11,𝑘12}∅^𝑗1+1_{2𝑘11+1,2𝑘12+1}×𝑥^𝑗+1_2𝑘1+1,𝑦^𝑗+1_2𝑘2+1. (20) Obviously, the calculation time complexity of the wavelet transform coefficient 𝛼¹_{𝑗,𝑘1,𝑘2}, 𝛼²_{𝑗,𝑘1,𝑘2}, and 𝛼³_{𝑗,𝑘1,𝑘2} is 𝑂(4^𝐽), which is decreased greatly than it in (6) which is 𝑂((1/3)4^2𝐽−1).

Substituting the wavelet transform efficient (20) into (14), we obtain 𝑢^𝐽𝑥,𝑦,𝑡_𝑛+1=𝑢^𝐽𝑥,𝑦,𝑡_𝑛+Δ𝑡2𝐹𝑡_𝑛,𝑥,𝑦,𝑢^𝐽𝑥,𝑦,𝑡_𝑛,𝑢^𝐽(1,0)𝑥,𝑦,𝑡_𝑛,𝑢^𝐽(0,1)𝑥,𝑦,𝑡_𝑛,𝑢^𝐽(2,0)𝑥,𝑦,𝑡_𝑛,𝑢^𝐽(1,1)𝑥,𝑦,𝑡_𝑛,𝑢^𝐽(0,2)𝑥,𝑦,𝑡_𝑛+𝐹𝑡_𝑛+1,𝑥,𝑦,𝑢^𝐽₀𝑥,𝑦,𝑡_𝑛+1,𝑢₀^𝐽(1,0)𝑥,𝑦,𝑡_𝑛+1,𝑢₀^𝐽(0,1)𝑥,𝑦,𝑡_𝑛+1,𝑢₀^𝐽(2,0)𝑥,𝑦,𝑡_𝑛+1,𝑢₀^𝐽(1,1)𝑥,𝑦,𝑡_𝑛+1,𝑢₀^𝐽(0,2)𝑥,𝑦,𝑡_𝑛+1 (21) and the derivative function 𝑢^{𝐽(𝑚,𝑛)}𝑥,𝑦=¹_𝑘01¹_=0𝑘02=0𝑢𝑥⁰_𝑘01,𝑦⁰_𝑘02∅_𝑘^0(𝑚,𝑛)_01,𝑘02+𝑥,𝑦_{𝐽−12𝑗=0𝑗−1𝑘112=0𝑗−1𝑘12=0}𝛼¹_{𝑗,𝑘11,𝑘12}∅^{𝑗+1(𝑚,𝑛)}_{2𝑘11+1,2𝑘12}𝑥,𝑦+𝛼²_{𝑗,𝑘11,𝑘12}∅^{𝑗+1(𝑚,𝑛)}_{2𝑘11,2𝑘12+1}𝑥,𝑦+𝛼³_{𝑗,𝑘11,𝑘12}∅^{𝑗+1(𝑚,𝑛)}_{2𝑘11+1,2𝑘12+1}𝑥,𝑦. (22) Obviously, the computation complexity is decreased greatly comparing with (12).

4. Numerical Experiences and Discussion

In this section, we take some images as examples to illustrate the efficiency of the novel wavelet interpolation operator based on HPM. The quasi-Shannon wavelet function [28–32] is employed to construct the novel interpolation operator, which is defined as 𝛿_Δ𝜎=sin(𝜋𝑥/Δ)−𝑥𝜋𝑥/Δexp²2𝜎², (23) where Δis the discrete space and 𝜎=𝑟Δ (𝑟 is a random constant) denotes the support range.

4.1. Adaptability of HPM-Based Wavelet Interpolation Operator

The first experiment image is a brain picture, which would be used to test the effectivness of the novel wavelet interpolation operator based on HPM. The result is shown in Figure 1.

Figure 1 

Image denoising results and the corresponding collocation points distribution at different time parameters 𝑡.

The image becomes more and more smooth with the image denoising processing, and the amount of the collocation points becomes less and less. Obviously, it is helpful to improve the efficiency of the algorithm. This shows that the operator constructed in this paper possesses the adaptability property.

4.2. Nonsensitivity to Time Step of the HPM-Based Wavelet Operator

The second image is solid circle with random noise. We try to illustrate the effectivness of the novel interpolation operator based on HPM comparing with the general operator. The results are shown in Figure 2.

Figure 2 

Comparison about the denoising results between the general wavelet interpolation operator and the HPM-base multilevel wavelet interpolation operator.

Using the novel HPM-based multilevel wavelet interpolation operator, the noise in images can be denoised rapidly. The denoising processing is slower with the original wavelet interpolation operator, and the image processing is less effective than that with the HPM-based multilevel wavelet interpolation operator on the same time step interval. One reason is because quasi-Shannon has no compact support property, and another is because the computational precision is so high that the noise is viewed as a part of the image. So, noises cannot be denoised rapidly. In addition, [32] has proved that HPM is not sensitive to the time step and is suitable for nonlinear equations.

4.3. Comparison between the Wavelet Interpolation Operator and the Finite Difference Operator

The third experiment image is a Metarhizium anisopliae picture, which would be used to test the effective of the novel wavelet interpolation operator based on HPM comparing with the finite difference operator. The result is shown in Figure 3. Although the time step of the difference method (𝜏=0.00003) is smaller than the wavelet method (𝜏=0.00005), the artifacts appeared in the denoising results with the difference method (Figures 3(b2) and 3(b3)), which is restricted by the wavelet interpolation method. It is well known that the basic function of the finite difference operator does not have the property of smoothness, and so taking smaller time step is the only way to restrict the artifacts in processed images. The wavelet function possesses merits of smoothness and compact support, which are helpful to avoid the artifacts. It is obvious that this can be used to improve the efficiency of the algorithm.

(a) HPM-based wavelet interpolation operator

(b) Finite difference operator

Figure 3 

The comparison about the smoothing result applying HAM.

5. Conclusions

We have introduced HPM into the construction of the dynamic wavelet interpolation operator. The experiments illustrate it is a better idea to consider the wavelet transform as a nonlinear processing. HPM-based wavelet interpolation operator possesses the adaptability and nonsensitivity to the time step. In addition, comparing to the difference method which is widely used in solving P-M model in recent years, the wavelet numerical method can restrict the artifacts appearing in the processed images. Obviously, all of these can be used to improve the efficiency of the algorithm. As a preprocessing step, it makes the thinning and linking of the edges unnecessary, it preserves the edge junctions, and it does not require a complicated comparison of images at different scales since shape and position are preserved at every single scale. We believe that this will prove useful in many tasks of image processing.

Acknowledgments

This work is supported by the National Key Technologies R & D Program of China under Grant no. 2012BAD35B02 and the National Natural Science Foundation of China under Grant no. 41171337.