Mathematical Problems in Engineering

Volume 2013 (2013), Article ID 890157, 22 pages

http://dx.doi.org/10.1155/2013/890157

Research Article

Image Denoising via Nonlinear Hybrid Diffusion

Xiaoping Ji, Dazhi Zhang, Zhichang Guo, and Boying Wu

Department of Mathematics, Harbin Institute of Technology, Harbin 150001, China

Received 3 December 2012; Accepted 15 February 2013

Academic Editor: Jocelyn Sabatier

Copyright © 2013 Xiaoping Ji et al. 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.

A nonlinear anisotropic hybrid diffusion equation is discussed for image denoising, which is a combination of mean curvature smoothing and Gaussian heat diffusion. First, we propose a new edge detection indicator, that is, the diffusivity function. Based on this diffusivity function, the new diffusion is nonlinear anisotropic and forward-backward. Unlike the Perona-Malik (PM) diffusion, the new forward-backward diffusion is adjustable and under control. Then, the existence, uniqueness, and long-time behavior of the new regularization equation of the model are established. Finally, using the explicit difference scheme (PM scheme) and implicit difference scheme (AOS scheme), we do numerical experiments for different images, respectively. Experimental results illustrate the effectiveness of the new model with respect to other known models.

1. Introduction

Image restoration and smoothing are important in problems ranging from medical diagnostic tests to defense applications such as target recognition. Over the past 20 years, the use of variational methods and nonlinear partial differential equations (PDEs) has significantly grown and evolved to address the image restoration problem. Let be the intensity of an image obtained from a noiseless image by adding Gaussian noise with zero mean, defined on a rectangle , and let represent the reconstructed image. The problem is to recover the restoration image , from the observed, noisy image , where the two are related by .

1.1. Nonlinear Diffusion

A large number of image restoration techniques are conveniently formulated using some nonlinear partial differential equations (PDEs) approach. The review article [1] provides a historical description of the use of PDEs in image processing. In [2], Perona and Malik developed an anisotropic diffusion scheme for image denoising. The basic idea of this nonlinear smoothing scheme was to smooth the image while preserving the edges in it. This was done by using equation

where is the noisy image and is the image to be smoothed and describes its evolution over time. The diffusivity controls the amount of diffusion. is also an edge indicator and a smooth nonincreasing function and has such properties as , , and , as . This ensures that strong edges are less blurred by the diffusion filter than noise and low-contrast details.

In [3], Iijima employs the following linear diffusivity function:

Because the model is linear isotropic diffusion, it cannot preserve the edge and some features. The PM diffusivity function [2] is usually

The PM diffusion is nonlinear anisotropic diffusion and can preserve the most features, especially edges in the image. Here are some of the previously employed diffusivity functions.

  Charbonnier diffusivity [4]:

  TV diffusivity [5]:

  Weickert diffusivity [6]:

Except the diffusivity functions, there are other diffusivity functions, such as BFB diffusivity [7] and FAB diffusivity [8, 9]. Well-posedness results are available for linear diffusivity, Charbonnier diffusivity, and TV diffusivity, since they result from convex potentials. For PM diffusivity, Weickert diffusivity, and BFB diffusivity, which can be related to nonconvex potentials, some well-posedness questions are open in the continuous setting [10, 11], while already a space-discretisation creates well-posed processes [12]. In practice, the ill-posedness results in a mild instability in the discrete problem. Regions of high gradients develop a “staircase” instability that involves dynamic coarsening of the steps as time evolves [13, 14].

To make the images more pleasing to the eye, it would be useful to reduce staircasing effect. Many models to reduce this effect have been proposed in the literature. A simple adjustment with practical applications is to include a short range mollifier in the nonlinear diffusion [15]. The new well-posed equation is given by

where is the Gaussian kernel, as described in Section 2. Existence and uniqueness of solutions to this modified Perona-Malik equation have been proved for initial data . Another way is to use a higher-order version of the Perona-Malik equation, examples of which are given in [16–18].

Some authors consider a new class of fractional-order anisotropic diffusion equations to remove the noise [19–27]. These proposed equations can be seen as generalizations of second-order and fourth-order anisotropic diffusion equations. Numerical results show that these methods can not only remove noise and eliminate the staircase effect efficiently in the nontextured region but also preserve the small details such as textures well in the textured region.

1.2. The TV Framework

The famous total variation method first proposed by Rudin et al. [28] consists in solving the following constrained minimization problem:

Here, the first constrain indicates that the noise has zero mean, and the second one uses a priori information that the standard deviation of the noise is . This problem is naturally linked to the unconstrained problem

Mathematically, this is reasonable, since it is natural to study solutions of this problem in the space of functions of bounded variation, , allowing for discontinuities which are necessary for edge reconstruction. The TV model has been studied extensively (see [29–32], et al.) and has proved to be an invaluable tool for preserving edges in image restoration problem. Given the success of TV-based diffusion, various modifications have been introduced. For instance, in [33], Strongand and Chan propose the Adaptive Total Variation model

