Image Denoising via Nonlinear Hybrid Diffusion

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 newdiffusion 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 themodel 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.


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  0 be the intensity of an image obtained from a noiseless image by adding Gaussian noise with zero mean, defined on a rectangle Ω ⊂ R 2 , and let  represent the reconstructed image.The problem is to recover the restoration image , from the observed, noisy image  0 , where the two are related by  0 =  + noise.

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 (0) = 1, () ≥ 0, and () → 0, 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 spacediscretisation 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 wellposed equation is given by   = div ( (     ∇  *      ) ∇) , in Ω × (0, ) ,  (0, ) =  0 , in Ω,   ⃗  = 0, on Ω × (0, ) , 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  0 ∈  2 (Ω).Another way is to use a higher-order version of the Perona-Malik equation, examples of which are given in [16][17][18].Some authors consider a new class of fractional-order anisotropic diffusion equations to remove the noise [19][20][21][22][23][24][25][26][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 min Mathematically, this is reasonable, since it is natural to study solutions of this problem in the space of functions of bounded variation, BV(Ω), allowing for discontinuities which are necessary for edge reconstruction.The TV model has been studied extensively (see [29][30][31][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 min 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][35][36][37][38][39][40][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 lim  → 0 ()→2, lim  → ∞ ()→1, 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, 1 ≤ () ≤ 2, 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.

Nonlinear Hybrid Diffusion Model
2.1.The New Model.In this paper, we propose the following model: ) , in Ω × (0, ) , where is the original image,  > 0,  and  > 0 are fixed constants, Ω is a bounded open domain of R  with the appropriate smooth boundary, and ⃗  denotes the unit outward normal to the boundary Ω.

Hybrid Diffusion.
As  → +∞, () → 1, the new divergence operator is changed as follows: div ( ∇ where the last term is the divergence operator of the mean curvature diffusion equation [42].However, when  = 0, () = 2, the original divergence operator is changed as follows: div ( ∇ 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.

The New Diffusivity Function.
Let where Now, we discuss the properties of () as follows.
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.

Forward-Backward Diffusion.
For the diffusivity function () it follows the new flux function F() 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: which implies that it is preferable to smooth more in the direction than in the -direction.
Remark 5. (1) From Proposition 4, the threshold value  0 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  ≥ 1, the part of the backward diffusion (  () < 0) is not evident; for  = 0.02, the flux function is similar to the PM flux..Using the similar skill in [15], the new model can be regularized by

The Modified Regularization Equation
where   () is the Gaussian kernel, namely, This small amount of linear filtering allows (|∇  * | 2 ) 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.

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 || = ∑  =1   ≤  all belong to  2 (Ω).  (Ω) is a Hilbert space with the norm The space  ∞ (0, ; ) consists of all functions  such that, for almost every  in (0, ),  belongs to .  ∞ (0, ; ) is a normed space with the norm We denote by  1 (Ω)  the dual of  1 (Ω).
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  ∈  2 (Ω), 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.
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 V ∈ Then by applying Proposition 7 on the linearized problem, we will prove that the problem (  ) has a unique solution   ∈  satisfying the estimates where  1 is the constant depending only on the constant ],   , and ‖ ‖  1 (Ω) .Choosing V =   in (  ), integrating over the interval (0, ), we arrive to the inequality which implies (46) From Proposition 1 and (48), noticing that |  ( 2 )| ≤ , we can deduce that Since ‖ ‖  1 (Ω) is small, letting /(8]) ≤ 1 yields ( 45) and (47).From ( 45)-( 47), we introduce the subspace  0 of  defined by By construction,  → () ≡   is a mapping from  0 into  0 .Moreover, one can prove that  0 is not empty, convex, and weakly compact in (0, ).
In order to use the Schauder fixed point theorem, we need to prove that the mapping  :  →   is weakly continuous from  0 into  0 .Let   be a sequence that converges weakly to some  in  0 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  0 to  = (); that is,  is weakly continuous.Consequently, thanks to Schauder's fixed-point theorem, there exists  ∈  0 such that  = () =   .The function   solves (32)- (34).Now, we turn to the proof of the uniqueness, following the idea in [44].Let  1 and  2 be two weak solutions of (32)- (34).For almost every  in [0, ] and  = 1, 2, we have where  4 is a constant that depends only on  1 , ], and   .From (58) and by using Young's inequality, we obtain from which we deduce Since  1 (0) =  2 (0) = , using Gronwall's inequality yields that is,  1 =  2 .

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.
with the constant  ≡ (Ω).Substituting (68) to (67) yields Multiplying this inequality by  2]/ and integrating over the interval (0, ) we arrive to the inequality Hence, we obtain the assertion of the theorem.
Next, let us build upon the maximum principle as follows.

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.
Next we we turn to the proof of the of the solution for the problem (84)-(86).Let  1 and  2 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 That is,  1 =  2 .
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 {  } ∞ =1 , denoted also by itself, such that, for all  ≥ 0, 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) inf ∈Ω  ≤  ≤ sup ∈Ω , which means that no new features are introduced in the image in process.(ii)  Ω , the mean of , is constant  Ω .
(iii) ∫ Ω ( −  Ω ) 2  tends to zero, which means that  converges in the  2 (Ω)-strong topology to the average of the initial data.

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

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.
The numerical algorithms for problems ( 12)-( 14) are given in the following: where 0 ≤  ≤ 1/4 for the numerical scheme to be stable.

The AOS Scheme.
Using the scheme in [47], (12) can be discretized as where N() 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 steadystate.

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 (128 × 128 pixels), a Lena image (300 × 300 pixels), and a boat image (512 × 512 pixels).For each image, a noisy observation is generated by adding the original image by Gaussian noise, standard deviation  ∈ {20, 35, 50}.
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  = 0.25.In the AOS scheme, besides the same influencing factor  with PM scheme, the time step  can be very large (in general,  = 2 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  = 50.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.
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  = 50 (see, e.g., Figures 5,8,and 11) and the visual resolution is quite respectable.But even for  = 20, the PSNR of our algorithm can be higher than that of PM and TV methods.
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

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.

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