Well-Posed Inhomogeneous Nonlinear Diffusion Scheme for Digital Image Denoising

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.

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 𝑢₀ be a noisy version of the true image 𝑢 with noise field 𝑛 of known variance 𝜎²:𝑢₀=𝑢+𝑛. (1.1) Our aim is to recover the noise-free image 𝑢 with edge preservation. To avoid the oversmoothing nature of the linear diffusion, Perona and Malik [2] introduced an edge indicator function 𝑔 for reducing the diffusion near edges via the following anisotropic diffusion scheme (ADS):𝜕𝑢(𝑥,𝑡)𝜕𝑡=div(𝑔(|∇𝑢(𝑥,𝑡)|)∇𝑢(𝑥,𝑡)),for𝑥∈Ω (1.2) with the initial condition 𝑢(𝑥,0)=𝑢₀(𝑥). The diffusion function 𝑔∶ℝ⁺→ℝ⁺ is chosen to be decreasing, such that 𝑔(0)=1,lim_𝑠→∞𝑔(𝑠)=0. If 𝑔(𝑠)≡1 in (1.2), then we recover the linear diffusion PDE. The diffusion functions 𝑔(⋅) proposed in [2] are𝑔₁𝑠(𝑠)=exp(−(𝐾)²),𝑔₂1(𝑠)=1+(𝑠/𝐾)², (1.3) where 𝐾>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 |∇𝑢|>𝐾, the PDE (1.2) turns into inverse diffusion, which is known to be ill-posed [3]. It is easily seen that the role of 𝑔 as an edge indicator function is important in reducing the noise and enhancing true edges. A more robust indicator function based on smoothed gradient, 𝑔(|∇𝐺_𝜌⋆𝑢|), where the Gaussian kernel of width 𝜌 is given by 𝐺_𝜌(𝑥)=(2𝜋𝜌²)⁻¹exp(−|𝑥|²/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 |∇𝑢| 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 𝑢₀ 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 Section 2, 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 function 𝑔 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 𝐿𝑒𝑛𝑎 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 𝑔 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 |∇𝑢|) shows that it is not reliable under noise, and loss of spatial coherence such as edge connectivity and continuity is high. Smoothed gradient (|∇𝐺_𝜌⋆𝑢|) 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.

(a)

(b)

(c)

(d)

Figure 1 

Image gradients based edge indicators. (a) Original image. (b) Noisy image. (c) 𝑔1(|∇𝑢|). (d) 𝑔1(|𝐺𝜎⋆𝑢|).

2.2. Edge Adaptive ADS

Instead of using only the edge indicator function 𝑔 alone to drive the diffusion process, we introduce a spatially adaptive term 𝛼 in the anisotropic scheme (𝛼-ADS):𝜕𝑢(𝑥,𝑡)𝜕𝑡=div(𝛼(𝑥,𝑡)𝑔(|∇𝑢(𝑥,𝑡)|)∇𝑢(𝑥,𝑡))−𝜆(𝑢−𝑢₀), (2.1) 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: 𝛼(𝑥,𝑡)∼0 on the edges,(2)in flat regions to allow smoothing: 𝛼(𝑥,𝑡)∼1 beyond the edges,(3)near edge pixels 𝛼(𝑥,𝑡) should vary smoothly between 0 and 1 so as to avoid spurious oscillations.

We choose a term of the form𝛼(𝑥,𝑡)=1−𝐸(𝑥,𝑡) (2.2) with 𝐸 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 image 𝑢 by a Gaussian 𝐺_𝜎,(2)estimating the second derivative ∇²(𝑢) using finite differences,(3)again using convolution for ∇²(𝑢) with a small Gaussian kernel,(4)thresholding the gradient of Step 1,(5)zero-crossings of Step 3 displayed if threshold of Step 4 is achieved.

Canny edge map 𝐶(𝑥,𝑡) is a binary output with edge pixels marked as 1 and nonedge pixels with 0. Hence we use𝐸(𝑥,𝑡)=𝐺_𝜌⋆𝐶(𝑥,𝑡), (2.3) where 𝐺_𝜌 smoothens the transition from detected edges to homogenous regions. Figure 2 shows a closeup of the edge map from the hat region of 𝐿𝑒𝑛𝑎 image given by a rectangle in Figure 1(b). Clearly the Canny edge detector-based smoothing parameter 𝐸(𝑥,𝑡) gives a good localization of edges (Figure 2(b)) and is devoid of noisy pixels. Compare this with the |∇𝑢| or the smoothed version |∇𝐺_𝜎⋆𝑢|-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.