in which they introduce a control factor, , which slows the diffusion at likely edges. The factor controls the speed of the diffusion and has demonstrated good results as it aids in noise reduction. It is also good at reconstructing edges, since the type of diffusion is the same as that of the original TV model.

The TV model is well posed, but TV-based denoising favors the piecewise constant solutions. Sometimes, this also causes “staircasing effect” [34–41], in which noisy smooth regions are processed into piecewise constant regions (see Figures 3–5). The blocky solution fails to satisfy visual impression and can develop false edges, which can mislead a human or computer into identifying erroneous features, not present in the true image.

Some authors consider another regularizing term to remove the noise [34], which is as follows:

where , and is monotonically decreasing. This model should reap the benefits of both isotropic and TV-based diffusions, as well as a combination of the two. However, it is difficult to study mathematically, since the lower semicontinuity of the functional is not readily evident. In [35], Chen et al. modify the model and propose a functional with variable exponent, , which is a combination of total variation based regularization and Gaussian smoothing.

From the models mentioned above, we can see that based on the PDE framework, the diffusivity functions affect the quality of the reconstructed image. In this paper, based on a new diffusivity function, we propose a new image denoising model which generalizes the approaches due to Perona and Malik [2], Chen et al. [35], and El-Fallah and Ford [42]. In the next section, we will describe this model more precisely. In Section 3, we prove the existence and uniqueness of the proposed model. The theorem can be proved by a similar argument developed in [15], but due to the presence of high degeneration and nonlinearity, more careful estimates are needed. We will give a modified proof in the sections. In the next two section, we first obtain some properties of the weak solution for the new model, and then using these properties, the long-time behavior of the proposed model is established. In Section 6, we describe an iterative scheme which converges to the solution. In the final section, we will give the numerical results which indicate the new model is able to preserve edges and denoise better than the existing methods, for instance, the TV model and PM model.

2. Nonlinear Hybrid Diffusion Model

2.1. The New Model

In this paper, we propose the following model:

where

is the original image, ,  and are fixed constants, is a bounded open domain of with the appropriate smooth boundary, and denotes the unit outward normal to the boundary .

2.2. Hybrid Diffusion

As , , the new divergence operator is changed as follows:

where the last term is the divergence operator of the mean curvature diffusion equation [42]. However, when , , the original divergence operator is changed as follows:

where the last term is the diffusion term of the heat equation.

Hence, the new model has a hybrid diffusion type which combined the mean curvature diffusion with the heat diffusion and has the following advantages.

(i)  Inside the regions where the magnitude of the gradient of is weak, the new model acts like the heat equation, resulting in isotropic smoothing. (ii)  Near the region's boundaries where the magnitude of the gradient is large, the new model acts like the mean curvature equation, resulting in anisotropic smoothing; the regularization is little and the edges are preserved.

2.3. The New Diffusivity Function

Let

where

Now, we discuss the properties of as follows.

Proposition 1. One has the following:

(a)   is a strict decreasing function, and   ,    for       (see Figure 1(a)); (b),   and       (see Figure 1(a)); (c);(d)   and   ; (e).

(a)

(b)

Figure 1: (a) Diffusivity function , for and the PM diffusivity for . (b) The Flux function for and PM flux function for .

Proof. By a direct calculation, we have

which implies (a)–(e).

Remark 2. In terms of the image processing, it is easy to see the following.

(1)From (a) and (b), the edge detection function is like that of the original Perona-Malik diffusion. (2)(c) implies that or as . (3)The diffusion coefficient is dependent on the exponent , which have the similar function with the fractional operator [19–27].

2.4. Forward-Backward Diffusion

For the diffusivity function it follows the new flux function which is defined by

where the variable stands for the norm of the gradient .

In the two-dimensional case, (32) can be replaced by [43]

where we have denoted by and the second derivatives of in the direction which is parallel to and in the orthogonal direction to , respectively:

Remark 3. From Proposition 1(d), we impose

which implies that it is preferable to smooth more in the -direction than in the -direction.

Proposition 4. There exists such that for , and for (see Figure 1(b)). Moreover,

Proof. By a direct calculation, we have

where

for . Then

for ; that is, . It is noticed that, for ,

Because of the continuousness of , there exists a unique point such that , and for , and for . For , that is, ,

From Proposition 4, the new model is of forward diffusion along isophotes (i.e., lines of constant grey value) and of forward-backward diffusion along flow lines (i.e., lines of maximal grey value variation).

Remark 5. (1) From Proposition 4, the threshold value about forward and backward diffusion is estimated as follows:

Therefore, plays the role of a control parameter separating forward from backward diffusion areas.

(2) From Figure 2, we can see, for , the part of the backward diffusion () is not evident; for , the flux function is similar to the PM flux.

