The Scientific World Journal

Volume 2014 (2014), Article ID 417486, 15 pages

http://dx.doi.org/10.1155/2014/417486

Research Article

HPM-Based Dynamic Sparse Grid Approach for Perona-Malik Equation

Shu-Li Mei¹ and De-Hai Zhu²

¹College of Information and Electrical Engineering, China Agricultural University, Beijing 100083, China

²Key Laboratory of Agricultural Information Acquisition Technology, Ministry of Agriculture, China Agricultural University, Beijing 100083, China

Received 24 January 2014; Accepted 2 March 2014; Published 23 June 2014

Academic Editors: D. Baleanu, H. Jafari, and C. M. Khalique

Copyright © 2014 Shu-Li Mei and De-Hai Zhu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The Perona-Malik equation is a famous image edge-preserved denoising model, which is represented as a nonlinear 2-dimension partial differential equation. Based on the homotopy perturbation method (HPM) and the multiscale interpolation theory, a dynamic sparse grid method for Perona-Malik was constructed in this paper. Compared with the traditional multiscale numerical techniques, the proposed method is independent of the basis function. In this method, a dynamic choice scheme of external grid points is proposed to eliminate the artifacts introduced by the partitioning technique. In order to decrease the calculation amount introduced by the change of the external grid points, the Newton interpolation technique is employed instead of the traditional Lagrange interpolation operator, and the condition number of the discretized matrix different equations is taken into account of the choice of the external grid points. Using the new numerical scheme, the time complexity of the sparse grid method for the image denoising is decreased to O(4J+2j) from O(43J), (). The experiment results show that the dynamic choice scheme of the external gird points can eliminate the boundary effect effectively and the efficiency can also be improved greatly comparing with the classical interval wavelets numerical methods.

1. Introduction

The nonlinear difference equation has been widely used in various fields in the past few decades such as the option pricing [1], stochastic analysis [2], hydrodynamics [3], and image processing [4]. Many powerful and efficient methods to find analytic solutions of nonlinear equation have drawn a lot of interest by a diverse group of scientists. These methods include the tanh-function method, the extended tanh-function method [5, 6], the sine-cosine method [7], the variational iteration method [8, 9], the homotopy perturbation method [10, 11], and Exp-function method [12].

As an excellent medical image processing model, the Perona-Malik model [4] has been widely used in image denoising in recent years. Perona-Malik model is a nonlinear 2-dimension partial differential equation in itself, which overcomes the drawback of the scale-space technique introduced by Witkin which 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 accurately the 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.

It is very difficult to find the exact analytical solution of the Perona-Malik model as it is a nonlinear partial differential equation. Conventional methods for numerical solutions of partial differential equations mostly fall into three classes: finite difference methods, finite element methods, and spectral methods. Briefly, the finite difference method consists in defining the different unknowns by their values on a discrete grid and in replacing differential operators by difference operators using neighboring points. In the finite element method, the equations are integrated against a set of linear independent test functions with small compact support, and the solution is considered as a linear combination of this set of test functions. In spectral methods, the unknown functions are developed along a basis of functions having global support. This development is truncated to a finite number of terms which satisfy a system of coupled ordinary differential equations in time. The advantage of using either of the first two numerical techniques is the simplicity in adapting to complex geometries, while the main advantage of spectral methods is the greater accuracy [13, 14].

If the solution of a partial differential equation is regular, any of the three above-mentioned numerical techniques can be applied successfully. It is obvious that most of the images are irregular. This makes the Perona-Malik equation particularly difficult to resolve numerically using the above-mentioned methods. Spectral methods are not easily implemented because the irregularity of the solution causes the loss of high accuracy. Moreover, the global support of the basis function induces the well-known Gibbs phenomenon, which appears as the artifacts in images. Wavelet analysis is a new numerical concept which allows one to represent a function in terms of a set of basis function, called wavelets, which are localized both in location and in scale. Up to now, the finite difference method is the primary numerical algorithm for Perona-Malik model, which can bring artifact into the images due to the nonsmoothness of the basis function of the finite difference method [15, 16] as has been said before. The 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 [17]. The artifacts in image can be eliminated with the wavelet numerical algorithm instead of the finite 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 finite difference method [18, 19]. This leads to a lower computational efficiency of wavelet transform in solving 2D nonlinear PDEs. Besides, most of the wavelet algorithms for solving partial differential equations can handle periodic boundary conditions easily. The treatment of general boundary conditions is still an open question especially in solving the nonlinear problems. Construction of the wavelet defined in the interval (interval wavelet) is another good choice to handle the boundary conditions [20, 21]. Compared to the interpolation wavelet, a linear mapping between the external collocation points and the interval ones was supplied in the interval wavelet. The choice of the external collocation points depends on the smoothness and the gradient near each collocation point of the solution of the PDEs. Besides, the condition number of the system of equations obtained by the wavelet collocation method should be taken into account.

