# Well-Posed Inhomogeneous Nonlinear Diffusion Scheme for Digital Image Denoising

Authors: V. B. Surya Prasath; Arindama Singh
Journal: Journal of Applied Mathematics
 (2010)

---

## Abstract

We study an inhomogeneous partial differential equation which includes a separate edge detection part to control smoothing in and around possible discontinuities, under the framework of
anisotropic diffusion. By incorporating edges found at multiple scales
via an adaptive edge detector-based indicator function, the proposed
scheme removes noise while respecting salient boundaries. We create
a smooth transition region around probable edges found and reduce
the diffusion rate near it by a gradient-based diffusion coefficient. In
contrast to the previous anisotropic diffusion schemes, we prove the
well-posedness of our scheme in the space of bounded variation. The
proposed scheme is general in the sense that it can be used with any
of the existing diffusion equations. Numerical simulations on noisy
images show the advantages of our scheme when compared to other
related schemes.


---

## Body

## 1. Introduction

Anisotropic diffusion-based image denoising [1] for gray images was started by Perona and Malik [2] in 1990. It uses a scale parameter-based adaptive image filtering using nonlinear diffusion. Let u0 be a noisy version of the true image u with noise field n of known variance σ2:(1.1)u0=u+n.
Our aim is to recover the noise-free image u with edge preservation. To avoid the oversmoothing nature of the linear diffusion, Perona and Malik [2] introduced an edge indicator function g for reducing the diffusion near edges via the following anisotropic diffusion scheme (ADS):(1.2)∂u(x,t)∂t=div(g(|∇u(x,t)|)∇u(x,t)),forx∈Ω
with the initial condition u(x,0)=u0(x). The diffusion function g:ℝ+→ℝ+ is chosen to be decreasing, such that g(0)=1,lims→∞g(s)=0. If g(s)≡1 in (1.2), then we recover the linear diffusion PDE. The diffusion functions g(·) proposed in [2] are(1.3)g1(s)=exp(-(sK)2),g2(s)=11+(s/K)2,
where K>0 is the so-called contrast parameter. Better numerical results are obtained in [2] using this class of anisotropic PDEs (1.2)-(1.3) instead of the linear diffusion. When |∇u|>K, the PDE (1.2) turns into inverse diffusion, which is known to be ill-posed [3]. It is easily seen that the role of g as an edge indicator function is important in reducing the noise and enhancing true edges. A more robust indicator function based on smoothed gradient, g(|∇Gρ⋆u|), where the Gaussian kernel of width ρ is given by Gρ(x)=(2πρ2)-1exp(-|x|2/2ρ2), can be used for edge identification. This presmoothing also alleviates the instability [4] associated with the ADS (1.2). However, edge localization is lost and noisy oscillations can remain in this method. Moreover, when the noise level is high, ADS can reduce the contrast of edges and reliance on the instable |∇u| alone can be a drawback.There are efforts to remedy the ill-posedness and to use better edge indicators in the diffusion process. Strong [5] used an adaptive parameter in total variation-based PDE using the relationship between energy minimization and PDEs [3]. Better results are obtained when compared to the classical total variation (TV) scheme of [6, 7], but the use of gradients alone and the noisy edge map provided by u0 can lead to noise amplification along edges. Recently, Ceccarelli et al. [8] devised well-posed ADS schemes by approximating the TV function. Kusnezow et al. [9] used a similar weight function of Strong [5] combined with linear diffusion for fast computations. Yu et al. [10] considered ADS in terms of kernel-based smoothing and proposed to use a diffusion function based on modified gradients. Barbu et al. [11] considered the variational-PDE problem in Sobolev space setting with different growth functions. Douiri et al. [12] use an edge-based weight in the regularization functional for diffuse optical tomography to control smoothing across edges. The exact edge information is used to compute the weight which in real images is not possible, since we do not know, a priori, the exact edge map of the true image. Indeed, this edge estimation problem is closely related to image denoising, and most of the edge detection schemes [13, 14] are based on this observation.Apart from staircasing artifacts in flat regions and noise amplification along edges, capturing multiscale edges are difficult in previous schemes. To circumvent the drawbacks a more stable and robust edge indicator function using the multiscale edge detectors as discussed in [13–15] can be used. Anisotropic diffusion scheme is usually employed as a preprocessing step, for noise removal and edge detection. Here, we reverse the trend and integrate the edge detection into the anisotropic scheme such as (1.2). The method can, in general, be used in any anisotropic PDE, and it provides a better coupling of edge detection and restoration of images under noise. We prove the well-posedness of the proposed scheme in the space of bounded variation BV(Ω) [16]. The analysis is based on monotone operators and approximation schemes. The discrete version is implemented and it satisfies the properties required for an edge preserving scale space. Preliminary numerical results were reported in [17].The rest of the paper is organized as follows. In Section2, we introduce the edge adaptive PDE scheme and in Section 3 its well-posedness is proved in the space of bounded variation functions. We show the numerical results on real and noisy images in Section 4. Finally, Section 5 concludes the paper.