Figure 2: Flux function for and PM flux function for .

(a) Noisy: 20, PSNR 22.08

(b) Original

(c) AOS: PSNR 39.00, MAE 1.86

(d) PMS: PSNR 40.70, MAE 1.62

(e) PM: PSNR 39.74, MAE 1.78

(f) TV: PSNR 37.08, MAE 2.52

Figure 3: Synthetic image () (a) Noisy image corrupted by Gaussian noise for . (b) Original image. (c) Our algorithm by AOS, , (11 steps). (d) Our algorithm by PMS, (62 steps). (e) PM algorithm, , (90 steps). (f) TV algorithm, (360 steps).

(a) Noisy: 35, PSNR 17.19

(b) Original

(c) AOS: PSNR 35.43, MAE 2.75

(d) PMS: PSNR 36.61, MAE 2.43

(e) PM: PSNR 34.81, MAE 2.76

(f) TV: PSNR 34.05, MAE 3.65

Figure 4: Synthetic image (). (a) Noisy image corrupted by Gaussian noise for . (b) Original image. (c) Our algorithm by AOS, , (11 steps). (d) Our algorithm by PMS, (115 steps). (e) PM algorithm, , (130 steps). (f) TV algorithm, (600 steps).

(a) Noisy: 50, PSNR 14.11

(b) Original

(c) AOS: PSNR 32.94, MAE 3.48

(d) PMS: PSNR 33.72, MAE 3.30

(e) PM: PSNR 31.90, MAE 3.66

(f) TV: PSNR 31.36, MAE 4.88

Figure 5: Synthetic image (). (a) Noisy image corrupted by Gaussian noise for . (b) Original image. (c) Our algorithm by AOS, , (11 steps). (d) Our algorithm by PMS, (115 steps). (e) PM algorithm, , (130 steps). (f) TV algorithm, (600 steps).

2.5. The Modified Regularization Equation

Using the similar skill in [15], the new model can be regularized by

where is the Gaussian kernel, namely,

This small amount of linear filtering allows to measure edges of in a more “global” sense, so that it is not easily affected by local discretization. It is noticed that equation (32)–(34) is forward diffusion. In [15], while use of the mollifier may seem to be counterproductive, since the original intention was to avoid the blurring caused by linear filtering, the results can be quite impressive and are in fact a great improvement over linear filtering. In the new model, the forward-backward diffusion under control by the factor , and therefore, we do not use this skill in numerical Experiments.

3. Existence and Uniqueness of Weak Solutions

In this section, we establish the existence and uniqueness of weak solutions of the proposed model (32)–(34) following the arguments in [15, 43, 44].

The standard notations are used throughout. We denote by , a positive integer, the set of all functions defined in such that and its distributional derivatives of order all belong to . is a Hilbert space with the norm

The space consists of all functions such that, for almost every in , belongs to . is a normed space with the norm

We denote by the dual of .

We introduce the solution space of the problem (32)–(34) as follows:

Obviously, is a Banach space with the norm

The solutions considered here are in the following weak sense.

Definition 6. A function is called a weak solution of the problem (32)–(34), if satisfies (32) and conditions (33) and (34) a.e. with derivatives of in the sense of distributions.

We will show the existence of weak solutions by the Schauder fixed point theorem. For this purpose, we need to discuss the corresponding linearized problem

Proposition 7. For any , the problem (40)–(42) admits a unique weak solution .

By classical theory, Proposition 7 can be proved by the Galerkin method (see [29, 31] for details).

Now, the theorem for the existence and uniqueness of weak solutions is stated as follow.

Theorem 8. Let and is appropriately small. Then the problem (32)–(34) admits one and only one weak solution such that , .

Proof. Firstly, we consider the proof of the existence, which is based on the Schauder fixed point argument. Let such that

We consider the following linear problem :

for all , a.e. . Since and satisfy (43), then and belong to and there exists a constant such that and a.e. . Since is decreasing and positive, it follows that a.e. in :

Then by applying Proposition 7 on the linearized problem, we will prove that the problem has a unique solution satisfying the estimates

where is the constant depending only on the constant , , and . Choosing in , integrating over the interval , we arrive to the inequality

which implies (46). Choosing in , integrating by parts yields

Integrating over the interval we arrive to that

From Proposition 1 and (48), noticing that , we can deduce that

Since is small, letting yields (45) and (47). From (45)–(47), we introduce the subspace of defined by

By construction, is a mapping from into . Moreover, one can prove that is not empty, convex, and weakly compact in .

In order to use the Schauder fixed point theorem, we need to prove that the mapping is weakly continuous from into . Let be a sequence that converges weakly to some in and let . We have to prove that converges weakly to . From (45)–(47), and classical results of compact inclusion in Sobolev spaces [45], we can extract from , respectively, from , a subsequence such that, for some , we have

