Variational Image Denoising Approach with Diffusion Porous Media Flow

Abstract

A novel PDE-based image denoising approach is proposed in this paper. One designs here a nonlinear filter for image noise reduction based on the diffusion flow generated by the porous media equation 𝜕𝑢/𝜕𝑡=Δ𝛽(𝑢), where 𝛽 is a nonlinear continuous function of the form 𝛽(𝑢)=𝜆𝑢^𝑚, 0<𝑚<1. With respect to standard 2D Gaussian smoothing and some nonlinear PDE-based filters, this one is more efficient to remove noise from degraded images and also to reduce “staircasing” effects and preserve the image edges.

1. Introduction

Image noise removal constitutes a very important process, often representing the first step of image analysis [1, 2]. Image denoising with feature preservation represents still a focus in the image processing domain, remaining a serious challenge for researchers. A robust denoising technique has to successfully reduce the amount of noise while preserving the edges and has no blurring effect on the processed image.

The most common noise results from the image acquisition system can be modeled as Gaussian random noise in most cases. Gaussian noise represents statistical noise having the probability density function equal to that of the normal distribution. Numerous image denoising models have been introduced in the past few decades. The conventional filters such as averaging filter, median filter, and 2D Gaussian filter are efficient in smoothing the noise, but also have the disadvantage of blurring image edges [1, 2]. For this reason, some nonlinear image filters, which produce more satisfactory noise reduction results and preserve better the integrity of edges and detail information, were introduced [3].

The nonlinear noise removal techniques based on Partial Differential Equations (PDEs) have been extensively studied in the last two decades [4–6]. The PDE models are the best candidates for a very efficient image denoising. We proposed several robust PDE-based models for noise removal and image restoration in our previous works [7–9]. In this paper we develop a novel PDE-based image denoising approach based on nonlinear diffusion.

The idea of using the diffusion in image processing arose from the use of the Gaussian filter in multiscale image analysis. Nonlinear diffusion methods reduce noise and enhance contours in images [8–11]. Numerous nonlinear PDE denoising approaches based on diffusion have been introduced since the early work of Perona and Malik in 1987 [10]. The anisotropic diffusion developed by them, also called Perona-Malik diffusion, was intended to smooth the image while preserving its edges. Many denoising schemes based on this influential work have been proposed in the last 25 years [10, 11]. More closely of our approach is perhaps Kacur and Mikula scheme [12] which we will briefly present later on.

In the next section of this paper we describe the porous media diffusion filtering model. Then, the explicit version of the fast diffusion filtering scheme is presented in Section 3. Our denoising technique is compared with Perona-Malik scheme and Kacur-Mikula scheme in the fourth section. Other method comparisons and the numerical experiments performed using the proposed approach are discussed in Section 5. The paper ends with a section of conclusions.

2. The Porous Media Diffusion Filtering

Numerous PDE image denoising and segmentation approaches are based on the following nonlinear diffusion equation [10, 11]: 𝜕𝑢(𝑎(()𝑢()(𝑥)𝜕𝑡−div∇𝑢))=0,in0,∞×Ω,0,𝑥=𝑔,𝑥∈Ω, (1) where 𝑔 is the original noised image and 𝑢=𝑢(𝑡,𝑥) is the filtered image. Here Ω⊂𝑅² is the image domain, 𝜕Ω is the image boundary, and 𝑎∶𝑅²→𝑅² is a given continuous and monotone function. Very often, 𝑎 is a gradient mapping of 𝑎=∇𝑗, where 𝑗∶𝑅²→𝑅is a convex continuous function.

To (1) boundary conditions of Dirichlet type: 𝑢=0 on (0,∞)×𝜕Ω or Neumann: 𝑎(∇𝑢)𝜐=0 on (0,∞)×𝜕Ω are taken. The special case, with 𝑗(𝑟)=|𝑟|_𝑑=(𝑟²₁+𝑟²₂)^1/2,  𝑟={𝑟₁,𝑟₂}, leads to the so-called total variation model. There is a large variety of nonlinear PDE models inspired by (1) which attempt to eliminate the staircasing effect or to recover images with edges and the total variation model is one of the most efficient. Here we propose a nonlinear filter based on the flow generated by the porous media equation: 𝜕𝑢(𝑢)()𝛽(𝑢)()𝜕𝑡−Δ𝛽=0in0,∞×Ω,=0in0,∞×𝜕Ω. (2)

