Mathematical Problems in Engineering

Mathematical Problems in Engineering / 2013 / Article
Special Issue

New Challenges in Fractional Systems

View this Special Issue

Research Article | Open Access

Volume 2013 |Article ID 890157 |

Xiaoping Ji, Dazhi Zhang, Zhichang Guo, Boying Wu, "Image Denoising via Nonlinear Hybrid Diffusion", Mathematical Problems in Engineering, vol. 2013, Article ID 890157, 22 pages, 2013.

Image Denoising via Nonlinear Hybrid Diffusion

Academic Editor: Jocelyn Sabatier
Received03 Dec 2012
Accepted15 Feb 2013
Published09 Apr 2013


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 [1618].

Some authors consider a new class of fractional-order anisotropic diffusion equations to remove the noise [1927]. 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 [2932], 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” [3441], in which noisy smooth regions are processed into piecewise constant regions (see Figures 35). 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) .

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 [1927].

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.

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 35 for the synthetic image, Figures 68 for the Lena image, and Figures 911 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.