To an image with pixels (), the Perona-Malik equation can be discretized into a system of ODEs with -dimension by the coupling technique of HPM [22–27] and the wavelet collocation method [28, 29]. The corresponding time complexity is about with the variational iterative method for the system of ODEs [30]. Obviously, it does not meet the requirement of the larger image processing. Partitioning technique is the effective measure to improve the efficiency of this problem. In other words, the image should be divided into several blocks before denoising to the images. In each of image blocks, the multiple programs can be executed simultaneously. This is similar to the finite element method to some extent. Obviously, if the size of the image blocks is adaptive to whole image, the algorithm efficiency can be improved furthermore. Our research focuses on the general frame of sparse grids and the dynamic choice scheme of the external grid points, which can be used to decrease the boundary effect of each image block, and so, we just talk about the even partitioning in this paper for simplification.

The sparse representation of functions via a linear combination of a small number of basic functions has recently received a lot of attention in several mathematical fields such as approximation theory as well as signal and image processing. The advantage of the sparse grid approach is that it can be extended to nonsmooth solutions by adaptive refinement methods; that is, it can capture the steep waves that appeared in the solution of the PDEs. The main objective of the paper is to present a dynamic choice scheme of the external grid points and a general sparse grid operator for solving the Perona-Malik equation. In other words, the dynamic sparse grid approach provides an adaptive choice scheme on both of the external and the internal grid points. In the presentation of the method, we try to be as general as possible, giving only the main philosophy of the method and leaving some freedom for further exploration of its applications. Both the boundary condition and the condition number are addressed in this work. The first is how to incorporate the dynamic choice scheme on external grid points with the interpolation wavelet basis to construct an effective algorithm of solving partial differential equation. The second is how to construct a stable, accurate, and efficient numerical algorithm for the image denoising model.

2. Construction of Dynamic Sparse Grid Operator

There are many ways to eliminate the boundary effect from the multiscale basis. A simple solution is the even 2-periodical extension of function , which is usually used in image analysis. Unfortunately, this extension generally produces discontinuities at the integers that are indicated by the large transform coefficients near the endpoints 0 and 1. Thus the constructed multiscale basis cannot exactly analyze the boundary behavior of a given function. To solve this problem, the popular method is using special boundary and interior scaling functions such as the interval wavelet to reduce the numerical problem at the boundaries. To the interpolation basis function, the common approach is to define the interpolation basis in the interval with the Lagrange multiplier. In fact, the Lagrange multiplier can be viewed as a map operator, which maps the external collocation points into the definition domain in the multiscale interpolation method. The choice of the amount of the external points relates to the smoothness and gradient near the boundary of the approximated function. In addition, another factor that we should take into account is the condition number of the system of ODEs obtained by the multiscale numerical method.

Obviously, the amount of the external collocation points should be different to different boundary conditions such as the smoothness, gradient near the boundary, and the condition number. In the partition technique about the image processing, the boundary conditions of the different image blocks are obviously different as the randomness of the image. In the representation, we try to give a dynamic choice scheme about the external collocation points to meet the requirement of the image partition technique, in which all above 3 factors are taken into account.

In the presentation of the method, we try to be as general as possible, giving only the main philosophy of the method and leaving some freedom for further exploration of its applications. We illustrate the method using two classical interpolation wavelets: Shannon wavelet and the autocorrelation function of Daubechies scaling functions. But we do not try to predict what wavelet is the best for our algorithm (it is simply impossible, due to the fact that some wavelets work better for some problems and worse for others).

2.1. Basis Functions with Interpolation Property

There are many wavelet functions which possess the interpolation property. The familiar interpolation wavelets family includes Shannon wavelet, Haar wavelet, and Faber-Schauder. Furthermore, it is easy to understand that the autocorrelation function of the orthogonal wavelet function also has the interpolation property. So, the autocorrelation function of the Daubechies scaling function is often employed to construct the wavelet collocation method.

The representation of Shannon wavelet [31, 32] is based upon approximating the Dirac delta function as a band-limited function and is given by

The Shannon wavelet possesses many excellent numerical properties such as interpolating, relative sparse, and orthogonal properties. A perceived disadvantage of (1) is that it tends to zero quite slowly as . A direct consequence of this is that there are a large number of grid points will contribute to the derivatives calculation of approximated function. For this reason Hoffman et al. [33] have suggested using the Shannon-Gabor wavelet as follows:

