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.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Recommended by Malgorzata Peszynska

Department of Mathematics, Indian Institute of Technology Madras, Chennai 600036, India

Received 21 August 2009; Revised 13 January 2010; Accepted 11 February 2010

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 : [figure omitted; refer to PDF] 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): [figure omitted; refer to PDF] with the initial condition u(x,0)=_u0 (x) . The diffusion function g:^...+ [arrow right]^...+ is chosen to be decreasing, such that g(0)=1, _{lim s[arrow right]∞} 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 [figure omitted; refer to PDF] 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 )-1} exp (-^|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 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 g 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) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

2.2. Edge Adaptive ADS

Instead of using only the edge indicator function g alone to drive the diffusion process, we introduce a spatially adaptive term α in the anisotropic scheme (α -ADS): [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] 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 image u 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 Step 1 ,

(5) zero-crossings of Step 3 displayed if threshold of Step 4 is achieved.

Canny edge map C(x,t) is a binary output with edge pixels marked as 1 and nonedge pixels with 0 . Hence we use [figure omitted; refer to PDF] 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 output C(x) . (b) Smooth edge indicator function E(x) with ρ=2 creates a band around C(x) (c) Diffusion weight αg(|∇u|) .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

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: [figure omitted; refer to PDF] with ^{[varphi][variant prime]} (s)=g(s)s . For example, for the diffusion functions in (1.3) we get the regularization functions as [figure omitted; refer to PDF] 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 [varphi]:...[arrow right][0,∞) satisfying the following.

(H1) [varphi] is a nondecreasing, convex, and even function with [varphi](0)=0 .

(H2) There exist a>0 and b≥0 such that as-b≤[varphi](s)≤as+b , for all s .

An important example in this class is the minimal surface function [varphi](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 [varphi](x,∇u)=α(x)[varphi](∇u) . We recall the definition of total variation of a function u∈^L1 (Ω) .

(a) Diffusion functions _g1 , _g2 given in (1.3) and TV function. (b) Corresponding regularization functions.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

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]

[figure omitted; refer to PDF]

(b) CL [4]

[figure omitted; refer to PDF]

(c) RO [6]

[figure omitted; refer to PDF]

(d) ST [5]

[figure omitted; refer to PDF]

(e) CD [8]

[figure omitted; refer to PDF]

(f) KH [9]

[figure omitted; refer to PDF]

(g) YW [10]

[figure omitted; refer to PDF]

(h) BB [11]

[figure omitted; refer to PDF]

(i) Our scheme

[figure omitted; refer to PDF]

Definition 3.1.

Let u∈^L1 (Ω) ; its total variation is defined as [figure omitted; refer to PDF]

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] [figure omitted; refer to PDF] Moreover, it satisfies ∂u/∂t∈^L2 (Ω×[0,T]) and u(x,0)=_u0 (x) .

The function space ^L2 (BV(Ω);[0,T]) is the set of all functions w:Ω×[0,T][arrow right]... 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): [figure omitted; refer to PDF]

Lemma 3.3.

Let the regularization function [varphi] satisfy the assumptions H1-H2 . If {_ui } converges to u in ^L1 (Ω) , then [figure omitted; refer to PDF]

Proof.

Let w∈^C01 (Ω) be such that |w(x)|≤α(x) for all x∈Ω . Then [figure omitted; refer to PDF] 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.

Let H be a Hilbert space, and let P(H) be its power set. Let ...9C;:H[arrow right]P(H) be a maximal monotone operator and _u0 ∈D(...9C;) , where the domain is D(...9C;)={x∈H:...9C;x≠∅} . Then there exists a unique function u(t):[0,∞)[arrow right]H such that [figure omitted; refer to PDF]

Lemma 3.5.

Consider the approximation PDE [figure omitted; refer to PDF] with [figure omitted; refer to PDF] 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, [figure omitted; refer to PDF]

Proof.

Since _u0 ∈BV(Ω)∩^L∞ (Ω) , existence of the sequence ^u0δ is guaranteed. Let _...9C;... (u)=λ(u-^u0δ )-div (_g... (x,|∇u|)∇u) . As _...9C;... 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 the ^L∞ 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 [figure omitted; refer to PDF] This implies that [figure omitted; refer to PDF] since the remaining integrals are all nonnegative. Let [figure omitted; refer to PDF] Then ...(t) is decreasing, nonnegative, and ...(0)=0 . Thus ...(t)=0 for all t . Subsequently, [figure omitted; refer to PDF] Letting s[arrow right]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 PDE ^u...δ satisfied the following inequality, for τ∈[0,T] : [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] and [varphi](x,s)≤_[varphi]... (x,s)≤[varphi](x,s)+... for any s .

From Lemma 3.5, it follows that ^u...δ is uniformly bounded and [figure omitted; refer to PDF] We have [figure omitted; refer to PDF]

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 [arrow right]0 .

Lemma 3.7.

There exists a subsequence {^{u_...nk δ} } such that as _...nk [arrow right]0 dddd [figure omitted; refer to PDF] 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 Lemma 3.6, ^u...δ and ∂^u...δ /∂t are uniformly bounded. When _...nk [arrow right]0 , for fixed δ>0 , we have a convergent subsequence ^{u_...nk δ} [arrow right]^uδ in ^L1 (Ω×[0,T]) for ^uδ ∈^L∞ (BV(Ω)∩^L∞ (Ω);[0,T]) , and (∂^{u_...nk δ} /∂t)...(∂^uδ /∂t) in ^L2 (Ω×[0,T]) .