The function 𝛽 can be multivalued but maximal monotone, that is, (𝛽(𝑟₁)−𝛽(𝑟₂))(𝑟₁−𝑟₂)≥0 for all 𝑟₁,𝑟₂≥0 and the range of 𝑟→𝑟+𝛽(𝑟) is all of 𝑅. This equation is well known as mathematical model for fluid diffusion in porous media and the typical example is 𝛽(𝑢)=𝜆|𝑢|^𝑚−1𝑢 where 𝜆>0 and 0<𝑚<∞. The case 𝑚>1 models low diffusion processes and 0<𝑚<1 fast diffusion. We note that the model (1) was already used in our previous works for restoring (inpainting) the grayscale images [7, 8].

Given the domain Ω⊂𝑅², we denote by 𝐿²(Ω), 𝑊^𝑘,𝑝(Ω), 𝑘=1,2, 1≤𝑝≤∞, the standard spaces of Lebesgue integrable functions and the Sobolev spaces on Ω. We set 𝐻^𝑘(Ω)=𝑊^𝑘,2(Ω),  𝐻¹₀(Ω)={𝑢∈𝐻¹(Ω);𝑢=0on𝜕Ω}, and 𝐻⁻¹(Ω) the dual space of 𝐻¹₀(Ω) with the norm denoting ||⋅||_𝐻⁻¹_(Ω). Consider here (2) with Dirichlet homogeneous boundary conditions, that is, 𝜕𝑢()(𝑢(()𝜕𝑡𝑡,𝑥−Δ𝛽𝑡,𝑥))=0in0,∞×Ω, (3) with 𝑢(0,𝑥)=𝑔(𝑥),𝑥∈Ω, and 𝛽(𝑢(𝑡,𝑥))=0 on (0,∞)×𝜕Ω, under the following assumption: 𝛽∶𝑅→𝑅 is monotonically increasing 𝛽(0)=0 and lim_𝑟→±∞𝛽(𝑟)=±∞.

Here Ω represents a bounded domain of 𝑅² with a sufficiently smooth boundary, 𝜕Ω. By strong solution 𝑢=𝑢(𝑡,𝑥) to (3) on [0,𝑇], we mean an absolutely continuous function 𝑢∶[0,𝑇]→𝐻⁻¹(Ω) such that the strong derivative (𝑑𝑢/𝑑𝑡)(𝑡)=lim_𝜀→0(𝑢(𝑡+𝜀)−𝑢(𝑡))/𝜀 in 𝐻⁻¹(Ω) exists a.e. on (0,𝑇),  𝛽(𝑢(𝑡,⋅))∈𝐻¹₀(Ω) a.e. 𝑡∈(0,𝑇) and 𝑑𝑢(𝑡)(𝑢((),𝑑𝑡−Δ𝛽𝑡,⋅))=0,a.e.𝑡∈0,𝑇 (4) where 𝑢(0)=𝑔, and the Laplace operator Δ is considered in the sense of distributions on domain Ω. It should be said that a strong solution to (4) is also a weak solution in the classical sense. Let us consider the nonlinear operator 𝐴∶𝐷(𝐴)⊂𝐻⁻¹(Ω)→𝐻⁻¹(Ω) defined by 𝐴𝑢=−Δ𝛽(𝑢),∀𝑢∈𝐷(𝐴) (5) with 𝐷(𝐴)={𝑢∈𝐿¹(Ω);𝛽(𝑢)∈𝐻¹₀(Ω)} and rewrite (4) as the infinite dimensional Cauchy problem: 𝑑𝑢(𝑡)(𝑡)().𝑑𝑡+𝐴𝑢=0,a.e.𝑡∈0,𝑇 (6)

We note [13, 14] that the nonlinear operator is maximal monotone (equivalently, m-accretive) in 𝐻⁻¹(Ω) and that it is just the subdifferential 𝜕𝜑 of the function 𝜑(𝑢)=∫_Ω𝜂(𝑢(𝑥))𝑑𝑥,if𝑢∈𝐿¹(Ω)∩𝐻⁻¹(Ω),ℎ(𝑢)∈𝐿¹(Ω)+∞,otherwise (7)

in the sense of convex analysis. Here, 𝜂(𝑟)=∫^𝑟₀𝛽(𝑠)𝑑𝑠. As a matter of fact, the distribution space 𝐻⁻¹(Ω) is the basic functional space to study problem (4). Then, we have the following.