where is the width parameter (or called window size). It has been proofed that (2) can improve the localized and asymptotic behavior of the Shannon scaling function. A consequence of this is that it ensures that derivatives at any one point are more dependent on the neighboring nodal values than on the nodal values further away from the point considered. However, the presence of the Gaussian window destroys the orthogonal properties possessed by the Shannon wavelet, effectively worsening the approximation to a Dirac delta function. In the following, the Shannon wavelet representation of the Dirac delta function is adopted, and it is shown that this representation ensures that the approach is identical to the weighted residual approach.

The autocorrelation functions of compactly supported scaling functions were first studied in the context of the Lagrange iterative interpolation scheme in [34]. Let be the autocorrelation function:

where is the scaling function which appears in the construction of compactly supported wavelet. The function is exactly the “fundamental function” of the symmetric iterative interpolation scheme introduced in [35]. Thus, there is a simple one-to-one correspondence between iterative interpolation schemes and compactly supported wavelet. In particular, the scaling function corresponding to Daubechies’s wavelet with two vanishing moments yields the scheme in [36]. In general, the scaling functions corresponding to Daubechies’s wavelets with vanishing moments lead to the iterative interpolation schemes which use the Lagrange polynomials of degree 2M. Additional variants of iterative interpolation schemes may be obtained using compactly supported wavelets described in [37].

2.2. Construction of Dynamic Interpolation Wavelet in Interval

According to the definition of the interval wavelet, the interval interpolation basis functions can be expressed as

where

where is the amount of the external collocation points; the amount of discrete points in the definition domain is ; and is the definition domain of the approximated function. Equations (4) and (5) illustrate that the interval wavelet is derived from the domain extension. The supplementary discrete points in the extended domain are called external points. The value of the approximated function at the external points can be obtained by Lagrange extrapolation method. Using the interval wavelet to approximate a function, the boundary effect can be left in the supplementary domain; that is, the boundary effect is eliminated in the definition domain.

According to (4) and (5), the interval wavelet approximant of the function can be expressed as

where is the given value at the discrete point . At the external points, can be obtained by extrapolation; that is,

So, the interval wavelet approximant of can be rewritten as

Let

Then,

where and correspond to the left and the right external points, respectively. They are obtained by Lagrange extrapolation using the internal collocation points near the boundary. So, the interval wavelet’s influence on the boundary effect can be attributed to Lagrange extrapolation. It should be pointed out that we did not care about the reliability of the extrapolation. The only function of the extrapolation is enlarging the definition domain of the given function which can avoid the boundary effect that occurred in the domain. Therefore, we can discuss the choice of by means of Lagrange inner- and extrapolation error polynomial as follows:

Equation (11) indicates that the approximation error is related to both the smoothness and the gradient of the original function near the boundary. Setting different can satisfy the error tolerance.

2.3. Adaptive Interval Interpolation Wavelet

The interval interpolation wavelet is often used to solve the diffusion PDEs with Neumann boundary conditions. The smoothness and gradient of the PDE’s solution usually vary with the time parameter. If the parameter is a constant, we have to take a bigger value in order to obtain result with higher calculation precision. But the bigger usually introduces the famous Gibbs phenomenon into the numerical solution, which usually makes the algorithm become invalid. In addition, the bigger will bring much more calculation. To keep higher numerical precision and save calculation, the best way is to design a procedure that can vary with the curve’s smoothness and gradient dynamically.

In this dynamic procedure, the error estimation equation (11) can be taken as the criterion about . But in most cases, we cannot know the smoothness and the derivative’s order of the original function. This can be solved by substituting the difference coefficient for the derivative. This is coincident with the Newton interpolation equation which is equivalent with Lagrange interpolation equation. In addition, the Lagrange interpolation algorithm has no inheritance which is the key feature of Newton interpolation. So, the basis function has to be calculated repeatedly as interpolation points are added into the calculation, which increases the computation complexity greatly. In contrast to the Lagrange method, the advantage of Newton interpolation method is that the basis function need not be recalculated as one point is added except only one more term which is needed to be added, which reduces the number of computed operations, especially the multiplication. So, it is convenient using the Newton interpolation method to construct the dynamic procedure.

2.3.1. Newton Interpolation

The expression of Newton interpolation can be written as

Substitute the Newton interpolation instead of the Lagrange interpolation into (29), which can be rewritten as

where