## 2. Proposed Scheme

### 2.1. Motivation

The diffusion functiong in ADS (1.2) uses the information provided by magnitude of the gradient to reduce diffusion near edges. Figure 1 shows the effect of this on Lena test image. Original image (Figure 1(a)) is corrupted with additive Gaussian noise of variance 20 and is shown in Figure 1(b). ADS (1.2) uses only the diffusion coefficient g as in (1.3). ADS restores the image in a piecewise constant manner but reduces overall contrast and creates staircasing artifacts near edges; see Figure 4(a). The magnitude of the gradient image (Figure 1(c), where whiter pixels represent maximal value of |∇u|) shows that it is not reliable under noise, and loss of spatial coherence such as edge connectivity and continuity is high. Smoothed gradient (|∇Gρ⋆u|) method of [4] reduces this loss but still leaves some spurious pixels due to noise (Figure 1(d)) and true locations of the edges are not preserved. This makes the result lacking clear boundaries; see Figure 4(b). Also this use of isotropic diffusion is against the very principle of ADS (1.2), which is anisotropic in nature.Image gradients based edge indicators. (a) Original image. (b) Noisy image. (c)g1(|∇u|). (d) g1(|Gσ⋆u|).
(a)(b)(c)(d)


### 2.2. Edge Adaptive ADS

Instead of using only the edge indicator functiong alone to drive the diffusion process, we introduce a spatially adaptive term α in the anisotropic scheme (α-ADS):(2.1)∂u(x,t)∂t=div(α(x,t)g(|∇u(x,t)|)∇u(x,t))-λ(u-u0),
where α provides a pixelwise edge characterization. As we will see in Section 3 and in the numerical experiments (Section 4) this avoids the stability problems of the gradient magnitude-based schemes and overcomes the localization error of smoothed gradient method. The requirements for α are as follows: (1)
at edge pixels to reduce the diffusion:α(x,t)~0 on the edges,(2)
in flat regions to allow smoothing:α(x,t)~1 beyond the edges,(3)
near edge pixelsα(x,t) should vary smoothly between 0 and 1 so as to avoid spurious oscillations.We choose a term of the form(2.2)α(x,t)=1-E(x,t)
with E as a robust and smooth edge indicator function.


### 2.3. Choice of Edge Indicator Function

