# Variational Image Denoising Approach with Diffusion Porous Media Flow

Authors: Tudor Barbu
Journal: Abstract and Applied Analysis
 (2013)

---

## 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∂u/∂t=Δβ(u),  where β is a nonlinear continuous function of the form β(u)=λum, 0<m<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.


---

## Body

## 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 Section3. 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]:
(1)∂u∂t-div(a(∇u))=0,in(0,∞)×Ω,u(0,x)=g(x),x∈Ω,
where g is the original noised image and u=u(t,x) is the filtered image. Here Ω⊂R2 is the image domain, ∂Ω is the image boundary, and a:R2→R2 is a given continuous and monotone function. Very often, a is a gradient mapping of a=∇j, where j:R2→Ris a convex continuous function.To (1) boundary conditions of Dirichlet type: u=0 on (0,∞)×∂Ω or Neumann: a(∇u)υ=0 on (0,∞)×∂Ω are taken. The special case, with j(r)=|r|d=(r12+r22)1/2,  r={r1,r2}, 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:
(2)∂u∂t-Δβ(u)=0in(0,∞)×Ω,β(u)=0in(0,∞)×∂Ω.The functionβ can be multivalued but maximal monotone, that is, (β(r1)-β(r2))(r1-r2)≥0 for all r1,r2≥0 and the range of r→r+β(r) is all of R. This equation is well known as mathematical model for fluid diffusion in porous media and the typical example is β(u)=λ|u|m-1u where λ>0 and 0<m<∞. The case m>1 models low diffusion processes and 0<m<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Ω⊂R2, we denote by L2(Ω), Wk,p(Ω), k=1,2, 1≤p≤∞, the standard spaces of Lebesgue integrable functions and the Sobolev spaces on Ω. We set Hk(Ω)=Wk,2(Ω),  H01(Ω)={u∈H1(Ω);u=0on∂Ω}, and H-1(Ω) the dual space of H01(Ω) with the norm denoting ||·||H-1(Ω). Consider here (2) with Dirichlet homogeneous boundary conditions, that is,
(3)∂u∂t(t,x)-Δβ(u(t,x))=0in(0,∞)×Ω,
with u(0,x)=g(x),x∈Ω, and β(u(t,x))=0 on (0,∞)×∂Ω, under the following assumption: β:R→R is monotonically increasing β(0)=0 and limr→±∞β(r)=±∞.HereΩ represents a bounded domain of R2 with a sufficiently smooth boundary, ∂Ω. By strong solution u=u(t,x) to (3) on [0,T], we mean an absolutely continuous function u:[0,T]→H-1(Ω) such that the strong derivative (du/dt)(t)=limε→0(u(t+ε)-u(t))/ε in H-1(Ω) exists a.e. on (0,T),  β(u(t,·))∈H01(Ω) a.e. t∈(0,T) and
(4)du(t)dt-Δβ(u(t,·))=0,a.e.t∈(0,T),
where u(0)=g, 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 A:D(A)⊂H-1(Ω)→H-1(Ω) defined by
(5)Au=-Δβ(u),∀u∈D(A)
with D(A)={u∈L1(Ω);β(u)∈H01(Ω)} and rewrite (4) as the infinite dimensional Cauchy problem:
(6)dudt(t)+Au(t)=0,a.e.t∈(0,T).We note [13, 14] that the nonlinear operator is maximal monotone (equivalently, m-accretive) in H-1(Ω) and that it is just the subdifferential ∂φ of the function (7)φ(u)={∫Ωη(u(x))dx,ifu∈L1(Ω)∩H-1(Ω),h(u)∈L1(Ω)+∞,otherwisein the sense of convex analysis. Here,η(r)=∫0rβ(s)ds. As a matter of fact, the distribution space H-1(Ω) is the basic functional space to study problem (4). Then, we have the following.Theorem 1.
Letg∈H-1(Ω). Then, (4) (equivalently, (6)) has a unique strong solution u:[0,∞)→H-1(Ω) satisfying
(8)tdudt∈L∞(0,T;H-1(Ω)),tβ(u)∈L∞(0,T;H01(Ω)).
Moreover, u is given by the exponential formula
(9)u(t)=limn→∞(I+tnA)-nguniformlyon[0,T].It should be recalled (see [13]) that this existence result extends to nonhomogeneous Cauchy problem du/dt+Au(t)=f(t),  t∈(0,T),u(0)=g.Theorem1 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,
(10)uk+1+hAuk+1=uk,u0=g,k=0,1,…,
where k=[t/h], is convergent to u(t). In the linear case, that is, β(u)=u, 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 flowu(t)=S(t)g is a continuous semigroup of contractions in H-1(Ω) in this case or in L1(Ω), in other situation, and the denoising effect of (6) is due to smoothing effect of S(t) on initial data, a property which is a characteristic to semigroups generated by the nonlinear operators A of the gradient type ∂φ, where φ is a convex and lower semicontinuous function. Since Theorem 1 applies to any initial data g∈H-1(Ω), which is a quite general space of distributions on Ω, the filter g→S(t)g can be used for smoothing very noised images g.It should be emphasized that, for filtering purpose, it suffices to computeu(t) at any time t>0 since the smoothing effect does not increase with t. On the contrary, since for t→0 we have [13]
(11)limt→∞u(t)=u∞inH-1(Ω),
where u∞ is a steady-state solution to A, that is, Au∞=0, one might expect that, for a large t, the filtered image could be far away from the original image g. This means that in the algorithm (10), the best filtered image can be obtained for k→∞,  kh≈ε, and ε arbitrarily small. If one wants to restore the image g via step descent algorithm, then (6) should be replaced by
(12)dudt(t)+Au(t)+λ(u(t)-g)=0,t≥0,λ>0
or, in its discrete version,
(13)uk+1+hAuk+1+λh(uk+1-g)=uk,k=0,1,…,
because, by (11), limt→∞u(t)=u∞, where Au∞+λ(u∞-g)=0, that is,(14)u∞=argmin{∫Ωη(u(x))dx+λ2∥u-g∥H-1(Ω)2}.This procedure was used in [7] for restoring and inpainting significantly degraded images g. We also have the generation formula limt→0(g-S(t)g)/t=AginH-1(Ω). This means that the operator-A is the infinitesimal generator of the semigroup S(t) [12]. The smoothing performance of the filter S(t) and its standard deviation is made precisely below.Theorem 2.
Letg∈L1(Ω) be such that β(g)∈H01(Ω). Then
(15)g-S(t)g=-tΔβ(g)+o(t)inH-1(Ω).
If j(g)∈L1(Ω), then
(16)∥g-S(t)g∥H-12(Ω)≤t∥j(g)∥L1(Ω),∀t≥0.Proof.
Formula (15) is immediate by the continuity of the function t→-Δβ(S(t)g) from [0,T] to H-1(Ω). As regards formula (16), we have g-S(t)g=-∫0tΔβ(S(τ)g)dτ,  forallt≥0, and this yields ∥g-S(t)g∥H-1(Ω)2≤t∫0t∥β(S(τ)g)∥H-1(Ω)2dτ,  forallt≥0, while by equation (d/dτ)S(τ)-Δβ(S(τ)g)=0,  forallτ>0, we have, by multiplying with β(S(τ)g) in H-1(Ω), ∫0t||β(S(τ)g)||H01(Ω)2dτ+∫Ωj(S(t)g)dx=∫Ωj(g)dx, for all t≥0.We also note that for (4) the maximum principle is also valid. More precisely, we have the following.Theorem 3.
Assume thatg≥0, a.e. in Ω. Then, u(t,x)≥0  a.e.  (t,x)∈(0,∞)×Ω.Proof.
We multiply (4) by (β(u))- and integrate on (0,t)×Ω to get u-≡0, as claimed.For denoising procedure, the filteru(t)=S(t)g, generated by (2) in the fast diffusion case, that is, for β(u)=λum,  0<m<1 (like the classical Gaussian filter), has the disadvantage that it is not localized as in this case the solution u=u(t,x) to (2) propagates with infinite speed. Moreover, u(t) has extinction in a finite time which depends on the norm of g. The situation is different in the low diffusion case β(u)=λum,  m>1, where the propagation is with finite speed, and so the denoising flow is localized.In (2) one might consider the Neumann boundary condition β(∇u)υ=0 on (0,∞)×∂Ω. In this case the operator A is given by (5) where D(A)={u∈L1(Ω),β(u)∈H1(Ω)} and Theorem 1 remains valid in the present situation; the basic space is H=(H1(Ω))′ instead of H-1(Ω).


## 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, uk+1=hΔ(uk-1)1/α+uk,  k=0,1,…,g=u0. We take g=gi,j≈g(ih,jh), where i∈[1,M] and j∈[1,N], and m=1/α. Thus, the finite difference scheme leads to the following iterative process:
(17)ui,jk+1=ui,jk+λ·((ui+1,jk-1)1/α+(ui-1,jk-1)1/α+(ui,j-1k-1)1/αui,jk+λ·+(ui,j+1k-1)1/α-4(ui,jk-1)1/α),u(i,1)=u(1,j)=0,∀i,j,
where λ=h/4l2,  ui,j0=gi,j, and k=1,…,K.Using the iterative scheme (17) with some properly selected parameters α, λ, and K, respectively, the initial noised image is successfully filtered in K steps. Therefore, the obtained uK represents the final image enhancement result. The choice of K 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 aK> 40 value, could produce a blurring effect on the processed image, while considering a very small K value, such as K < 5, could result in an unsatisfactory image denoising result. A great K 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 u∈L1(Ω)∩H-1(Ω) and η is the potential function corresponding to nonlinear diffusivity function β:R→R. The function u→||u-g||H-1(Ω)2 is a penalty term which forces the restored image u=u(x) to stay close to the initial image g=g(x).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:
(18)min{∫Ωj(∇u)dx+λ2∫Ω∥u-g∥2dx;u∈H01(Ω)}.The fact that the distance fromu to g is taken in the norm ||·||H-1(Ω) which is considerably weaker than the L2(Ω)-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 u to minimization problem satisfies the nonlinear elliptic equation:
(19)-Δβ(u)+λ(u-g)=0inΩ
with boundary conditions β(u)∈H01(Ω) (or Neumann flux boundary conditions as the basic space is (H1(Ω))′, the dual of H1(Ω), instead of H01(Ω)). As seen earlier, the minimization problem can be replaced by the evolution equations (4) or (13)).This means that∇β(u)=β′(u)∇u∈L2(Ω) 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 A0:D(A0)⊂L2(Ω)→L2(Ω) the operator A0u=-Δu,  D(A0)=H01(Ω)∩H2(Ω) and recall that A0 is an isometry from H01(Ω) to H0-1(Ω) for u=A01/2y, we can rewrite the minimization problem as
(20)Min{λ2∫Ωh(A01/2y)dxMin+λ2∥y-A0-1/2g∥L2(Ω)2;y∈D(A01/2)},
where A01/2 is the square root of A0. The Euler-Lagrange optimality condition for variational problem (20) can be expressed as
(21)A01/2β(A01/2y)+λ(y-A0-1/2g)=0.
Keeping in mind that A01/2 is defined by
(22)∫ΩA01/2uA01/2vdx=-∫ΩΔu(x)v(x)dx,u∈D(A0),v∈D(A01/2),
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
(23)β(r)=sgnr={1ifr>0,-1ifr<0,r∈R,[-1,1]ifr=0.
Indeed, in this case, (21) reduces to
(24)A1/2(A1/2|A1/2y|)+λ(y-A1/2g)=0on{x∈Ω;A1/2y(x)≠0}.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 Figures1 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:
(25)∂∂tb(x,u)-div(g(|∇Gσ*β(x,u)|)∇β(x,u))=f(u0-u)in(0,T)×Ω
with Neumann boundary condition ∇β(x,u)·υ=0 on (0,T)×Ω.Gaussian noise removal results produced by our technique and other approaches.
(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 modelNoise removal results provided by our model and other denoising techniques.
(a)
Original Baboon image
(b)
Gaussian noise
(c)
2D Gaussian filtering
(d)
Average filtering
(e)
Perona-Malik scheme
(f)
Our PDE based ModelHereb(x,·) and β(x,·) are continuous and monotonically increasing, g and f are Lipschitzian, and Gσ is a smoothing kernel. In the special case f=0,  g(r)=r, and b, β independent of x problem (25) reduces to (2) with Neumann boundary conditions and β replaced by β·b-1. In general, (25) can be rewritten as
(26)∂y∂t-g(|∇Gσ*β(x,u)|)Δβ(x,u)-∇g(|∇Gσ*β(x,u)|)∇β(x,u)=f(u0-b-1(u))in(0,T)×Ω
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 N=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 imageg and its restored version u, having the [M×N] size, the error value ∑x=1M∑y=1N(u(x,y)-g(x,y))2 is calculated. One can see the values of this norm of the image error parameter in the next tables.In Figure1(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 Figure1(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 Table1. 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.
Denoising technique
Gaussian filter
Average filter
Perona-Malik filter
Proposed AD filter
Norm of the error
6.40
×
1
0
3
6.05
×
1
0
3
6.84
×
1
0
3
5.15
×
1
0
3Another image denoising method comparison example is displayed in Figure2. 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.
Denoising technique
Gaussian filter
Average filter
Perona-Malik filter
Proposed AD filter
Norm of the error
7.20
×
1
0
3
6.14
×
1
0
3
6.55
×
1
0
3
5.36
×
1
0
3Our 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 K 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 asspeckle [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.



---

*Source: 856876-2013-01-30.xml*