(a)

(b)

(c)

Figure 2 

Canny edge detector-based indicator function. (a) Canny binary output 𝐶(𝑥). (b) Smooth edge indicator function 𝐸(𝑥) with 𝜌=2 creates a band around 𝐶(𝑥) (c) Diffusion weight 𝛼𝑔(|∇𝑢|).

One can use any of the mentioned edge detectors into the spatially adaptive term 𝛼(𝑥,𝑡) 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 Ω⊂ℝ² 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:min_{𝑢∈BV(Ω)}𝜆{ℰ(𝑢)=2_Ω|𝑢−𝑢₀|²𝑑𝑥+_Ω𝛼𝜙(|∇𝑢|)𝑑𝑥}, (3.1) with 𝜙^′(𝑠)=𝑔(𝑠)𝑠. For example, for the diffusion functions in (1.3) we get the regularization functions as 𝜙₁𝑠(𝑠)=1−exp(−(𝐾)²),𝜙₂𝑠(𝑠)=log(1+(𝐾)²). (3.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 𝑢=𝑢₀=𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡; 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 exist 𝑎>0 and 𝑏≥0 such that 𝑎𝑠−𝑏≤𝜙(𝑠)≤𝑎𝑠+𝑏, for all 𝑠.

An important example in this class is the minimal surface function 𝜙(𝑠)=√1+𝑠², 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 𝜙(𝑥,∇𝑢)=𝛼(𝑥)𝜙(∇𝑢). We recall the definition of total variation of a function 𝑢∈𝐿¹(Ω).

(a)

(b)

Figure 3 

(a) Diffusion functions 𝑔1, 𝑔2 given in (1.3) and TV function. (b) Corresponding regularization functions.

(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 scheme

Figure 4 

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.

Definition 3.1. Let 𝑢∈𝐿¹(Ω); its total variation is defined as _Ω|∇𝑢|∶=sup_|𝑤|≤1{_Ω𝑢div𝑤𝑑𝑥∶𝑤∈𝐶¹₀(Ω)}. (3.3)

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. 𝑢∈𝐿²(BV(Ω);[0,𝑇]) is a weak solution of the PDE (2.1) if for every 𝑣∈𝐿²(BV(Ω);[0,𝑇]), and 𝜏∈[0,𝑇]^𝜏_0Ω𝜙(𝑥,∇𝑢)𝑑𝑥𝑑𝑡−𝜆^𝜏_0Ω(𝑢−𝑢₀≤)(𝑣−𝑢)𝑑𝑥𝑑𝑡_𝜏0Ω𝜕𝑢𝜕𝑡(𝑣−𝑢)𝑑𝑥𝑑𝑡+^𝜏_0Ω𝜙(𝑥,∇𝑣)𝑑𝑥𝑑𝑡,a.e.in𝜏. (3.4) Moreover, it satisfies 𝜕𝑢/𝜕𝑡∈𝐿²(Ω×[0,𝑇]) and 𝑢(𝑥,0)=𝑢₀(𝑥).

The function space 𝐿²(BV(Ω);[0,𝑇]) is the set of all functions 𝑤∶Ω×[0,𝑇]→ℝ such that, 𝑤(⋅,𝑡)∈BV(Ω) for each 𝑡 and 𝑤(𝑥,⋅)∈𝐿²([0,𝑇]) for all 𝑥. We consider the following regularized linear growth function (see [7]) in the approximation problem of (2.1):𝜙_𝜖1(𝑥,𝑠)∶=2𝜖𝜕𝜙(𝑥,𝑠)|𝜕𝑠^𝑠=𝜖𝑠²−𝜖2𝜕𝜙(𝑥,𝑠)|𝜕𝑠^𝑠=𝜖1+𝜙(𝑥,𝜖),if𝑠<𝜖,𝜙(𝑥,𝑠),if𝜖≤𝑠≤𝜖,𝜖2𝜕𝜙(𝑥,𝑠)|𝜕𝑠^𝑠=1/𝜖𝑠²−12𝜖𝜕𝜙(𝑥,𝑠)|𝜕𝑠^𝑠=𝜖1+𝜙(𝑥,1/𝜖),if𝑠>𝜖. (3.5)

Lemma 3.3. Let the regularization function 𝜙 satisfy the assumptions 𝐻1−𝐻2. If {𝑢_𝑖} converges to 𝑢 in 𝐿¹(Ω), then _Ω𝜙(𝑥,∇𝑢)𝑑𝑥≤liminf_𝑖→∞Ω𝜙(𝑥,∇𝑢_𝑖)𝑑𝑥. (3.6)

Proof. Let 𝑤∈𝐶¹₀(Ω) be such that |𝑤(𝑥)|≤𝛼(𝑥) for all 𝑥∈Ω. Then _Ω𝑢div𝑤𝑑𝑥=lim_𝑖→∞Ω𝑢_𝑖div𝑤𝑑𝑥≤liminf_𝑖→∞Ω𝛼|∇𝑢_𝑖|𝑑𝑥, (3.7) where the inequality follows from Definition 3.1 of total variation of 𝑢. Taking supremum over 𝑤 gives the result.

Note that 𝛼 defined in (2.2) is a smooth function, that is, 𝛼(𝑥)∈𝐶^∞(Ω;(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 𝛼∈𝐶⁰(Ω;(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. Let 𝐻 be 𝑎 Hilbert space, and let 𝑃(𝐻) be its power set. Let 𝒜∶𝐻→𝑃(𝐻) be a maximal monotone operator and 𝑢₀∈𝐷(𝒜), where the domain is 𝐷(𝒜)={𝑥∈𝐻∶𝒜𝑥≠∅}. Then there exists a unique function 𝑢(𝑡)∶[0,∞)→𝐻 such that 0∈𝜕𝑢𝜕𝑡+𝒜(𝑢),𝑢(0)=𝑢₀. (3.8)

Lemma 3.5. Consider the approximation PDE −𝜖Δ𝑢+𝜕𝑢𝜕𝑡−div(𝑔_𝜖(𝑥,|∇𝑢|)∇𝑢)+𝜆(𝑢−𝑢^𝛿₀)=0onΩ×[0,𝑇],𝑢=𝑢^𝛿₀onΩ×{0} (3.9) with 𝑔_𝜖1(𝑥,𝜉)=𝜉𝜕𝜙_𝜖(𝑥,𝜉)𝜕𝜉,𝑢^𝛿₀∈𝐶^∞(𝑢Ω),_𝛿0→𝑢₀in𝐿²(Ω)with‖𝑢^𝛿₀‖_𝐿∞(Ω)≤𝐶‖𝑢₀‖_𝐿∞(Ω). (3.10) There exists a unique solution 𝑢^𝛿_𝜖∈𝐿^∞(𝑊^1,2(Ω);[0,𝑇]) for the approximation PDE. The solution 𝑢^𝛿_𝜖 is bounded in 𝐿^∞(Ω) by the initial value 𝑢₀, that is, ‖𝑢^𝛿_𝜖‖_{𝐿∞(Ω×[0,𝑇])}≤𝐶‖𝑢₀‖_𝐿∞(Ω). (3.11)

Proof. Since 𝑢₀∈BV(Ω)∩𝐿^∞(Ω), existence of the sequence 𝑢^𝛿₀ is guaranteed. Let 𝒜_𝜖(𝑢)=𝜆(𝑢−𝑢^𝛿₀)−div(𝑔_𝜖(𝑥,|∇𝑢|)∇𝑢). As 𝒜_𝜖 is the derivative of a convex lower semicontinuous functional, it is a maximal monotone operator. The existence of 𝑢^𝛿_𝜖 follows from Theorem 3.4.
To get the 𝐿^∞ bound for the solution we proceed as follows. Let 𝐵=‖𝑢^𝛿₀‖_𝐿∞(Ω) and 𝑓⁺=max(𝑓,0). Multiply the approximation PDE by the term 𝑒^−𝑠𝑡(𝑒^−𝑠𝑡𝑢^𝛿_𝜖−𝐵)⁺. Integrating the PDE, we obtain𝜖_Ω∇𝑢^𝛿_𝜖𝑒^−𝑠𝑡⋅𝑒^−𝑠𝑡∇𝑢^𝛿_𝜖𝑑𝑥+_Ω𝜕𝑢^𝛿_𝜖𝑒𝜕𝑡^−𝑠𝑡(𝑒^−𝑠𝑡𝑢^𝛿_𝜖−𝐵)⁺+𝑑𝑥_Ω[𝑔_𝜖(𝑥,|∇𝑢^𝛿_𝜖|)∇𝑢^𝛿_𝜖]𝑒^−𝑠𝑡⋅𝑒^−𝑠𝑡∇𝑢^𝛿_𝜖𝑑𝑥+𝜆_Ω(𝑢^𝛿_𝜖−𝑢^𝛿₀)𝑒^−𝑠𝑡(𝑒^−𝑠𝑡𝑢^𝛿_𝜖−𝐵)⁺𝑑𝑥=0. (3.12)This implies that _Ω𝜕𝑢^𝛿_𝜖𝑒𝜕𝑡^−𝑠𝑡(𝑒^−𝑠𝑡𝑢^𝛿_𝜖−𝐵)⁺𝑑𝑥≤0 (3.13) since the remaining integrals are all nonnegative. Let ℐ(𝑡)=_Ω|(𝑒^−𝑠𝑡𝑢^𝛿_𝜖−𝐵)⁺|²𝑑𝑥. (3.14) Then ℐ(𝑡) is decreasing, nonnegative, and ℐ(0)=0. Thus ℐ(𝑡)=0 for all 𝑡. Subsequently, 𝑢^𝛿_𝜀(𝑡)≤𝐵𝑒^𝑠𝑡a.e.onΩfor𝑠,𝑡>0. (3.15) Letting 𝑠→0, we get 𝑢^𝛿_𝜖(𝑡)≤‖𝑢^𝛿₀‖_𝐿∞(Ω). Similarly multiplying the term 𝑒^−𝑠𝑡(−𝐵−𝑒^−𝑠𝑡𝑢^𝛿_𝜖)⁺ with the approximation PDE and integrating, we obtain 𝑢^𝛿_𝜖(𝑡)≥−‖𝑢^𝛿₀‖_𝐿∞(Ω).

Lemma 3.6. The weak solution of the approximation PDE 𝑢^𝛿_𝜖 satisfied the following inequality, for 𝜏∈[0,𝑇]: ^𝜏_0Ω(𝜕𝑢^𝛿_𝜖)𝜕𝑡²𝜖𝑑𝑥𝑑𝑡+2_Ω|∇𝑢^𝛿_𝜖(𝑡)|²𝑑𝑥+_Ω𝜙(𝑥,∇𝑢^𝛿_𝜖≤𝜖(𝑡))𝑑𝑥2_Ω|∇𝑢^𝛿₀|²𝑑𝑥+_Ω𝜙(𝑥,∇𝑢^𝛿₀)𝑑𝑥+𝜖,a.ein𝜏. (3.16) Moreover, the weak solutions 𝑢^𝛿_𝜖 and 𝜕𝑢^𝛿_𝜖/𝜕𝑡 are uniformly bounded in ϵ with respect to the 𝐿^∞(Ω×[0,𝑇]) and 𝐿²(Ω×[0,𝑇]) norms, respectively.

Proof. The proof of the inequality follows from the identity ^𝜏_0Ω(𝜕𝑢^𝛿_𝜖)𝜕𝑡²𝜖𝑑𝑥𝑑𝑡+2_Ω|∇𝑢^𝛿_𝜖(𝑡)|²𝑑𝑥+_Ω𝜙_𝜖(𝑥,∇𝑢^𝛿_𝜖𝜆(𝑡))𝑑𝑥+2_Ω(𝑢^𝛿_𝜖−𝑢^𝛿₀)²=𝜀𝑑𝑥2_Ω|∇𝑢^𝛿₀|²𝑑𝑥+_Ω𝜙_𝜖(𝑥,∇𝑢^𝛿₀)𝑑𝑥 (3.17)and 𝜙(𝑥,𝑠)≤𝜙_𝜖(𝑥,𝑠)≤𝜙(𝑥,𝑠)+𝜖 for any 𝑠.
From Lemma 3.5, it follows that 𝑢^𝛿_𝜖 is uniformly bounded and ^𝑇_0Ω(𝜕𝑢^𝛿_𝜖)𝜕𝑡²𝑑𝑥𝑑𝑡+_Ω𝜙(𝑥,∇𝑢^𝛿_𝜖𝜖)𝑑𝑥≤2_Ω|∇𝑢^𝛿₀|²𝑑𝑥+_Ω𝜙(𝑥,∇𝑢^𝛿₀)𝑑𝑥+𝜖. (3.18) We have ^𝑇_0Ω(|𝜕𝑢^𝛿_𝜖𝜕𝑡|+𝜙(∇𝑢^𝛿_𝜖))𝑑𝑥𝑑𝑡≤𝐶‖𝑢^𝛿₀‖_𝐿∞(Ω). (3.19)

Let {𝜖_𝑛} be a sequence converging to 0. Consider the corresponding weak solutions of the approximation PDE {𝑢^𝛿_𝜖𝑛}. In the following lemma, we prove the convergence of a subsequence of weak solutions {𝑢^𝛿_𝜖𝑛𝑘} as 𝜖_𝑛𝑘→0.

Lemma 3.7. There exists a subsequence {𝑢^𝛿_𝜖𝑛𝑘} such that as 𝜖_𝑛𝑘→0 dddd 𝑢^𝛿_𝜖𝑛𝑘→𝑢^𝛿stronglyin𝐿¹(Ω×[0,𝑇]),a.e.inΩ×[0,𝑇],𝜕𝑢^𝛿_𝜖𝑛𝑘⇀𝜕𝑡𝜕𝑢^𝛿𝜕𝑡weaklyin𝐿²(Ω×[0,𝑇]) (3.20) where 𝑢^𝛿 is a weak solution of the PDE (2.1) with 𝑢^𝛿(0)=𝑢^𝛿₀. Moreover ‖𝑢^𝛿‖_𝐿^∞_{(Ω×[0,𝑇])}≤𝐶‖𝑢₀‖_𝐿∞(Ω).

Proof. By Lemma 3.6, 𝑢^𝛿_𝜖 and 𝜕𝑢^𝛿_𝜖/𝜕𝑡 are uniformly bounded. When 𝜖_𝑛𝑘→0, for fixed 𝛿>0, we have a convergent subsequence 𝑢^𝛿_𝜖𝑛𝑘→𝑢^𝛿 in 𝐿¹(Ω×[0,𝑇]) for 𝑢^𝛿∈𝐿^∞(BV(Ω)∩𝐿^∞(Ω);[0,𝑇]), and (𝜕𝑢^𝛿_𝜖𝑛𝑘/𝜕𝑡)⇀(𝜕𝑢^𝛿/𝜕𝑡) in 𝐿²(Ω×[0,𝑇]).
Since 𝑢^𝛿_𝜖𝑛𝑘 is the weak solution of the approximation PDE, it satisfies the corresponding weak solution formulation in Definition (3.4). That is, for every 𝑣∈𝐿²(BV(Ω);[0,𝑇]) and 𝜏∈[0,𝑇], ^𝜏_0Ω𝜙_𝜖𝑛𝑘(𝑥,∇𝑢^𝛿_𝜖𝑛𝑘)𝑑𝑥𝑑𝑡−𝜆^𝜏_0Ω(𝑢^𝛿_𝜖𝑛𝑘−𝑢^𝛿₀)(𝑣−𝑢^𝛿_𝜖𝑛𝑘≤𝜖)𝑑𝑥𝑑𝑡_𝑛𝑘2^𝜏_0Ω|∇𝑣|²𝑑𝑥𝑑𝑡+^𝜏_0Ω𝜙_𝜖𝑛𝑘+(𝑥,∇𝑣)𝑑𝑥𝑑𝑡_𝜏0Ω𝜕𝑢^𝛿_𝜖𝑛𝑘𝜕𝑡(𝑣−𝑢^𝛿_𝜖𝑛𝑘)𝑑𝑥𝑑𝑡a.e.in𝜏. (3.21) When 𝜖_𝑛𝑘→0, we obtain, using Lemma 3.3, ^𝜏_0Ω𝜙(𝑥,∇𝑢^𝛿)𝑑𝑥𝑑𝑡−𝜆^𝜏_0Ω(𝑢^𝛿−𝑢^𝛿₀)(𝑣−𝑢^𝛿≤)𝑑𝑥𝑑𝑡_𝜏0Ω𝜙(𝑥,∇𝑣)𝑑𝑥𝑑𝑡+^𝜏_0Ω𝜕𝑢^𝛿𝜕𝑡(𝑣−𝑢^𝛿)𝑑𝑥𝑑𝑡. (3.22) Thus, 𝑢^𝛿 is a weak solution of the PDE (2.1) with 𝑢^𝛿(𝑥,0)=𝑢^𝛿₀(𝑥).

Now, we state and prove our main result.

Theorem 3.8. If 𝑢₀∈BV(Ω)∩𝐿^∞(Ω), then there exists a unique weak solution 𝑢∈𝐿^∞(BV(Ω)∩𝐿^∞(Ω);[0,𝑇]) of (2.1).

Proof. From Lemma 3.7 we get 𝑢^𝛿 as a weak solution of the PDE (2.1) with 𝑢^𝛿(𝑥,0)=𝑢^𝛿₀(𝑥). Using Lemma 3.7 for the sequence {𝑢^𝛿} as 𝛿→0, we get 𝑢^𝛿→𝑢stronglyin𝐿¹(Ω×[0,𝑇]),a.e.inΩ×[0,𝑇],𝜕𝑢^𝛿⇀𝜕𝑡𝜕𝑢𝜕𝑡weaklyin𝐿²(Ω×[0,𝑇]). (3.23) Finally we let 𝛿→0 in the following inequality: ^𝜏_0Ω𝜙(𝑥,∇𝑢^𝛿)𝑑𝑥𝑑𝑡−𝜆^𝜏_0Ω(𝑢^𝛿−𝑢^𝛿₀)(𝑣−𝑢^𝛿≤)𝑑𝑥𝑑𝑡_𝜏0Ω𝜙(𝑥,∇𝑣)𝑑𝑥𝑑𝑡+^𝜏_0Ω𝜕𝑢^𝛿𝜕𝑡(𝑣−𝑢^𝛿)𝑑𝑥𝑑𝑡 (3.24) to obtain ^𝜏_0Ω𝜙(𝑥,∇𝑢)𝑑𝑥𝑑𝑡−𝜆^𝜏_0Ω(𝑢−𝑢₀≤)(𝑣−𝑢)𝑑𝑥𝑑𝑡_𝜏0Ω𝜙(𝑥,∇𝑣)𝑑𝑥𝑑𝑡+^𝜏_0Ω𝜕𝑢𝜕𝑡(𝑣−𝑢)𝑑𝑥𝑑𝑡. (3.25) We see that 𝑢∈𝐿^∞(BV(Ω)∩𝐿^∞(Ω);[0,𝑇]) is a weak solution of (2.1). Uniqueness of the solution follows from the weak solution inequality (3.4) above.

As a consequence of Theorem 3.8, we also obtain a bound for the solution ‖𝑢‖_𝐿^∞≤𝐶‖𝑢₀‖_𝐿∞, 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 𝛼∈𝐶^∞(Ω) (since 𝐺_𝜌∈𝐶^∞(Ω)) 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 ℎ as the grid size, and 𝑈^𝑡_𝑖𝑗 as the intensity value 𝑢(𝑖,𝑗) at iteration 𝑡. 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 𝑈^𝑡+1=[1−𝜏𝐴(𝑈^𝑡)]⁻¹𝑈^𝑡, (4.1) where 𝜏 is the time step, 𝐴(𝑈^𝑡)=[𝑎_𝑖𝑗(𝑈^𝑡)], and𝑎_𝑖𝑗(𝑈^𝑡𝛾)∶=^𝑡𝑖+𝛾^𝑡_𝑗2ℎ²,𝑗∈𝒩_𝑖,−_{𝑘∈𝒩𝑖}𝛾^𝑡_𝑖+𝛾^𝑡_𝑘2ℎ²,𝑗=𝑖,0,otherwise (4.2) with 𝛾_𝑖=𝛼_𝑖𝑔_𝑖 and 𝒩_𝑖 is the discrete neighborhood of the pixel 𝑖. For n-D images the semi-implicit scheme is written as𝑈^𝑡+1=[1−𝜏^𝑛_𝑙=1𝐴_𝑙(𝑈^𝑡)]⁻¹𝑈^𝑡. (4.3) The matrix 𝐴_𝑙=(𝑎^𝑙_𝑖𝑗)_𝑖𝑗 corresponds to the derivatives along the 𝑙th 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 𝐿𝑒𝑛𝑎 image given in Figure 1(b). We show results of applying our scheme (2.1) with regularized TV function 𝑔(𝑠)=(𝜖²+𝑠²)^−1/2, 𝜖=10⁻⁶ 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. Figure 5(a) shows the effect of the same against different noise levels (5≤𝜎≤25) for our scheme. For an 𝑚×𝑛 image 𝑢, the PSNR value of an estimated image 𝑢 is given by PSNR=10log₁₀(𝑚×𝑛_𝑥∑_∈Ω(𝑢(𝑥)−𝑢(𝑥))²)(𝑑𝐵). (4.4) 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 𝐿𝑒𝑛𝑎 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)

(b)

Figure 5 

(a) Noise level 𝜎 versus PSNR (dB) values. (b) Profile line taken across 𝐿𝑒𝑛𝑎 at pixel position 250.

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.

Acknowledgment

The authors thank the anonymous reviewers for their comments which improved the content and presentation of the paper and the handling editor Professor Malgorzata Peszynska for her patience during the review period.