# Image Denoising via Nonlinear Hybrid Diffusion

Authors: Xiaoping Ji; Dazhi Zhang; Zhichang Guo; Boying Wu
Journal: Mathematical Problems in Engineering
 (2013)

---

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


---

## Body

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


### 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
(1)ut=div(c(|∇u|)∇u),u(0,x)=f,
where f is the noisy image and u is the image to be smoothed and ut describes its evolution over time. The diffusivity c(|∇u|) controls the amount of diffusion. c(s) is also an edge indicator and a smooth nonincreasing function and has such properties as c(0)=1, c(s)≥0, and c(s)→0, as s→∞. 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:
(2)c(|x|)=1.Because the model is linear isotropic diffusion, it cannot preserve the edge and some features. The PM diffusivity function [2] is usually
(3)cPM(|x|)=11+|x|2/K2.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]:
(4)c(|x|)=11+|x|2/K2.TV diffusivity [5]:
(5)c(|x|)=1|x|.Weickert diffusivity [6]:
(6)c(|x|)={1,|x|=01-exp(-3.31488|x|8/k8),|x|>0.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
(7)∂u∂t=div(c(|∇Gσ*u|)∇u),inΩ×(0,T),u(0,x)=u0,inΩ,∂u∂n→=0,on∂Ω×(0,T),
where Gσ 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 u0∈L2(Ω). 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:
(8)min∫Ω‍|∇u|,∫Ω‍udx=∫Ω‍u0,∫Ω‍|u-u0|22dx=η2.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
(9)minu∈BV(Ω)E(u)=∫Ω‍|∇u|+λ∥u-u0∥L2(Ω)2.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–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
(10)minu∈BV(Ω)Eg(u)=∫Ω‍g(x)|∇u|+λ∥u-u0∥L2(Ω)2,
in which they introduce a control factor, g(x), which slows the diffusion at likely edges. The factor g(x) 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:
(11)infuE(u)=∫Ω‍|∇u|p(|∇u|)dx+λ∥f-u∥L2(Ω)2,
where lims→0p(s)→2,lims→∞p(s)→1, and p 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≤p(x)≤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.


## 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
(1)ut=div(c(|∇u|)∇u),u(0,x)=f,
where f is the noisy image and u is the image to be smoothed and ut describes its evolution over time. The diffusivity c(|∇u|) controls the amount of diffusion. c(s) is also an edge indicator and a smooth nonincreasing function and has such properties as c(0)=1, c(s)≥0, and c(s)→0, as s→∞. 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:
(2)c(|x|)=1.Because the model is linear isotropic diffusion, it cannot preserve the edge and some features. The PM diffusivity function [2] is usually
(3)cPM(|x|)=11+|x|2/K2.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]:
(4)c(|x|)=11+|x|2/K2.TV diffusivity [5]:
(5)c(|x|)=1|x|.Weickert diffusivity [6]:
(6)c(|x|)={1,|x|=01-exp(-3.31488|x|8/k8),|x|>0.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
(7)∂u∂t=div(c(|∇Gσ*u|)∇u),inΩ×(0,T),u(0,x)=u0,inΩ,∂u∂n→=0,on∂Ω×(0,T),
where Gσ 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 u0∈L2(Ω). 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:
(8)min∫Ω‍|∇u|,∫Ω‍udx=∫Ω‍u0,∫Ω‍|u-u0|22dx=η2.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
(9)minu∈BV(Ω)E(u)=∫Ω‍|∇u|+λ∥u-u0∥L2(Ω)2.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–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
(10)minu∈BV(Ω)Eg(u)=∫Ω‍g(x)|∇u|+λ∥u-u0∥L2(Ω)2,
in which they introduce a control factor, g(x), which slows the diffusion at likely edges. The factor g(x) 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:
(11)infuE(u)=∫Ω‍|∇u|p(|∇u|)dx+λ∥f-u∥L2(Ω)2,
where lims→0p(s)→2,lims→∞p(s)→1, and p 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≤p(x)≤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.


## 2. Nonlinear Hybrid Diffusion Model

### 2.1. The New Model

In this paper, we propose the following model:(12)∂u∂t=div(∇u(|∇u|2+1)(2-p(|∇u|2))/2),inΩ×(0,T),(13)u(0,x)=f,inΩ,(14)∂u∂n→=0,on∂Ω×(0,T),
where
(15)p(|∇u|2)=1+11+k|∇u|2,k>0,σ>0,f is the original image, k>0,  σ and T>0 are fixed constants, Ω is a bounded open domain of ℝN with the appropriate smooth boundary, and n→ denotes the unit outward normal to the boundary ∂Ω.


### 2.2. Hybrid Diffusion

Ass→+∞, p(s)→1, the new divergence operator is changed as follows:
(16)div(∇u(|∇u|2+1)(2-p(|∇u|))/2)⟶div(∇u|∇u|2+1),
where the last term is the divergence operator of the mean curvature diffusion equation [42]. However, when s=0, p(s)=2, the original divergence operator is changed as follows:
(17)div(∇u(|∇u|2+1)(2-p(|∇u|))/2)⟶Δu,
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 ofu 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(18)C(s)=(1+s)(p(s)-2)/2,
where
(19)p(s)=1+11+ks.Now, we discuss the properties ofC(s) as follows.Proposition 1.
One has the following:(a)
C:[0,+∞)→(0,+∞)   is a strict decreasing function, and    0≤C(s)≤1,    for   s∈ℝ¯+   (see Figure1(a));(b)
C(0)=1,   and   lims→∞C(s)=0    (see Figure1(a));(c)
lims→+∞sC′(s)C(s)=-1/2;(d)
lims→+∞1+sC(s)=1   and   lims→+∞sC(s2)=1;(e)
lims→+∞sC′(s2)=k/2.(a) Diffusivity functionC(s), for  k=0.01,0.1,10 and the PM diffusivity for k=5. (b) The Flux function ℱ(s) for  k=0.01,0.1,10 and PM flux function for k=5.
(a)
(b)Proof.
By a direct calculation, we have(20)C(s)=exp{-ksln(1+s)2(1+ks)},C′(s)=-kC(s)s(ks+1)+(1+s)ln(1+s)2(1+s)(1+ks)2<0,
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 functionC(s) is like that of the original Perona-Malik diffusion.(2)
(c) implies thatC(s)≈1/s or C(s)≈1/1+s as s→+∞.(3)
The diffusion coefficientC(s) is dependent on the exponent p(s), which have the similar function with the fractional operator [19–27].


### 2.4. Forward-Backward Diffusion

For the diffusivity functionC(s) it follows the new flux function ℱ(s) which is defined by
(21)ℱ(s):=sC(s2),
where the variable s stands for the norm of the gradient |∇u|.In the two-dimensional case,  (32) can be replaced by [43]
(22)∂tu=div(C(|∇u|2)∇u)=C(|∇u|2)uTT+ℱ′(∇u)uNN,
where we have denoted by uNN and uTT the second derivatives of u in the direction N(x,y)=∇u/|∇u| which is parallel to ∇u and T(x,y) in the orthogonal direction to N(x,y), respectively:
(23)uNN=ux2uxx+uy2uyy+2uxuyuxy|∇u|2,uTT=ux2uyy+uy2uxx-2uxuyuxy|∇u|2.Remark 3.
From Proposition1(d), we impose
(24)lims→∞ℱ′(s)C(s2)=0,
which implies that it is preferable to smooth more in the T-direction than in the N-direction.Proposition 4.
There existss0∈(0,+∞) such that ℱ′(s)≥0 for |s|≤s0, and ℱ′(s)<0 for |s|>s0 (see Figure 1(b)). Moreover,
(25)s0>ek-1.Proof.
By a direct calculation, we have(26)ℱ′(s)=C(s2)+2s2C′(s2)=C(s2)h(s2)(1+s2)(1+ks2)2,
where
(27)h(s)=(k2+k)s2+(2k+1)s+1-ks(1+s)ln(1+s),
for s>0. Then
(28)h′(s)=(2k2+k)s+(2k+1)-k(2s+1)ln(1+s),h′′(s)=k((2k-1)-2ln(1+s)+11+s)<2k(k-ln(1+s))<0,
for k<ln(1+s); that is, s>ek-1. It is noticed that, for s=ek-1,
(29)h(s)=ks2+(2k+1-k2)s+1>0,lims→+∞h(s)=-∞.Because of the continuousness ofg(s), there exists a unique point s0∈(0,+∞) such that h(s0)=0, and h(s)≥0 for ek-1<s≤s0, and h(s)<0 for s>s0. For k>ln(1+s), that is,  0<s<ek-1,
(30)h(s)=(k2+k)s2+(2k+1)s+1-ks(1+s)ln(1+s)>ks2+(2k+1-k2)s+1>0.From Proposition4, 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 Proposition4, the threshold value s0 about forward and backward diffusion is estimated as follows:
(31)s0>ek-1.
Therefore, k plays the role of a control parameter separating forward from backward diffusion areas.
(2) From Figure2, we can see, for  k≥1, the part of the backward diffusion (F′(s)<0) is not evident; for k=0.02, the flux function is similar to the PM flux.Figure 2
Flux functionℱ(s) for  k=0.01,0.02,0.1,1,10 and PM flux function for k=5.Synthetic image (128×128) (a) Noisy image corrupted by Gaussian noise for σ=20. (b) Original image. (c) Our algorithm by AOS, k=0.02, τ=2 (11 steps). (d) Our algorithm by PMS, τ=0.25 (62 steps). (e) PM algorithm, K=5, τ=0.25 (90 steps). (f) TV algorithm, τ=0.1 (360 steps).
(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.52Synthetic image (128×128). (a) Noisy image corrupted by Gaussian noise for σ=35. (b) Original image. (c) Our algorithm by AOS, k=0.02, τ=2 (11 steps). (d) Our algorithm by PMS, τ=0.25 (115 steps). (e) PM algorithm, K=7, τ=0.25 (130 steps). (f) TV algorithm, τ=0.1 (600 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.65Synthetic image (128×128). (a) Noisy image corrupted by Gaussian noise for σ=50. (b) Original image. (c) Our algorithm by AOS, k=0.02, τ=2 (11 steps). (d) Our algorithm by PMS, τ=0.25 (115 steps). (e) PM algorithm, K=7, τ=0.25 (130 steps). (f) TV algorithm, τ=0.1 (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


### 2.5. The Modified Regularization Equation

Using the similar skill in [15], the new model can be regularized by
(32)∂u∂t=div(C(|∇Gσ*u|2)∇u),inΩ×(0,T),(33)u(0,x)=f,inΩ,(34)∂u∂n→=0,on∂Ω×(0,T),
where Gσ(x) is the Gaussian kernel, namely,
(35)Gσ(x)=1(4πσ)N/2exp(-|x|24σ).This small amount of linear filtering allowsC(|∇Gσ*u|2) to measure edges of u 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 k, and therefore, we do not use this skill in numerical Experiments.


## 2.1. The New Model

In this paper, we propose the following model:(12)∂u∂t=div(∇u(|∇u|2+1)(2-p(|∇u|2))/2),inΩ×(0,T),(13)u(0,x)=f,inΩ,(14)∂u∂n→=0,on∂Ω×(0,T),
where
(15)p(|∇u|2)=1+11+k|∇u|2,k>0,σ>0,f is the original image, k>0,  σ and T>0 are fixed constants, Ω is a bounded open domain of ℝN with the appropriate smooth boundary, and n→ denotes the unit outward normal to the boundary ∂Ω.


## 2.2. Hybrid Diffusion

Ass→+∞, p(s)→1, the new divergence operator is changed as follows:
(16)div(∇u(|∇u|2+1)(2-p(|∇u|))/2)⟶div(∇u|∇u|2+1),
where the last term is the divergence operator of the mean curvature diffusion equation [42]. However, when s=0, p(s)=2, the original divergence operator is changed as follows:
(17)div(∇u(|∇u|2+1)(2-p(|∇u|))/2)⟶Δu,
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 ofu 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(18)C(s)=(1+s)(p(s)-2)/2,
where
(19)p(s)=1+11+ks.Now, we discuss the properties ofC(s) as follows.Proposition 1.
One has the following:(a)
C:[0,+∞)→(0,+∞)   is a strict decreasing function, and    0≤C(s)≤1,    for   s∈ℝ¯+   (see Figure1(a));(b)
C(0)=1,   and   lims→∞C(s)=0    (see Figure1(a));(c)
lims→+∞sC′(s)C(s)=-1/2;(d)
lims→+∞1+sC(s)=1   and   lims→+∞sC(s2)=1;(e)
lims→+∞sC′(s2)=k/2.(a) Diffusivity functionC(s), for  k=0.01,0.1,10 and the PM diffusivity for k=5. (b) The Flux function ℱ(s) for  k=0.01,0.1,10 and PM flux function for k=5.
(a)
(b)Proof.
By a direct calculation, we have(20)C(s)=exp{-ksln(1+s)2(1+ks)},C′(s)=-kC(s)s(ks+1)+(1+s)ln(1+s)2(1+s)(1+ks)2<0,
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 functionC(s) is like that of the original Perona-Malik diffusion.(2)
(c) implies thatC(s)≈1/s or C(s)≈1/1+s as s→+∞.(3)
The diffusion coefficientC(s) is dependent on the exponent p(s), which have the similar function with the fractional operator [19–27].


## 2.4. Forward-Backward Diffusion

For the diffusivity functionC(s) it follows the new flux function ℱ(s) which is defined by
(21)ℱ(s):=sC(s2),
where the variable s stands for the norm of the gradient |∇u|.In the two-dimensional case,  (32) can be replaced by [43]
(22)∂tu=div(C(|∇u|2)∇u)=C(|∇u|2)uTT+ℱ′(∇u)uNN,
where we have denoted by uNN and uTT the second derivatives of u in the direction N(x,y)=∇u/|∇u| which is parallel to ∇u and T(x,y) in the orthogonal direction to N(x,y), respectively:
(23)uNN=ux2uxx+uy2uyy+2uxuyuxy|∇u|2,uTT=ux2uyy+uy2uxx-2uxuyuxy|∇u|2.Remark 3.
From Proposition1(d), we impose
(24)lims→∞ℱ′(s)C(s2)=0,
which implies that it is preferable to smooth more in the T-direction than in the N-direction.Proposition 4.
There existss0∈(0,+∞) such that ℱ′(s)≥0 for |s|≤s0, and ℱ′(s)<0 for |s|>s0 (see Figure 1(b)). Moreover,
(25)s0>ek-1.Proof.
By a direct calculation, we have(26)ℱ′(s)=C(s2)+2s2C′(s2)=C(s2)h(s2)(1+s2)(1+ks2)2,
where
(27)h(s)=(k2+k)s2+(2k+1)s+1-ks(1+s)ln(1+s),
for s>0. Then
(28)h′(s)=(2k2+k)s+(2k+1)-k(2s+1)ln(1+s),h′′(s)=k((2k-1)-2ln(1+s)+11+s)<2k(k-ln(1+s))<0,
for k<ln(1+s); that is, s>ek-1. It is noticed that, for s=ek-1,
(29)h(s)=ks2+(2k+1-k2)s+1>0,lims→+∞h(s)=-∞.Because of the continuousness ofg(s), there exists a unique point s0∈(0,+∞) such that h(s0)=0, and h(s)≥0 for ek-1<s≤s0, and h(s)<0 for s>s0. For k>ln(1+s), that is,  0<s<ek-1,
(30)h(s)=(k2+k)s2+(2k+1)s+1-ks(1+s)ln(1+s)>ks2+(2k+1-k2)s+1>0.From Proposition4, 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 Proposition4, the threshold value s0 about forward and backward diffusion is estimated as follows:
(31)s0>ek-1.
Therefore, k plays the role of a control parameter separating forward from backward diffusion areas.
(2) From Figure2, we can see, for  k≥1, the part of the backward diffusion (F′(s)<0) is not evident; for k=0.02, the flux function is similar to the PM flux.Figure 2
Flux functionℱ(s) for  k=0.01,0.02,0.1,1,10 and PM flux function for k=5.Synthetic image (128×128) (a) Noisy image corrupted by Gaussian noise for σ=20. (b) Original image. (c) Our algorithm by AOS, k=0.02, τ=2 (11 steps). (d) Our algorithm by PMS, τ=0.25 (62 steps). (e) PM algorithm, K=5, τ=0.25 (90 steps). (f) TV algorithm, τ=0.1 (360 steps).
(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.52Synthetic image (128×128). (a) Noisy image corrupted by Gaussian noise for σ=35. (b) Original image. (c) Our algorithm by AOS, k=0.02, τ=2 (11 steps). (d) Our algorithm by PMS, τ=0.25 (115 steps). (e) PM algorithm, K=7, τ=0.25 (130 steps). (f) TV algorithm, τ=0.1 (600 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.65Synthetic image (128×128). (a) Noisy image corrupted by Gaussian noise for σ=50. (b) Original image. (c) Our algorithm by AOS, k=0.02, τ=2 (11 steps). (d) Our algorithm by PMS, τ=0.25 (115 steps). (e) PM algorithm, K=7, τ=0.25 (130 steps). (f) TV algorithm, τ=0.1 (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


## 2.5. The Modified Regularization Equation

Using the similar skill in [15], the new model can be regularized by
(32)∂u∂t=div(C(|∇Gσ*u|2)∇u),inΩ×(0,T),(33)u(0,x)=f,inΩ,(34)∂u∂n→=0,on∂Ω×(0,T),
where Gσ(x) is the Gaussian kernel, namely,
(35)Gσ(x)=1(4πσ)N/2exp(-|x|24σ).This small amount of linear filtering allowsC(|∇Gσ*u|2) to measure edges of u 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 k, 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 byHk(Ω), k a positive integer, the set of all functions u defined in Ω such that u and its distributional derivatives ∂mu/∂xm of order |m|=∑j=1N‍mj≤k all belong to L2(Ω). Hk(Ω) is a Hilbert space with the norm
(36)∥u∥Hk(Ω)=(∑|m|≤k‍∫Ω‍|∂mu∂xm|2dx)1/2.
The space L∞(0,T;X) consists of all functions u such that, for almost every t in (0,T), u belongs to X. L∞(0,T;X) is a normed space with the norm
(37)∥u∥L∞(0,T;X)=esssup0≤t≤T∥u(·,t)∥X.
We denote by H1(Ω)′ the dual of H1(Ω).We introduce the solution spaceW of the problem (32)–(34) as follows:
(38)W(0,T)={w∈L∞(0,T;H1(Ω));dwdt∈L2(QT)}.
Obviously, W is a Banach space with the norm
(39)∥w∥W=∥w∥L∞(0,T;H1(Ω))+∥∂w∂t∥L2(QT).
The solutions considered here are in the following weak sense.Definition 6.
A functionu is called a weak solution of the problem (32)–(34), if u∈Wsatisfies (32) and conditions (33) and (34) a.e. with derivatives of u 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(40)∂u∂t=div(C(|∇Gσ*w|2)∇u),inΩ×(0,T),(41)u(0,x)=f,inΩ,(42)∂u∂n→=0,on∂Ω×(0,T).Proposition 7.
For anyf∈L2(Ω), the problem (40)–(42) admits a unique weak solution u∈W.By classical theory, Proposition7 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.
Letf∈H1(Ω) and ∥f∥H1(Ω) is appropriately small. Then the problem (32)–(34) admits one and only one weak solution u(x,t) such that ∂u/∂t∈L2(QT), u∈L∞(0,T;H1(Ω)∩u∈C([0,T],L2(Ω)).Proof.
Firstly, we consider the proof of the existence, which is based on the Schauder fixed point argument. Letw∈W such that
(43)∥w∥L∞(0,T;L2(Ω))≤∥f∥L2(Ω),∥∂w∂t∥L2(QT)≤∥f∥H1(Ω).
We consider the following linear problem (Pw):
(Pw)〈du(t)dt,v〉H1(Ω)′×H1(Ω)+∫Ω‍C(|∇Gσ*w|2)∇u(t)∇v(t)dx=0
for all v∈H1(Ω), a.e. t∈[0,T]. Since w and ∂w/∂t satisfy (43), then |∇Gσ*w| and |∇Gσ*(∂w/∂t)| belong to L∞((0,T);C∞(Ω)) and there exists a constant M=M(Gσ,∥f∥H1(Ω)) such that |∇Gσ*w|≤M and |∇Gσ*(∂w/∂t)|≤M a.e. t,forallx∈Ω. Since C(s) is decreasing and positive, it follows that a.e. in (0,T)×Ω:
(44)1≥C(|∇Gσ*w|2)≥C(M2)=ν.
Then by applying Proposition 7 on the linearized problem, we will prove that the problem (Pw)  has a unique solution uw∈W satisfying the estimates
(45)∥uw∥L∞(0,T;H1(Ω))≤c1,(46)∥uw∥L∞(0,T;L2(Ω))≤∥f∥L2(Ω),(47)∥∂uw∂t∥L2(QT)≤∥f∥H1(Ω),
where c1 is the constant depending only on the constant ν, Gσ, and ∥f∥H1(Ω). Choosing v=uw in (Pw), integrating over the interval (0,t), we arrive to the inequality
(48)12∫Ω‍uw2dx+ν∫0t‍∫Ω‍|∇uw|2dxds≤12∫Ω‍f2dx,
which implies (46). Choosing v=∂uw/∂t in (Pw), integrating by parts yields
(49)∫Ω‍(∂uw∂t)2dx+12∫Ω‍C(|∇Gσ*w|2)∂|∇uw|2∂tdx=0.
Integrating over the interval (0,t) we arrive to that
(50)∫0t‍∫Ω‍(∂uw∂s)2dxds+12∫Ω‍C(|∇Gσ*w|2)|∇uw|2dx=12∫Ω‍C(|∇Gσ*w|2)|∇f|2dx+12∫0t‍∫Ω‍C′(|∇Gσ*w|2)×∇(Gσ*w)∇(Gσ*∂w∂s)|∇uw|2dxds.
From Proposition 1 and (48), noticing that |C′(s2)s|≤k, we can deduce that
(51)∫0t‍∫Ω‍(∂uw∂s)2dxds+ν2∫Ω‍|∇u|2dx≤∫Ω‍|∇f|2dx+Mk4ν∫Ω‍f2dxds.
Since ∥f∥H1(Ω) is small, letting Mk/(8ν)≤1 yields (45) and (47). From (45)–(47), we introduce the subspace W0 of W defined by
(52)W0={∥dwdt∥L2(QT)w∈W(0,T),w(0)=fa∥w∥L∞(0,T;H1(Ω))≤c1,a∥w∥L∞(0,T;L2(Ω))≤∥f∥L2(Ω),a∥dwdt∥L2(QT)≤∥f∥H1(Ω)}.
By construction, w→S(w)≡uw is a mapping from W0 into W0. Moreover, one can prove that W0 is not empty, convex, and weakly compact in W(0,T).
In order to use the Schauder fixed point theorem, we need to prove that the mappingS:w→uw is weakly continuous from W0 into W0. Let wj be a sequence that converges weakly to some w in W0 and let uj=uwj. We have to prove that S(wj)=uj converges weakly to S(w)=uw. From (45)–(47), and classical results of compact inclusion in Sobolev spaces [45], we can extract from wj,  respectively, from uj, a subsequence such that, for some u, we have
(53)dujdt⇀dudt,weaklyinL2(QT),uj→u,inL∞(0,T;L2(Ω)),∂uj∂xk⇀∂u∂xk,weakly*inL∞(0,T;L2(Ω)),wj⟶w,inL∞(0,T;L2(Ω)),∂Gσ∂xk*wj⟶∂Gσ∂xk*w,inL2(Ω),a.e.on(0,T)×Ω,C(|∇Gσ*wj|2)⟶C(|∇Gσ*w|2),inL2(0,T;L2(Ω)),uj(0)⟶u(0),inL2(Ω).The above convergence allows us to pass to the limit in the problem (Pwj) and obtain u=uw=S(w). Moreover, since the solution is unique, the whole sequence uj=S(wj) converges weakly in W0 to u=S(w); that is, S is weakly continuous. Consequently, thanks to Schauder's fixed-point theorem, there exists w∈W0 such that w=S(w)=uw. The function uw solves (32)–(34).
Now, we turn to the proof of the uniqueness, following the idea in [44]. Let u1 and u2 be two weak solutions of (32)–(34). For almost every t in [0,T] and i=1,2, we have
(54)ddt(u1-u2)(t)-div(α1(t)∇(u1-u2)(t))=div((α1-α2)(t)∇u2(t)),(u1-u2)(0,x)=0,inΩ,∂(u1-u2)∂n→=0,on∂Ω×(0,T),
in the distribution sense, where
(55)αi=C(|∇Gσ*ui|2).Then multiplying the above equality by (u1-u2), integrating over Ω, and using the Neumann boundary condition, we get a.e. t∈[0,T],
(56)12ddt∫Ω‍(u1(t)-u2(t))2dx+∫Ω‍α1|∇u1(t)-∇u2(t)|2dx=-∫Ω‍(α1-α2)∇u2(t)·(∇u1(t)-∇u2(t))dx.
Since C(s) is decreasing and positive, it follows that a.e. in (0,T)×Ω, αi≥ν, which implies from (56),
(57)12ddt∫Ω‍(u1(t)-u2(t))2dx+ν∫Ω‍|∇u1(t)-∇u2(t)|2dx≤∫Ω‍|α1-α2||∇u2(t)·(∇u1(t)-∇u2(t))|dx.
Moveover, since C(s), Gσ1, and Gσ are smooth, we have
(58)|α1(t)-α2(t)|L∞(Ω)≤C4|u1(t)-u2(t)|L2(Ω),
where C4 is a constant that depends only on g1, ν, and Gσ. From (58) and by using Young's inequality, we obtain
(59)12ddt∫Ω‍(u1(t)-u2(t))2dx+ν∫Ω‍|∇u1(t)-∇u2(t)|2dx≤12νC42∫Ω‍(u1(t)-u2(t))2dx∫Ω‍|∇u2(t)|2dx≤c+ν2∫‍|∇(u1-u2)(t)|2dx,
from which we deduce
(60)ddt∫Ω‍(u1(t)-u2(t))2dx≤1νC42∫Ω‍(u1(t)-u2(t))2dx∫Ω‍|∇u2(t)|2dx.
Since u1(0)=u2(0)=f, using Gronwall's inequality yields
(61)∫Ω‍(u1(t)-u2(t))2dx≤0;
that is, u1=u2.Remark 9.
Letu be the weak solution of problem (32)–(34) obtained in the proof of Theorem 8. Then from the proof we get that u∈L∞(ℝ+;H1(Ω)),∂u/∂t∈L2(Q∞), where Q∞=Ω×ℝ+.


## 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.
Assumeu  is the weak solutions of problem (32)–(34) with the initial data f. Then
(62)∫Ω‍(u-f)dx=0,∥u(·,t)-fΩ∥L2(Ω)≤e-νt/μ∥f-fΩ∥L2(Ω),a.e.t∈[0,+∞), where fΩ=(1|Ω|)∫Ω‍fdx, and |Ω| is Lebesgue measure of Ω.Proof.
Letu   be the solutions for problem (32)–(34) with the initial data f. For almost every t in [0,T], we have
(63)∂u∂t=div(C(|∇Gσ*u|2)∇u),inΩ×(0,T),(64)u(0,x)=f,inΩ,(65)∂u∂n→=0,on∂Ω×(0,T),
in the distribution sense. Integrating over the interval (0,t) and using the Neumann boundary condition yield
(66)∫Ω‍(u-f)dx=0.
Then, multiplying the above equality (63) by (u-fΩ), and integrating over Ω, and integrating by parts yield
(67)12ddt∫Ω‍(u-fΩ)2dx+∫Ω‍C(|∇Gσ*u|2)|∇u|2dx=0.
Using the following Poincaré-Wirtinger inequality [46, page 148], we have
(68)∥u-1|Ω|∫Ω‍udx∥L2(Ω)2=∥u-fΩ∥L2(Ω)2≤μ∫Ω‍|∇u|2dx,
with the constant μ≡μ(Ω). Substituting (68) to (67) yields
(69)ddt∫Ω‍(u-fΩ)2dx≤-2νμ∫Ω‍(u-fΩ)2dx.
Multiplying this inequality bye2νt/μ and integrating over the interval (0,t) we arrive to the inequality
(70)∫Ω‍(u-fΩ)2dx≤e-2νt/μ∫Ω‍(f-fΩ)2dx.
Hence, we obtain the assertion of the theorem.Next, let us build upon the maximum principle as follows.Theorem 11.
Letu be the weak solutions of problem (32)–(34) with the initial data f and f∈L∞(Ω). Then
(71)infx∈Ωf≤u≤supx∈Ωf.Proof.
LetI:=supx∈Ωf, and  J:=infx∈Ωf. Multiply (32) by (u-I)+, where
(72)(u-I)+={u-I,ifu-M>0,0,otherwise,
and integrate over Ω to get
(73)12ddt∫Ω‍(u(t)-I)+2dx+∫Ω‍C(|∇Gσ*ui|2)|∇(u(t)-I)+|2dx=0.
Then
(74)12ddt∫Ω‍(u(t)-I)+2dx≤0.
Therefore, (1/2)(d/dt)∫Ω‍(u(t)-I)+2dx is decreasing in t, and since
(75)∫Ω‍(u(t)-I)+2dx≥0,∫Ω‍(u(t)-I)+2dx|t=0=0,
we have that
(76)∫Ω‍(u(t)-I)+2dx=0,∀t∈[0,+∞),
and so
(77)u(t)≤supx∈Ωfa.e.onΩ,∀t>0.
Multiplying (32) by (u-J)-, a similar argument yields that u≥J for all t∈[0,+∞). Equation (71) is followed directly.


## 5. Behavior ast→∞

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.
Letu be the weak solutions of problem (32)–(34) with the initial data f∈L∞(Ω)∩H1(Ω), and limn→∞tn=+∞. Then(i)
for allτ≥0,  {u(tn+τ,·)}n=1∞→fΩ in L2(Ω),(ii)
there exists a subsequence of{tn}n=1∞, denoted also by itself, such that for all  T≥0, {u(·,tn+·)}n=1∞ converges to fΩ weakly in L2((0,T),H1(Ω)) and strongly in L2(QT), and {(∂u(·,tn+·))/∂t}n=1∞ converges to 0 weakly in L2(QT).Proof.
For allT≥0, since {u(tn+·,·)}n=1∞ are uniform bounded in L∞((0,T);H1(Ω))∩L∞(QT), and ∂u/∂t is bounded in L2(QT), and then there exist g(t,x) and the subsequence of {tn}n=1∞ which is independent on T, and denoted by {tnj}j=1∞, such that
(78)u(tnj+·,·)⇀gweaklyinL∞((0,T);H1(Ω))∩L∞(QT),u(tnj+·,·)⟶gstronglyinL2(QT),∂u(tnj+·,·)∂t⇀∂g∂tweaklyinL2(QT),
where g∈L∞((0,T);H1(Ω))∩L∞(QT). From Theorem 10, we have
(79)∫Ω‍(u(tnj+t,x)-fΩ)2dx≤e-2ν(tnj+t)∫Ω‍(f-fΩ)2dx.
Letting j→∞, we can obtain the first part of Lemma 12. For (79), integrate over (0,T) to get
(80)∫0T‍∫Ω‍(u(tnj+t,x)-fΩ)2dx≤2νe-2νtnj(1-e-2νT)∫Ω‍(f-fΩ)2dx.
Letting j→∞, we can obtain that
(81)u(tnj+·,·)⟶fΩstronglyinL2(QT).Noticing that  u(tnj+·,·)→gstronglyinL2(QT), we have
(82)g(t,x)=fΩ,a.e.(x,t)∈QT.
Since (∂u(tnj+·,·))/∂t⇀∂g/∂t  weakly in  L2(QT),  and therefore, we have
(83)∂u(tnj+·,·)∂t⇀0weaklyinL2(QT).
Hence, we obtain the remaining part of the lemma.Now let us consider the following problem:(84)div(C(|∇Gσ*u|2)∇u)=0,inΩ,(85)∂u∂n→=0,on∂Ω,(86)∫Ω‍(u-f)=0.Theorem 13.
Assumef∈L1(Ω). Then the problem (84)–(86) admits one and only one weak solution u∈H1(Ω) such that
(87)∫Ω‍C(|∇Gσ*u|2)∇u∇v=0,∀v∈C∞(Ω¯),∫Ω‍(u-f)=0.Proof.
It is clearly thatu=fΩ 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 u1 and u2 be two weak solutions of (84)–(86). Multiplying (84) by u, integrating over Ω, and using the Neumann boundary condition, we get
(88)∫Ω‍C(|∇Gσ*u|2)|∇u|2=0.
Using the following Poincaré-Wirtinger inequality, we have
(89)∥u-1|Ω|∫Ω‍udx∥L2(Ω)2=∥u-fΩ∥L2(Ω)2≤μ∫Ω‍|∇u|2dx,
with the constant μ≡μ(Ω). Substituting (86) and (89) to (88) yields
(90)∫Ω‍(u-fΩ)2dx=0.
Then
(91)∫Ω‍(u1-u2)2dx≤∫Ω‍(u1-fΩ)2dx+∫Ω‍(u2-fΩ)2dx=0.
That is, u1=u2.Theorem 14.
letu be the weak solution of problem (32)–(34), Then when t→∞, u tends to be steady-state solution fΩ, that is, the solution for Problem (84)–(86).Proof.
Let(92)un(·,·)=u(tn+·,·).
Then un is the weak solutions of problem (32)–(34) with the initial data u(tn,·). From Lemma 12, we obtain there exists a subsequence of {un}n=1∞, denoted also by itself, such that, for all T≥0,
(93)un⇀fΩweaklyinL∞((0,T);H1(Ω))∩L∞(QT),un⟶fΩstronglyinL2(QT),∂un∂t⇀0weaklyinL2(QT),
which implies that as t→∞, u tends to be steady-state solution fΩ, which is the unique solution for problem (84)–(86).Remark 15.
From Theorems10, 11, and 14, we can observe the following:  (i)
infx∈Ωf≤u≤supx∈Ωf, which means that no new features are introduced in the image in process.(ii)
uΩ, the mean of u, is constant fΩ.(iii)
∫Ω‍(u-fΩ)2dx tends to zero, which means that u converges in the L2(Ω)-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.
Letf∈H1(Ω). The sequence {un}n=1 defined by solving the iterative scheme
(94)∂un+1∂t=div(C(|∇Gσ*un|2)∇un+1),in(0,T)×Ω,un+1(0,x)=f,inΩ,∂un+1∂n→=0,on∂(0,T)×Ω
converges in C([0,T],L2(Ω)) to the strong solution of (32)–(34).Proof.
We denote byαn=C(|∇Gσ*un|2).  By Proposition 7, problem (94) has a unique solution un+1. It is clear that
(95)αn≥C(|∇Gσ*f|L∞(Ω)2)a.e.in(0,T)×Ω.
Now we verify that the sequence {u}n=1∞ converges in C([0,T],L2(Ω)) to u, the strong solution of (32)–(34).
As in Section3, from the estimate (60), we have
(96)ddt∫Ω‍(un+1(t)-u(t))2dx≤1νC42∫Ω‍(un(t)-u(t))2dx∫Ω‍|∇u(t)|2dx.
Moreover, we have
(97)∫Ω‍(f-u)2dx≤C0∀t∈[0,T],
where C0 is a constant which only depends on ∥f∥H1(Ω). Then Gronwall's inequality yields, for any t∈[0,T]:
(98)∫Ω‍(u1(t)-u(t))2dx≤C0∫0t‍a(s)ds,
where
(99)a(s)=C42∫Ω‍|∇u(s)|2dx.
By (96) and (98), we can deduce
(100)ddt(∥u2(t)-u(t)∥L2(Ω)2)≤C0a(t)∫0t‍a(s)ds,
and thus,
(101)∥u2(t)-u(t)∥L2(Ω)2≤C012(∫0t‍a(s)ds)2.
Finally, we obtain by iterating
(102)∥un+1(t)-u(t)∥L2(Ω)2≤C01(n+1)!(∫0T‍a(s)ds)n+1,
which implies that the sequence {u}n=1∞ converges in C([0,T],L2(Ω)) 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 h=1 and the time step by τ. Thus, we have
(103)∇N(ui,j)=ui-1,j-ui,j,∇S(ui,j)=ui+1,j-ui,j,∇E(ui,j)=ui,j+1-ui,j,∇W(ui,j)=ui,j-1-ui,j.
The numerical algorithms for problems (12)–(14) are given in the following:
(104)ui,jn+1=ui,jn+τ(C(|∇N(ui,jn)|2)·∇N(ui,jn)+C(|∇S(ui,jn)|2)·∇S(ui,jn)+C(|∇E(ui,jn)|2)·∇E(ui,jn)+C(|∇W(ui,jn)|2)·∇W(ui,jn)),
where 0≤τ≤1/4 for the numerical scheme to be stable.


### 7.2. The AOS Scheme

Using the scheme in [47], (12) can be discretized as
(105)un+1=1m∑l=1m‍[I-mτAl(uk)]-1un,
where A(un)=[aij(un)],
(106)aij(un):={Cin+Cjn2h2[j∈𝒩(i)]-∑n∈𝒩(i)‍Cin+CNn2h2(j=i),0(else),Cin∶=C[12∑p,q∈𝒩(i)‍(upn-uqn2h)],
where 𝒩(i) is the set of the two neighbors of pixel i (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 (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}.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(107)PSNR=10log10(2552MN∥uO-u∥22),MAE=∥uO-u∥1MN,
where uO and u 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 factork and the time step τ=0.25. In the AOS scheme, besides the same influencing factor k 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 t→∞, 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 Figures3–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 K is necessary.Lenna image (300×300). (a) Noisy image corrupted by Gaussian noise for σ=20. (b) Original image. (c) Our algorithm by AOS, k=1, τ=2 (4 steps). (d) Our algorithm by PMS, τ=0.25 (47 steps). (e) PM algorithm, K=5, τ=0.25 (55 steps). (f) TV algorithm, τ=0.1 (182 steps).
(a)
Noisy:σ=20, PSNR=22.09
(b)
Original
(c)
AOS:PSNR=29.72, 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.38Lenna image (300×300). (a) Noisy image corrupted by Gaussian noise for σ=35. (b) Original image. (c) Our algorithm by AOS, k=0.02, τ=2 (8 steps). (d) Our algorithm by PMS, τ=0.25 (96 steps). (e) PM algorithm, K=7, τ=0.25 (84 steps). (f) TV algorithm, τ=0.1 (400 steps).
(a)
Noisy:σ=35, PSNR =17.21
(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.23Lenna image (300×300). (a) Noisy image corrupted by Gaussian noise for σ=50. (b) Original image. (c) Our algorithm by AOS, k=0.02, τ=2 (12 steps). (d) Our algorithm by PMS, τ=0.25 (154 steps). (e) PM algorithm, K=9, τ=0.25 (108 steps). (f) TV algorithm, τ=0.1 (650 steps).
(a)
Noisy:σ=50, PSNR =14.12
(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.93Boat image (512×512). (a) Noisy image corrupted by Gaussian noise for σ=20. (b) Original image. (c) Our algorithm by AOS, k=0.2, τ=2 (4 steps). (d) Our algorithm by PMS, τ=0.25 (47 steps). (e) PM algorithm, K=5, τ=0.25 (55 steps). (f) TV algorithm, τ=0.1 (194 steps).
(a)
Noisy:σ=20, PSNR =22.11
(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.97Boat image (512×512). (a) Noisy image corrupted by Gaussian noise for σ=35. (b) Original image. (c) Our algorithm by AOS, k=0.2, τ=2 (9 steps). (d) Our algorithm by PMS, τ=0.25 (110 steps). (e) PM algorithm, K=7, τ=0.25 (90 steps). (f) TV algorithm, τ=0.1 (450 steps).
(a)
Noisy:σ=35, PSNR =17.23
(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.98Boat image (512×512). (a) Noisy image corrupted by Gaussian noise for σ=50. (b) Original image. (c) Our algorithm by AOS, k=0.2, τ=2 (15 steps). (d) Our algorithm by PMS, τ=0.25 (170 steps). (e) PM algorithm, K=9, τ=0.25 (115 steps). (f) TV algorithm, τ=0.1 (710 steps).
(a)
Noisy:σ=50, PSNR =14.14
(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.98Figures3, 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
PSNR, MAE, and CPU time (seconds) of all methods.
PSNR
MAE
CPU time(s)
The synthetic image (128×128)
σ
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 (300×300)
σ
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 (512×512)
σ
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.17Table1 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).


## 7.1. The PM Scheme

To discretize (12), the finite difference scheme in [2] is used. Denote the space step by h=1 and the time step by τ. Thus, we have
(103)∇N(ui,j)=ui-1,j-ui,j,∇S(ui,j)=ui+1,j-ui,j,∇E(ui,j)=ui,j+1-ui,j,∇W(ui,j)=ui,j-1-ui,j.
The numerical algorithms for problems (12)–(14) are given in the following:
(104)ui,jn+1=ui,jn+τ(C(|∇N(ui,jn)|2)·∇N(ui,jn)+C(|∇S(ui,jn)|2)·∇S(ui,jn)+C(|∇E(ui,jn)|2)·∇E(ui,jn)+C(|∇W(ui,jn)|2)·∇W(ui,jn)),
where 0≤τ≤1/4 for the numerical scheme to be stable.


## 7.2. The AOS Scheme

Using the scheme in [47], (12) can be discretized as
(105)un+1=1m∑l=1m‍[I-mτAl(uk)]-1un,
where A(un)=[aij(un)],
(106)aij(un):={Cin+Cjn2h2[j∈𝒩(i)]-∑n∈𝒩(i)‍Cin+CNn2h2(j=i),0(else),Cin∶=C[12∑p,q∈𝒩(i)‍(upn-uqn2h)],
where 𝒩(i) is the set of the two neighbors of pixel i (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 (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}.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(107)PSNR=10log10(2552MN∥uO-u∥22),MAE=∥uO-u∥1MN,
where uO and u 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 factork and the time step τ=0.25. In the AOS scheme, besides the same influencing factor k 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 t→∞, 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 Figures3–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 K is necessary.Lenna image (300×300). (a) Noisy image corrupted by Gaussian noise for σ=20. (b) Original image. (c) Our algorithm by AOS, k=1, τ=2 (4 steps). (d) Our algorithm by PMS, τ=0.25 (47 steps). (e) PM algorithm, K=5, τ=0.25 (55 steps). (f) TV algorithm, τ=0.1 (182 steps).
(a)
Noisy:σ=20, PSNR=22.09
(b)
Original
(c)
AOS:PSNR=29.72, 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.38Lenna image (300×300). (a) Noisy image corrupted by Gaussian noise for σ=35. (b) Original image. (c) Our algorithm by AOS, k=0.02, τ=2 (8 steps). (d) Our algorithm by PMS, τ=0.25 (96 steps). (e) PM algorithm, K=7, τ=0.25 (84 steps). (f) TV algorithm, τ=0.1 (400 steps).
(a)
Noisy:σ=35, PSNR =17.21
(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.23Lenna image (300×300). (a) Noisy image corrupted by Gaussian noise for σ=50. (b) Original image. (c) Our algorithm by AOS, k=0.02, τ=2 (12 steps). (d) Our algorithm by PMS, τ=0.25 (154 steps). (e) PM algorithm, K=9, τ=0.25 (108 steps). (f) TV algorithm, τ=0.1 (650 steps).
(a)
Noisy:σ=50, PSNR =14.12
(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.93Boat image (512×512). (a) Noisy image corrupted by Gaussian noise for σ=20. (b) Original image. (c) Our algorithm by AOS, k=0.2, τ=2 (4 steps). (d) Our algorithm by PMS, τ=0.25 (47 steps). (e) PM algorithm, K=5, τ=0.25 (55 steps). (f) TV algorithm, τ=0.1 (194 steps).
(a)
Noisy:σ=20, PSNR =22.11
(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.97Boat image (512×512). (a) Noisy image corrupted by Gaussian noise for σ=35. (b) Original image. (c) Our algorithm by AOS, k=0.2, τ=2 (9 steps). (d) Our algorithm by PMS, τ=0.25 (110 steps). (e) PM algorithm, K=7, τ=0.25 (90 steps). (f) TV algorithm, τ=0.1 (450 steps).
(a)
Noisy:σ=35, PSNR =17.23
(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.98Boat image (512×512). (a) Noisy image corrupted by Gaussian noise for σ=50. (b) Original image. (c) Our algorithm by AOS, k=0.2, τ=2 (15 steps). (d) Our algorithm by PMS, τ=0.25 (170 steps). (e) PM algorithm, K=9, τ=0.25 (115 steps). (f) TV algorithm, τ=0.1 (710 steps).
(a)
Noisy:σ=50, PSNR =14.14
(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.98Figures3, 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
PSNR, MAE, and CPU time (seconds) of all methods.
PSNR
MAE
CPU time(s)
The synthetic image (128×128)
σ
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 (300×300)
σ
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 (512×512)
σ
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.17Table1 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 functionC(s), 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.



---

*Source: 890157-2013-04-09.xml*