2.3.2. Relation between the Newton Interpolation Error and the Choice of

It is well known that the Newton interpolation is equivalent to the Lagrange interpolation. The corresponding error estimation can be expressed as

And the simplest criterion to terminate the dynamic choice on is ( is the absolute error tolerance). Obviously, it is difficult to define which should meet the precision requirement of all approximated curves. In fact, the difference coefficient can be used directly as the criterion; that is,

As mentioned above, once the curves with lower order smoothness are approximated by higher order polynomial expression, the errors will become bigger on the contrary. In fact, even if the is infinite, the computational precision cannot be satisfied except increasing computational complexity. To avoid this, we design the termination procedure of dynamic choice about as follows:

 if , then , else, if or , then , else, if or , then ,

2.3.3. L and the Condition Number of the System of Algebraic Equations

In the field of numerical analysis, the condition number of a function with respect to an argument measures how much the output value of the function can change for a small change in the input argument. This is used to measure how sensitive a function is to changes or errors in the input and how many errors in the output result from an error in the input. There is no doubt that the choice of can change the condition number of the system of algebraic equations discretized by the wavelet interpolation operator or the finite difference method. Therefore, the choice of should take the condition number into account. In fact, if the condition number , then you may lose up to digits of accuracy on top of what would be lost to the numerical method due to loss of precision from arithmetic methods [34]. According to the general rule of thumb, the choice should follow the rule as follows:

2.3.4. Relation between and Computation Complexity

The computational complexity of interpolation calculation is not proportional to the increasing points. The former is mainly up to the computation amount of and the derivative operations. Obviously, according to (2), the increase in computational complexity is when the number of extension points increases by 1. But the computational complexity of adaptively increasing collocation points is related to the different wavelet functions. For the wavelet with compact support property such as Daubechies wavelet and Shannon wavelet, the value of is impossible to be infinite. For Haar wavelet which has no smoothness property, L can be taken as 0 at most since it need not be extended. For Faber-Schauder wavelet, L can be taken as 1 at most. For Daubechies wavelet, L can be taken as different values according to the order of vanishing moments, but it must be finite. For the wavelets without compact support property, L can take value dynamically, such as Shannon wavelet. The computational complexity of increasing points is mainly up to the wavelet function of itself.

3. Construction of the Multilevel Interpolation Operator Based on the Interval Wavelet

Let the definition domain of the image be ; the discretization points can be defined as , where is a scale parameter and and are position parameters. So,

In addition, 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 discretized as follows:

where and are constants, which denote the wavelet scale number and the maximum of the scale number, respectively. , , and are the wavelet coefficients at the points . According to the interpolation wavelet transform theory, the wavelet coefficients can be written as

where denotes the multilevel interpolation operator. In order to obtain the multilevel interpolation operator, it is necessary to express the wavelet coefficients as a weighted sum of in all of the collocation points in the J level. Therefore, we should give the definition of the restriction operator as follows:

Using the restriction operator, ,  and   can be rewritten as

Introducing the extension operators C1, C2, and C3, and substituting (22) into (20), the wavelet coefficients can be rewritten as

and are similar with . From above equation, the extension operator can be obtained as

and can be obtained with the same method. Therefore, the calculation time complexity of the wavelet transform coefficients , , and is .

Substituting , , and and C1, C2, and C3 into (2), the multilevel wavelet interpolation operator can be obtained as

Then, (10) can be rewritten as

Substituting (26) into (13), the multilevel wavelet discretization scheme of Perona-Malik model can be obtained.

The purpose of constructing the multilevel sparse grid approach 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 operators, HPM would be employed to construct a novel interpolation operator in this section.

For convenience, and its derivative in (5) should be rewritten as

respectively, where

The value of at is denoted by , and

is denoted by . And then, a linear homotopy function can be constructed as

It is easy to identify the homotopy parameter as

According to the perturbation theory, the solution of (31) can be expressed as the power series expansion of :

Substituting (29) into (27) and rearranging based on powers of -terms, we have

According to HPM, we obtain the wavelet coefficients at as follows:

Obviously, the calculation time complexity of the wavelet transform coefficients is , which is decreased greatly than in (11) which is .

Substituting the wavelet transform coefficient (35) into (29), we obtain

And the derivative function

Obviously, the computation complexity is decreased greatly comparing with (26).

4. Numerical Experiences and Discussion