Various edge detectors [13–15] exist in digital image processing literature. These detectors are based on a principle of gradient maxima. Among the available edge detectors [14] such as Sobel's, Prewitt's, Laplacian of Gaussian (LOG), and Marr-Hildreth's Zero-crossing that use maxima points of gradient to decide about the edge pixels, Canny's edge detector [13, 15] has been proved to be the best in terms of spatial coherency and edge continuity. The general method of Canny involves the following steps: (1)
convolution of the given imageu by a Gaussian Gσ,(2)
estimating the second derivative∇2(u) using finite differences,(3)
again using convolution for∇2(u) with a small Gaussian kernel,(4)
thresholding the gradient of Step1,(5)
zero-crossings of Step3 displayed if threshold of Step 4 is achieved.Canny edge mapC(x,t) is a binary output with edge pixels marked as 1 and nonedge pixels with 0. Hence we use(2.3)E(x,t)=Gρ⋆C(x,t),
where Gρ smoothens the transition from detected edges to homogenous regions. Figure 2 shows a closeup of the edge map from the hat region of Lena image given by a rectangle in Figure 1(b). Clearly the Canny edge detector-based smoothing parameter E(x,t) gives a good localization of edges (Figure 2(b)) and is devoid of noisy pixels. Compare this with the |∇u| or the smoothed version |∇Gσ⋆u|-based edge indicators (Figures 1(c) and 1(d)) where spatial coherency is missing and edge continuity is lost. Also, in the diffusion process based on these measures, the noisy pixels found in homogenous regions are propagated resulting in staircase artifacts.Canny edge detector-based indicator function. (a) Canny binary outputC(x). (b) Smooth edge indicator function E(x) with ρ=2 creates a band around C(x) (c) Diffusion weight αg(|∇u|).
(a)(b)(c)One can use any of the mentioned edge detectors into the spatially adaptive termα(x,t) in (2.2). We choose Canny detector because of its superiority and its computational efficiency. This trade-off between optimal performance and computation time can be decided according to the problem at hand.


## 2.1. Motivation

The diffusion functiong in ADS (1.2) uses the information provided by magnitude of the gradient to reduce diffusion near edges. Figure 1 shows the effect of this on Lena test image. Original image (Figure 1(a)) is corrupted with additive Gaussian noise of variance 20 and is shown in Figure 1(b). ADS (1.2) uses only the diffusion coefficient g as in (1.3). ADS restores the image in a piecewise constant manner but reduces overall contrast and creates staircasing artifacts near edges; see Figure 4(a). The magnitude of the gradient image (Figure 1(c), where whiter pixels represent maximal value of |∇u|) shows that it is not reliable under noise, and loss of spatial coherence such as edge connectivity and continuity is high. Smoothed gradient (|∇Gρ⋆u|) method of [4] reduces this loss but still leaves some spurious pixels due to noise (Figure 1(d)) and true locations of the edges are not preserved. This makes the result lacking clear boundaries; see Figure 4(b). Also this use of isotropic diffusion is against the very principle of ADS (1.2), which is anisotropic in nature.Image gradients based edge indicators. (a) Original image. (b) Noisy image. (c)g1(|∇u|). (d) g1(|Gσ⋆u|).
(a)(b)(c)(d)


## 2.2. Edge Adaptive ADS

Instead of using only the edge indicator functiong alone to drive the diffusion process, we introduce a spatially adaptive term α in the anisotropic scheme (α-ADS):(2.1)∂u(x,t)∂t=div(α(x,t)g(|∇u(x,t)|)∇u(x,t))-λ(u-u0),
where α provides a pixelwise edge characterization. As we will see in Section 3 and in the numerical experiments (Section 4) this avoids the stability problems of the gradient magnitude-based schemes and overcomes the localization error of smoothed gradient method. The requirements for α are as follows: (1)
at edge pixels to reduce the diffusion:α(x,t)~0 on the edges,(2)
in flat regions to allow smoothing:α(x,t)~1 beyond the edges,(3)
near edge pixelsα(x,t) should vary smoothly between 0 and 1 so as to avoid spurious oscillations.We choose a term of the form(2.2)α(x,t)=1-E(x,t)
with E as a robust and smooth edge indicator function.


## 2.3. Choice of Edge Indicator Function