The above convergence allows us to pass to the limit in the problem () and obtain . Moreover, since the solution is unique, the whole sequence converges weakly in to ; that is, is weakly continuous. Consequently, thanks to Schauder's fixed-point theorem, there exists such that . The function solves (32)–(34).

Now, we turn to the proof of the uniqueness, following the idea in [44]. Let and be two weak solutions of (32)–(34). For almost every in and , we have

in the distribution sense, where

Then multiplying the above equality by , integrating over , and using the Neumann boundary condition, we get a.e. ,

Since is decreasing and positive, it follows that a.e. in , , which implies from (56),

Moveover, since , , and are smooth, we have

where is a constant that depends only on , , and . From (58) and by using Young's inequality, we obtain

from which we deduce

Since , using Gronwall's inequality yields

that is, .

Remark 9. Let be the weak solution of problem (32)–(34) obtained in the proof of Theorem 8. Then from the proof we get that , where .

4. Some Properties of Weak Solution

In this section, we first investigate the continuity with respect to initial data of the weak solution for (32)–(34), and then investigate the stability of the weak solution and the maximum principle. According to the uniqueness proof in Theorem 8, we obtain the following theorem.

Theorem 10. Assume is the weak solutions of problem (32)–(34) with the initial data . Then

, where , and is Lebesgue measure of .

Proof. Let be the solutions for problem (32)–(34) with the initial data . For almost every in , we have

in the distribution sense. Integrating over the interval and using the Neumann boundary condition yield

Then, multiplying the above equality (63) by , and integrating over , and integrating by parts yield

Using the following Poincaré-Wirtinger inequality [46, page 148], we have

with the constant . Substituting (68) to (67) yields

Multiplying this inequality by and integrating over the interval we arrive to the inequality

Hence, we obtain the assertion of the theorem.

Next, let us build upon the maximum principle as follows.

Theorem 11. Let be the weak solutions of problem (32)–(34) with the initial data and . Then

Proof. Let , and . Multiply (32) by , where

and integrate over to get

Then

Therefore, is decreasing in , and since

we have that

and so

Multiplying (32) by , a similar argument yields that for all . Equation (71) is followed directly.

5. Behavior as

In this section, we investigate the asymptotic behavior of the weak solution as time tends to infinity and obtain the equilibrium weak solution.

Lemma 12. Let be the weak solutions of problem (32)–(34) with the initial data , and . Then

(i)for all , in ,(ii)there exists a subsequence of , denoted also by itself, such that for all , converges to weakly in and strongly in , and converges to weakly in .

Proof. For all , since are uniform bounded in ;, and is bounded in , and then there exist and the subsequence of which is independent on , and denoted by , such that

where . From Theorem 10, we have

Letting , we can obtain the first part of Lemma 12. For (79), integrate over to get

Letting , we can obtain that

Noticing that , we have

Since weakly in , and therefore, we have

Hence, we obtain the remaining part of the lemma.

Now let us consider the following problem:

Theorem 13. Assume . Then the problem (84)–(86) admits one and only one weak solution such that

Proof. It is clearly that is one solution for the problem (84)–(86).

Next we we turn to the proof of the uniqueness of the solution for the problem (84)–(86). Let and be two weak solutions of (84)–(86). Multiplying (84) by , integrating over , and using the Neumann boundary condition, we get

Using the following Poincaré-Wirtinger inequality, we have

with the constant . Substituting (86) and (89) to (88) yields

Then

That is, .

Theorem 14. let be the weak solution of problem (32)–(34), Then when , tends to be steady-state solution , that is, the solution for Problem (84)–(86).

Proof. Let

Then is the weak solutions of problem (32)–(34) with the initial data . From Lemma 12, we obtain there exists a subsequence of , denoted also by itself, such that, for all ,

which implies that as , tends to be steady-state solution , which is the unique solution for problem (84)–(86).

Remark 15. From Theorems 10, 11, and 14, we can observe the following:

(i), which means that no new features are introduced in the image in process.(ii), the mean of , is constant .(iii) tends to zero, which means that converges in the -strong topology to the average of the initial data.

6. Convergent Iterative Scheme

A convergent iterative scheme for (32) is given in this section.

Theorem 16. Let . The sequence defined by solving the iterative scheme

converges in to the strong solution of (32)–(34).

Proof. We denote by . By Proposition 7, problem (94) has a unique solution . It is clear that

Now we verify that the sequence converges in to , the strong solution of (32)–(34).

As in Section 3, from the estimate (60), we have

Moreover, we have

where is a constant which only depends on . Then Gronwall's inequality yields, for any :

where

By (96) and (98), we can deduce

and thus,

Finally, we obtain by iterating

which implies that the sequence converges in to the strong solution of (32)–(34).

7. Numerical Implementation

We present in this section some numerical examples illustrating the capability of our model. We also compare it with the known models (PM and TV). In the next two sections, two numerical discrete schemes, the PM scheme (PMS) and the AOS scheme, will be proposed.

