Image Denoising via Nonlinear Hybrid Diffusion

Abstract

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

1. Introduction

Image restoration and smoothing are important in problems ranging from medical diagnostic tests to defense applications such as target recognition. Over the past 20 years, the use of variational methods and nonlinear partial differential equations (PDEs) has significantly grown and evolved to address the image restoration problem. Let 𝑢₀ be the intensity of an image obtained from a noiseless image by adding Gaussian noise with zero mean, defined on a rectangle Ω⊂ℝ², and let 𝑢 represent the reconstructed image. The problem is to recover the restoration image 𝑢, from the observed, noisy image 𝑢₀, where the two are related by 𝑢₀=𝑢+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 𝑢_𝑡(𝑐),𝑢()=div(|∇𝑢|)∇𝑢0,𝑥=𝑓,(1) 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: 𝑐(|𝑥|)=1.(2)

Because the model is linear isotropic diffusion, it cannot preserve the edge and some features. The PM diffusivity function [2] is usually 𝑐_PM𝑥=1(||)1+|𝑥|²/𝐾².(3)

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]: 𝑐𝑥=1(||)1+|𝑥|²/𝐾².(4)  TV diffusivity [5]: 𝑐𝑥=1(||)|𝑥|.(5)  Weickert diffusivity [6]: 𝑐𝑥=|𝑥|(||)1,=01−exp−3.31488|𝑥|⁸/𝑘⁸,|𝑥|>0.(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 𝜕𝑢𝑐𝜕𝑡=div∇𝐺_𝜎(),𝑢()∗𝑢∇𝑢,inΩ×0,𝑇0,𝑥=𝑢₀,inΩ,𝜕𝑢𝜕𝑛=0,on𝜕Ω×(0,𝑇),(7) where 𝐺_𝜎 is the Gaussian kernel, as described in Section 2. Existence and uniqueness of solutions to this modified Perona-Malik equation have been proved for initial data 𝑢₀∈𝐿²(Ω). Another way is to use a higher-order version of the Perona-Malik equation, examples of which are given in [16–18].

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

1.2. The TV Framework

The famous total variation method first proposed by Rudin et al. [28] consists in solving the following constrained minimization problem: min_Ω|∇𝑢_Ω|,𝑢𝑑𝑥=_Ω𝑢₀,_Ω𝑢−𝑢₀²₂𝑑𝑥=𝜂².(8)

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_{𝑢∈BV(Ω)}𝐸(𝑢)=_Ω|∇𝑢|+𝜆𝑢−𝑢₀²_𝐿2(Ω).(9)

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 min_{𝑢∈BV(Ω)}𝐸_𝑔(𝑢)=_Ω𝑔(𝑥)|∇𝑢|+𝜆𝑢−𝑢₀²_𝐿2(Ω),(10) in which they introduce a control factor, 𝑔(𝑥), which slows the diffusion at likely edges. The factor 𝑔(𝑥) controls the speed of the diffusion and has demonstrated good results as it aids in noise reduction. It is also good at reconstructing edges, since the type of diffusion is the same as that of the original TV model.

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

Some authors consider another regularizing term to remove the noise [34], which is as follows: inf_𝑢𝐸(𝑢)=_Ω|∇𝑢|^{𝑝(|∇𝑢|)}𝑑𝑥+𝜆𝑓−𝑢²_𝐿2(Ω),(11) 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.

2. Nonlinear Hybrid Diffusion Model

2.1. The New Model

In this paper, we propose the following model: 𝜕𝑢𝜕𝑡=div∇𝑢|∇𝑢|²+1^{(2−𝑝(|∇𝑢|2))/2}()𝑢(),inΩ×0,𝑇,(12)0,𝑥=𝑓,inΩ,(13)𝜕𝑢𝜕𝑛=0,on𝜕Ω×(0,𝑇),(14) where 𝑝|∇𝑢|²1=1+1+𝑘|∇𝑢|²,𝑘>0,𝜎>0,(15)𝑓 is the original image, 𝑘>0,  𝜎 and 𝑇>0 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 𝑠→+∞, 𝑝(𝑠)→1, the new divergence operator is changed as follows: div∇𝑢|∇𝑢|²+1^{(2−𝑝(|∇𝑢|))/2}⟶div∇𝑢|∇𝑢|²+1,(16) 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∇𝑢|∇𝑢|²+1^{(2−𝑝(|∇𝑢|))/2}⟶Δ𝑢,(17) 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 𝐶(𝑠)=(1+𝑠)^{(𝑝(𝑠)−2)/2},(18) where 𝑝(𝑠)1=1+.1+𝑘𝑠(19)

Now, we discuss the properties of 𝐶(𝑠) as follows.

Proposition 1. One has the following:(a)𝐶∶[0,+∞)→(0,+∞)   is a strict decreasing function, and   0≤𝐶(𝑠)≤1,    for   𝑠∈ℝ₊     (see Figure 1(a)); (b)𝐶(0)=1,   and   lim_𝑠→∞𝐶(𝑠)=0     (see Figure 1(a)); (c)lim_𝑠→+∞𝑠𝐶^′(𝑠)𝐶(𝑠)=−1/2;(d)lim_𝑠→+∞√1+𝑠𝐶(𝑠)=1   and   lim_𝑠→+∞𝑠𝐶(𝑠²)=1; (e)lim_𝑠→+∞𝑠𝐶^′(𝑠²)=𝑘/2.

(a)

(b)

Figure 1 

(a) Diffusivity function 𝐶(𝑠), for 𝑘=0.01,0.1,10 and the PM diffusivity for 𝑘=5. (b) The Flux function ℱ(𝑠) for 𝑘=0.01,0.1,10 and PM flux function for 𝑘=5.

Proof. By a direct calculation, we have 𝐶(𝑠)−()=exp𝑘𝑠ln1+𝑠2(),𝐶1+𝑘𝑠_′(𝑠)(𝑠)𝑠()+()()=−𝑘𝐶𝑘𝑠+11+𝑠ln1+𝑠2(1+𝑠)(1+𝑘𝑠)²<0,(20) 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 𝐶(𝑠)≈1/√𝑠 or 𝐶(𝑠)≈1/√1+𝑠 as 𝑠→+∞. (3)The diffusion coefficient 𝐶(𝑠) is dependent on the exponent 𝑝(𝑠), which have the similar function with the fractional operator [19–27].

2.4. Forward-Backward Diffusion

For the diffusivity function 𝐶(𝑠) it follows the new flux function ℱ(𝑠) which is defined by ℱ(𝑠)∶=𝑠𝐶𝑠²,(21) where the variable 𝑠 stands for the norm of the gradient |∇𝑢|.

In the two-dimensional case, (32) can be replaced by [43] 𝜕_𝑡𝑢=div𝐶|∇𝑢|²||∇𝑢=𝐶∇𝑢²𝑢_𝑇𝑇+ℱ^′(∇𝑢)𝑢_𝑁𝑁,(22) where we have denoted by 𝑢_𝑁𝑁 and 𝑢_𝑇𝑇 the second derivatives of 𝑢 in the direction 𝑁(𝑥,𝑦)=∇𝑢/|∇𝑢| which is parallel to ∇𝑢 and 𝑇(𝑥,𝑦) in the orthogonal direction to 𝑁(𝑥,𝑦), respectively: 𝑢_𝑁𝑁=𝑢²_𝑥𝑢_𝑥𝑥+𝑢²_𝑦𝑢_𝑦𝑦+2𝑢_𝑥𝑢_𝑦𝑢_𝑥𝑦|∇𝑢|²,𝑢_𝑇𝑇=𝑢²_𝑥𝑢_𝑦𝑦+𝑢²_𝑦𝑢_𝑥𝑥−2𝑢_𝑥𝑢_𝑦𝑢_𝑥𝑦|∇𝑢|².(23)

Remark 3. From Proposition 1(d), we impose lim_𝑠→∞ℱ^′(𝑠)𝐶𝑠²=0,(24) which implies that it is preferable to smooth more in the 𝑇-direction than in the 𝑁-direction.

Proposition 4. There exists 𝑠₀∈(0,+∞) such that ℱ^′(𝑠)≥0 for |𝑠|≤𝑠₀, and ℱ^′(𝑠)<0 for |𝑠|>𝑠₀ (see Figure 1(b)). Moreover, 𝑠₀>𝑒^𝑘−1.(25)

Proof. By a direct calculation, we have ℱ^′(𝑠)=𝐶𝑠²+2𝑠²𝐶^′𝑠²=𝐶𝑠²ℎ𝑠²1+𝑠²1+𝑘𝑠²²,(26) where ℎ(𝑠)=𝑘²+𝑘𝑠²+(2𝑘+1)𝑠+1−𝑘𝑠(1+𝑠)ln(1+𝑠),(27) for 𝑠>0. Then ℎ^′(𝑠)=2𝑘²()()(),ℎ+𝑘𝑠+2𝑘+1−𝑘2𝑠+1ln1+𝑠_′′(𝑠)()()+1=𝑘2𝑘−1−2ln1+𝑠((1+𝑠<2𝑘𝑘−ln1+𝑠))<0,(28) for 𝑘<ln(1+𝑠); that is, 𝑠>𝑒^𝑘−1. It is noticed that, for 𝑠=𝑒^𝑘−1, ℎ(𝑠)=𝑘𝑠²+2𝑘+1−𝑘²𝑠+1>0,lim_𝑠→+∞ℎ(𝑠)=−∞.(29)

Because of the continuousness of 𝑔(𝑠), there exists a unique point 𝑠₀∈(0,+∞) such that ℎ(𝑠₀)=0, and ℎ(𝑠)≥0 for 𝑒^𝑘−1<𝑠≤𝑠₀, and ℎ(𝑠)<0 for 𝑠>𝑠₀. For 𝑘>ln(1+𝑠), that is, 0<𝑠<𝑒^𝑘−1, ℎ(𝑠)=𝑘²+𝑘𝑠²+(2𝑘+1)𝑠+1−𝑘𝑠(1+𝑠)ln(1+𝑠)>𝑘𝑠²+2𝑘+1−𝑘²𝑠+1>0.(30)

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: 𝑠₀>𝑒^𝑘−1.(31) 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.

Figure 2 

Flux function ℱ(𝑠) for 𝑘=0.01,0.02,0.1,1,10 and PM flux function for 𝑘=5.

(a) Noisy: 𝜎= 20, PSNR = 22.08

(b) Original

(c) AOS: PSNR = 39.00, MAE = 1.86

(d) PMS: PSNR = 40.70, MAE = 1.62

(e) PM: PSNR = 39.74, MAE = 1.78

(f) TV: PSNR = 37.08, MAE = 2.52

Figure 3 

Synthetic image (128×128) (a) Noisy image corrupted by Gaussian noise for 𝜎=20. (b) Original image. (c) Our algorithm by AOS, 𝑘=0.02, 𝜏=2 (11 steps). (d) Our algorithm by PMS, 𝜏=0.25 (62 steps). (e) PM algorithm, 𝐾=5, 𝜏=0.25 (90 steps). (f) TV algorithm, 𝜏=0.1 (360 steps).

(a) Noisy: 𝜎= 35, PSNR = 17.19

(b) Original

(c) AOS: PSNR = 35.43, MAE = 2.75

(d) PMS: PSNR = 36.61, MAE = 2.43

(e) PM: PSNR = 34.81, MAE = 2.76

(f) TV: PSNR = 34.05, MAE = 3.65

Figure 4 

Synthetic image (128×128). (a) Noisy image corrupted by Gaussian noise for 𝜎=35. (b) Original image. (c) Our algorithm by AOS, 𝑘=0.02, 𝜏=2 (11 steps). (d) Our algorithm by PMS, 𝜏=0.25 (115 steps). (e) PM algorithm, 𝐾=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

Figure 5 

Synthetic image (128×128). (a) Noisy image corrupted by Gaussian noise for 𝜎=50. (b) Original image. (c) Our algorithm by AOS, 𝑘=0.02, 𝜏=2 (11 steps). (d) Our algorithm by PMS, 𝜏=0.25 (115 steps). (e) PM algorithm, 𝐾=7, 𝜏=0.25 (130 steps). (f) TV algorithm, 𝜏=0.1 (600 steps).

2.5. The Modified Regularization Equation

Using the similar skill in [15], the new model can be regularized by 𝜕𝑢𝐶𝜕𝑡=div∇𝐺_𝜎∗𝑢²(),𝑢()∇𝑢,inΩ×0,𝑇(32)0,𝑥=𝑓,inΩ,(33)𝜕𝑢𝜕𝑛=0,on𝜕Ω×(0,𝑇),(34) where 𝐺_𝜎(𝑥) is the Gaussian kernel, namely, 𝐺_𝜎(𝑥)=1(4𝜋𝜎)^𝑁/2−|𝑥|exp².4𝜎(35)

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 |𝑚|=∑^𝑁_𝑗=1𝑚_𝑗≤𝑘 all belong to 𝐿²(Ω). 𝐻^𝑘(Ω) is a Hilbert space with the norm ‖𝑢‖_𝐻^𝑘_(Ω)=_{|𝑚|≤𝑘Ω}𝜕^𝑚𝑢𝜕𝑥^𝑚2𝑑𝑥^1/2.(36) The space 𝐿^∞(0,𝑇;𝑋) consists of all functions 𝑢 such that, for almost every 𝑡 in (0,𝑇), 𝑢 belongs to 𝑋. 𝐿^∞(0,𝑇;𝑋) is a normed space with the norm ‖𝑢‖_𝐿^∞_{(0,𝑇;𝑋)}=esssup_{0≤𝑡≤𝑇}‖𝑢(⋅,𝑡)‖_𝑋.(37) We denote by 𝐻¹(Ω)^′ the dual of 𝐻¹(Ω).

We introduce the solution space 𝑊 of the problem (32)–(34) as follows: 𝑊(0,𝑇)=𝑤∈𝐿^∞0,𝑇;𝐻¹(Ω);𝑑𝑤𝑑𝑡∈𝐿²𝑄_𝑇.(38) Obviously, 𝑊 is a Banach space with the norm ‖𝑤‖_𝑊=‖𝑤‖_𝐿^∞_(0,𝑇;𝐻¹_(Ω))+𝜕𝑤𝜕𝑡_𝐿²_(𝑄𝑇).(39) 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 𝜕𝑢𝐶𝜕𝑡=div∇𝐺_𝜎∗𝑤²(),𝑢()∇𝑢,inΩ×0,𝑇(40)0,𝑥=𝑓,inΩ,(41)𝜕𝑢𝜕𝑛=0,on𝜕Ω×(0,𝑇).(42)

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 𝜕𝑢/𝜕𝑡∈𝐿²(𝑄_𝑇), 𝑢∈𝐿^∞(0,𝑇;𝐻¹(Ω)∩𝑢∈𝐶([0,𝑇],𝐿²(Ω)).

Proof. Firstly, we consider the proof of the existence, which is based on the Schauder fixed point argument. Let 𝑤∈𝑊 such that ‖𝑤‖_𝐿^∞_(0,𝑇;𝐿²_(Ω))≤𝑓_𝐿2(Ω),𝜕𝑤𝜕𝑡_𝐿²_(𝑄𝑇)≤𝑓_𝐻1(Ω).(43) We consider the following linear problem 𝑃_𝑤: 𝑑𝑢(𝑡)𝑑𝑡,𝑣_{𝐻1(Ω)′×𝐻1(Ω)}+_Ω𝐶∇𝐺_𝜎∗𝑤²(𝑡)(𝑡)𝑃∇𝑢∇𝑣𝑑𝑥=0^𝑤 for all 𝑣∈𝐻¹(Ω), a.e. 𝑡∈[0,𝑇]. Since 𝑤 and 𝜕𝑤/𝜕𝑡 satisfy (43), then |∇𝐺_𝜎∗𝑤| and |∇𝐺_𝜎∗(𝜕𝑤/𝜕𝑡)| belong to 𝐿^∞((0,𝑇);𝐶^∞(Ω)) and there exists a constant 𝑀=𝑀(𝐺_𝜎,‖𝑓‖_𝐻¹_(Ω)) such that |∇𝐺_𝜎∗𝑤|≤𝑀 and |∇𝐺_𝜎∗(𝜕𝑤/𝜕𝑡)|≤𝑀 a.e. 𝑡,forall𝑥∈Ω. Since 𝐶(𝑠) is decreasing and positive, it follows that a.e. in (0,𝑇)×Ω: 1≥𝐶∇𝐺_𝜎∗𝑤²≥𝐶𝑀²=𝜈.(44) Then by applying Proposition 7 on the linearized problem, we will prove that the problem 𝑃_𝑤 has a unique solution 𝑢_𝑤∈𝑊 satisfying the estimates 𝑢_{𝑤𝐿∞(0,𝑇;𝐻1(Ω))}≤𝑐₁,𝑢(45)_{𝑤𝐿∞(0,𝑇;𝐿2(Ω))}≤𝑓_𝐿2(Ω),(46)𝜕𝑢_𝑤𝜕𝑡_𝐿²_(𝑄𝑇)≤𝑓_𝐻1(Ω),(47) where 𝑐₁ is the constant depending only on the constant 𝜈, 𝐺_𝜎, and ‖𝑓‖_𝐻¹_(Ω). Choosing 𝑣=𝑢_𝑤 in 𝑃_𝑤, integrating over the interval (0,𝑡), we arrive to the inequality 12_Ω𝑢²_𝑤𝑑𝑥+𝜈^𝑡_0Ω∇𝑢_𝑤²1𝑑𝑥𝑑𝑠≤2_Ω𝑓²𝑑𝑥,(48) which implies (46). Choosing 𝑣=𝜕𝑢_𝑤/𝜕𝑡 in 𝑃_𝑤, integrating by parts yields _Ω𝜕𝑢_𝑤𝜕𝑡²1𝑑𝑥+2_Ω𝐶∇𝐺_𝜎∗𝑤²𝜕∇𝑢_𝑤²𝜕𝑡𝑑𝑥=0.(49) Integrating over the interval (0,𝑡) we arrive to that ^𝑡_0Ω𝜕𝑢_𝑤𝜕𝑠²1𝑑𝑥𝑑𝑠+2_Ω𝐶∇𝐺_𝜎∗𝑤²∇𝑢_𝑤²=1𝑑𝑥2_Ω𝐶∇𝐺_𝜎∗𝑤²∇𝑓²+1𝑑𝑥2^𝑡_0Ω𝐶^′∇𝐺_𝜎∗𝑤²×∇𝐺_𝜎∗𝑤∇𝐺_𝜎∗𝜕𝑤𝜕𝑠∇𝑢_𝑤²𝑑𝑥𝑑𝑠.(50) From Proposition 1 and (48), noticing that |𝐶^′(𝑠²)𝑠|≤𝑘, we can deduce that ^𝑡_0Ω𝜕𝑢_𝑤𝜕𝑠²𝜈𝑑𝑥𝑑𝑠+2_Ω|∇𝑢|²≤𝑑𝑥_Ω∇𝑓²𝑑𝑥+𝑀𝑘4𝜈_Ω𝑓²𝑑𝑥𝑑𝑠.(51) Since ‖𝑓‖_𝐻¹_(Ω) is small, letting 𝑀𝑘/(8𝜈)≤1 yields (45) and (47). From (45)–(47), we introduce the subspace 𝑊₀ of 𝑊 defined by 𝑊₀=()(0)‖𝑤‖𝑤∈𝑊0,𝑇,𝑤=𝑓_{𝐿∞(0,𝑇;𝐻1(Ω))}≤𝑐₁,‖𝑤‖_𝐿^∞_(0,𝑇;𝐿²_(Ω))≤𝑓_𝐿2(Ω),𝑑𝑤𝑑𝑡_𝐿²_(𝑄𝑇)≤𝑓_𝐻1(Ω).(52) By construction, 𝑤→𝑆(𝑤)≡𝑢_𝑤 is a mapping from 𝑊₀ into 𝑊₀. Moreover, one can prove that 𝑊₀ 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 𝑊₀ 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 𝑑𝑢_𝑗⇀𝑑𝑡𝑑𝑢𝑑𝑡,weaklyin𝐿²𝑄_𝑇,𝑢_𝑗→𝑢,in𝐿^∞0,𝑇;𝐿²(Ω),𝜕𝑢_𝑗𝜕𝑥_𝑘⇀𝜕𝑢𝜕𝑥_𝑘,weakly∗in𝐿^∞0,𝑇;𝐿²(Ω),𝑤_𝑗⟶𝑤,in𝐿^∞0,𝑇;𝐿²(Ω),𝜕𝐺_𝜎𝜕𝑥_𝑘∗𝑤_𝑗⟶𝜕𝐺_𝜎𝜕𝑥_𝑘∗𝑤,in𝐿²(Ω)()𝐶,a.e.on0,𝑇×Ω,∇𝐺_𝜎∗𝑤_𝑗²⟶𝐶∇𝐺_𝜎∗𝑤²,in𝐿²0,𝑇;𝐿²(Ω),𝑢_𝑗(0)⟶𝑢(0),in𝐿²(Ω).(53)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 [0,𝑇] and 𝑖=1,2, we have 𝑑𝑢𝑑𝑡¹−𝑢₂(𝑡)−div𝛼₁(𝑡)∇𝑢₁−𝑢₂(𝑡)=div𝛼₁−𝛼₂(𝑡)∇𝑢₂(𝑡),𝑢₁−𝑢₂()𝜕𝑢0,𝑥=0,inΩ,₁−𝑢₂𝜕𝑛=0,on𝜕Ω×(0,𝑇),(54) in the distribution sense, where 𝛼_𝑖=𝐶∇𝐺_𝜎∗𝑢_𝑖².(55)Then multiplying the above equality by (𝑢₁−𝑢₂), integrating over Ω, and using the Neumann boundary condition, we get a.e. 𝑡∈[0,𝑇], 12𝑑𝑑𝑡_Ω𝑢₁(𝑡)−𝑢₂(𝑡)²𝑑𝑥+_Ω𝛼₁∇𝑢₁(𝑡)−∇𝑢₂(𝑡)²𝑑𝑥=−_Ω𝛼₁−𝛼₂∇𝑢₂(𝑡)⋅∇𝑢₁(𝑡)−∇𝑢₂(𝑡)𝑑𝑥.(56) Since 𝐶(𝑠) is decreasing and positive, it follows that a.e. in (0,𝑇)×Ω, 𝛼_𝑖≥𝜈, which implies from (56), 12𝑑𝑑𝑡_Ω𝑢₁(𝑡)−𝑢₂(𝑡)²𝑑𝑥+𝜈_Ω∇𝑢₁(𝑡)−∇𝑢₂(𝑡)²≤𝑑𝑥_Ω𝛼₁−𝛼₂∇𝑢₂(𝑡)⋅∇𝑢₁(𝑡)−∇𝑢₂(𝑡)𝑑𝑥.(57) Moveover, since 𝐶(𝑠), 𝐺_𝜎1, and 𝐺_𝜎 are smooth, we have 𝛼₁(𝑡)−𝛼₂(𝑡)_𝐿∞(Ω)≤𝐶₄𝑢₁(𝑡)−𝑢₂(𝑡)_𝐿2(Ω),(58) where 𝐶₄ is a constant that depends only on 𝑔₁, 𝜈, and 𝐺_𝜎. From (58) and by using Young's inequality, we obtain 12𝑑𝑑𝑡_Ω𝑢₁(𝑡)−𝑢₂(𝑡)²𝑑𝑥+𝜈_Ω∇𝑢₁(𝑡)−∇𝑢₂(𝑡)²≤1𝑑𝑥𝐶2𝜈²⁴_Ω𝑢₁(𝑡)−𝑢₂(𝑡)²𝑑𝑥_Ω∇𝑢₂(𝑡)²+𝜈𝑑𝑥2∇𝑢₁−𝑢₂(𝑡)²𝑑𝑥,(59) from which we deduce 𝑑𝑑𝑡_Ω𝑢₁(𝑡)−𝑢₂(𝑡)²≤1𝑑𝑥𝜈𝐶²_4Ω𝑢₁(𝑡)−𝑢₂(𝑡)²𝑑𝑥_Ω∇𝑢₂(𝑡)²𝑑𝑥.(60) Since 𝑢₁(0)=𝑢₂(0)=𝑓, using Gronwall's inequality yields _Ω𝑢₁(𝑡)−𝑢₂(𝑡)²𝑑𝑥≤0;(61) 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 𝑢()^Ω𝑢−𝑓𝑑𝑥=0,⋅,𝑡−𝑓_Ω𝐿2(Ω)≤𝑒^{−𝜈𝑡/𝜇}𝑓−𝑓_Ω𝐿2(Ω),(62)𝑎.𝑒.𝑡∈[0,+∞), where 𝑓_Ω=1∫|Ω|^Ω𝑓𝑑𝑥, and |Ω| is Lebesgue measure of Ω.

Proof. Let 𝑢 be the solutions for problem (32)–(34) with the initial data 𝑓. For almost every 𝑡 in [0,𝑇], we have 𝜕𝑢𝐶𝜕𝑡=div∇𝐺_𝜎∗𝑢²(),𝑢()∇𝑢,inΩ×0,𝑇(63)0,𝑥=𝑓,inΩ,(64)𝜕𝑢𝜕𝑛=0,on𝜕Ω×(0,𝑇),(65) in the distribution sense. Integrating over the interval (0,𝑡) and using the Neumann boundary condition yield _Ω𝑢−𝑓𝑑𝑥=0.(66) Then, multiplying the above equality (63) by (𝑢−𝑓_Ω), and integrating over Ω, and integrating by parts yield 12𝑑𝑑𝑡_Ω𝑢−𝑓_Ω²𝑑𝑥+_Ω𝐶∇𝐺_𝜎∗𝑢²|∇𝑢|²𝑑𝑥=0.(67) Using the following Poincaré-Wirtinger inequality [46, page 148], we have 1𝑢−|Ω|_Ω𝑢𝑑𝑥²_𝐿2(Ω)=𝑢−𝑓_Ω²_𝐿2(Ω)≤𝜇_Ω|∇𝑢|²𝑑𝑥,(68) with the constant 𝜇≡𝜇(Ω). Substituting (68) to (67) yields 𝑑𝑑𝑡_Ω𝑢−𝑓_Ω²𝑑𝑥≤−2𝜈𝜇_Ω𝑢−𝑓_Ω²𝑑𝑥.(69)
Multiplying this inequality by 𝑒^{2𝜈𝑡/𝜇} and integrating over the interval (0,𝑡) we arrive to the inequality _Ω𝑢−𝑓_Ω²𝑑𝑥≤𝑒^{−2𝜈𝑡/𝜇}_Ω𝑓−𝑓_Ω²𝑑𝑥.(70)
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 inf_𝑥∈Ω𝑓≤𝑢≤sup_𝑥∈Ω𝑓.(71)

Proof. Let 𝐼∶=sup_𝑥∈Ω𝑓, and 𝐽∶=inf_𝑥∈Ω𝑓. Multiply (32) by (𝑢−𝐼)₊, where (𝑢−𝐼)₊=𝑢−𝐼,if𝑢−𝑀>0,0,otherwise,(72) and integrate over Ω to get 12𝑑𝑑𝑡_Ω(𝑢(𝑡)−𝐼)²₊+𝑑𝑥_Ω𝐶∇𝐺_𝜎∗𝑢_𝑖²∇(𝑢(𝑡)−𝐼)₊²𝑑𝑥=0.(73) Then 12𝑑𝑑𝑡_Ω(𝑢(𝑡)−𝐼)²₊𝑑𝑥≤0.(74) Therefore, (1/2)(𝑑/𝑑𝑡)∫_Ω(𝑢(𝑡)−𝐼)²₊𝑑𝑥 is decreasing in 𝑡, and since _Ω(𝑢(𝑡)−𝐼)²₊𝑑𝑥≥0,_Ω(𝑢(𝑡)−𝐼)²₊𝑑𝑥_𝑡=0=0,(75) we have that _Ω(𝑢(𝑡)−𝐼)²₊𝑑𝑥=0,∀𝑡∈[0,+∞),(76) and so 𝑢(𝑡)≤sup_𝑥∈Ω𝑓a.e.onΩ,∀𝑡>0.(77) Multiplying (32) by (𝑢−𝐽)₋, a similar argument yields that 𝑢≥𝐽 for all 𝑡∈[0,+∞). 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 lim_𝑛→∞𝑡_𝑛=+∞. Then(i)for all 𝜏≥0, {𝑢(𝑡_𝑛+𝜏,⋅)}^∞_𝑛=1→𝑓_Ω in 𝐿²(Ω),(ii)there exists a subsequence of {𝑡_𝑛}^∞_𝑛=1, denoted also by itself, such that for all 𝑇≥0, {𝑢(⋅,𝑡_𝑛+⋅)}^∞_𝑛=1 converges to 𝑓_Ω weakly in 𝐿²((0,𝑇),𝐻¹(Ω)) and strongly in 𝐿²(𝑄_𝑇), and {(𝜕𝑢(⋅,𝑡_𝑛+⋅))/𝜕𝑡}^∞_𝑛=1 converges to 0 weakly in 𝐿²(𝑄_𝑇).

Proof. For all 𝑇≥0, since {𝑢(𝑡_𝑛+⋅,⋅)}^∞_𝑛=1 are uniform bounded in 𝐿^∞((0,𝑇);𝐻¹(Ω))∩𝐿^∞(𝑄_𝑇), and 𝜕𝑢/𝜕𝑡 is bounded in 𝐿²(𝑄_𝑇), and then there exist 𝑔(𝑡,𝑥) and the subsequence of {𝑡_𝑛}^∞_𝑛=1 which is independent on 𝑇, and denoted by {𝑡_𝑛𝑗}^∞_𝑗=1, such that 𝑢𝑡_𝑛𝑗𝐿+⋅,⋅⇀𝑔weaklyin_∞(0,𝑇);𝐻¹(Ω)∩𝐿^∞𝑄_𝑇,𝑢𝑡_𝑛𝑗+⋅,⋅⟶𝑔stronglyin𝐿²𝑄_𝑇,𝜕𝑢𝑡_𝑛𝑗+⋅,⋅⇀𝜕𝑡𝜕𝑔𝜕𝑡weaklyin𝐿²𝑄_𝑇,(78) where 𝑔∈𝐿^∞((0,𝑇);𝐻¹(Ω))∩𝐿^∞(𝑄_𝑇). From Theorem 10, we have _Ω𝑢𝑡_𝑛𝑗+𝑡,𝑥−𝑓_Ω²𝑑𝑥≤𝑒^{−2𝜈(𝑡𝑛𝑗+𝑡)}_Ω𝑓−𝑓_Ω²𝑑𝑥.(79) Letting 𝑗→∞, we can obtain the first part of Lemma 12. For (79), integrate over (0,𝑇) to get ^𝑇_0Ω𝑢𝑡_𝑛𝑗+𝑡,𝑥−𝑓_Ω²𝑑𝑥≤2𝜈𝑒^{−2𝜈𝑡𝑛𝑗}1−𝑒^−2𝜈𝑇_Ω𝑓−𝑓_Ω²𝑑𝑥.(80) Letting 𝑗→∞, we can obtain that 𝑢𝑡_𝑛𝑗+⋅,⋅⟶𝑓_Ωstronglyin𝐿²𝑄_𝑇.(81)Noticing that 𝑢(𝑡_𝑛𝑗+⋅,⋅)→𝑔stronglyin𝐿²(𝑄_𝑇), we have 𝑔(𝑡,𝑥)=𝑓_Ω,a.e.(𝑥,𝑡)∈𝑄_𝑇.(82) Since (𝜕𝑢(𝑡_𝑛𝑗+⋅,⋅))/𝜕𝑡⇀𝜕𝑔/𝜕𝑡 weakly in 𝐿²(𝑄_𝑇), and therefore, we have 𝜕𝑢𝑡_𝑛𝑗+⋅,⋅𝜕𝑡⇀0weaklyin𝐿²𝑄_𝑇.(83) Hence, we obtain the remaining part of the lemma.

Now let us consider the following problem: div𝐶∇𝐺_𝜎∗𝑢²∇𝑢=0,inΩ,(84)𝜕𝑢𝜕𝑛=0,on𝜕Ω,(85)_Ω𝑢−𝑓=0.(86)

Theorem 13. Assume 𝑓∈𝐿¹(Ω). Then the problem (84)–(86) admits one and only one weak solution 𝑢∈𝐻¹(Ω) such that _Ω𝐶∇𝐺_𝜎∗𝑢²∇𝑢∇𝑣=0,∀𝑣∈𝐶^∞Ω,_Ω𝑢−𝑓=0.(87)

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 _Ω𝐶∇𝐺_𝜎∗𝑢²|∇𝑢|²=0.(88) Using the following Poincaré-Wirtinger inequality, we have 1𝑢−|Ω|_Ω𝑢𝑑𝑥²_𝐿2(Ω)=𝑢−𝑓_Ω²_𝐿2(Ω)≤𝜇_Ω|∇𝑢|²𝑑𝑥,(89) with the constant 𝜇≡𝜇(Ω). Substituting (86) and (89) to (88) yields _Ω𝑢−𝑓_Ω²𝑑𝑥=0.(90) Then _Ω𝑢₁−𝑢₂²𝑑𝑥≤_Ω𝑢₁−𝑓_Ω²𝑑𝑥+_Ω𝑢₂−𝑓_Ω²𝑑𝑥=0.(91) 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 𝑢^𝑛(⋅,⋅)=𝑢𝑡_𝑛+⋅,⋅.(92) 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, 𝑢^𝑛⇀𝑓_Ωweaklyin𝐿^∞(0,𝑇);𝐻¹(Ω)∩𝐿^∞𝑄_𝑇,𝑢^𝑛⟶𝑓_Ωstronglyin𝐿²𝑄_𝑇,𝜕𝑢^𝑛𝜕𝑡⇀0weaklyin𝐿²𝑄_𝑇,(93) 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)∫_Ω(𝑢−𝑓_Ω)²𝑑𝑥 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 {𝑢^𝑛}_𝑛=1 defined by solving the iterative scheme 𝜕𝑢^𝑛+1𝐶𝜕𝑡=div∇𝐺_𝜎∗𝑢^𝑛2∇𝑢^𝑛+1()𝑢,𝑖𝑛0,𝑇×Ω,_𝑛+1(0,𝑥)=𝑓,𝑖𝑛Ω,𝜕𝑢^𝑛+1𝜕𝑛=0,𝑜𝑛𝜕(0,𝑇)×Ω(94) converges in 𝐶([0,𝑇],𝐿²(Ω)) to the strong solution of (32)–(34).

Proof. We denote by 𝛼^𝑛=𝐶(|∇𝐺_𝜎∗𝑢^𝑛|²). By Proposition 7, problem (94) has a unique solution 𝑢^𝑛+1. It is clear that 𝛼^𝑛≥𝐶∇𝐺_𝜎∗𝑓²_𝐿∞(Ω)a.e.in(0,𝑇)×Ω.(95) Now we verify that the sequence {𝑢}^∞_𝑛=1 converges in 𝐶([0,𝑇],𝐿²(Ω)) to 𝑢, the strong solution of (32)–(34).
As in Section 3, from the estimate (60), we have 𝑑𝑑𝑡_Ω𝑢^𝑛+1(𝑡)−𝑢(𝑡)²≤1𝑑𝑥𝜈𝐶²_4Ω𝑢^𝑛(𝑡)−𝑢(𝑡)²𝑑𝑥_Ω|∇𝑢(𝑡)|²𝑑𝑥.(96) Moreover, we have _Ω𝑓−𝑢²𝑑𝑥≤𝐶₀∀𝑡∈[0,𝑇],(97) where 𝐶₀ is a constant which only depends on ‖𝑓‖_𝐻¹_(Ω). Then Gronwall's inequality yields, for any 𝑡∈[0,𝑇]: _Ω𝑢¹(𝑡)−𝑢(𝑡)²𝑑𝑥≤𝐶₀^𝑡₀𝑎(𝑠)𝑑𝑠,(98) where 𝑎(𝑠)=𝐶²_4Ω|∇𝑢(𝑠)|²𝑑𝑥.(99)
By (96) and (98), we can deduce 𝑑𝑢𝑑𝑡²(𝑡)−𝑢(𝑡)²_𝐿2(Ω)≤𝐶₀𝑎(𝑡)^𝑡₀𝑎(𝑠)𝑑𝑠,(100) and thus, 𝑢²(𝑡)−𝑢(𝑡)²_𝐿2(Ω)≤𝐶₀12^𝑡₀𝑎(𝑠)𝑑𝑠².(101) Finally, we obtain by iterating 𝑢^𝑛+1(𝑡)−𝑢(𝑡)²_𝐿2(Ω)≤𝐶₀1(𝑛+1)!^𝑇₀𝑎(𝑠)𝑑𝑠^𝑛+1,(102) which implies that the sequence {𝑢}^∞_𝑛=1 converges in 𝐶([0,𝑇],𝐿²(Ω)) 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 ℎ=1 and the time step by 𝜏. Thus, we have ∇_𝑁𝑢_𝑖,𝑗=𝑢_{𝑖−1,𝑗}−𝑢_𝑖,𝑗,∇_𝑆𝑢_𝑖,𝑗=𝑢_𝑖+1,𝑗−𝑢_𝑖,𝑗,∇_𝐸𝑢_𝑖,𝑗=𝑢_𝑖,𝑗+1−𝑢_𝑖,𝑗,∇_𝑊𝑢_𝑖,𝑗=𝑢_{𝑖,𝑗−1}−𝑢_𝑖,𝑗.(103) The numerical algorithms for problems (12)–(14) are given in the following: 𝑢^𝑛+1_𝑖,𝑗=𝑢^𝑛_𝑖,𝑗+𝜏𝐶∇_𝑁𝑢^𝑛_𝑖,𝑗²⋅∇_𝑁𝑢^𝑛_𝑖,𝑗+𝐶∇_𝑆𝑢^𝑛_𝑖,𝑗²⋅∇_𝑆𝑢^𝑛_𝑖,𝑗+𝐶∇_𝐸𝑢^𝑛_𝑖,𝑗²⋅∇_𝐸𝑢^𝑛_𝑖,𝑗+𝐶∇_𝑊𝑢^𝑛_𝑖,𝑗²⋅∇_𝑊𝑢^𝑛_𝑖,𝑗,(104) 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 𝑢^𝑛+1=1𝑚^𝑚_𝑙=1𝐼−𝑚𝜏𝐴_𝑙𝑢^𝑘−1𝑢^𝑛,(105) where 𝐴(𝑢^𝑛)=[𝑎_𝑖𝑗(𝑢^𝑛)], 𝑎_𝑖𝑗𝑢^𝑛𝐶∶=^𝑛𝑖+𝐶^𝑛_𝑗2ℎ²𝑗∈𝒩(𝑖)−∑_{(𝑖)𝑛∈𝒩}𝐶^𝑛_𝑖+𝐶^𝑛_𝑁2ℎ²,0(),𝐶𝑗=𝑖else_𝑛𝑖1∶=𝐶2_{(𝑖)𝑝,𝑞∈𝒩}𝑢^𝑛_𝑝−𝑢^𝑛_𝑞,2ℎ(106) 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 (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 PSNR=10log₁₀255²MN𝑢_𝑂−𝑢²₂,𝑢MAE=^𝑂−𝑢₁,MN(107) 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.

(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.38

Figure 6 

Lenna image (300×300). (a) Noisy image corrupted by Gaussian noise for 𝜎=20. (b) Original image. (c) Our algorithm by AOS, 𝑘=1, 𝜏=2 (4 steps). (d) Our algorithm by PMS, 𝜏=0.25 (47 steps). (e) PM algorithm, 𝐾=5, 𝜏=0.25 (55 steps). (f) TV algorithm, 𝜏=0.1 (182 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.23

Figure 7 

Lenna image (300×300). (a) Noisy image corrupted by Gaussian noise for 𝜎=35. (b) Original image. (c) Our algorithm by AOS, 𝑘=0.02, 𝜏=2 (8 steps). (d) Our algorithm by PMS, 𝜏=0.25 (96 steps). (e) PM algorithm, 𝐾=7, 𝜏=0.25 (84 steps). (f) TV algorithm, 𝜏=0.1 (400 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.93

Figure 8 

Lenna image (300×300). (a) Noisy image corrupted by Gaussian noise for 𝜎=50. (b) Original image. (c) Our algorithm by AOS, 𝑘=0.02, 𝜏=2 (12 steps). (d) Our algorithm by PMS, 𝜏=0.25 (154 steps). (e) PM algorithm, 𝐾=9, 𝜏=0.25 (108 steps). (f) TV algorithm, 𝜏=0.1 (650 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.97

Figure 9 

Boat image (512×512). (a) Noisy image corrupted by Gaussian noise for 𝜎=20. (b) Original image. (c) Our algorithm by AOS, 𝑘=0.2, 𝜏=2 (4 steps). (d) Our algorithm by PMS, 𝜏=0.25 (47 steps). (e) PM algorithm, 𝐾=5, 𝜏=0.25 (55 steps). (f) TV algorithm, 𝜏=0.1 (194 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.98

Figure 10 

Boat image (512×512). (a) Noisy image corrupted by Gaussian noise for 𝜎=35. (b) Original image. (c) Our algorithm by AOS, 𝑘=0.2, 𝜏=2 (9 steps). (d) Our algorithm by PMS, 𝜏=0.25 (110 steps). (e) PM algorithm, 𝐾=7, 𝜏=0.25 (90 steps). (f) TV algorithm, 𝜏=0.1 (450 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.98

Figure 11 

Boat image (512×512). (a) Noisy image corrupted by Gaussian noise for 𝜎=50. (b) Original image. (c) Our algorithm by AOS, 𝑘=0.2, 𝜏=2 (15 steps). (d) Our algorithm by PMS, 𝜏=0.25 (170 steps). (e) PM algorithm, 𝐾=9, 𝜏=0.25 (115 steps). (f) TV algorithm, 𝜏=0.1 (710 steps).

Figures 3, 4, and 5 illustrate the proposed model is able to reconstruct sharp edges and nonuniform regions while avoiding staircasing. TV-based diffusion reconstructs sharp edges, but the staircasing effect is clear evidence. PM-based diffusion also reconstructs sharp edges but creates isolated black and white speckles in the denoise image. The proposed model reconstructs sharp edges as effectively as PM-based diffusion and recovers smooth regions as effectively as pure isotropic diffusion (in particular, without staircasing). The denoising performance results are tabulated in Table 1 where the best PSNR and MAE value is shown in boldface. The PSNR improvement brought by our approach can be quite high particularly for 𝜎=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.

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

8. Conclusions

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

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

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

Acknowledgments

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