Various edge detectors [13–15] exist in digital image processing literature. These detectors are based on a principle of gradient maxima. Among the available edge detectors [14] such as Sobel's, Prewitt's, Laplacian of Gaussian (LOG), and Marr-Hildreth's Zero-crossing that use maxima points of gradient to decide about the edge pixels, Canny's edge detector [13, 15] has been proved to be the best in terms of spatial coherency and edge continuity. The general method of Canny involves the following steps: (1)
convolution of the given imageu by a Gaussian Gσ,(2)
estimating the second derivative∇2(u) using finite differences,(3)
again using convolution for∇2(u) with a small Gaussian kernel,(4)
thresholding the gradient of Step1,(5)
zero-crossings of Step3 displayed if threshold of Step 4 is achieved.Canny edge mapC(x,t) is a binary output with edge pixels marked as 1 and nonedge pixels with 0. Hence we use(2.3)E(x,t)=Gρ⋆C(x,t),
where Gρ smoothens the transition from detected edges to homogenous regions. Figure 2 shows a closeup of the edge map from the hat region of Lena image given by a rectangle in Figure 1(b). Clearly the Canny edge detector-based smoothing parameter E(x,t) gives a good localization of edges (Figure 2(b)) and is devoid of noisy pixels. Compare this with the |∇u| or the smoothed version |∇Gσ⋆u|-based edge indicators (Figures 1(c) and 1(d)) where spatial coherency is missing and edge continuity is lost. Also, in the diffusion process based on these measures, the noisy pixels found in homogenous regions are propagated resulting in staircase artifacts.Canny edge detector-based indicator function. (a) Canny binary outputC(x). (b) Smooth edge indicator function E(x) with ρ=2 creates a band around C(x) (c) Diffusion weight αg(|∇u|).
(a)(b)(c)One can use any of the mentioned edge detectors into the spatially adaptive termα(x,t) in (2.2). We choose Canny detector because of its superiority and its computational efficiency. This trade-off between optimal performance and computation time can be decided according to the problem at hand.


## 3. Existence and Uniqueness