7.1. The PM Scheme

To discretize (12), the finite difference scheme in [2] is used. Denote the space step by and the time step by . Thus, we have

The numerical algorithms for problems (12)–(14) are given in the following:

where for the numerical scheme to be stable.

7.2. The AOS Scheme

Using the scheme in [47], (12) can be discretized as

where ,

where is the set of the two neighbors of pixel (boundary pixels have only one neighbor).

AOS schemes with large time steps still reveal average grey value invariance, stability based on extremum principle, Lyapunov functionals, and convergence to a constant steady-state.

7.3. Comparison with Other Methods

For comparison purposes, some very classical noise removal algorithms from the literature are considered, such as the PM Algorithm [2] (see (1)–(3)) and the TV algorithm [28] (see (9)).

The denoising algorithms were tested on three images: a synthetic image ( pixels), a Lena image ( pixels), and a boat image ( pixels). For each image, a noisy observation is generated by adding the original image by Gaussian noise, standard deviation .

Peak-signal-to-noise-ratio (PSNR) and the mean absolute-deviation error (MAE) are used to measure the quality of the restoration results. They are defined as

where and are the original image and the restored image, respectively. The stopping criterion of all methods is set to achieve the maximal PSNR or the best MAE.

For fair comparison, the parameters of PM and TV were tweaked manually to reach their best performance level. In the PM scheme, there are two parameters: the influencing factor and the time step . In the AOS scheme, besides the same influencing factor with PM scheme, the time step can be very large (in general, for the maximal PSNR). Notice that the parameters of our method are very stable with respect to the image. From these experiments, we also observe that the PSNR reaches a maximum rapidly and decreases rapidly. So the steady-state solution is arrived when , but the time evolution may be stopped earlier to achieve an optimal tradeoff between noise removal and edge preservation (the time when the largest PSNR achieves).

The results are depicted in Figures 3–5 for the synthetic image, Figures 6–8 for the Lena image, and Figures 9–11 for the boat image. Our methods do a good job at restoring faint geometrical structures of the images even for high values of ; see for instance the results on the boat image for . Our algorithm performs among the best and even outperforms its competitors most of the time both visually and quantitatively as revealed by the PSNR and MAE values. For TV method, the number of iterations necessary to satisfy the stopping rule rapidly increases when increases. For PM method, the appropriate parameter is necessary.

(a) Noisy: ,

(b) Original

(c) AOS: , MAE 6.08

(d) PMS: PSNR 29.58, MAE 6.15

(e) PM: PSNR 28.81, MAE 6.46

(f) TV: PSNR 29.15, MAE 6.38

Figure 6: Lenna image (). (a) Noisy image corrupted by Gaussian noise for . (b) Original image. (c) Our algorithm by AOS, , (4 steps). (d) Our algorithm by PMS, (47 steps). (e) PM algorithm, , (55 steps). (f) TV algorithm, (182 steps).

(a) Noisy: , PSNR

(b) Original

(c) AOS: PSNR 27.45, MAE 7.77

(d) PMS: PSNR 27.03, MAE 8.08

(e) PM: PSNR 25.55, MAE 8.98

(f) TV: PSNR 26.87, MAE 8.23

Figure 7: Lenna image (). (a) Noisy image corrupted by Gaussian noise for . (b) Original image. (c) Our algorithm by AOS, , (8 steps). (d) Our algorithm by PMS, (96 steps). (e) PM algorithm, , (84 steps). (f) TV algorithm, (400 steps).

(a) Noisy: , PSNR

(b) Original

(c) AOS: PSNR 26.03, MAE 9.11

(d) PMS: PSNR 25.60, MAE 9.46

(e) PM: PSNR 23.81, MAE 10.77

(f) TV: PSNR 24.43, MAE 9.93

Figure 8: Lenna image (). (a) Noisy image corrupted by Gaussian noise for . (b) Original image. (c) Our algorithm by AOS, , (12 steps). (d) Our algorithm by PMS, (154 steps). (e) PM algorithm, , (108 steps). (f) TV algorithm, (650 steps).

(a) Noisy: , PSNR

(b) Original

(c) AOS: PSNR 29.93, MAE 5.73

(d) PMS: PSNR 29.75, MAE 5.76

(e) PM: PSNR 28.64, MAE 6.22

(f) TV: PSNR 29.13, MAE 5.97

Figure 9: Boat image (). (a) Noisy image corrupted by Gaussian noise for . (b) Original image. (c) Our algorithm by AOS, , (4 steps). (d) Our algorithm by PMS, (47 steps). (e) PM algorithm, , (55 steps). (f) TV algorithm, (194 steps).

(a) Noisy: , PSNR

(b) Original

(c) AOS: PSNR 27.59, MAE 7.41

(d) PMS: PSNR 27.38, MAE 7.45

(e) PM: PSNR 25.85, MAE 8.22