Theorem 1. Let 𝑔∈𝐻⁻¹(Ω). Then, (4) (equivalently, (6)) has a unique strong solution 𝑢∶[0,∞)→𝐻⁻¹(Ω) satisfying 𝑡𝑑𝑢𝑑𝑡∈𝐿^∞0,𝑇;𝐻⁻¹(Ω),𝑡𝛽(𝑢)∈𝐿^∞0,𝑇;𝐻¹₀(Ω). (8) Moreover, 𝑢 is given by the exponential formula 𝑢(𝑡)=lim_𝑛→∞𝑡𝐼+𝑛𝐴^−𝑛𝑔uniformly𝑜𝑛[0,𝑇]. (9)

It should be recalled (see [13]) that this existence result extends to nonhomogeneous Cauchy problem 𝑑𝑢/𝑑𝑡+𝐴𝑢(𝑡)=𝑓(𝑡),  𝑡∈(0,𝑇),𝑢(0)=𝑔.

Theorem 1 is a consequence of the standard existence theory for the infinite dimensional Cauchy problem associated with maximal monotone (m-accretive) operators in Hilbert spaces [13, 14]. It should be emphasized that Theorem 1 applies to (1) as well as to other nonlinear diffusion equations of monotone type arising in the filtering theory. In particular, (9) amounts to saying that the finite difference scheme, 𝑢_𝑘+1+ℎ𝐴𝑢_𝑘+1=𝑢_𝑘,𝑢₀=𝑔,𝑘=0,1,…, (10) where 𝑘=[𝑡/ℎ], is convergent to 𝑢(𝑡). In the linear case, that is, 𝛽(𝑢)=𝑢, this scheme is at origin of high-order Sobolev gradient method developed in [15, 16]. Formula (9) is known in the literature as the Crandall-Liggett exponential formula.

The flow 𝑢(𝑡)=𝑆(𝑡)𝑔 is a continuous semigroup of contractions in 𝐻⁻¹(Ω) in this case or in 𝐿¹(Ω), in other situation, and the denoising effect of (6) is due to smoothing effect of 𝑆(𝑡) on initial data, a property which is a characteristic to semigroups generated by the nonlinear operators 𝐴 of the gradient type 𝜕𝜑, where 𝜑 is a convex and lower semicontinuous function. Since Theorem 1 applies to any initial data 𝑔∈𝐻⁻¹(Ω), which is a quite general space of distributions on Ω, the filter 𝑔→𝑆(𝑡)𝑔 can be used for smoothing very noised images 𝑔.

It should be emphasized that, for filtering purpose, it suffices to compute 𝑢(𝑡) at any time 𝑡>0 since the smoothing effect does not increase with 𝑡. On the contrary, since for 𝑡→0 we have [13] lim_𝑡→∞𝑢(𝑡)=𝑢_∞in𝐻⁻¹(Ω), (11) where 𝑢_∞ is a steady-state solution to A, that is, 𝐴𝑢_∞=0, one might expect that, for a large t, the filtered image could be far away from the original image 𝑔. This means that in the algorithm (10), the best filtered image can be obtained for 𝑘→∞,  𝑘ℎ≈𝜀, and 𝜀 arbitrarily small. If one wants to restore the image 𝑔 via step descent algorithm, then (6) should be replaced by 𝑑𝑢(𝑡)(𝑡)𝑢(𝑡)𝑑𝑡+𝐴𝑢+𝜆−𝑔=0,𝑡≥0,𝜆>0 (12) or, in its discrete version, 𝑢_𝑘+1+ℎ𝐴𝑢_𝑘+1+𝜆ℎ𝑢_𝑘+1−𝑔=𝑢_𝑘,𝑘=0,1,…, (13) because, by (11), lim_𝑡→∞𝑢(𝑡)=𝑢_∞, where 𝐴𝑢_∞+𝜆(𝑢_∞−𝑔)=0, that is,𝑢_∞=argmin_Ω𝜂(𝑢(𝑥𝜆))𝑑𝑥+2𝑢−𝑔²_𝐻−1(Ω). (14)

This procedure was used in [7] for restoring and inpainting significantly degraded images 𝑔. We also have the generation formula lim_𝑡→0(𝑔−𝑆(𝑡)𝑔)/𝑡=𝐴𝑔in𝐻⁻¹(Ω). This means that the operator-𝐴 is the infinitesimal generator of the semigroup 𝑆(𝑡) [12]. The smoothing performance of the filter 𝑆(𝑡) and its standard deviation is made precisely below.