LetΩ⊂ℝ2 be a bounded domain. The proposed spatially adaptive version α-ADS (2.1) can be considered as a gradient descent form of a regularization functional [3]. The proposed scheme (2.1) is equivalent to a descent method for solving the following energy minimization problem:(3.1)minu∈BV(Ω){ℰ(u)=λ2∫Ω|u-u0|2dx+∫Ωαϕ(|∇u|)dx},
with ϕ′(s)=g(s)s. For example, for the diffusion functions in (1.3) we get the regularization functions as (3.2)ϕ1(s)=1-exp(-(sK)2),ϕ2(s)=log(1+(sK)2).
Since these two functions are nonconvex, the minimization problem (3.1) and ADS (1.2) do not have a unique solution. In these nonconvex regularization function cases the corresponding energy minimization problem (3.1) does not have a solution at all, except when u=u0=constant; see [3] for an analysis between the PDE (1.2) and minimization problem (3.1). To obtain a well-posed scheme we restrict ourselves to the class of linear growth functions for ϕ:ℝ→[0,∞) satisfying the following.(H1)
ϕ is a nondecreasing, convex, and even function with ϕ(0)=0.(H2)
There exista>0 and b≥0 such that as-b≤ϕ(s)≤as+b, for all s.An important example in this class is the minimal surface functionϕ(s)=1+s2, which is related to the total variation (TV) regularization [6, 7]; see Figure 3. We assume the preliminary results about the space of functions of bounded variation BV(Ω); see [16]. Throughout this paper, we use the notation ϕ(x,∇u)=α(x)ϕ(∇u). We recall the definition of total variation of a function u∈L1(Ω).(a) Diffusion functionsg1, g2 given in (1.3) and TV function. (b) Corresponding regularization functions.
(a)(b)Results of (a) Perona and Malik scheme [2], (b) Catté et al. scheme [4], (c) Rudin et al. scheme [6], (d) Strong scheme [5], (e) Ceccarelli et al. scheme [8], (f) Kusnezow et al. scheme [9], (g) Yu et al. scheme [10], (h) Barbu et al. scheme [11], (i) Our α-ADS scheme.
(a)
PM [2](b)
CL [4](c)
RO [6](d)
ST [5](e)
CD [8](f)
KH [9](g)
YW [10](h)
BB [11](i)
Our schemeDefinition 3.1.
Letu∈L1(Ω); its total variation is defined as
(3.3)∫Ω|∇u|:=sup|w|≤1{∫Ωudivwdx:w∈C01(Ω)}.The total variation is in fact a Radon measure; see [1, page 50] for more details. Making use of the properties of the convex functions of measures from [18] we will prove that the α-ADS (2.1) has a unique weak solution in the following sense.Definition 3.2.
u∈L2(BV(Ω);[0,T]) is a weak solution of the PDE (2.1) if for every v∈L2(BV(Ω);[0,T]), and τ∈[0,T](3.4)∫0τ∫Ωϕ(x,∇u)dxdt-λ∫0τ∫Ω(u-u0)(v-u)dxdt≤∫0τ∫Ω∂u∂t(v-u)dxdt+∫0τ∫Ωϕ(x,∇v)dxdt,a.e.inτ.
Moreover, it satisfies ∂u/∂t∈L2(Ω×[0,T]) and u(x,0)=u0(x).The function spaceL2(BV(Ω);[0,T]) is the set of all functions w:Ω×[0,T]→ℝ such that, w(·,t)∈BV(Ω) for each t and w(x,·)∈L2([0,T]) for all x. We consider the following regularized linear growth function (see [7]) in the approximation problem of (2.1):(3.5)ϕϵ(x,s):={12ϵ∂ϕ(x,s)∂s|s=ϵs2-ϵ2∂ϕ(x,s)∂s|s=ϵ+ϕ(x,ϵ),ifs<ϵ,ϕ(x,s),ifϵ≤s≤1ϵ,ϵ2∂ϕ(x,s)∂s|s=1/ϵs2-12ϵ∂ϕ(x,s)∂s|s=ϵ+ϕ(x,1/ϵ),ifs>1ϵ.Lemma 3.3.
Let the regularization functionϕ satisfy the assumptions H1-H2. If {ui} converges to u in L1(Ω), then
(3.6)∫Ωϕ(x,∇u)dx≤lim infi→∞∫Ωϕ(x,∇ui)dx.Proof.
Letw∈C01(Ω) be such that |w(x)|≤α(x) for all x∈Ω. Then
(3.7)∫Ωudivwdx=limi→∞∫Ωuidivwdx≤lim infi→∞∫Ωα|∇ui|dx,
where the inequality follows from Definition 3.1 of total variation of u. Taking supremum over w gives the result.Note thatα defined in (2.2) is a smooth function, that is, α(x)∈C∞(Ω;(0,∞)) due to the Gaussian presmoothing; see (2.3). The above lemma can also be proved using a general assumption about the adaptive parameter such as α∈C0(Ω;(0,∞)); see also [19, Theorem 1]. We next need the following result from nonlinear semigroup theory; see Theorem 3.1 in [20].Theorem 3.4.
LetH be a Hilbert space, and let P(H) be its power set. Let 𝒜:H→P(H) be a maximal monotone operator and u0∈D(𝒜), where the domain is D(𝒜)={x∈H:𝒜x≠∅}. Then there exists a unique function u(t):[0,∞)→H such that
(3.8)0∈∂u∂t+𝒜(u),u(0)=u0.Lemma 3.5.
Consider the approximation PDE(3.9)-ϵΔu+∂u∂t-div(gϵ(x,|∇u|)∇u)+λ(u-u0δ)=0onΩ×[0,T],u=u0δonΩ×{0}
with
(3.10)gϵ(x,ξ)=1ξ∂ϕϵ(x,ξ)∂ξ,u0δ∈C∞(Ω̅),u0δ→u0inL2(Ω)with∥u0δ∥L∞(Ω)≤C∥u0∥L∞(Ω).
There exists a unique solution uϵδ∈L∞(W1,2(Ω);[0,T]) for the approximation PDE. The solution uϵδ is bounded in L∞(Ω) by the initial value u0, that is,
(3.11)∥uϵδ∥L∞(Ω×[0,T])≤C∥u0∥L∞(Ω).Proof.
Sinceu0∈BV(Ω)∩L∞(Ω), existence of the sequence u0δ is guaranteed. Let 𝒜ϵ(u)=λ(u-u0δ)-div(gϵ(x,|∇u|)∇u). As 𝒜ϵ is the derivative of a convex lower semicontinuous functional, it is a maximal monotone operator. The existence of uϵδ follows from Theorem 3.4.
To get theL∞ bound for the solution we proceed as follows. Let B=∥u0δ∥L∞(Ω) and f+=max(f,0). Multiply the approximation PDE by the term e-st(e-stuϵδ-B)+. Integrating the PDE, we obtain(3.12)ϵ∫Ω∇uϵδe-st·e-st∇uϵδdx+∫Ω∂uϵδ∂te-st(e-stuϵδ-B)+dx+∫Ω[gϵ(x,|∇uϵδ|)∇uϵδ]e-st·e-st∇uϵδdx+λ∫Ω(uϵδ-u0δ)e-st(e-stuϵδ-B)+dx=0.This implies that
(3.13)∫Ω∂uϵδ∂te-st(e-stuϵδ-B)+dx≤0
since the remaining integrals are all nonnegative. Let
(3.14)ℐ(t)=∫Ω|(e-stuϵδ-B)+|2dx.
Then ℐ(t) is decreasing, nonnegative, and ℐ(0)=0. Thus ℐ(t)=0 for all t. Subsequently,
(3.15)uεδ(t)≤Besta.e.onΩfors,t>0.
Letting s→0, we get uϵδ(t)≤∥u0δ∥L∞(Ω). Similarly multiplying the term e-st(-B-e-stuϵδ)+ with the approximation PDE and integrating, we obtain uϵδ(t)≥-∥u0δ∥L∞(Ω).Lemma 3.6.
The weak solution of the approximation PDEuϵδ satisfied the following inequality, for τ∈[0,T]:
(3.16)∫0τ∫Ω(∂uϵδ∂t)2dxdt+ϵ2∫Ω|∇uϵδ(t)|2dx+∫Ωϕ(x,∇uϵδ(t))dx≤ϵ2∫Ω|∇u0δ|2‍dx+∫Ωϕ(x,∇u0δ)dx+ϵ,a.einτ.
Moreover, the weak solutions uϵδ and ∂uϵδ/∂t are uniformly bounded in ϵ with respect to the L∞(Ω×[0,T]) and L2(Ω×[0,T]) norms, respectively.Proof.
The proof of the inequality follows from the identity(3.17)∫0τ∫Ω(∂uϵδ∂t)2dxdt+ϵ2∫Ω|∇uϵδ(t)|2dx+∫Ωϕϵ(x,∇uϵδ(t))dx+λ2∫Ω(uϵδ-u0δ)2dx=ε2∫Ω‍|∇u0δ|2dx+∫Ωϕϵ(x,∇u0δ)dxand ϕ(x,s)≤ϕϵ(x,s)≤ϕ(x,s)+ϵ for any s.
From Lemma3.5, it follows that uϵδ is uniformly bounded and
(3.18)∫0T∫Ω(∂uϵδ∂t)2dxdt+∫Ωϕ(x,∇uϵδ)dx≤ϵ2∫Ω|∇u0δ|2dx+∫Ωϕ(x,∇u0δ)dx+ϵ.
We have
(3.19)∫0T∫Ω(|∂uϵδ∂t|+ϕ(∇uϵδ))dxdt≤C∥u0δ∥L∞(Ω).Let{ϵn} be a sequence converging to 0. Consider the corresponding weak solutions of the approximation PDE {uϵnδ}. In the following lemma, we prove the convergence of a subsequence of weak solutions {uϵnkδ} as ϵnk→0.Lemma 3.7.
There exists a subsequence{uϵnkδ} such that as ϵnk→0 dddd
(3.20)uϵnkδ→uδstronglyinL1(Ω×[0,T]),a.e.inΩ×[0,T],∂uϵnkδ∂t⇀∂uδ∂tweaklyinL2(Ω×[0,T])
where uδ is a weak solution of the PDE (2.1) with uδ(0)=u0δ. Moreover ∥uδ∥L∞(Ω×[0,T])≤C∥u0∥L∞(Ω).Proof.
By Lemma3.6, uϵδ and ∂uϵδ/∂t are uniformly bounded. When ϵnk→0, for fixed δ>0, we have a convergent subsequence uϵnkδ→uδ in L1(Ω×[0,T]) for uδ∈L∞(BV(Ω)∩L∞(Ω);[0,T]), and (∂uϵnkδ/∂t)⇀(∂uδ/∂t) in L2(Ω×[0,T]).
Sinceuϵnkδ is the weak solution of the approximation PDE, it satisfies the corresponding weak solution formulation in Definition (3.4). That is, for every v∈L2(BV(Ω);[0,T]) and τ∈[0,T],
(3.21)∫0τ∫Ωϕϵnk(x,∇uϵnkδ)dxdt-λ∫0τ∫Ω(uϵnkδ-u0δ)(v-uϵnkδ)dxdt≤ϵnk2∫0τ∫Ω|∇v|2dxdt+∫0τ∫Ωϕϵnk(x,∇v)dxdt+∫0τ∫Ω∂uϵnkδ∂t(v-uϵnkδ)dxdta.e.inτ.
When ϵnk→0, we obtain, using Lemma 3.3,
(3.22)∫0τ∫Ωϕ(x,∇uδ)dxdt-λ∫0τ∫Ω(uδ-u0δ)(v-uδ)dxdt≤∫0τ∫Ωϕ(x,∇v)dxdt+∫0τ∫Ω∂uδ∂t(v-uδ)dxdt.
Thus, uδ is a weak solution of the PDE (2.1) with uδ(x,0)=u0δ(x).Now, we state and prove our main result.Theorem 3.8.
Ifu0∈BV(Ω)∩L∞(Ω), then there exists a unique weak solution u∈L∞(BV(Ω)∩L∞(Ω);[0,T]) of (2.1).Proof.
From Lemma3.7 we get uδ as a weak solution of the PDE (2.1) with uδ(x,0)=u0δ(x). Using Lemma 3.7 for the sequence {uδ} as δ→0, we get
(3.23)uδ→ustronglyinL1(Ω×[0,T]),a.e.inΩ×[0,T],∂uδ∂t⇀∂u∂tweaklyinL2(Ω×[0,T]).
Finally we let δ→0 in the following inequality:
(3.24)∫0τ∫Ωϕ(x,∇uδ)dxdt-λ∫0τ∫Ω(uδ-u0δ)(v-uδ)dxdt≤∫0τ∫Ωϕ(x,∇v)dxdt+∫0τ∫Ω∂uδ∂t(v-uδ)dxdt
to obtain
(3.25)∫0τ∫Ωϕ(x,∇u)dxdt-λ∫0τ∫Ω(u-u0)(v-u)dxdt≤∫0τ∫Ωϕ(x,∇v)dxdt+∫0τ∫Ω∂u∂t(v-u)dxdt.
We see that u∈L∞(BV(Ω)∩L∞(Ω);[0,T]) is a weak solution of (2.1). Uniqueness of the solution follows from the weak solution inequality (3.4) above.As a consequence of Theorem3.8, we also obtain a bound for the solution ∥u∥L∞≤C∥u0∥L∞, that is, the maximum principle holds for the solution of (2.1); this guarantees that no new structures are created and α-ADS (2.1) satisfies the scale space axiom. Note that though α∈C∞(Ω) (since Gρ∈C∞(Ω)) in (2.3), if we use the original ADS functions (1.3) we cannot obtain the well-posedness of the PDE (2.1), since the corresponding ϕ functions are not convex.


## 4. Numerical Experiments

We use finite differences to discretize the proposed edge adaptive PDE (2.1) for images normalized to be in the range [0,1]. We take h as the grid size, and Uijt as the intensity value u(i,j) at iteration t. Instead of the classical explicit scheme, which severely restricts the step size, we make use of the unconditionally stable semi-implicit scheme. In 1D with matrix-vector notation it reads as (4.1)Ut+1=[1-τA(Ut)]-1Ut,
where τ is the time step, A(Ut)=[aij(Ut)], and(4.2)aij(Ut):={γit+γjt2h2,j∈𝒩i,-∑k∈𝒩iγit+γkt2h2,j=i,0,otherwise
with γi=αigi and 𝒩i is the discrete neighborhood of the pixel i. For n-D images the semi-implicit scheme is written as(4.3)Ut+1=[1-τ∑l=1nAl(Ut)]-1Ut.
The matrix Al=(aijl)ij corresponds to the derivatives along the lth coordinate axis.We use the Canny edge detector with its default settings for (2.3) in MATLAB7.4, and the computations are done in a desktop computer with Pentium IV, 2.45 GHz processor. It took nearly a minute for 100 iterations of our scheme, since the computation of Canny edge detector is done at each iteration. Further reduction in execution time can be achieved if we are to use other faster edge detectors. The numerical example, given in Figure 4, shows the results obtained for the noisy Lena image given in Figure 1(b). We show results of applying our scheme (2.1) with regularized TV function g(s)=(ϵ2+s2)-1/2, ϵ=10-6 in Figure 4(i) and also compare with the following schemes: (I)
classical schemes: ADS (1.2) of Perona and Malik [2] (Figure 4(a)), Catté et al.'s smoothed gradient scheme [4] (Figure 4(b)), and Rudin et al.'s total variation scheme [6] (Figure 4(c)),(II)
adaptive schemes: Strong's adaptive total variation scheme [5] (Figure 4(d)), Ceccarelli et al.'s approximate TV scheme [8] (Figure 4(e)), Kusnezow et al.'s adaptive linear diffusion scheme [9] (Figure 4(f)), Yu et al.'s kernel-based ADS [10] (Figure 4(g)), and Barbu et al.'s variational PDE scheme [11] (Figure 4(h)).We use the peak signal-to-noise ratio (PSNR) for comparing our scheme with other schemes. The parameters are tuned to get the best possible PSNR value for each of the scheme compared. Figure5(a) shows the effect of the same against different noise levels (5≤σ≤25) for our scheme. For an m×n image u, the PSNR value of an estimated image û is given by (4.4)PSNR=10log10(m×n∑x∈Ω(u(x)-û(x))2)(dB).
We see that our scheme produces a stable PSNR value as the noise level increases and attains highest values among other related schemes as well. To show the strong smoothing along flat regions and edge preservation, Figure 5(b) shows a signal line taken across Lena image from Figure 4(i), whose position is indicated in Figure 1(a) by a line. As can be seen our scheme can remove small scale texture details as the edge detector in (2.2) is not sufficient to capture them. Incorporation of textural measures into α points at further improvements in this direction.(a) Noise levelσ versus PSNR (dB) values. (b) Profile line taken across Lena at pixel position 250.
(a)(b)


## 5. Conclusions

A class of well-posed inhomogeneous diffusion schemes for image denoising is studied in this paper. By integrating multiscale edge detectors into the divergence term of anisotropic diffusion PDE we obtain edge preserving restoration of noisy images. Unlike the classical anisotropic diffusion schemes, which use only gradients to find edge pixels, we have proposed to include an adaptive parameter computed from Canny edge detector. Using an approximation scheme and the theory of monotone operators, well-posedness of the proposed inhomogeneous PDE is proved in the space of functions of bounded variation. Numerical examples illustrate that the scheme performs well under noise. Comparison with related schemes is undertaken to show that the expected improvement can really be achieved. Similar analysis for higher-order diffusion schemes and comparing various edge detectors for performance evaluation provide future directions.



---

*Source: 763847-2010-04-27.xml*