(f) TV: PSNR 26.50, MAE 7.98

Figure 10: Boat image (). (a) Noisy image corrupted by Gaussian noise for . (b) Original image. (c) Our algorithm by AOS, , (9 steps). (d) Our algorithm by PMS, (110 steps). (e) PM algorithm, , (90 steps). (f) TV algorithm, (450 steps).

(a) Noisy: , PSNR

(b) Original

(c) AOS: PSNR 26.22, MAE 8.59

(d) PMS: PSNR 26.02, MAE 8.68

(e) PM: PSNR 24.17, MAE 9.93

(f) TV: PSNR 25.63, MAE 8.98

Figure 11: Boat image (). (a) Noisy image corrupted by Gaussian noise for . (b) Original image. (c) Our algorithm by AOS, , (15 steps). (d) Our algorithm by PMS, (170 steps). (e) PM algorithm, , (115 steps). (f) TV algorithm, (710 steps).

Figures 3, 4, and 5 illustrate the proposed model is able to reconstruct sharp edges and nonuniform regions while avoiding staircasing. TV-based diffusion reconstructs sharp edges, but the staircasing effect is clear evidence. PM-based diffusion also reconstructs sharp edges but creates isolated black and white speckles in the denoise image. The proposed model reconstructs sharp edges as effectively as PM-based diffusion and recovers smooth regions as effectively as pure isotropic diffusion (in particular, without staircasing). The denoising performance results are tabulated in Table 1 where the best PSNR and MAE value is shown in boldface. The PSNR improvement brought by our approach can be quite high particularly for (see, e.g., Figures 5, 8, and 11) and the visual resolution is quite respectable. But even for , the PSNR of our algorithm can be higher than that of PM and TV methods.

Table 1: PSNR, MAE, and CPU time (seconds) of all methods.

PSNR

MAE
CPU time(s)

The synthetic image ()

20
35
50

20
35
55

20
35
55

AOS
39.00
35.43
32.94
AOS
1.86
2.75
3.48
AOS
0.43
0.41
0.54

PMS
40.70
36.61
33.72
PMS
1.65
2.43
3.30
PMS
1.33
2.65
3.91

PM
39.74
34.81
31.90
PM
1.78
2.76
3.66
PM
0.70
1.04
1.82

TV
37.08
34.05
31.36
TV
2.52
3.65
4.88
TV
4.41
7.57
11.15

The Lena image ()

20
35
50

20
35
50

20
35
50

AOS
29.72
27.44
26.03
AOS
6.08
7.77
9.11
AOS
0.96
1.94
2.13

PMS
29.58
27.03
25.60
PMS
6.15
8.08
9.46
PMS
5.05
14.80
18.77

PM
28.81
25.55
23.81
PM
6.46
8.98
10.77
PM
2.30
4.78
12.24

TV
29.15
26.87
24.43
TV
6.38
8.23
9.93
TV
14.32
32.68
49.26

The boat image ()

20
35
50

20
35
50

20
35
50

AOS
29.93
27.59
26.22
AOS
5.73
7.41
8.59
AOS
4.02
6.19
9.58

PMS
29.75
27.38
26.02
PMS
5.76
7.45
8.68
PMS
19.76
39.00
60.44

PM
28.64
25.85
24.17
PM
6.22
8.22
9.93
PM
9.24
14.90
18.90

TV
29.13
26.50
25.63
TV
5.97
7.98
8.98
TV
49.51
109.15
172.17

Table 1 summarizes the computational times for all algorithms. From [47], we know the AOS is a high efficient algorithm. It is less than twice the typical effort needed for an explicit scheme, a rather low price for gaining absolute stability. Moreover, the new algorithm by AOS scheme performs high PSNRs on real images (Figures 6, 7, 8, 9, 10 and 11).

8. Conclusions

This work proposes quite an original, efficient method for noise removal. Noise removal is a difficult problem that arises in various applications relevant to active imaging system.

The main ingredients of our method are as follows. (1) Dependent on the diffusivity function , the new model is hybrid diffusion which is combination of mean curvature smoothing and Gaussian heat diffusion. (2) The new diffusion is forward-backward diffusion, but the backward diffusion is under control and the restored image does not create new features. (3) There are less parameters in the new model and the resultant algorithm is insensitive to these parameters. (4) The new model can be performed by AOS scheme, which is very efficient.

Our experimental results demonstrate that the new algorithm is very efficient and the quality of restored images by our method is quite well.

Acknowledgments

The authors would like to express their sincere thanks to the referees for their valuable suggestions in the revision of the paper which contributed greatly to this work. This work was partially supported by the Fundamental Research Funds for the Central Universities (Grants HIT.NSRIF.2011003 and HIT.NSRIF.2012065), the National Science Foundation of China (Grant no. 11271100), the Aerospace Supported Fund, China, under Contract no. 2011-HT-HGD-06, China Postdoctoral Science Foundation funded project, Grant no. 2012M510933, and also the 985 project of Harbin Institute of Technology.