Theorem 2. Let 𝑔∈𝐿¹(Ω) be such that 𝛽(𝑔)∈𝐻¹₀(Ω). Then 𝑔−𝑆(𝑡)𝑔=−𝑡Δ𝛽𝑔+𝑜(𝑡)in𝐻⁻¹(Ω). (15) If 𝑗(𝑔)∈𝐿¹(Ω), then 𝑔−𝑆(𝑡)𝑔²_𝐻−1(Ω)≤𝑡𝑗𝑔_𝐿1(Ω),∀𝑡≥0. (16)

Proof. Formula (15) is immediate by the continuity of the function 𝑡→−Δ𝛽(𝑆(𝑡)𝑔) from [0,𝑇] to 𝐻⁻¹(Ω). As regards formula (16), we have 𝑔−𝑆(𝑡)𝑔=−∫^𝑡₀Δ𝛽(𝑆(𝜏)𝑔)𝑑𝜏,  forall𝑡≥0, and this yields ‖𝑔−𝑆(𝑡)𝑔‖²_𝐻−1(Ω)≤𝑡∫^𝑡₀‖𝛽(𝑆(𝜏)𝑔)‖²_𝐻−1(Ω)𝑑𝜏,  forall𝑡≥0, while by equation (𝑑/𝑑𝜏)𝑆(𝜏)−Δ𝛽(𝑆(𝜏)𝑔)=0,  forall𝜏>0, we have, by multiplying with 𝛽(𝑆(𝜏)𝑔) in 𝐻⁻¹(Ω), ∫^𝑡₀||𝛽(𝑆(𝜏)𝑔)||²_𝐻10(Ω)𝑑𝜏+∫_Ω𝑗(𝑆(𝑡)𝑔)𝑑𝑥=∫_Ω𝑗(𝑔)𝑑𝑥, for all 𝑡≥0.

We also note that for (4) the maximum principle is also valid. More precisely, we have the following.

Theorem 3. Assume that 𝑔≥0, a.e. in Ω. Then, 𝑢(𝑡,𝑥)≥0  𝑎.𝑒.  (𝑡,𝑥)∈(0,∞)×Ω.

Proof. We multiply (4) by (𝛽(𝑢))⁻ and integrate on (0,𝑡)×Ω to get 𝑢⁻≡0, as claimed.

For denoising procedure, the filter 𝑢(𝑡)=𝑆(𝑡)𝑔, generated by (2) in the fast diffusion case, that is, for 𝛽(𝑢)=𝜆𝑢^𝑚,  0<𝑚<1 (like the classical Gaussian filter), has the disadvantage that it is not localized as in this case the solution 𝑢=𝑢(𝑡,𝑥) to (2) propagates with infinite speed. Moreover, 𝑢(𝑡) has extinction in a finite time which depends on the norm of 𝑔. The situation is different in the low diffusion case 𝛽(𝑢)=𝜆𝑢^𝑚,  𝑚>1, where the propagation is with finite speed, and so the denoising flow is localized.

In (2) one might consider the Neumann boundary condition 𝛽(∇𝑢)𝜐=0 on (0,∞)×𝜕Ω. In this case the operator 𝐴 is given by (5) where 𝐷(𝐴)={𝑢∈𝐿¹(Ω),𝛽(𝑢)∈𝐻¹(Ω)} and Theorem 1 remains valid in the present situation; the basic space is 𝐻=(𝐻¹(Ω))^′ instead of 𝐻⁻¹(Ω).

3. The Explicit Scheme of the PDE Algorithm

Now, let us consider the explicit version of the scheme (10), obtained by the fast diffusion model. Namely, 𝑢^𝑘+1=ℎΔ(𝑢^𝑘−1)^1/𝛼+𝑢^𝑘,  𝑘=0,1,…,𝑔=𝑢⁰. We take 𝑔=𝑔_𝑖,𝑗≈𝑔(𝑖ℎ,𝑗ℎ), where 𝑖∈[1,𝑀] and 𝑗∈[1,𝑁], and 𝑚=1/𝛼. Thus, the finite difference scheme leads to the following iterative process: 𝑢^𝑘+1_𝑖,𝑗=𝑢^𝑘_𝑖,𝑗+𝜆⋅𝑢^𝑘−1_𝑖+1,𝑗^1/𝛼+𝑢^𝑘−1_{𝑖−1,𝑗}^1/𝛼+𝑢^𝑘−1_{𝑖,𝑗−1}^1/𝛼+𝑢^𝑘−1_𝑖,𝑗+1^1/𝛼−4𝑢^𝑘−1_𝑖,𝑗^1/𝛼,𝑢(𝑖,1)=𝑢1,𝑗=0,∀𝑖,𝑗, (17) where 𝜆=ℎ/4𝑙²,  𝑢⁰_𝑖,𝑗=𝑔_𝑖,𝑗, and 𝑘=1,…,𝐾.