In this section, we take some images as examples to illustrate the efficiency of the dynamic interval wavelet interpolation operator based on HPM in partitioning technique on the image processing. In fact, the partitioning technique is a scheme to divide the image into several subimages in the multiscale wavelet numerical method to improve the efficiency. The dynamic interval wavelet provides an adaptive choice scheme for the external collocation points to eliminate the boundary effect of the subimages. Perona-Malik equation is employed as the denoising model, which is an anisotropic diffusion image denoising model that was proposed by Perona and Malik. It has been widely used in various image processing fields. It can be represented as the nonlinear partial differential equations:

where denotes pixel position, is the time parameter, is the 2D image being processed, is the image after diffusion processing, and is the initial value. denotes the divergence operator, denotes the gradient operator, and denotes the diffusion coefficient, which is nonnegative decreasing function of the image gradient modulus. It is usually taken as

or

where is a constant.

Two different medical images are taken as examples to test the characteristic of different interpolation wavelets, which is showed in Figure 1. One is the human brain (Figure 1(a)), which has so clear contour that that image cannot be represented as a continuous function near the contour. The Gibbs phenomenon is possible to be introduced into the image near the boundary. So, this can be used to test the advantages of the multiscale wavelet approximation comparing with the difference operator. Another one is the image of the locust coelom, which has many microgrooves without clear boundary. This image is used to test the characteristic of different interpolation wavelets, which is showed in Figure 1(b).

(a) Image of human brain

(b) Image of locust coelom

Figure 1: Original images.

4.1. Comparison between the Sparse Grid Approach and the Finite Different Method

It has been mentioned above that the brain image is used to test the difference between the sparse grid approach and the finite difference method and the difference between different wavelet functions which are taken as the basis functions in the sparse grid approach. In this experiment, all the results are obtained by solving the Perona-Malik equation with different methods, which have been showed in Figure 2.

(a) Original image

(b) Autocorrelation function of Daubechies scaling function

(c) Shannon scaling function

(d) Finite difference method

Figure 2: Comparison between different numerical methods for image denoising (time step , terminal time ).

Two interpolation wavelet scaling functions are employed to test the dynamic sparse grid approach for image denoising proposed in this paper. The Shannon wavelet possesses the smoothness and or the orthogonality but has no compact support property. Daubechies scaling function possess almost all the excellent properties in numerical algorithm such as smoothness, orthogonality, and compact support property. But what we utilized in this research is the autocorrelative function of the Daubechies scaling function, which keeps the better edge preserving property although it loses the orthogonality. It can be easily observed from Figure 2(c) that the evident artifacts appeared in the denoised image obtained by the Shannon sparse grid approach. That is, the Gibbs phenomenon has appeared in the Shannon scaling function representation of the image near the boundary. In contrast to the Shannon wavelet, the denoised image (Figure 2(b)) obtained by the Daubechies wavelet sparse grid approach has clear boundary. It is easy to understand that the compact support property of the wavelet scaling function is helpful to eliminate the Gibbs phenomenon and so to improve the numerical performances of the wavelet numerical methods.

Comparing with the sparse grid approaches, the finite difference method utilizes the difference operator to approximate the derivative in Perona-Malik equation, which decreases the value of the derivative to some extent. Therefore, the edge of the brain contour is smoothed in denoised images; this is showed in Figure 2(d). It should be noticed that the edge of the denoised image obtained by the Shannon wavelet sparse grid approach is more clear than that obtained by the finite difference method, in despite of the appearing artifacts.

4.2. Comparison between the Dynamic Interval Wavelet and the Static Interval Wavelet

For convenience of comparison, we call the interval interpolation wavelet constructed by the Lagrange interpolator as static interval wavelet and the interval interpolation wavelet based on the Newton interpolator as the dynamic interval wavelet. The difference between two interval wavelets above is the choice of the parameter . The value of is constant to the static interval wavelet, and it varies with both of the boundary condition and the condition number of the system ODEs obtained from the sparse grid approach.

The purpose of constructing of the interval wavelet is to eliminate the boundary effect in the partitioning technique on the image denoising process. In this section, the image of locust coelom (300 * 300 pixels) is taken as example to compare the difference between the dynamic and static interval wavelets. According to the partitioning technique, the image is divided evenly into 25 parts for simplification. So, the size of each image block is 60 * 60 pixels (Figure 3). According to the sparse grid approach based on HPM, the calculation amount decreases from (300 * 300)³ to 25 * (60 * 60)³. It has been mentioned that there are many ways to eliminate the boundary effect such as the extension method and the interval wavelet method. There is no doubt that the interval wavelet method is more efficient than the extension method. According to the interval interpolation wavelet based on the Lagrange interpolator, the amount of the external collocation points is a constant. With increase of , the calculation amount will increase correspondingly.