References

V. Caselles, J.-M. Morel, G. Sapiro, and A. Tannenbaum, “Introduction to the special issue on partial differential equations and geometry-driven diffusion in image processing and analysis,” IEEE Transactions on Image Processing, vol. 7, no. 3, pp. 269–273, 1998.

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.

T. Iijima, “Basic theory on normalization of pattern (in case of typical onedimensional pattern),” Bulletin of the Electrotechnical Laboratory, vol. 26, pp. 368–388, 1962 (Japanese).

P. Charbonnier, L. Blanc-Feraud, G. Aubert, and M. Barlaud, “Two deterministic half-quadratic regularization algorithms for computed imaging,” in Proceedings of IEEE International Conference on Image Processing, vol. 2, pp. 168–172, IEEE Computer Society Press, Austin, Tex, USA, November 1994.

J. Weickert, Anisotropic Diffusion in Image Processing, European Consortium for Mathematics in Industry, B. G. Teubner, Stuttgart, Germany, 1998.

F. Andreu, C. Ballester, V. Caselles, and J. M. Mazón, “Minimizing total variation flow,” Differential and Integral Equations, vol. 14, no. 3, pp. 321–360, 2001.

S. L. Keeling and R. Stollberger, “Nonlinear anisotropic diffusion filtering for multiscale edge enhancement,” Inverse Problems, vol. 18, no. 1, pp. 175–190, 2002.

G. Gilboa, N. Sochen, and Y. Y. Zeevi, “Forward-and-backward diffusion processes for adaptive image enhancement and denoising,” IEEE Transactions on Image Processing, vol. 11, no. 7, pp. 689–703, 2002.

B. Smolka, “Combined forward and backward anisotropic diffusion filtering of color images,” in Pattern Recognition, L. van Gool, Ed., vol. 2449 of Lecture Notes in Computer Science, pp. 314–320, Springer, Berlin, Germany, 2002.

S. Kichenassamy, “The Perona-Malik paradox,” SIAM Journal on Applied Mathematics, vol. 57, no. 5, pp. 1328–1342, 1997.

B. Kawohl and N. Kutev, “Maximum and comparison principle for one-dimensional anisotropic diffusion,” Mathematische Annalen, vol. 311, no. 1, pp. 107–123, 1998.

J. Weickert and B. Benhamouda, “A semidiscrete nonlinear scale-space theory and its relation to the Perona-Malik paradox,” in Advances in Computer Vision, F. Solina, W. G. Kropatsch, R. Klette, and R. Bajcsy, Eds., pp. 1–10, Springer, Wien, Austria, 1997.

S. Esedoglu, An analysis of the Perona-Malik Scheme [Ph.D. thesis], Courant Institute of Mathematical Sciences, New York, NY, USA, 2000.

S. Kichenassamy, “The Perona-Malik paradox,” SIAM Journal on Applied Mathematics, vol. 57, no. 5, pp. 1328–1342, 1997.

L. Alvarez, P.-L. Lions, and J.-M. Morel, “Image selective smoothing and edge detection by nonlinear diffusion. II,” SIAM Journal on Numerical Analysis, vol. 29, no. 3, pp. 845–866, 1992.

J. Tumblin and G. Turk, “LCIS: a boundary hierarchy for detail-preserving contrast reduction,” in Proceedings of the SIGGRAPH Annual Conference on Computer Graphics, pp. 83–90, Los Angeles, Calif, USA, August 1999.

G. W. Wei, “Generalized Perona-Malik equation for image restoration,” IEEE Signal Processing Letters, vol. 6, no. 7, pp. 165–167, 1999.

Y. L. You and M. Kaveh, “Fourth-order partial differential equations for noise removal,” IEEE Transactions on Image Processing, vol. 9, no. 10, pp. 1723–1730, 2000.

B. Mathieu, P. Melchior, A. Oustaloup, and C. Ceyral, “Fractional differentiation for edge detection,” Signal Processing, vol. 83, no. 11, pp. 2421–2432, 2003.

R. W. Ibrahim, “On generalized Srivastava-Owa fractional operators in the unit disk,” Advances in Difference Equations, vol. 2011, pp. 1–12, 2011.

