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

Authors: Shu-Li Mei
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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 fromO((1/3)42J−1) to O(4J), where J 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.


---

## Body

## 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:(1)∂u(x,y,t)∂t=div(c(|∇u|)∇u),u(x,y,0)=f(x,y),
where (x,y) denotes pixel position, t is the time parameter, f(x,y) is the 2D image being processed, u(x,y,t) is the image after diffusion processing, and u(x,y,0) is the initial value. div denotes the divergence operator, and the ∇u denotes the gradient operator; c(|∇u|) denotes the diffusion coefficient, which is a nonnegative decreasing function of the image gradient modulus. It is usually taken as
(2)c(|∇u|)=11+(∇u/k)2
or
(3)c(|∇u|)=exp[-(∇uk)2],
where k is a constant.


### 2.2. Wavelet Numerical Discretization Scheme of Perona-Malik Model

Let the definition domain of the image be(xmin,xmax)×(ymin,ymax); the discretization points can be defined as (xk1j,yk2j), where j is a scale parameter and k1 and k2 are position parameters. So,
(4)xk1j=xmin+k1xmax-xmin2j,yk1j=ymin+k2ymax-ymin2j,j,k1,k2∈ℤ.In addition,wk1,k2j(m,n)(x,y) denotes the multiscale wavelet function and the corresponding mth and nth derivatives with respect to x and y, respectively. The level set function  ϕ(x,y,t) and the corresponding derivative function can be descretized as follows:
(5)uJ(m,n)(x,y,t)=∑k01=01∑k02=01u(xk010,yk020)wk01,k020(m,n)(x,y)+∑j=0J-1∑k11=02j-1∑k12=02j-1[αj,k11,k121(t)w2k11+1,2k12j+1(m,n)(x,y)+αj,k11,k122(t)w2k11,2k12+1j+1(m,n)(x,y)+αj,k11,k123(t)w2k11+1,2k12+1j+1(m,n)(x,y)],
where j and J are constants, which denote the wavelet scale number and the maximum of the scale number, respectively.  αj,k11,k121, αj,k11,k122,  and αj,k11,k123are the wavelet coefficients at the point (xk1j,yk2j). According to the interpolation wavelet transform theory, the wavelet coefficients can be written as
(6)αj,k1,k21=u(xj+1,2k1+1,yj+1,2k2)-Iju(xj+1,2k1+1,yj+1,2k2),αj,k1,k22=u(xj+1,2k1,yj+1,2k2+1)-Iju(xj+1,2k1,yj+1,2k2+1),αj,k1,k23=u(xj+1,2k1+1,yj+1,2k2+1)-Iju(xj+1,2k1+1,yj+1,2k2+1),
where Ij  denotes the multilevel interpolation operator. In order to obtain the multilevel interpolation operator, it is necessary to express the wavelet coefficients  αj,k1,k21, αj,k1,k22, and  αj,k1,k23 as a weighted sum of u in all of the collocation points in the J-level. Therefore, we should give the definition of the restriction operator as follows:
(7)Rk1,k2,m1,m2l,l,j,j={1,xk1l=xm1j,yk2l=ym2j0,otherwise.
Using the restriction operator, u(x2k1+1j+1,y2k2j+1),u(x2k1j+1,y2k2+1j+1),  and   u(x2k1+1j+1,y2k2+1j+1) can be rewritten as
(8)u(x2k1+1j+1,y2k2j+1)=∑n1=02J∑n2=02JR2k1+1,2k2,n1,n2j+1,j+1,J,Ju(xn1J,yn2J),u(x2k1j+1,y2k2+1j+1)=∑n1=02J∑n2=02JR2k1,2k2+1,n1,n2j+1,j+1,J,Ju(xn1J,yn2J),u(x2k1+1j+1,y2k2+1j+1)=∑n1=02J∑n2=02JR2k1+1,2k2+1,n1,n2j+1,j+1,J,Ju(xn1J,yn2J).
Introducing the extension operators C1, C2, and C3 and substituting (8) into (6), the wavelet coefficients can be rewritten as
(9)αj,k1,k21=∑n1=02J∑n2=02JR2k1+1,2k2,n1,n2j+1,j+1,J,Ju(xn1J,yn2J)-[∑n1=02J∑n2=02J∑k01=02j0∑k02=02j0Rk01,k02,n1,n2j0,j0,J,Ju(xn1J,yn2J)×ϕk01,k02j0(x2k1+1j+1,y2k2j+1)+∑j1=j0j-1∑n1=02J∑n2=02J∑k11=02j1∑k12=02j1(C1k11,k12,n1,n2j1,j1,J,Jϕ2k11+1,2k12j1+1×(x2k1+1j+1,y2k2j+1)u(xn1J,yn2J)+C2k11,k12,n1,n2j1,j1,J,Jϕ2k11,2k12+1j1+1×(x2k1+1j+1,y2k2j+1)u(xn1J,yn2J)+C3k11,k12,n1,n2j1,j1,J,J×ϕ2k11+1,2k12+1j1+1(x2k1+1j+1,y2k2j+1)×u(xn1J,yn2J))∑k02=02j0Rk01,k02,n1,n2j0,j0,J,J]=∑n1=02J∑n2=02JC1k1,k2,n1,n2j,j,J,Ju(xn1J,yn2J).αj,k1,k22and αj,k1,k23 are similar to αj,k1,k21. From the above equation, the extension operator can be obtained as
(10)C1k1,k2,n1,n2j,j,J,J=R2k1+1,2k2,n1,n2j+1,j+1,J,J-[∑k01=02j0∑k02=02j0Rk01,k02,n1,n2j0,j0,J,Ju(xn1J,yn2J)×ϕk01,k02j0(x2k1+1j+1,y2k2j+1)+∑j1=j0j-1∑n2=02J∑k11=02j1(C1k11,k12,n1,n2j1,j1,J,Jϕ2k11+1,2k12j1+1×(x2k1+1j+1,y2k2j+1)u(xn1J,yn2J)+C2k11,k12,n1,n2j1,j1,J,Jϕ2k11,2k12+1j1+1×(x2k1+1j+1,y2k2j+1)u(xn1J,yn2J)+C3k11,k12,n1,n2j1,j1,J,Jϕ2k11+1,2k12+1j1+1×(x2k1+1j+1,y2k2j+1)×u(xn1J,yn2J))-∑k01=02j0∑k02=02j0Rk01,k02,n1,n2j0,j0,J,J].C2 and C3 can be obtained with the same method. Therefore, the calculation time complexity of the wavelet transform coefficient αj,k11,k121, αj,k11,k122, and αj,k11,k223is O((1/3)42J-1).Substitutingαj,k11,k121, αj,k11,k122,  and αj,k11,k223  and C1, C2 and C3 into (5), the multilevel wavelet interpolation operator can be obtained as
(11)In1,n2(x,y)=∑k01=02j0∑k02=02j0Rk01,k02,n1,n2j0,j0,J,Jϕk01,k02j0(x,y)+∑j=j0J-1∑k1=02j∑k2=02j(C1k1,k2,n1,n2j,j,J,Jϕ2k1+1,2k2j+1×(x,y)u(xn1J,yn2J)+C2k1,k2,n1,n2j,j,J,Jϕ2k1,2k2+1j+1×(x,y)u(xn1J,yn2J)+C3k1,k2,n1,n2j,j,J,Jϕ2k1+1,2k2+1j+1×(x,y)u(xn1J,yn2J)).
Then, (5) can be rewritten as
(12)uJ(m,n)(x,y,t)=∑n12J∑n22JIn1,n2(x,y)u(xn1J,yn2J).
Substituting (12) into (1), the multilevel wavelet discretization scheme of Perona-Malik model can be obtained.


## 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:(1)∂u(x,y,t)∂t=div(c(|∇u|)∇u),u(x,y,0)=f(x,y),
where (x,y) denotes pixel position, t is the time parameter, f(x,y) is the 2D image being processed, u(x,y,t) is the image after diffusion processing, and u(x,y,0) is the initial value. div denotes the divergence operator, and the ∇u denotes the gradient operator; c(|∇u|) denotes the diffusion coefficient, which is a nonnegative decreasing function of the image gradient modulus. It is usually taken as
(2)c(|∇u|)=11+(∇u/k)2
or
(3)c(|∇u|)=exp[-(∇uk)2],
where k is a constant.


## 2.2. Wavelet Numerical Discretization Scheme of Perona-Malik Model

Let the definition domain of the image be(xmin,xmax)×(ymin,ymax); the discretization points can be defined as (xk1j,yk2j), where j is a scale parameter and k1 and k2 are position parameters. So,
(4)xk1j=xmin+k1xmax-xmin2j,yk1j=ymin+k2ymax-ymin2j,j,k1,k2∈ℤ.In addition,wk1,k2j(m,n)(x,y) denotes the multiscale wavelet function and the corresponding mth and nth derivatives with respect to x and y, respectively. The level set function  ϕ(x,y,t) and the corresponding derivative function can be descretized as follows:
(5)uJ(m,n)(x,y,t)=∑k01=01∑k02=01u(xk010,yk020)wk01,k020(m,n)(x,y)+∑j=0J-1∑k11=02j-1∑k12=02j-1[αj,k11,k121(t)w2k11+1,2k12j+1(m,n)(x,y)+αj,k11,k122(t)w2k11,2k12+1j+1(m,n)(x,y)+αj,k11,k123(t)w2k11+1,2k12+1j+1(m,n)(x,y)],
where j and J are constants, which denote the wavelet scale number and the maximum of the scale number, respectively.  αj,k11,k121, αj,k11,k122,  and αj,k11,k123are the wavelet coefficients at the point (xk1j,yk2j). According to the interpolation wavelet transform theory, the wavelet coefficients can be written as
(6)αj,k1,k21=u(xj+1,2k1+1,yj+1,2k2)-Iju(xj+1,2k1+1,yj+1,2k2),αj,k1,k22=u(xj+1,2k1,yj+1,2k2+1)-Iju(xj+1,2k1,yj+1,2k2+1),αj,k1,k23=u(xj+1,2k1+1,yj+1,2k2+1)-Iju(xj+1,2k1+1,yj+1,2k2+1),
where Ij  denotes the multilevel interpolation operator. In order to obtain the multilevel interpolation operator, it is necessary to express the wavelet coefficients  αj,k1,k21, αj,k1,k22, and  αj,k1,k23 as a weighted sum of u in all of the collocation points in the J-level. Therefore, we should give the definition of the restriction operator as follows:
(7)Rk1,k2,m1,m2l,l,j,j={1,xk1l=xm1j,yk2l=ym2j0,otherwise.
Using the restriction operator, u(x2k1+1j+1,y2k2j+1),u(x2k1j+1,y2k2+1j+1),  and   u(x2k1+1j+1,y2k2+1j+1) can be rewritten as
(8)u(x2k1+1j+1,y2k2j+1)=∑n1=02J∑n2=02JR2k1+1,2k2,n1,n2j+1,j+1,J,Ju(xn1J,yn2J),u(x2k1j+1,y2k2+1j+1)=∑n1=02J∑n2=02JR2k1,2k2+1,n1,n2j+1,j+1,J,Ju(xn1J,yn2J),u(x2k1+1j+1,y2k2+1j+1)=∑n1=02J∑n2=02JR2k1+1,2k2+1,n1,n2j+1,j+1,J,Ju(xn1J,yn2J).
Introducing the extension operators C1, C2, and C3 and substituting (8) into (6), the wavelet coefficients can be rewritten as
(9)αj,k1,k21=∑n1=02J∑n2=02JR2k1+1,2k2,n1,n2j+1,j+1,J,Ju(xn1J,yn2J)-[∑n1=02J∑n2=02J∑k01=02j0∑k02=02j0Rk01,k02,n1,n2j0,j0,J,Ju(xn1J,yn2J)×ϕk01,k02j0(x2k1+1j+1,y2k2j+1)+∑j1=j0j-1∑n1=02J∑n2=02J∑k11=02j1∑k12=02j1(C1k11,k12,n1,n2j1,j1,J,Jϕ2k11+1,2k12j1+1×(x2k1+1j+1,y2k2j+1)u(xn1J,yn2J)+C2k11,k12,n1,n2j1,j1,J,Jϕ2k11,2k12+1j1+1×(x2k1+1j+1,y2k2j+1)u(xn1J,yn2J)+C3k11,k12,n1,n2j1,j1,J,J×ϕ2k11+1,2k12+1j1+1(x2k1+1j+1,y2k2j+1)×u(xn1J,yn2J))∑k02=02j0Rk01,k02,n1,n2j0,j0,J,J]=∑n1=02J∑n2=02JC1k1,k2,n1,n2j,j,J,Ju(xn1J,yn2J).αj,k1,k22and αj,k1,k23 are similar to αj,k1,k21. From the above equation, the extension operator can be obtained as
(10)C1k1,k2,n1,n2j,j,J,J=R2k1+1,2k2,n1,n2j+1,j+1,J,J-[∑k01=02j0∑k02=02j0Rk01,k02,n1,n2j0,j0,J,Ju(xn1J,yn2J)×ϕk01,k02j0(x2k1+1j+1,y2k2j+1)+∑j1=j0j-1∑n2=02J∑k11=02j1(C1k11,k12,n1,n2j1,j1,J,Jϕ2k11+1,2k12j1+1×(x2k1+1j+1,y2k2j+1)u(xn1J,yn2J)+C2k11,k12,n1,n2j1,j1,J,Jϕ2k11,2k12+1j1+1×(x2k1+1j+1,y2k2j+1)u(xn1J,yn2J)+C3k11,k12,n1,n2j1,j1,J,Jϕ2k11+1,2k12+1j1+1×(x2k1+1j+1,y2k2j+1)×u(xn1J,yn2J))-∑k01=02j0∑k02=02j0Rk01,k02,n1,n2j0,j0,J,J].C2 and C3 can be obtained with the same method. Therefore, the calculation time complexity of the wavelet transform coefficient αj,k11,k121, αj,k11,k122, and αj,k11,k223is O((1/3)42J-1).Substitutingαj,k11,k121, αj,k11,k122,  and αj,k11,k223  and C1, C2 and C3 into (5), the multilevel wavelet interpolation operator can be obtained as
(11)In1,n2(x,y)=∑k01=02j0∑k02=02j0Rk01,k02,n1,n2j0,j0,J,Jϕk01,k02j0(x,y)+∑j=j0J-1∑k1=02j∑k2=02j(C1k1,k2,n1,n2j,j,J,Jϕ2k1+1,2k2j+1×(x,y)u(xn1J,yn2J)+C2k1,k2,n1,n2j,j,J,Jϕ2k1,2k2+1j+1×(x,y)u(xn1J,yn2J)+C3k1,k2,n1,n2j,j,J,Jϕ2k1+1,2k2+1j+1×(x,y)u(xn1J,yn2J)).
Then, (5) can be rewritten as
(12)uJ(m,n)(x,y,t)=∑n12J∑n22JIn1,n2(x,y)u(xn1J,yn2J).
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 parametert. In other words, the wavelet coefficient should vary with the time parameter t. 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,u and its derivative in (1) should be rewritten as
(13)∂u∂t=F(t,x,y,u,∂u∂x,∂u∂y,∂2u∂x2,∂2u∂x∂y,∂2u∂y2),(t>0),u(x,y,0)=u0(x,y),duJ(x,y,t)dt=F[t,x,y,uJ(x,y,t),uJ(1,0)(x,y,t),uJ(0,1)(x,y,t),uJ(2,0)(x,y,t),uJ(1,1)(x,y,t),uJ(0,2)(x,y,t)],
respectively, where
(14)uJ(x,y,t)=∑k01=01∑k02=01u(xk010,yk020)∅k01,k020(x,y)+∑j=0J-1∑k11=02j-1∑k12=02j-1[αj,k11,k121∅2k11+1,2k12j+1(x,y)+αj,k11,k122∅2k11,2k12+1j+1(x,y)+αj,k11,k123∅2k11+1,2k12+1j+1(x,y)].
The value of uJ(x,y,tn) at tn is denoted as un, and
(15)F[tn,x,y,uJ(x,y,tn),uJ(1,0)(x,y,tn),uJ(0,1)(x,y,tn),uJ(2,0)(x,y,tn),uJ(1,1)(x,y,tn),uJ(0,2)(x,y,tn)]
is denoted as Fn. And then a linear homotopy function can be constructed as
(16)uJ(x,y,t)=(1-ε)Fn+εFn+1.
It is easy to identify the homotopy parameter as
(17)ε(t)=t-tntn+1-tn,t∈[tn,tn+1]∴ε∈[0,1].
According to the perturbation theory, the solution of (16) can be expressed as the power series expansion of ε(18)uJ=u0J+εu1J+ε2u2J+⋯.
Substituting (18) into (16) and rearranging based on powers of  ε-terms, we have
(19)ε0∶u0J=Fnε1∶u1J=Fn+1-Fn⋮
According to HPM, we obtain the wavelet coefficient αj,k1,k21(tn+1),    αj,k1,k22(tn+1),   and  αj,k1,k23(tn+1) at tn  as follows:
(20)αj,k1,k21=u(x2k1+1j+1,y2k2j+1)-Iju(x2k1+1j+1,y2k2j+1)=u(x2k1+1j+1,y2k2j+1)-[∑k01=01∑k02=01u(xk010,yk020)∅k01,k020(x2k1+1j+1,y2k2j+1)+∑j1=0j-1∑k11=02j1∑k12=02j1(αj1,k11,k121∅2k11+1,2k12j1+1×(x2k1+1j+1,y2k2j+1)+αj1,k11,k122∅2k11,2k12+1j1+1×(x2k1+1j+1,y2k2j+1)+αj1,k11,k123∅2k11+1,2k12+1j1+1×(x2k1+1j+1,y2k2j+1))∑k02=01u(xk010,yk020)],αj,k1,k22=u(x2k1j+1,y2k2+1j+1)-Iju(x2k1j+1,y2k2+1j+1)=u(x2k1j+1,y2k2+1j+1)-[∑k01=01∑k02=01u(xk010,yk020)∅k01,k020×(x2k1j+1,y2k2+1j+1)+∑j1=0j-1∑k11=02j1∑k12=02j1(αj1,k11,k121∅2k11+1,2k12j1+1×(x2k1j+1,y2k2+1j+1)+αj1,k11,k122∅2k11,2k12+1j1+1×(x2k1j+1,y2k2+1j+1)+αj1,k11,k123∅2k11+1,2k12+1j1+1×(x2k1j+1,y2k2+1j+1))∑k01=01∑k02=01u(xk010,yk020)],αj,k1,k23=u(x2k1+1j+1,y2k2+1j+1)-Iju(x2k1+1j+1,y2k2+1j+1)=u(x2k1+1j+1,y2k2+1j+1)-[∑k01=01∑k02=01u(xk010,yk020)∅k01,k020×(x2k1+1j+1,y2k2+1j+1)+∑j1=0j-1∑k11=02j1∑k12=02j1(αj1,k11,k121∅2k11+1,2k12j1+1×(x2k1+1j+1,y2k2+1j+1)+αj1,k11,k122∅2k11,2k12+1j1+1×(x2k1+1j+1,y2k2+1j+1)+αj1,k11,k123∅2k11+1,2k12+1j1+1×(x2k1+1j+1,y2k2+1j+1))∑k01=01∑k02=01u(xk010,yk020)].
Obviously, the calculation time complexity of the wavelet transform coefficient αj,k1,k21, αj,k1,k22, and αj,k1,k23 is O(4J), which is decreased greatly than it in (6) which is  O((1/3)42J-1).Substituting the wavelet transform efficient (20) into (14), we obtain
(21)uJ(x,y,tn+1)=uJ(x,y,tn)+Δt2[F(tn,x,y,uJ(x,y,tn),uJ(1,0)(x,y,tn),uJ(0,1)(x,y,tn),uJ(2,0)(x,y,tn),uJ(1,1)(x,y,tn),uJ(0,2)(x,y,tn))+F(tn+1,x,y,u0J(x,y,tn+1),u0J(1,0)(x,y,tn+1),u0J(0,1)(x,y,tn+1),u0J(2,0)(x,y,tn+1),u0J(1,1)(x,y,tn+1),u0J(0,2)(x,y,tn+1))]
and the derivative function
(22)uJ(m,n)(x,y)=∑k01=01∑k02=01u(xk010,yk020)∅k01,k020(m,n)(x,y)+∑j=0J-1∑k11=02j-1∑k12=02j-1[αj,k11,k121∅2k11+1,2k12j+1(m,n)(x,y)+αj,k11,k122∅2k11,2k12+1j+1(m,n)(x,y)+αj,k11,k123∅2k11+1,2k12+1j+1(m,n)(x,y)].
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
(23)δΔσ=sin(πx/Δ)πx/Δexp(-x22σ2),
where Δis the discrete space and σ=rΔ (r 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 Figure1.Figure 1
Image denoising results and the corresponding collocation points distribution at different time parameterst.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 Figure2.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 aMetarhizium 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.The comparison about the smoothing result applying HAM.
(a)
HPM-based wavelet interpolation operator
(b)
Finite difference operator


## 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 Figure1.Figure 1
Image denoising results and the corresponding collocation points distribution at different time parameterst.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 Figure2.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 aMetarhizium 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.The comparison about the smoothing result applying HAM.
(a)
HPM-based wavelet interpolation operator
(b)
Finite difference operator


## 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.



---

*Source: 309418-2013-09-25.xml*