Using the iterative scheme (17) with some properly selected parameters 𝛼, 𝜆, and 𝐾, respectively, the initial noised image is successfully filtered in 𝐾 steps. Therefore, the obtained 𝑢^𝐾 represents the final image enhancement result. The choice of 𝐾 in our simulation was dictated by the numerous tests we have performed in specific examples.

A proper selection of the parameter values is very important and cannot be a priori defined. It turns out that the selection of a very large number of iterations, for example, using a K > 40 value, could produce a blurring effect on the processed image, while considering a very small 𝐾 value, such as K < 5, could result in an unsatisfactory image denoising result. A great 𝐾 value increases also the computational complexity of this filtering process, producing a much higher computation time.

Also, using a large enough 𝜆 value, such as 𝜆>5, could increase the degradation of the noised image. A very small 𝜆 parameter, such as 𝜆<0.1, produces no visible smoothing results. Obviously, the parameter 𝛼 must satisfy the condition 1/𝛼∈(0,1), for a successfully noise removal.

4. Comparison with Other Noise Removal Techniques

The denoising model discussed above is generated by the minimization problem given by (14), where 𝑢∈𝐿¹(Ω)∩𝐻⁻¹(Ω) and 𝜂 is the potential function corresponding to nonlinear diffusivity function 𝛽∶𝑅→𝑅. The function 𝑢→||𝑢−𝑔||²_𝐻−1(Ω) is a penalty term which forces the restored image 𝑢=𝑢(𝑥) to stay close to the initial image 𝑔=𝑔(𝑥).

We compared our nonlinear diffusion-based noise removal technique with Perona-Malik denoising scheme [5, 10, 17] and other more general techniques (see [12]). These approaches based on (1) reduces to the minimization problem: min_Ω𝑗()𝜆∇𝑢𝑑𝑥+2_Ω𝑢−𝑔²𝑑𝑥;𝑢∈𝐻¹₀(Ω). (18)

The fact that the distance from 𝑢 to 𝑔 is taken in the norm ||⋅||_𝐻⁻¹_(Ω) which is considerably weaker than the 𝐿²(Ω)-norm considered in the Perona-Malik model as well as in most of the denoising models [9–11] has the advantage that it allows to work with very degraded initial images which practically are not represented by Lebesgue integrable functions but by distributions. However, it should be said that our model has a considerable better smoothing effect than that proposed by Perona and Malik. Indeed, as seen above, the solution 𝑢 to minimization problem satisfies the nonlinear elliptic equation: −Δ𝛽(𝑢)+𝜆𝑢−𝑔=0inΩ (19) with boundary conditions 𝛽(𝑢)∈𝐻¹₀(Ω) (or Neumann flux boundary conditions as the basic space is (𝐻¹(Ω))^′, the dual of 𝐻¹(Ω), instead of 𝐻¹₀(Ω)). As seen earlier, the minimization problem can be replaced by the evolution equations (4) or (13)).

This means that ∇𝛽(𝑢)=𝛽^′(𝑢)∇𝑢∈𝐿²(Ω) which indicates a smoothing effect comparable with the standard one, but for a suitable choice of 𝛽 one can avoid the “staircasing” effect which is common in denoising procedure with high smoothing effect.