Since ^{u_...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] , [figure omitted; refer to PDF] When _...nk [arrow right]0 , we obtain, using Lemma 3.3, [figure omitted; refer to PDF] 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.

If _u0 ∈BV(Ω)∩^L∞ (Ω) , then there exists a unique weak solution u∈^L∞ (BV(Ω)∩^L∞ (Ω);[0,T]) of (2.1).

Proof.

From Lemma 3.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 δ[arrow right]0 , we get [figure omitted; refer to PDF] Finally we let δ[arrow right]0 in the following inequality: [figure omitted; refer to PDF] to obtain [figure omitted; refer to PDF] 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 Theorem 3.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 [varphi] 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 [figure omitted; refer to PDF] where τ is the time step, A(^Ut )=[_aij (^Ut )] , and [figure omitted; refer to PDF] with _γi =_αigi and _...A9;i is the discrete neighborhood of the pixel i . For n-D images the semi-implicit scheme is written as [figure omitted; refer to PDF] The matrix _Al =(^aijl_)ij corresponds to the derivatives along the l 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 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 , ...=1^0-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. Figure 5(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 [figure omitted; refer to PDF] 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) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

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.

References

[1] G. Aubert, P. Kornprobst Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations , vol. 147, of Applied Mathematical Sciences, pp. xxxii+377, Springer, New York, NY, USA, 2006., 2nd.

[2] P. Perona, J. Malik, "Scale-space and edge detection using anisotropic diffusion," IEEE Transactions on Pattern Analysis and Machine Intelligence , vol. 12, no. 7, pp. 629-639, 1990.

[3] Y.-L. You, W. Xu, A. Tannenbaum, M. Kaveh, "Behavioral analysis of anisotropic diffusion in image processing," IEEE Transactions on Image Processing , vol. 5, no. 11, pp. 1539-1553, 1996., [email protected]

[4] F. Catté, P.-L. Lions, J.-M. Morel, T. Coll, "Image selective smoothing and edge detection by nonlinear diffusion," SIAM Journal on Numerical Analysis , vol. 29, no. 1, pp. 182-193, 1992.

[5] D. Strong ftp://ftp.math.ucla.edu/pub/camreport/cam97-38.ps.gz Adaptive total variation minimizing image restoration , Ph.D. thesis, UCLA Mathematics Department, Los Angeles, Calif, USA, 1997.

[6] L. I. Rudin, S. Osher, E. Fatemi, "Nonlinear total variation based noise removal algorithms," Physica D , vol. 60, no. 1-4, pp. 259-268, 1992.

[7] A. Chambolle, P.-L. Lions, "Image recovery via total variation minimization and related problems," Numerische Mathematik , vol. 76, no. 2, pp. 167-188, 1997.

[8] M. Ceccarelli, V. De Simone, A. Murli, "Well-posed anisotropic diffusion for image denoising," IEE Proceedings: Vision, Image and Signal Processing , vol. 149, no. 4, pp. 244-252, 2002.

[9] W. Kusnezow, W. Horn, R. P. Würtz, "Fast image processing with constraints by solving linear PDEs," Electronic Letters in Computer Vision and Image Analysis , vol. 6, no. 2, pp. 22-35, 2007.

[10] J. Yu, Y. Wang, Y. Shen, "Noise reduction and edge detection via kernel anisotropic diffusion," Pattern Recognition Letters , vol. 29, no. 10, pp. 1496-1503, 2008., [email protected]

[11] T. Barbu, V. Barbu, V. Biga, D. Coca, "A PDE variational approach to image denoising and restoration," Nonlinear Analysis: Real World Applications , vol. 10, no. 3, pp. 1351-1361, 2009.

[12] A. Douiri, M. Schweiger, J. Riley, S. R. Arridge, "Anisotropic diffusion regularization methods for diffuse optical tomography using edge prior information," Measurement Science and Technology , vol. 18, no. 1, pp. 87-95, 2007., [email protected]

[13] M. Basu, "Gaussian-based edge-detection methods--a survey," IEEE Transactions on Systems, Man and Cybernetics Part C , vol. 32, no. 3, pp. 252-260, 2002., [email protected]

[14] R. C. Gonzelez, R. E. Woods Digital Image Processing , Pearson, Upper Saddle River, NJ, USA, 2002.

[15] J. F. Canny, "A computational approach to edge detection," IEEE Transactions on Pattern Analysis and Machine Intelligence , vol. 8, no. 6, pp. 679-698, 1986.

[16] E. Giusti Minimal Surfaces and Functions of Bounded Variation , vol. 80, of Monographs in Mathematics, pp. xii+240, Birkhäuser, Basel, Switzerland, 1984.

[17] V. B. S. Prasath, A. Singh, "Edge detectors based anisotropic diffusion for enhancement of digital images," in Proceedings of the 6th Indian Conference on Computer Vision, Graphics and Image Processing (ICVGIP '08), pp. 33-38, IEEE Computer Society Press, Bhubaneswar, India, [email protected]; [email protected], December 2008.

[18] F. Demengel, R. Temam, "Convex functions of a measure and applications," Indiana University Mathematics Journal , vol. 33, no. 5, pp. 673-709, 1984.

[19] L. C. Evans, R. F. Gariepy Measure Theory and Fine Properties of Functions , of Studies in Advanced Mathematics, pp. viii+268, CRC Press, Boca Raton, Fla, USA, 1992.

[20] H. Brézis Opérateurs Maximaux Monotones et Semi-Groupes de Contractions dans les Espaces de Hilbert , of North-Holland Mathematics Studies, no. 5, pp. vi+183, North-Holland, Amsterdam, The Netherlands, 1973.