Y. Pu, “Fractional calculus approach to texture of digital image,” in Proceedings of the 8th International Conference on Signal Processing (ICSP '06), vol. 2, pp. 15–20, November 2006.

Y. Zhang, Y. Pu, and J. Zhou, “Construction of fractional differential masks based on Riemann-Liouville definition,” Journal of Computational Information Systems, vol. 6, no. 10, pp. 3191–3200, 2010.

D. Chen, Y. Chen, and D. Xue, “Digital fractional order Savitzky-Golay differentiator,” IEEE Transaction on Circuits and Systems, vol. 11, no. 58, pp. 758–762, 2011.

J. Hu, Y. Pu, and J. Zhou, “Fractional integral denoising algorithm and implementation of fractional integral filter,” Journal of Computational Information Systems, vol. 7, no. 3, pp. 729–736, 2011.

Y. Liu, Y. Pu, and J. Zhou, “Design of image denoising filter based on fractional integral,” Journal of Computational Information Systems, vol. 6, no. 9, pp. 2839–2847, 2010.

J. Hu, Y. Pu, and J. Zhou, “A novel image denoising algorithm based on Riemann-Liouville definition,” Journal of Computers, vol. 6, no. 7, pp. 1332–1338, 2011.

H. A. Jalab and R. W. Ibrahim, “Denoising algorithm based on generalized fractional integral operator with two parameters,” Discrete Dynamics in Nature and Society, vol. 2012, Article ID 529849, 14 pages, 2012.

L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D: Nonlinear Phenomena, vol. 60, no. 1–4, pp. 259–268, 1992.

R. Acar and C. R. Vogel, “Analysis of bounded variation penalty methods for ill-posed problems,” Inverse Problems, vol. 10, no. 6, pp. 1217–1229, 1994.

P. Blomgren, T. F. Chan, P. Mulet, and C. K. Wong, “Total variation image restoration: numerical methods and extensions,” in Proceedings of the International Conference on Image Processing, pp. 384–387, October 1997.

V. Caselles, J.-M. Morel, and C. Sbert, “An axiomatic approach to image interpolation,” IEEE Transactions on Image Processing, vol. 7, no. 3, pp. 376–386, 1998.

A. Chambolle and P.-L. Lions, “Image recovery via total variation minimization and related problems,” Numerische Mathematik, vol. 76, no. 2, pp. 167–188, 1997.

D. M. Strongand and T. F. Chan, “Spatially and scale adaptive total variation based regularization and anisotropic diffusion in image processing,” Tech. Rep. CAM96-46(UCLA), University of California, Los Angeles, Calif, USA, 1996.

P. Blomgren, T. F. Chan, P. Mulet, and C. K. Wong, “Total variation image restoration: numerical methods and extensions,” in Proceedings of the International Conference on Image Processing, pp. 384–387, October 1997.

Y. Chen, S. Levine, and M. Rao, “Variable exponent, linear growth functionals in image restoration,” SIAM Journal on Applied Mathematics, vol. 66, no. 4, pp. 1383–1406, 2006.

A. Chambolle and P.-L. Lions, “Image recovery via total variation minimization and related problems,” Numerische Mathematik, vol. 76, no. 2, pp. 167–188, 1997.

D. C. Dobson and F. Santosa, “Recovery of blocky images from noisy and blurred data,” SIAM Journal on Applied Mathematics, vol. 56, no. 4, pp. 1181–1198, 1996.

M. Nikolova, “Weakly constrained minimization: application to the estimation of images and signals involving constant regions,” Journal of Mathematical Imaging and Vision, vol. 21, no. 2, pp. 155–175, 2004.

W. Ring, “Structural properties of solutions to total variation regularization problems,” Mathematical Modelling and Numerical Analysis, vol. 34, no. 4, pp. 799–810, 2000.

R. T. Whitaker and S. M. Pizer, “A mulit-scale approach to nonuniform diffusion,” Computer Vision, Graphics, and Image Processing: Image Understand, vol. 57, no. 1, pp. 99–110, 1993.

Y. L. You, W. Xu, A. Tannenbaum, and M. Kaveh, “Behavioral analysis of anisotropic diffusion in image processing,” IEEE Transactions on Image Processing, vol. 5, no. 11, pp. 1539–1553, 1996.

A. I. El-Fallah and G. E. Ford, “On mean curvature diffusion in nonlinear image filtering,” Pattern Recognition Letters, vol. 19, no. 5-6, pp. 433–437, 1998.

G. Aubert and P. Kornprobst, Mathematical Problems in Image Processing: PDEs and the Calculus of Variations, vol. 147 of Applied Mathematical Sciences, Springer, New York, NY, USA, 2002.

L. C. Evans, Partial Differential Equations, vol. 19 of Graduate Studies in Mathematics, American Mathematical Society, Providence, RI, USA, 1998.

R. A. Adams, Sobolev Spaces, vol. 65 of Pure and Applied Mathematics Series of Monographs and Textbooks, Academic Press, New York, NY, USA, 1975.

L. Ambrosio, N. Fusco, and D. Pallara, Functions of Bounded Variation and Free Discontinuity Problems, Oxford Mathematical Monographs, Oxford University Press, Oxford, UK, 2000.

J. Weickert, B. M. Ter Haar Romeny, and M. A. Viergever, “Efficient and reliable schemes for nonlinear diffusion filtering,” IEEE Transactions on Image Processing, vol. 7, no. 3, pp. 398–410, 1998.