A reformulation of problem (14) and, implicitly, of the dynamic model (3) allows a closer comparison of this model with that inspired by the Perona-Malik classical denoising technique [10]. If we denote by 𝐴₀∶𝐷(𝐴₀)⊂𝐿²(Ω)→𝐿²(Ω) the operator 𝐴₀𝑢=−Δ𝑢,  𝐷(𝐴₀)=𝐻¹₀(Ω)∩𝐻²(Ω) and recall that 𝐴₀ is an isometry from 𝐻¹₀(Ω) to 𝐻₀⁻¹(Ω) for 𝑢=𝐴₀^1/2𝑦, we can rewrite the minimization problem as Min_Ωℎ𝐴₀^1/2𝑦+𝜆𝑑𝑥2𝑦−𝐴₀^−1/2𝑔²_𝐿2(Ω);𝑦∈𝐷𝐴₀^1/2, (20) where 𝐴₀^1/2 is the square root of 𝐴₀. The Euler-Lagrange optimality condition for variational problem (20) can be expressed as 𝐴₀^1/2𝛽𝐴₀^1/2𝑦+𝜆𝑦−𝐴₀^−1/2𝑔=0. (21) Keeping in mind that 𝐴₀^1/2 is defined by _Ω𝐴₀^1/2𝑢𝐴₀^1/2𝑣𝑑𝑥=−_ΩΔ𝑢(𝑥)𝑣(𝑥)𝑑𝑥,𝑢∈𝐷𝐴₀,𝑣∈𝐷𝐴₀^1/2, (22) we may interpret (21) (resp., (20)) in an appropriate sense (1). From this perspective, the denoising model (19) is close to the “total variation model” (1) if 𝛽 is taken as 𝛽(𝑟)[]=sgn𝑟=1if𝑟>0,−1if𝑟<0,𝑟∈𝑅,−1,1if𝑟=0. (23) Indeed, in this case, (21) reduces to 𝐴^1/2𝐴^1/2𝐴^1/2𝑦+𝜆𝑦−𝐴^1/2𝑔=0on𝑥∈Ω;𝐴^1/2𝑦(𝑥)≠0. (24)

However, the present model is more convenient that the “total variation model,” which is constructed in a nonenergetic space (the space of function with bounded variation) and so hard to treat from the computational point of view. As a matter of fact, by the regularization necessary to construct a viable numerical scheme, the “total variation model” loses most of the theoretical advantages regarding the sharp edge detection and elimination of the staircasing effect.

The simulation developed in the next section confirms this important advantage (see the results depicted in Figures 1 and 2). We must emphasize also that there are numerous empirical denoising schemes obtained by modifying the original Perona-Malik model [5, 10, 11, 17–19], most of them not well posed from the mathematical point of view. The Kačur and Mikula denoising model [12] is based on the boundary value problem: 𝜕𝑏()𝑔𝜕𝑡𝑥,𝑢−div∇𝐺_𝜎(()𝑢∗𝛽𝑥,𝑢)∇𝛽𝑥,𝑢=𝑓₀−𝑢in(0,𝑇)×Ω (25) with Neumann boundary condition ∇𝛽(𝑥,𝑢)⋅𝜐=0 on (0,𝑇)×Ω.

(a) Original Lena image

(b) Image affected by Gaussian noise

(c) Classic Gaussian filtering result

(d) Average filtering result

(e) Perona-Malik denoising scheme

(f) Denoising result of the proposed model

Figure 1 

Gaussian noise removal results produced by our technique and other approaches.

(a) Original Baboon image

(b) Gaussian noise

(c) 2D Gaussian filtering

(d) Average filtering

(e) Perona-Malik scheme

(f) Our PDE based Model

Figure 2 

Noise removal results provided by our model and other denoising techniques.

Here 𝑏(𝑥,⋅) and 𝛽(𝑥,⋅) are continuous and monotonically increasing, 𝑔 and 𝑓 are Lipschitzian, and 𝐺_𝜎 is a smoothing kernel. In the special case 𝑓=0,  𝑔(𝑟)=𝑟, and b, 𝛽 independent of 𝑥 problem (25) reduces to (2) with Neumann boundary conditions and 𝛽 replaced by 𝛽⋅𝑏⁻¹. In general, (25) can be rewritten as 𝜕𝑦𝜕𝑡−𝑔∇𝐺_𝜎∗𝛽(𝑥,𝑢)Δ𝛽(𝑥,𝑢)−∇𝑔∇𝐺_𝜎(()𝑢∗𝛽𝑥,𝑢)∇𝛽𝑥,𝑢=𝑓₀−𝑏⁻¹(𝑢)in(0,𝑇)×Ω (26) with appropriate Neumann boundary conditions. It should be said that in this general form the Cauchy problem for (25) or (26) is not well posed and so the finite difference scheme (10) is not convergent. The model considered here is mathematically well posed. It should be mentioned however that in [12] a semi-implicit scheme is designed which is convergent to a weak solution to (25).