(a) Static interval wavelet ()

(b) Static interval wavelet ()

(c) Static interval wavelet ()

(d) Dynamic interval wavelet

Figure 3: Comparison between the dynamic interval wavelet and the static interval wavelet.

L is taken as 1, 2, and 3, respectively, in the experiments. It is easy to be observed from Figures 3(a)–3(c) that there are more collocation points near the boundary of each of block images in all 3 cases. In fact, the adaptive increase of the collocation points can also eliminate the boundary effect. Therefore, there are no artifacts appearing in the denoised images in the first two cases. But the increase of the collocation points can increase the calculation amount greatly. According to the definition of the interval interpolation wavelet based on the Lagrange operator, the increase of can improve the smoothness and the precision of the approximated function near the boundary. This is helpful to decrease the boundary effect in theory. In contrast to the theory, the collocation points in the whole image domain increased so much that the artifacts appeared in the denoised subimages when comparing to other two cases (Figure 3(c)). As a matter of fact, this is caused by the increase of the condition number of the system of ODEs obtained by the sparse grid approach. That is, the increase of the value of can induce the condition number change greatly; this is showed in Table 1. It has been pointed out in Section 2 that if the condition number , then you may lose up to digits of accuracy on top of what would be lost to the numerical method due to loss of precision from arithmetic methods. This also illustrates that the condition number must be taken into account in the dynamic interval wavelet. Figure 3(d) is the result obtained by the dynamic interval sparse grid approach. The distribution of the collocation points in Figure 3(d) is just correlative with the image content itself and is not correlative with the partitioning scheme of the image anymore. The amount of the wavelet collocation points also decreased accordingly.

Table 1: Condition number of each image block at different times.

Block number
Condition number

1
9.6829
7.7193
6.7935
6.8257
5.9730

2
13.7539
12.8756
11.8760
11.8703
11.3319

3
9.2841
8.8832
8.1376
7.8737
6.8713

4
9.6829
7.7193
6.7935
6.8257
5.9730

5
1.6816
1.6931
1.7153
1.7378
1.8923

6
13.9357
12.9657
11.8891
11.8757
10.3356

7
13.7543
12.8757
11.8765
10.8701
9.3318

8
13.7556
12.8781
11.8776
11.8734
11.3329

9
43.9354
32.9663
31.8882
28.8749
21.3346

10
3.2389
2.9137
2.2917
1.9365
1.2919

11
8.5692
6.9416
6.6971
5.2875
4.3907

12
31.9823
29.3266
28.8330
25.6736
20.1976

13
15.1617
14.7818
13.8967
12.5738
12.3428

14
12.9614
12.1973
10.9887
10.1725
9.3189

15
4.6593
3.7938
2.5476
2.3789
1.9916

16
14.6835
13.9864
13.1838
12.8897
11.3156

17
12.8113
12.1031
10.1763
9.5627
8.2474

18
14.9834
13.6523
12.7719
12.1658
11.7829

19
13.6689
12.1791
11.1782
10.5977
8.9664

20
125.4782
110.3379
98.9073
89.7761
80.2749

21
1.6816
1.6931
1.7153
1.7378
1.8923

22
2.4589
2.6623
3.7662
3.8955
4.5811

23
43.2983
39.6744
36.7943
32.1079
28.6179

24
211.5877
198.7219
180.7089
86.9125
81.8510

25
185.7428
170.5897
160.0987
52.1757
83.0421

4.3. Adaptability of the Wavelet Collocation Points

In this research, the dynamic interpolation operator was viewed as a nonlinear problem, and so HPM is employed in construction of the dynamic interpolation wavelet defined in interval. This is helpful to improve the efficiency of the multilevel wavelet interpolation operators. In this section, the autocorrelation function of the Daubechies scaling function is employed to construct the dynamic interval wavelet. The brain image is taken as example to test the precision and efficiency of the HPM-based dynamic interval wavelet proposed in this research. The experiment results were showed in Figure 4. It is easy to be observed that the noise pixels of the brain images were smoothed completely and the edges of the brain contour were preserved perfectly. With the increase of the iteration times, more and more trivial objects such as the noise pixels are being smoothed, and more areas in the brain image are becoming smoother. Accordingly, the amount of the wavelet collocation points should be smaller and smaller. This has been illustrated in Figures 4(a) and 4(b). In this experiment, the time step ; the definition domain of the parameter is [0, 0.001]. The experiment results show that the amount of the wavelet collocation points decreases from 23488 to 19413 with the parameter increasing from 0.0005 to 0.001. According to the finite difference method, the amount of the collocation points should be 90000, which is greater than the sparse grid approach, evidently. This illustrates that the dynamic interval sparse grid approach proposed in this paper is more efficient than the finite difference method.