5. Numerical Experiments and Method Comparisons

The PDE-based image denoising technique proposed here has been tested on various image datasets, satisfactory filtering results being obtained. We have filtered hundreds of grayscale images affected by various levels of Gaussian noise using the described approach. An important advantage of our filtering technique is that it performs a robust noise reduction while preserving the image edges [20, 21].

We have obtained the best denoising results for the following set of parameters of the diffusion model provided by formula (18): 𝛼=2, which corresponds to the physical model of diffusion in plasma, 𝜆=1.5 and 𝑁=20. Some grayscale image reduction examples based on these parameter values are provided in Figures 1 and 2, respectively.

Numerous method comparisons have also been performed. The denoising performance of our method has been compared with performances of other noise removal techniques, such as the 2D Gaussian filter, the averaging filter, and Perona-Malik anisotropic diffusion scheme. From the performed numerical experiments it is obvious that the approach introduced here works better than these well-known filtering algorithms, as one can observe in the figures.

In order to assess the performance levels of each image denoising technique, one computes the norm of the error image. Thus, for any initial image 𝑔 and its restored version u, having the [𝑀×𝑁] size, the error value ∑^𝑀_𝑥=1∑^𝑁_𝑦=1(𝑢(𝑥,𝑦)−𝑔(𝑥,𝑦))² is calculated. One can see the values of this norm of the image error parameter in the next tables.

In Figure 1(a), the standard [512×512] image of Lena is displayed in the grayscale form. Then, it is corrupted by an amount of Gaussian noise characterized by parameters 0.2 (mean) and 0.02 (variance), the noising result being displayed in Figure 1(b).

In Figure 1(c) the image denoising result produced by the classic [3×3] Gaussian 2D filter kernel is displayed, while the noise reduction obtained with an [3×3] averaging filter kernel is represented in Figure 1(d). The noise removal produced by the Perona-Malik approach is displayed in Figure 1(e), while the denoising result provided by the nonlinear PDE model proposed here is represented in Figure 1(f).

The values of the norm of the error image corresponding to all these image filtering methods are displayed in Table 1. One can see that our anisotropic diffusion-based noise reduction approach performs better than the other filters, minimizing the respective error.

Table 1 

Values of norm of the error image parameter for various noise reduction approaches.

Another image denoising method comparison example is displayed in Figure 2. The same noise removal approaches (standard Gaussian, average filtering, Perona-Malik, and our technique) are applied on the standard Baboon image corrupted by the same amount of Gaussian noise, their results being represented in Figures 2(c)–2(f). The corresponding values of the norm of the error image parameter are displayed in Table 2. One can observe the minimum error value obtained by our proposed AD-based noise removal algorithm.

Table 2 

Values of norm of the error image parameter for various noise reduction approaches.

Our filtering algorithm increases the image quality not only by performing a robust Gaussian noise reduction, but also by enhancing the grayscale image contrast. Also, the proposed denoising technique runs quite fast, a [512×512] digital image being processed in less than 1 s. The time complexity decreases proportionally with the size of 𝐾 parameter.

6. Conclusions

A novel PDE-based image noise reduction technique has been proposed in this paper. We have designed a robust nonlinear image filter based on the flow generated by the porous media equation.

The provided nonlinear diffusion-based denoising technique performs successfully for grayscale images corrupted by 2D Gaussian noise, producing also an improved image contrast. Our method is also an edge preserving noise removal technique, which represents an important advantage.

It has been compared with some influential anisotropic diffusion approaches, like that proposed by Perona and Malik [5, 10], the obtained comparison results being very encouraging for us. Also, it provides better smoothing results than many other denoising techniques, while executing almost as fast as them, given its quite low-time complexity.

We have also tested this PDE-based model on other types of image noise and obtained mixed filtering results. While our technique performs somewhat well for images affected by some types of noise such as speckle [9], it does not succeed in removing properly some other types of noise, like salt and pepper. Therefore, designing some robust filters for other types of image noise will represent the focus of our future work in the image processing domain.

Because of its edge-preserving and contrast-enhancing character, the technique described in this paper can be successfully applied in image analysis domains requiring image object emphasizing and detection.

Acknowledgments

This work was mainly supported by the project PN-II-ID-PCE-2011-3-0027-160/5.10.2011 financed by UEFSCDI, Romania. It was supported also by the Institute of Computer Science, The Romanian Academy, Iasi, Romania.