(a)

(b)

Figure 4: Adaptability of the multiscale sparse grid approach in image denoising (time step ).

5. Conclusions

The dynamic interval wavelet and the corresponding numerical method proposed in this paper are intrinsically an adaptive choice scheme on the external collocation points. In partitioning technique about the image processing, the dynamic sparse grid approach can be used to eliminate the boundary effect and improve the algorithm efficiency. In this method, the wavelet interpolation operator is constructed based on the homotopy perturbation method, which can decrease the calculation amount greatly. In addition, comparing with the finite difference method, the dynamic interval sparse grid approach can preserve the object edge more clear, especially in the case that the edge is sharper. For simplification, the image is divided evenly into several parts according to the partitioning scheme in the experiments. It is obvious that the partitioning scheme can be adaptive, which can improve the efficiency furthermore.

Conflict of Interests

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

Acknowledgments

This work is supported by the National Natural Science Foundation of China (no. 41171337) and the National Key Technology R&D Program of China (no. 2012BAD35B02).

References

H.-H. Yan, “Adaptive wavelet precise integration method for nonlinear Black-Scholes model based on variational iteration method,” Abstract and Applied Analysis, vol. 2013, Article ID 735919, 6 pages, 2013.

L.-W. Liu, “Interval wavelet numerical method on fokker-planck equations for nonlinear random system,” Advances in Mathematical Physics, vol. 2013, Article ID 651357, 7 pages, 2013.

R.-Y. Xing, “Wavelet-based homotopy analysis method for nonlinear matrix system and its application in burgers equation,” Mathematical Problems in Engineering, vol. 2013, Article ID 982810, 7 pages, 2013.

P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 12, no. 7, pp. 629–639, 1990.

A.-M. Wazwaz, “The tanh-coth method for solitons and kink solutions for nonlinear parabolic equations,” Applied Mathematics and Computation, vol. 188, no. 2, pp. 1467–1475, 2007.

A.-M. Wazwaz, “The tan h method: solitons and periodic solutions for the Dodd-Bullough-Mikhailov and the Tzitzeica-Dodd-Bullough equations,” Chaos, Solitons & Fractals, vol. 25, no. 1, pp. 55–63, 2005.

A.-M. Wazwaz, “The sine-cosine method for obtaining solutions with compact and noncompact structures,” Applied Mathematics and Computation, vol. 159, no. 2, pp. 559–576, 2004.

H. Jafari, H. Tajadodi, and D. Baleanu, “A modified variational iteration method for solving fractional Riccati differential equation by Adomian polynomials,” Fractional Calculus and Applied Analysis, vol. 16, no. 1, pp. 109–122, 2013.

H. Jafari, M. Saeidy, and D. Baleanu, “The variational iteration method for solving n-th order fuzzy differential equations,” Central European Journal of Physics, vol. 10, no. 1, pp. 76–85, 2012.

H. Jafari, M. Alipour, and H. Tajadodi, “Convergence of homotopy perturbation method for solving integral equations,” Thai Journal of Mathematics, vol. 8, no. 3, pp. 511–520, 2010.

H. Jafari, C. Chun, S. Seifi, and M. Saeidy, “Analytical solution for nonlinear gas dynamic equation by homotopy analysis method,” Applications and Applied Mathematics, vol. 4, no. 1, pp. 149–154, 2009.

H. Jafari, N. Kadkhoda, and C. M. Khalique, “Exact solutions of equation using lie symmetry approach along with the simplest equation and Exp-function methods,” Abstract and Applied Analysis, vol. 2012, Article ID 350287, 7 pages, 2012.

O. V. Vasilyev and S. Paolucci, “A dynamically adaptive multilevel wavelet collocation method for solving partial differential equations in a finite domain,” Journal of Computational Physics, vol. 125, no. 2, pp. 498–512, 1996.

H. Schaeffer, R. Caflisch, C. D. Hauck, and S. Osher, “Sparse dynamics for partial differential equations,” Proceedings of the National Academy of Sciences of the United States of America, vol. 110, no. 17, pp. 6634–6639, 2013.

C. Lopez-Molina, B. de Baets, J. Cerron, M. Galar, and H. Bustince, “A generalization of the Perona-Malik anisotropic diffusion method using restricted dissimilarity functions,” International Journal of Computational Intelligence Systems, vol. 6, no. 1, pp. 14–28, 2013.

I. Daubechies and G. Teschke, “Variational image restoration by means of wavelets: simultaneous decomposition, deblurring, and denoising,” Applied and Computational Harmonic Analysis, vol. 19, no. 1, pp. 1–16, 2005.

I. Daubechies and G. Teschke, “Wavelet based image decomposition by variational functionals,” in Wavelet Applications in Industrial Processing, vol. 5266 of Proceedings of SPIE, pp. 94–105, The International Society for Optical Engineering, Providence, RI, USA, February 2004.

C. Cattani, “Shannon wavelets theory,” Mathematical Problems in Engineering, vol. 2008, Article ID 164808, 24 pages, 2008.

C. Cattani, “Harmonic wavelets towards the solution of nonlinear PDE,” Computers and Mathematics with Applications, vol. 50, no. 8-9, pp. 1191–1210, 2005.

S.-L. Mei, H.-L. Lv, and Q. Ma, “Construction of interval wavelet based on restricted variational principle and its application for solving differential equations,” Mathematical Problems in Engineering, vol. 2008, Article ID 629253, 14 pages, 2008.

S.-L. Mei, Q.-S. Lu, S.-W. Zhang, and L. Jin, “Adaptive interval wavelet precise integration method for partial differential equations,” Applied Mathematics and Mechanics, vol. 26, no. 3, pp. 364–371, 2005.

J.-H. He, “Asymptotic methods for solitary solutions and compactons,” Abstract and Applied Analysis, vol. 2012, Article ID 916793, 130 pages, 2012.

J.-H. He, “New interpretation of homotopy perturbation method,” International Journal of Modern Physics B, vol. 20, no. 18, pp. 2561–2568, 2006.

A. A. Elbeleze, A. Kiliçman, and B. M. Taib, “Homotopy perturbation method for fractional black-scholes european option pricing equations using sumudu transform,” Mathematical Problems in Engineering, vol. 2013, Article ID 524852, 7 pages, 2013.

A. R. Ghotbi, H. Bararnia, G. Domairry, and A. Barari, “Investigation of a powerful analytical method into natural convection boundary layer flow,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 5, pp. 2222–2228, 2009.

S. H. Hosein Nia, A. N. Ranjbar, D. D. Ganji, H. Soltani, and J. Ghasemi, “Maintaining the stability of nonlinear differential equations by the enhancement of HPM,” Physics Letters A, vol. 372, no. 16, pp. 2855–2861, 2008.

M. Esmaeilpour and D. D. Ganji, “Application of He's homotopy perturbation method to boundary layer flow and convection heat transfer over a flat plate,” Physics Letters A, vol. 372, no. 1, pp. 33–38, 2007.

S.-L. Pang, “Wavelet numerical method for nonlinear random system,” Transactions of the Chinese Society of Agricultural Machinery, vol. 38, no. 3, pp. 168–170, 2007.

E. Quak and N. Weyrich, “Decomposition and reconstruction algorithms for spline wavelets on a bounded interval,” Applied and Computational Harmonic Analysis, vol. 1, no. 3, pp. 217–231, 1994.

S.-L. Mei and S.-W. Zhang, “Coupling technique of variational iteration and homotopy perturbation methods for nonlinear matrix differential equations,” Computers and Mathematics with Applications, vol. 54, no. 7-8, pp. 1092–1100, 2007.

C. Cattani, “Connection coefficients of Shannon wavelets,” Mathematical Modelling and Analysis, vol. 11, no. 2, pp. 117–132, 2006.

C. Cattani, “Shannon wavelets for the solution of integrodifferential equations,” Mathematical Problems in Engineering, vol. 2010, Article ID 408418, 22 pages, 2010.

D. K. Hoffman, G. W. Wei, D. S. Zhang, and D. J. Kouri, “Shannon-Gabor wavelet distributed approximating functional,” Chemical Physics Letters, vol. 287, no. 1-2, pp. 119–124, 1998.

S. Dubuc, “Interpolation through an iterative scheme,” Journal of Mathematical Analysis and Applications, vol. 114, no. 1, pp. 185–204, 1986.

G. Deslauriers and S. Dubuc, “Symmetric iterative interpolation processes,” Constructive Approximation, vol. 5, no. 1, pp. 49–68, 1989.

I. Daubechies, “Orthonormal bases of compactly supported wavelets,” Communications on Pure and Applied Mathematics, vol. 41, no. 7, pp. 909–996, 1988.

E. Ward Cheney and D. R. Kincaid, Numerical Mathematics and Computing, Thomson Higher Education, Belmont, Calif, USA, 6th edition, 2008.