On image denoising methods

On image denoising methods Antoni Buades ∗ † Bartomeu Coll ∗ Jean Michel Morel † Abstract The search for efficient image denoising methods still is a valid challenge, at the crossing of functional analysis and statistics. In spite of the sophistication of the recently proposed methods, most algorithms have not yet attained a desirable level of applicability. All show an outstanding performance when the image model corresponds to the algorithm assumptions, but fail in general and create artifacts or remove image fine structures. The main focus of this paper is, first, to define a general mathematical and experimental methodology to compare and classify classical image denoising algorithms, second, to propose an algorithm (Non Local Means) addressing the preservation of structure in a digital image. The mathematical analysis is based on the analysis of the “method noise”, defined as the difference between a digital image and its denoised version. The NL-means algorithm is also proven to be asymptotically optimal under a generic statistical image model. The denoising performance of all considered methods are compared in four ways ; mathematical: asymptotic order of magnitude of the method noise under regularity assumptions; perceptual-mathematical: the algorithms artifacts and their explanation as a violation of the image model; quantitative experimental: by tables of L2 distances of the denoised version to the original image. The most powerful evaluation method seems, however, to be the visualization of the method noise on natural images. The more this method noise looks like a real white noise, the better the method. 1 Introduction 1.1 Digital images and noise The need for efficient image restoration methods has grown with the massive production of digital images and movies of all kinds, often taken in poor conditions. No matter how good cameras are, an image improvement is always desirable to extend their range of action. A digital image is generally encoded as a matrix of grey level or color values. In the case of a movie, this matrix has three dimensions, the third one corresponding to time. Each pair (i, u(i)) where u(i) is the value at i is called pixel, for “picture element”. In the case of grey level images, i is a point on a 2D grid and u(i) is a real value. In the case of classical color images, u(i) is a triplet of values for the red, green and blue components. All of what we shall say applies identically to movies, 3D images and color or multispectral images. For a sake of simplicity in notation and display of experiments, we shall here be contented with rectangular 2D grey-level images. The two main limitations in image accuracy are categorized as blur and noise. Blur is intrinsic to image acquisition systems, as digital images have a finite number of samples and must respect the Shannon-Nyquist sampling conditions [34]. The second main image perturbation is noise. Each one of the pixel values u(i) is the result of a light intensity measurement, usually made by a CCD matrix coupled with a light focusing system. Each captor of the CCD is roughly a square ∗ Universitat † Centre de les Illes Balears, Ctra. Valldemossa Km. 7.5, 07122 Palma de Mallorca, Spain de Mathématiques et Leurs Applications. ENS Cachan 61, Av du Président Wilson 94235 Cachan, France 1 in which the number of incoming photons is being counted for a fixed period corresponding to the obturation time. When the light source is constant, the number of photons received by each pixel fluctuates around its average√in accordance with the central limit theorem. In other terms one can expect fluctuations of order n for n incoming photons. In addition, each captor, if not adequately cooled, receives heat spurious photons. The resulting perturbation is usually called “obscurity noise”. In a first rough approximation one can write v(i) = u(i) + n(i), where v(i) is the observed value, u(i) would be the “true” value at pixel i, namely the one which would be observed by averaging the photon counting on a long period of time, and n(i) is the noise perturbation. As indicated, the amount of noise is signal-dependent, that is n(i) is larger when u(i) is larger. In noise models, the normalized values of n(i) and n(j) at different pixels are assumed to be independent random variables and one talks about “white noise”. 1.2 Signal and noise ratios A good quality photograph (for visual inspection) has about 256 grey level values, where 0 represents black and 255 represents white. Measuring the amount of noise by its standard deviation, σ(n), one can define the signal noise ratio (SNR) as SN R = σ(u) , σ(n) where σ(u) denotes the empirical standard deviation of u(i), ! 12 Ã 1 X 2 (u(i) − u) σ(u) = |I| i P 1 and u = |I| i∈I u(i) is the average grey level value. The standard deviation of the noise can also be obtained as an empirical measurement or formally computed when the noise model and parameters are known. A good quality image has a standard deviation of about 60. The best way to test the effect of noise on a standard digital image is to add a gaussian white noise, in which case n(i) are i.i.d. gaussian real variables. When σ(n) = 3, no visible alteration is usually observed. Thus, a 60 3 ≃ 20 signal to noise ratio is nearly invisible. Surprisingly enough, one can add white noise up to 21 ratio and still see everything in a picture ! This fact is illustrated in Figure 1 and constitutes a major enigma of vision. It justifies the many attempts to define convincing denoising algorithms. As we shall see, the results have been rather deceptive. Denoising algorithms see no difference between small details and noise, and therefore remove them. In many cases, they create new distortions and the researchers are so much used to them as to have created a taxonomy of denoising artifacts: “ringing”, “blur”, “staircase effect”, “checkerboard effect”. This fact is not quite a surprise. Indeed, to the best of our knowledge, all denoising algorithms are based on • a noise model • a generic image smoothness model, local or global. In experimental settings, the noise model is perfectly precise. So the weak point of the algorithms is the inadequacy of the image model. All of the methods assume that the noise is oscillatory, and that the image is smooth, or piecewise smooth. So they try to separate the smooth or patchy part (the image) from the oscillatory one. Actually, many fine structures in images are as oscillatory as noise is; conversely, white noise has low frequencies and therefore smooth components. Thus a separation method based on smoothness arguments only is hazardous. 2 Figure 1: A digital image with standard deviation 55, the same with noise added (standard deviation 3), the signal noise ratio being therefore equal to 18, and the same with signal noise ratio slightly larger than 2. In this second image, no alteration is visible. In the third, a conspicuous noise with standard deviation 25 has been added but, surprisingly enough with a two to 1 signal ratio, all details of the original image still are visible. 1.3 The “method noise” All denoising methods depend on a filtering parameter h. This parameter measures the degree of filtering applied to the image. For most methods, the parameter h depends on an estimation of the noise variance σ 2 . One can define the result of a denoising method Dh as a decomposition of any image v as v = Dh v + n(Dh , v), (1) where 1. Dh v is more smooth than v 2. n(Dh , v) is the noise guessed by the method. Now, it is not enough to smooth v to ensure that n(Dh , v) will look like a noise. The more recent methods are actually not contented with a smoothing, but try to recover lost information in n(Dh , v) [21], [27]. So the focus is on n(Dh , v). Following the same idea but applying the denoising method to the original image, we define the method noise. Definition 1.3.1 (Method noise) Let u be an image and Dh a denoising operator depending on h. Then we define the method noise of u as the image difference n(Dh , u) = u − Dh (u). (2) This method noise should be as similar to a white noise as possible. In addition, since we would like the original image u not to be altered by denoising methods, the method noise should be as small as possible for the functions with the right regularity. According to the preceding discussion, four criteria can and will be taken into account in the comparison of denoising methods: • a display of typical artifacts in denoised images. • a formal computation of the method noise on smooth images, evaluating how small it is in accordance with image local smoothness. 3 • a comparative display of the method noise of each method on real images with σ = 2.5. We mentioned that a noise standard deviation smaller than 3 is subliminal and it is expected that most digitization methods allow themselves this kind of noise. • a classical comparison method based on noise simulation : it consists of taking a good quality image, add gaussian white noise to it with known σ and then compute the best image recovered from the noisy one by each method. A table of L2 distances from the restored to the original can be established. The L2 distance does not provide a good quality assessment. However, it reflects well the relative performances of algorithms. On top of this, in two cases, a proof of asymptotic recovery of the image can be obtained by statistical arguments. 1.4 Which methods to compare ? We had to make a selection of the denoising methods we wished to compare. Here a difficulty arises, as most original methods have caused an abundant literature proposing many improvements. So we tried to get the best available version, but keeping the simple and genuine character of the original method : no hybrid method. So we shall analyze : 1. the Gaussian smoothing model (Gabor [19]), where the smoothness of u is measured by the R Dirichlet integral |Du|2 . 2. the anisotropic filtering model (Perona-Malik [30], Alvarez et al. [1]). 3. the Rudin-Osher-Fatemi [15] total variation model. 4. the Yaroslavsky ([41], [40]) neighborhood filters and an elegant variant, the SUSAN filter (Smith and Brady) [35]. 5. the Wiener local empirical filter as implemented by Yaroslavsky [40]. 6. the translation invariant wavelet thresholding [7], a simple and performing variant of the wavelet thresholding [11]. 7. DUDE, the discrete universal denoiser [26]. 8. The non local means (NL-means) algorithm. This last algorithm may deserve some more precise mention here because of the simple closed formula defining it, that is, Z (Ga ∗|u(x+.)−u(y+.)|2 )(0) 1 h2 e− u(y) dy, N L(u)(x) = C(x) R (Ga ∗|u(x+.)−u(z+.)|2 )(0) h2 dz is the where Ga is a Gaussian kernel of standard deviation a, C(x) = e− normalizing factor and h acts as a filtering parameter. This amounts to say that the denoised value at x is a mean of the values of all pixels whose gaussian neighborhood looks like the neighborhood of x. 4 1.5 What is left We do not draw into comparison the hybrid methods, in particular the total variation + wavelets ([6], [12], [20]). Such methods are significant improvements of the simple methods but are impossible to draw into a benchmark : their efficiency depends a lot upon the choice of wavelet dictionaries and the kind of image. Second, we do not draw into the comparison the method introduced recently by Y. Meyer [24], whose aim it is to decompose the image into a BV part and a texture part (the so called u + v methods), and even into three terms, namely u + v + w where u is the BV part, v is the “texture” ∞ part belonging to the dual space of BV , denoted by G, and w belongs to the Besov space Ḃ−1,∞ , a space characterized by the fact that the wavelet coefficients have a uniform bound. G is proposed by Y. Meyer as the right space to model oscillatory patterns such as textures. The main focus of this method is not denoising, yet. Because of the different and more ambitious scopes of the Meyer method, [2, 38, 29], which makes it parameter and implementation-dependent, we could not draw it into the discussion. Last but not least, let us mention the bandlets [17] and curvelets [36] transforms for image analysis. These methods also are separation methods between the geometric part and the oscillatory part of the image and intend to find an accurate and compressed version of the geometric part. Incidentally, they may be considered as denoising methods in geometric images, as the oscillatory part then contains part of the noise. Those methods are closely related to the total variation method and to the wavelet thresholding and we shall be contented with those simpler representative. 1.6 Plan of the paper Section 2 computes formally the method noise for the best elementary local smoothing methods, namely gaussian smoothing, anisotropic smoothing (mean curvature motion), total variation minimization and the neighborhood filters. For all of them we prove or recall the asymptotic expansion of the filter at smooth points of the image and therefore obtain a formal expression of the method noise. This expression permits to characterize places where the filter performs well and where it fails. In section 3, we treat the Wiener-like methods, which proceed by a soft or hard threshold on frequency or space-frequency coefficients. We examine in turn the Wiener-Fourier filter, then the Yaroslavsky Local adaptive filters, based on DCT and the wavelet thresholding method. Of course the gaussian smoothing belongs to both classes of filters. We finally describe the universal denoiser DUDE, but we cannot draw it into the comparison as its direct application to grey level images is unpractical so far (we discuss its feasibility). In section 5, we introduce the Non Local means (NL-means) filter. This method is not easily classified in the preceding terminology, since it can work adaptively in a local or non local way. We first give a proof that this algorithm is asymptotically consistent (it gives back the conditional expectation of each pixel value given an observed neighborhood) under the assumption that the image is a fairly general stationary random process. The works of Efros and Leung [13] and Levina [18] have shown that this assumption is sound for images having enough samples in each texture patch. In section 6, we compare all algorithms from several points of view, do a performance classification and explain why the NL-means algorithm shares the consistency properties of most of the aforementioned algorithms. 2 Local smoothing filters The original image u is defined in a bounded domain Ω ⊂ R2 , and denoted by u(x) for x = (x, y) ∈ R2 . This continuous image is usually interpreted as the Shannon interpolation of a discrete grid of samples, [34] and is therefore analytic. The distance between two consecutive samples will be denoted by ε. The noise itself is a discrete phenomenon on the sampling grid. According to the usual screen and printing visualization practice, we do not interpolate the noise samples ni as a band limited function, 5 but rather as a piecewise constant function, constant on each pixel i and equal to ni . 1 We write |x| = (x2 + y 2 ) 2 and x1 .x2 = x1 x2 + y1 y2 their scalar product and denote the derivatives ∂2u of u by ux = ∂u uy = ∂u uxy = ∂x∂y . The gradient of u is written Du = (ux , uy ) and the ∂x , ∂y , Laplacian of u ∆u = uxx + uyy . 2.1 Gaussian smoothing By Riesz theorem, image isotropic linear filtering boils down to a convolution of the image by a linear symmetric kernel. The smoothing requirement is usually expressed by the positivity of the kernel. |x|2 − 4h2 1 . In that The paradigm of such kernels is of course the gaussian kernel x → Gh (x) = (4πh 2) e case, Gh has standard deviation h and it is easily seen that Theorem 2.1.1 (Gabor 1960) The image method noise of the convolution with a gaussian kernel Gh is u − Gh ∗ u = −h2 ∆u + o(h2 ). A similar result is actually valid for any positive radial kernel with bounded variance, so one can keep the gaussian example without loss of generality. The preceding estimate is valid if h is small enough. On the other hand, the noise reduction properties depend upon the fact that the neighborhood involved in the smoothing is large enough, so that the noise gets reduced by averaging. So in the following we assume that h = kε, where k stands for the number of samples of the function u and of the noise in an interval of length h, k must be much larger than 1 to assure the noise reduction. The effect of a gaussian smoothing on the noise can be evaluated at a reference pixel i = 0. At this pixel, X XZ Gh (x)n(x)dx = ε2 Gh (i)ni Gh ∗ n(0) = i∈I Pi i∈I where we recall that n(x) is been interpolated as a piecewise function, the Pi square pixels centered in i have size ε2 and Gh (i) denotes the mean value of the function Gh on the pixel i. Denoting by V ar(X) the variance of a random variable X, the additivity of variances of independent centered random variables yields Z X ε2 σ 2 V ar(Gh ∗ n(0)) = ε4 Gh (i)2 σ 2 ≃ σ 2 ε2 Gh (x)2 dx = . 8πh2 i So we have proved Theorem 2.1.2 Let n(x) be the a piecewise constant white noise, with n(x) = ni on each square pixel i. Assume that the ni are i.i.d. with zero mean and variance σ 2 . Then the “noise residue” after a gaussian convolution of n by Gh satisfies V ar(Gh ∗ n(0)) ≃ ε2 σ 2 . 8πh2 In other terms, the standard deviation of the noise, which can be identified with the noise amplitude, is multiplied by h√ε8π . Theorems 2.1.1 and 2.1.2 traduce the delicate equilibrium between noise reduction and image destruction by any linear smoothing. Denoising does not alter the image at points where it is smooth at a scale h much larger than the sampling scale ε. The first theorem tells us, that the method noise of the gaussian denoising method is zero in harmonic parts of the image. A Gaussian convolution is optimal on harmonic functions, and performs instead poorly on singular parts of u, namely edges or texture, where the Laplacian of the image is large. See Figure 3. 6 2.2 Anisotropic filters and curvature motion The anisotropic filter (AF ) attempts to avoid the blurring effect of the gaussian by convolving the image u at x only in the direction orthogonal to Du(x). The idea of such filter goes back to Perona and Malik [30] and actually again to Gabor [19]. Set Z Du(x)⊥ AFh u(x) = Gh (t)u(x + t )dt, |Du(x)| x2 1 e− 2h2 is the one-dimensional for x such that Du(x) 6= 0 and where (x, y)⊥ = (−y, x) and Gh = √2πh Gauss function with variance h2 . If one assumes that the original image u is twice continuously differentiable (C 2 ) at x, it is easily shown by a second order Taylor expansion that Theorem 2.2.1 The image method noise of an anisotropic filter AFh is Du⊥ Du⊥ 1 1 , ) = − h2 |Du|curv(u)(x), u(x) − AFh u(x) ≃ − h2 D2 u( 2 |Du| |Du| 2 where the relation holds when Du(x) 6= 0. By curv(u)(x), we denote the curvature, i.e. the signed inverse of the radius of curvature of the level line passing by x. When Du(x) 6= 0, this means curv(u) = uxx u2y − 2uxy ux uy + uyy u2x 3 (u2x + u2y ) 2 . This method noise is zero wherever u behaves locally like a one variable function, u(x, y) = f (ax + by + c). In such a case, the level line of u is locally the straight line with equation ax + by + c = 0 and the gradient of f may instead be very large. In other terms, with anisotropic filtering, an edge can be maintained. On the other hand, we have to evaluate the gaussian noise reduction. This is easily done by a one-dimensional adaptation of Theorem 2.1.2. Notice that the noise on a grid is√ not isotropic ; so the gaussian average when Du is parallel to one coordinate axis is made roughly on 2 more samples than the gaussian average in the diagonal direction. Theorem 2.2.2 By anisotropic gaussian smoothing, when ε is small enough with respect to h, the noise residue satisfies ε Var (AFh (n)) ≤ √ σ2 2πh In other terms, the standard deviation of the noise n is multiplied by a factor at most equal to ε ( √2πh )1/2 , this maximal value being attained in the diagonals. ⊥ (x ) Proof: Let D be the line x + t Du |Du(x)| passing by x, parameterized by t ∈ R and denote by Pi , i ∈ I the pixels which meet D, n(i) the noise value, constant on pixel Pi , and εi the length of the ⊥ (x ) intersection of D ∩ Pi . Denote by g(i) the average of Gh (x + t Du |Du(x)| ) on D ∩ Pi . We One then has AFh n(x) ≃ X εi n(i)g(i). i The n(i) are i.i.d. with standard variation σ and therefore X X V ar(AFh (n)) = ε2i σ 2 g(i)2 ≤ σ 2 max(εi ) εi g(i)2 . i i 7 This yields Var (AFh (n)) ≤ √ 2εσ 2 Z Gh (t)2 dt = √ ε σ2 . 2πh ¤ There are many versions of AFh , all yielding an asymptotic estimate equivalent to the one in Theorem 2.2.1 : the famous median filter [15], an inf-sup filter on segments centered at x [9], and the clever numerical implementation of the mean curvature equation in [3]. So all of those filters have in common the good preservation of edges, but perform poorly on flat regions and are worse there than a gaussian blur. This is patent in Figure 3. 2.3 Total variation The Total variation minimization was introduced by Rudin, Osher and Fatemi [32, 33]. The original image u is supposed to derive from a simple geometric description, a set of connected sets, the objects, along with their smooth contours, or edges. The image is smooth inside the objects but with jumps across the boundaries. The functional space modelling these properties is BV (Ω), the space of R integrable functions with finite total variation T VΩ (u) = |Du|, where Du is assumed to be a Radon measure. Given a noisy image v(x), the above mentioned authors proposed to recover the original image u(x) as the solution of the constrained minimization problem arg min T VΩ (u) u (3) subject to the noise constraints Z Z |u(x) − v(x)|2 dx = σ 2 . (u(x) − v(x))dx = 0 and Ω Ω The solution u must be as regular as possible in the sense of the total variation, while the difference v(x) − u(x) is treated as an error, with a prescribed energy. The constraints prescribe the right mean and variance to u − v, but do not ensure that it be similar to a noise (see a thorough discussion in [24]). The preceding problem is naturally linked to the unconstrained problem Z |v(x) − u(x)|2 dx (4) arg min T VΩ (u) + λ u Ω for a given Lagrange multiplier λ. The above functional is strictly convex and lower semicontinuous with respect to the weak-star topology of BV. Therefore the minimum exists, is unique and computable (see (e.g.) [4].) The parameter λ controls the trade off between the regularity and fidelity terms. As λ gets smaller the weight of the regularity term increases. Therefore λ is related to the degree of filtering of the solution of the minimization problem. Let us denote by TVFλ (v) the solution of problem (4) for a given value of λ. The Euler Lagrange equation associated with the minimization problem is given by 1 curv(u)(x) = 0 (u(x) − v(x)) − 2λ (see [32]). Thus, Theorem 2.3.1 The image method noise of the total variation minimization (4) is 1 curv(TVFλ (u))(x). 2λ As in the anisotropic case, straight edges are maintained because of their small curvature. However, details and texture can be over smoothed if λ is too small, as is shown in Figure 3. u(x) − TVFλ (u)(x) = − 8 2.4 Iterated Total Variation refinement In the original TV model the removed noise, v(x) − u(x), is treated as an error and is no longer studied. In practice, some structures and texture are present in this error. Many recent works have been introduced trying to avoid this effect. The simplest one, due to Burger et al. [27], consists of an iteration of the original model. The new algorithm is as follows. 1. First solve the original TV model u1 = arg min u∈BV ½Z |∇u(x)|dx + λ Z ¾ (v(x) − u(x))2 dx . 2. Perform a correction step to obtain ½Z ¾ Z 2 |∇u(x)|dx + λ (v(x) + n1 (x) − u(x)) dx , u2 = arg min u∈BV where n1 is the noise estimated by the first step, v(x) = u1 (x) + n1 (x). The correction step adds this first estimate of the noise to the original image. 3. Iterate : compute uk+1 as a minimizer of the modified total variation minimization, ½Z ¾ Z 2 uk+1 = arg min |∇u(x)|dx + λ (v(x) + nk (x) − u(x)) dx , u∈BV where v + nk = uk+1 + nk+1 . Some results are presented in [27] which clarify the nature of the above sequence: • uk converges monotonically in L2 to v, the noisy image, as k → ∞. • uk approaches the noisy free image monotonically in the Bregman distance associated with the BV seminorm, at least until kuk̄ − uk ≤ σ 2 , where u is the original image and σ is the standard deviation of the added noise. These two results indicate how to stop the sequence and choose uk̄ . It is enough to proceed iteratively until the result gets noisier or the distance kuk̄ − uk2 gets smaller than σ 2 . The new solution has more details preserved, as Figure 3 shows. One can attempt an explanation of this improvement by computing the method noise. The above iterated denoising strategy being quite general, one can make the computations for a linear denoising operator T as well. In that case, the method noise after one iteration is u − T (u + n1 ) = = n1 − (T (u + n1 ) − T (u)) n1 − T (n1 ). In the linear case, this amounts to say that the first estimated noise n1 is filtered again and its smooth components added back to the original. 9 2.5 Neighborhood filters The previous filters are based on a notion of spatial neighborhood or proximity. Neighborhood filters instead take into account grey level values to define neighboring pixels. In the simplest and more extreme case, the denoised value at pixel i is an average of values at pixels which have a grey level value close to u(i). The grey level neighborhood is therefore B(i, h) = {j ∈ I | u(i) − h < u(j) < u(i) + h} . This is a fully non-local algorithm, since pixels belonging to the whole image are used for the estimation at pixel i. This algorithm can be written in a more continuous form Z |u(y)−u(x)|2 1 h2 NFh u(x) = dy, u(y)e− C(x) Ω |u(y)−u(x)|2 R h2 dy is the normalization factor. where Ω ⊂ R2 is an open and bounded set, and C(x) = Ω e− The first question to address here is the consistency of such a filter, namely how close the denoised version is to the original when u is smooth. Lemma 2.5.1 Suppose u is Lipschitz in Ω and h > 0, then C(x) ≥ O(h2 ). Given x, y ∈ Ω, by the Mean Value Theorem, |u(x) − u(y)| ≤ K|x − y| for some real |u(y)−u(x)|2 |u(y)−u(x)|2 R R 2 h2 h2 constant K. Then, C(x) = Ω e− dy ≥ B(x,h) e− dy ≥ e−K O(h2 ). ¤ Proof: Proposition 2.5.2 (Method noise estimate). Suppose u is a Lipschitz bounded √ function on Ω, where Ω is an open and bounded domain of R2 . Then |u(x) − NFh u(x)| = O(h − log h), for h small, 0 < h < 1, x ∈ Ω. Proof: Let x be a point of Ω and for a given B and h, B, h ∈ R, consider the set Dh = {y ∈ Ω | |u(y) − u(x)| ≤ Bh}. Then |u(x) − NFh u(x)| ≤ 1 C Z Dh e− |u( y)−u(x)|2 h2 |u(y) − u(x)|dy + 1 C Z e− y)−u(x)|2 |u( |u(y) − u(x)|dy. h2 c Dh By one hand, considering that |u(y)−u(x)| ≤ Bh for y ∈ Dh and R e− y)−u(x)|2 |u( dy ≤ C(x) one y )−u(x)|2 2 h2 ≤ e−B sees that the first Rterm is bounded by Bh. On the other hand, considering that e− for y ∈ / Dh , that Dc |u(y) − u(x)|dy is bounded, and C ≥ O(h2 ) one deduces that the second term h Dh h2 |u( 2 has an order O(h−2 e−B ). Finally, choosing B such that B 2 = −3 log h yields p 2 |u(x) − NFh u(x)| ≤ Bh + O(h−2 e−B ) = O(h − log h) + O(h) √ and so the method noise has order O(h − log h). ¤ The Yaroslavsky ([40, 41]) neighborhood filters consider mixed neighborhoods B(i, h) ∩ Bρ (i), where Bρ (i) is a ball of center i and radius ρ. So the method takes an average of the values of pixels which are both close in grey level and spatial distance. The SUSAN filter [35] proposes the following elegant close formula for such filters, Z |u(y)−u(x)|2 1 h2 dy, u(y)e− YNFh,ρ (x) = C(x) Bρ (x) 10 |u(y)−u(x)|2 R h2 dy is the normalization factor. where C(x) = Bρ (x) e− First, we study the method noise in the 1D case. In that case, u denotes a one dimensional signal. Theorem 2.5.3 Suppose u ∈ C 2 ((a, b)), a, b ∈ R. Then, for h, ρ ∈ R+ and h small enough u(s) − YNFρ,h (s) ≃ − where h2 ′′ h u (s) f ( |u′ (s)|), 2 ρ 1 3 − 35 t2 . 1 − 13 t2 f (t) = Proof: Let s ∈ (a, b) and h, ρ ∈ R+ . Then u(s) − YNFρ,h (s) = − R h −h 1 (u(s+t)−u(s))2 − ρ2 e dt Z h −h − (u(s + t) − u(s))e (u(s+t)−u(s))2 ρ2 dt. 2 We take the Taylor expansion of u(s+t) and the exponential function e−x , then we integrate obtaining that 5 ′2 ′′ h3 u′′ u u − 3h 5ρ 2 3 , u(s) − YNFρ,h (s) ≃ − 3 ′2 2h − 2h3ρu2 for h small enough. The method noise follows from the above expression. ¤ The previous result shows that the neighborhood filtering method noise can be locally written as a weighted Laplacian. The weighting function f gives the sign and magnitude of the method noise. The function f takes positive values, where the method noise behaves as a Laplacian, and negative values, where the method noise behaves as an inverse Laplacian. The zeros and the discontinuity points of f represent the singular points where the behavior of the method changes. The magnitude of this change is much larger near the discontinuities of f producing a shock or staircase effect. Figure 2 illustrates this fact and the asymptotical behavior of the filter. Next we give the analogous result for 2D images. Theorem 2.5.4 Suppose u ∈ C 2 (Ω), Ω ⊂ R2 . Then, for h, ρ ∈ R+ and h small enough µ ¶ h2 h Du Du h 2 u(x) − YNFρ,h (x) ≃ − ∆u f ( |Du|) + D u( , ) g( |Du|) 8 ρ |Du| |Du| ρ where f (t) = 1 − 21 t2 , 1 − 41 t2 g(t) = −t2 . 1 − 41 t2 Z (u(x + t) − u(x))e Proof: Let x ∈ Ω and h, ρ ∈ R+ . Then u(x) − YNFρ,h (x) = − 1 (u( x +t)−u(x))2 R − ρ2 dt e Bh (0) − Bh (0) (u( x+t)−u(t))2 ρ2 dt. 2 We take the Taylor expansion of u(x + t), and the exponential function e−y . Then, we take polar coordinates and integrate, obtaining u(x)−YNFρ,h (x) ≃ πh4 8 ∆u − πh6 2 16ρ2 (ux uxx + u2y uxx + u2x uyy + u2x uxx ) − πh2 − 11 h4 π 2 4ρ2 (ux πh6 2 8ρ2 (ux uxx + u2y ) + 2ux uy uxy + u2y uyy ) , Figure 2: One dimensional neighborhood filter experience. Top: Original signal, u(s) = sin 4πs for ′ | s ∈ (0, 1). Middle left: Weighting function f ( h|u ρ ) for the sinus signal and h = 0.0038, ρ = 0.2. Middle center: signal after 5 iterations with the same parameters. The discontinuity points of f produces a shock effect on the filtered signal. Middle right: Asymptotically filtered signal. Bottom left: Weighting function for h = 0.0038, ρ = 4.0. Bottom center: signal after 5 iterations. The signal is filtered by a gaussian smoothing, since f is positive, and therefore the method noise behaves like a Laplacian. Bottom right: Asymptotically filtered signal. for h small enough. Grouping the terms of above expression we get the desired result. ¤ In order to make clear the behavior of the previous method noise, we comment the different cases depending on the value of t = hρ |Du|. 2 • If t is near zero the method noise is equal to − h8 ∆u, the Laplacian operator, and therefore the neighborhood filter behaves like a Gaussian smoothing. √ 2 Du Du , |Du| ). Therefore, the structures with • If t is near 2 the method noise is equal to − h8 D2 u( |Du| Du Du 2 a large value of D u( |Du| , |Du| ) will not be preserved. According to Haralick, edge points are Du Du , |Du| ) is zero ; on both sides of the edges it can be instead large and the points where D2 u( |Du| method actually enhances the edges by making the image flat on both sides. This is a typical shock filter effect [28]. √ 2 • If t ∈ (0, 2) and |f (t)| ≃ |g(t)| the method noise equals to − h8 |f (t)||Du|curvu(x), the curvature operator, and therefore the filter behaves like an anisotropic filter. ´ ³ 2 Du Du , |Du| ) . • When t is much larger than 2 the method noise equals to − h8 ∆u + 2D2 u( |Du| The method is very instable near t = 2 where the functions f (t) and g(t) have an asymptotical discontinuity. This instability can deal to unwanted shock effects and artifacts, as Figure 3 shows. 12 Figure 3: Denoising experience on a natural image. From left to right and from top to bottom: noisy image (standard deviation 20), gaussian convolution, anisotropic filter, total variation minimization, total variation iteration and Yaroslavsky neighborhood filter. 3 Frequency domain filters Let u be the original image defined on the grid I . The image is supposed to be modified by the addition of a signal independent white noise N . N is a random process where N (i) are i.i.d, zero mean and have constant variance σ 2 . The resulting noisy process depends on the random noise component, and therefore is modeled as a random field V , V (i) = u(i) + N (i). (5) Given a noise observation n(i), v(i) denotes the observed noisy image, v(i) = u(i) + n(i). (6) Let B = {gα }α∈A be an orthogonal basis of R|I| . The noisy process is transformed as VB (α) = uB (α) + NB (α), (7) where VB (α) = hV, gα i, uB (α) = hu, gα i, NB (α) = hN, gα i are the scalar products of V, u and N with gα ∈ B. Noise coefficients NB (α) remain uncorrelated and zero mean, but the variances are multiplied by kgα k2 , E[NB (α)NB (β)] = X gα (m)gβ (n)E[N (m)N (n)] m,n∈I = hgα , gβ iσ 2 = σ 2 kgα k2 δ[α − β]. 13 Frequency domain filters are applied independently to every transform coefficient VB (α) and then the solution is estimated by the inverse transform of the new coefficients. Noisy coefficients VB (α) are modified to a(α)VB (α). This is a non linear algorithm because a(α) depends on the value VB (α). The inverse transform yields the estimate X Û = DV = a(α) VB (α) gα . (8) α∈A D is also called a diagonal operator. Let us look for the frequency domain filter D which minimizes a certain estimation error. This error is based on the square euclidean distance, and it is averaged over the noise distribution. Definition 3.0.5 Let u be the original image, N a white noise and V = u + N . Let D be a frequency domain filter. Define the risk of D as r(D, u) = E{||u − DV ||2 } (9) where the expectation is taken over the noise distribution. The following theorem, which is easily proved, gives the diagonal operator Dinf that minimizes the risk, Dinf = arg min r(D, u). D Theorem 3.0.6 The operator Dinf which minimizes the risk is given by the family {a(α)}α a(α) = and the corresponding risk is rinf (u) = X s∈S |uB (α)|2 , |uB (α)|2 + kgα k2 σ 2 kgα k4 |uB (α)|2 σ 2 . |uB (α)|2 + kgα k2 σ 2 (10) (11) The previous optimal operator attenuates all noisy coefficients in order to minimize the risk. If one restricts a(α) to be 0 or 1, one gets a projection operator. In that case, a subset of coefficients is kept and the rest gets cancelled. The projection operator that minimizes the risk r(D, u) is obtained by the following family of {a(α)}α , ½ 1 |uB (α)|2 ≥ ||gα ||2 σ 2 a(α) = 0 |uB (α)|2 < ||gα ||2 σ 2 and the corresponding risk is rp (u) = X ||gα ||2 min(|uB (α)|2 , ||gα ||2 σ 2 ). Note that both filters are ideal operators because they depend on the coefficients uB (α) of the original image, which are not known. We call, as classical, Fourier Wiener Filter the optimal operator (10) where B is a Fourier Basis. This is an ideal filter, since it uses the (unknown) Fourier transform of the original image. By the use of the Fourier basis (see Figure 4), global image characteristics may prevail over local ones and create spurious periodic patterns. To avoid this effect, the basis must take into account more local features, as the wavelet and local DCT transforms do. The search for the ideal basis associated with each image still is open. At the moment, the way seems to be a dictionary of basis instead of one single basis, [21]. 14 Figure 4: Fourier Wiener filter experiment. Top Left: Degraded image by an additive white noise of σ = 15. Top Right: Fourier Wiener filter solution. Down: Zoom on three different zones of the solution. The image is filtered as a whole and therefore a uniform texture is spread all over the image. 3.1 Local adaptive filters in Transform Domain The local adaptive filters have been introduced by L. Yaroslavsky [41, 42]. In this case, the noisy image is analyzed in a moving window and in each position of the window its spectrum is computed and modified. Finally, an inverse transform is used to estimate only the signal value in the central pixel of the window. Let i ∈ I be a pixel and W = W (i) a window centered in i. Then the DCT transform of W is computed and modified. The original image coefficients of W, uB,W (α) are estimated and the optimal attenuation of Theorem 3.0.6 is applied. Finally, only the center pixel of the restored window is used. This method is called the Empirical Wiener Filter. In order to approximate uB,W (α) one can take averages on the additive noise model, that is, E|VB,W (α)|2 = |uB,W (α)|2 + σ 2 ||gα ||2 . Denoting by µ = σ||gα ||, the unknown original coefficients can be written as |uB,W (α)|2 = E|VB,W (α)|2 − µ2 . The observed coefficients |vB,W (α)|2 are used to approximate E|VB,W (α)|2 , and the estimated original coefficients are replaced in the optimal attenuation, leading to the family {a(α)}α , where ¾ ½ |vB,W (α)|2 − µ2 . a(α) = max 0, |vB,W (α)|2 Denote by EW Fµ (i) the filter given by the previous family of coefficients. The method noise of the EW Fµ (i) is easily computed, as proved in the following theorem. 15 Theorem 3.1.1 Let u be an image defined in a grid I and let i ∈ I be a pixel. Let W = W (i) be a window centered in the pixel i. Then the method noise of the EW Fµ (i) is given by u(i) − EW Fµ (i) = X {α||vB,W (α)|<µ} vB,W (α) gα (i) + X {α||vB,W (α)|>=µ} µ2 vB,W (α) gα (i). |vB,W (α)|2 The presence of an edge in the window W will produce a great amount of large coefficients and as a consequence, the cancelation of these coefficients will produce oscillations. Then, spurious cosines will also appear in the image under the form of chessboard patterns, see Figure 5. 3.2 Wavelet thresholding Let B = {gα }α∈A be an orthonormal basis of wavelets [23]. Let us discuss two procedures modifying the noisy coefficients, called wavelet thresholding methods (D. Donoho et al. [11]). The first procedure is a projection operator which approximates the ideal projection (10). It is called hard thresholding, and cancels coefficients smaller than a certain threshold µ, ½ 1 |vB (α)| > µ a(α) = 0 |vB (α)| ≤ µ Let us denote this operator by HWTµ (v). This procedure is based on the idea that the image is represented with large wavelet coefficients, which are kept, whereas the noise is distributed across small coefficients, which are cancelled. The performance of the method depends on the capacity of approximating u by a small set of large coefficients. Wavelets are for example an adapted representation for smooth functions. Theorem 3.2.1 Let u be an image defined in a grid I. The method noise of a hard thresholding HWTµ (u) is X uB (α)gα u − HWTµ (u) = {α||uB (α)|<µ} Unfortunately, edges lead to a great amount of wavelet coefficients lower than the threshold, but not zero. The cancellation of these wavelet coefficients causes small oscillations near the edges, i.e. a Gibbs-like phenomenon. Spurious wavelets can also be seen in the restored image due to the cancellation of small coefficients : see Figure 5. D. Donoho [10] showed that these effects can be partially avoided with the use of a soft thresholding, ( vB (α)−sgn(vB (α))µ |vB (α)| ≥ µ vB (α) a(α) = 0 |vB (α)| < µ which will be denoted by SWT µ(v). The continuity of the soft thresholding operator better preserves the structure of the wavelet coefficients, reducing the oscillations near discontinuities. Note that a soft thresholding attenuates all coefficients in order to reduce the noise, as an ideal operator does. As we shall see at the end of this paper, the L2 norm of the method noise is lessened when replacing the hard by a soft threshold. See Figures 5 and 14 for a comparison of the both method noises. Theorem 3.2.2 Let u be an image defined in a grid I. The method noise of a soft thresholding SWTµ (u) is X X sgn(uB (α)) gα uB (α)gα + µ u − SWTµ (u) = {α||uB (α)|>µ} {α||uB (α)|<µ} 16 A simple example can show how to fix the threshold µ. Suppose the original image u is zero, then vB (α) = nB (α), and therefore the threshold µ must be taken over the maximum of noise coefficients to ensure their suppression and the recovery of the original image. It can be pshown that the maximum amplitude of a white noise has a high probability of being smallerp than σ 2 log |I|. It can be proved that the risk of a wavelet thresholding with the threshold µ = σ 2 log |I| is near the risk rp of the optimal projection, see [11, 23]. p Theorem 3.2.3 The risk rt (u) of a hard or soft thresholding with the threshold µ = σ 2 log |I| is such that for all |I| ≥ 4 rt (u) ≤ (2 log |I| + 1)(σ 2 + rp (u)). (12) The factor 2 log |I| is optimal among all the diagonal operators in B, lim inf sup |I|−>∞ D∈DB u∈R|I| E{||u − DV ||2 } 1 = 1. σ 2 + rp (u) 2log|I| (13) In practice the optimal threshold µ is very high and cancels too many coefficients not produced by the noise. A threshold lower than the optimal is used in the experiments and produces much better results, see Figure 5. For a hard thresholding the threshold is fixed to 3 ∗ σ. For a soft thresholding this threshold still is too high ; it is better fixed at 32 σ. 3.3 Translation invariant wavelet thresholding R. Coifman and D. Donoho [7] improved the wavelet thresholding methods by averaging the estimation of all translations of the degraded signal. Calling v p (i) the translated signal v(i − p), the wavelet coefficients of the original and translated signals can be very different, and they are not related by a simple translation or permutation, vBp (α) = hv(n − p), gα (n)i = hv(n), gα (n + p)i. The vectors gα (n + p) are not in general in the basis B = {gα }α∈A , and therefore the estimation of the translated signal is not related to the estimation of v. This new algorithm yields an estimate ûp for every translated v p of the original image, X ûp = Dv p = a(α)vBp (α)gα . (14) α∈A The translation invariant thresholding is obtained averaging all these estimators after a translation in the inverse sense, 1 X p û (i + p), |I| (15) p∈I and will be denoted by TIWT (v). The Gibbs effect is considerably reduced by the translation invariant wavelet thresholding, (see Figure 5), because the average of different estimations of the image reduces the oscillations. This is therefore the version we shall use in the comparison section. 4 The Universal Denoiser The recent work by Ordentlich et al. [26] has led to the proposition of a “universal denoiser” for digital images. The authors assume that the noise model is fully known, namely the probability transition 17 Figure 5: Denoising experiment on a natural image. From left to right and from top to bottom: noisy image (standard deviation 20), Fourier Wiener filter (ideal filter), the DCT empirical Wiener filter, the wavelet hard thresholding, the soft wavelet thresholding and the translation invariant wavelet hard thresholding. matrix Π(a, b) , where a, b ∈ A, the finite alphabet of all possible values for the image. In order to fix ideas, we shall assume as in the rest of this paper that the noise is additive gaussian, in which case (a−b)2 1 one simply has Π(a, b) = √2πσ e− 2σ2 for the probability of observing b when the real value was a. The authors also fix an error cost Λ(a, b) which, to fix ideas, we can take to be a quadratic function Λ(a, b) = (a − b)2 , namely the cost of mistaking a for b. The authors fix a neighborhood shape, say, a square discrete window deprived of its center i, Ñi = Ni \ {i} around each pixel i. Then the question is : once the image has been observed in the window Ñi , what is the best estimate we can make from the observation of the full image ? Compute, for each possible value b of u(i) the number of windows Nj in the image such the restrictions of u to Ñj and Ñi coincide and the observed value at the pixel j is b. This number is called m(b, Ni ) and the line vector (m(b, Ni ))b∈A is denoted by m(Ni ). Then the authors propose to compute the denoised value of u at i as ũ(i) = arg min m(Ni )Π−1 (Λb ⊗ Πu(i) ), b∈A where w ⊗ v = (w(b)v(b)) denotes the vector obtained by multiplying each component of u by each component of v, u(i) is the observed value at i, and we denote by Xa the a-column of a matrix X. The authors prove that this denoiser is universal in the sense “of asymptotically achieving, without access to any information on the statistics of the clean signal, the same performance as the best denoiser that does have access to this information”. This asymptote is achieved when the size of the image tends to infinity. The authors present an implementation valid for binary images with an impulse noise, with excellent results. The reason of these limitations in implementation are clear : first the 18 matrix Π is of very low dimension and invertible for impulse noise. If instead we consider as above a gaussian noise, then the application of Π−1 amounts to deconvolve a signal by a gaussian, which is a rather ill-conditioned method. All the same, it is doable, while the computation of m certainly is not for a large alphabet, like the one involved in grey tone images (256 values). Even supposing that the learning window Ni has the minimal possible size of 9, the number of possible such windows is about 2569 which turns out to be much larger than the number of observable windows in an image (whose typical size amounts to 106 pixels). Thus, it is fully unrealistic to hope to observe the number m(Ni ) in an image, particularly if it is corrupted by noise. In the next section, we propose a more realistic method to learn the image model from the image itself. 5 Non local means algorithm (NL-means) The local smoothing methods and the frequency domain filters aim at a noise reduction and at a reconstruction of the main geometrical configurations but not at the preservation of the fine structure, details and texture. Due to the regularity assumptions on the original image of previous methods, details and fine structures are smoothed out because they behave in all functional aspects as noise. The NL-means algorithm we shall now discuss tries to take advantage of the high degree of redundancy of any natural image. By this, we simply mean that every small window in a natural image has many similar windows in the same image. This fact is patent for windows close by, at one pixel distance and in that case we go back to a local regularity assumption. Now in a very general sense inspired by the neighborhood filters, one can define as “neighborhood of a pixel i” any set of pixels j in the image such that a window around j looks like a window around i. All pixels in that neighborhood can be used for predicting the value at i, as was first shown in [13]. The fact that such a self-similarity exists is a regularity assumption, actually more general and more accurate than all regularity assumptions we have considered in section 2. It also generalizes a periodicity assumption of the image. Let v be the noisy image observation and let x ∈ Ω. The NL-means algorithm estimates the value of x as an average of the values of all the pixels whose gaussian neighborhood looks like the neighborhood of x, Z (Ga ∗|v(x+.)−v(y+.)|2 )(0) 1 h2 N L(v)(x) = e− v(y) dy, C(x) R (Ga ∗|v(x+.)−v(z+.)|2 )(0) h2 dz is the where Ga is a Gaussian kernel with standard deviation a, C(x) = e− normalizing factor and h acts as a filtering parameter. Since we are considering images defined on a discrete grid I, we shall give a discrete description of the NL-means algorithm and some consistency results. 5.1 Description Given a discrete noisy image v = {v(i) | i ∈ I}, the estimated value N L(v)(i) is computed as a weighted average of all the pixels in the image, X N L(v)(i) = w(i, j)v(j), j∈I where the family of weights {w(i, j)}j depends on the similarity between the pixels i and j, and satisfy P the usual conditions 0 ≤ w(i, j) ≤ 1 and j w(i, j) = 1. In order to compute the similarity between the image pixels, define a neighborhood system on I. Definition 5.1.1 (Neighborhoods) A neighborhood system on I is a family N = {Ni }i∈I of subsets of I such that for all i ∈ I, 19 (i) i ∈ Ni (ii) j ∈ Ni ⇒ i ∈ Nj The subset Ni is called the neighborhood or the similarity window of i. The subset Ñi will denote Ni \{i}. The similarity windows can have different sizes and shapes to better adapt to the image. For simplicity we will use square windows of fixed size. The restriction of v to a neighborhood Ni will be denoted by v(Ni ), v(Ni ) = (v(j), j ∈ Ni ). The similarity between two pixels i and j will depend on the similarity of the intensity gray level vectors v(Ni ) and v(Nj ). The pixels with a similar grey level neighborhood to v(Ni ) will have larger weights in the average, see Figure 6. Figure 6: q1 and q2 will have a large weight because their similarity windows are similar to that of p. On the other side the weight w(p,q3) will be much small because intensity grey values in the similarity windows are much different. In order to compute the similarity of the intensity gray level vectors v(Ni ) and v(Nj ), one can compute a gaussian weighted Euclidean distance, kv(Ni ) − v(Nj )k22,a . Efros and Leung [13] showed that the L2 distance is a reliable measure for the comparison of image windows in a texture patch. Now, this measure is so much the more adapted to any additive white noise as such a noise alters the distance between windows in a uniform way. Indeed, E||v(Ni ) − v(Nj )||22,a = ||u(Ni ) − u(Nj )||22,a + 2σ 2 where u and v are respectively the original and noisy images and σ 2 is the noise variance. This equality shows that, in expectation, the Euclidean distance preserves the order of similarity between pixels. So the most similar pixels to i in v also are expected to be the most similar pixels of i in u. The weights associated with the quadratic distances are defined by 2 w(i, j) = j )||2,a 1 − ||v(Ni )−v(N h2 e , Z(i) 20 ||v(Ni )−v(Nj )||2 2,a P h2 where Z(i) is the normalizing factor Z(i) = j e− and the parameter h controls the decay of the exponential function and therefore the decay of the weights as a function of the Euclidean distances. 5.2 A consistency theorem for NL-means The NL-means algorithm is intuitively consistent under stationarity conditions, saying that one can find many samples of every image detail. In fact, we shall be assuming that the image is a fairly general stationary random process. Under these assumptions, for every pixel i, the NL-means algorithm converges to the conditional expectation of i knowing its neighborhood. In the case of an additive or multiplicative white noise model, this expectation is in fact the solution to a minimization problem. Let X and Y denote two random vectors with values on Rp and R respectively. Let fX , fY denote the probability distribution functions of X, Y and let fXY denote the joint probability distribution function of X and Y . Let us recall briefly the definition of the conditional expectation. Definition 5.2.1 i) Define the probability distribution function of Y conditioned to X as, ( fXY (x,y) fX (x) > 0 fX (x) f (y | x) = 0 otherwise for all x ∈ Rp and y ∈ R. ii) Define the conditional expectation of Y given {X = x} as the expectation with respect to the conditional distribution f (y | x) Z E[Y | X = x] = y f (y | x) dy, for all x ∈ Rp . Remark 5.2.2 The conditional expectation is a function of X and therefore a new random variable g(X) which is denoted by E[Y | X]. Let now V be a random field and N a neighborhood system on I. Let Z denote the sequence of random variables Zi = {Yi , Xi }i∈I where Yi = V (i) is real valued and Xi = V (Ñi ) is Rp valued. Recall that Ñi = Ni \{i}. Let us restrict Z to the n first elements {Yi , Xi }ni=1 . Let us define the function rn (x), rn (x) = Rn (x)/fˆn (x) where n Xi − x 1 X K( ), fˆn (x) = p nh i=1 h Rn (x) = n Xi − x 1 X φ(Yi )K( ), p nh i=1 h (16) (17) φ is an integrable real valued function, K is a nonnegative kernel and x ∈ Rp Let X and Y be distributed as X1 and Y1 . Under this form the NL-means algorithm can be seen as an instance for the exponential operator of the Nadaraya-Watson estimator [25, 39]. This is an estimator of the conditional expectation r(x) = E[φ(Y ) | X = x]. Some definitions are needed for the statement of the main result. Definition 5.2.3 A stochastic process {Zt | t = 1, 2, . . .}, with Zt defined on some probability space (Ω, A, P), is said to be (strict-sense) stationary if for any finite partition {t1 , t2 , · · · , tn } the joint distributions Ft1 ,t2 ,···,tn (x1 , x2 , · · · , xn ) are the same as the joint distributions Ft1 +τ,t2 +τ,···,tn +τ (x1 , x2 , · · · , xn ) for any τ ∈ N . 21 In the case of images, this stationary condition amounts to say that as the size of the image grows, we are able to find in the image many similar patches for all the details of the image. This is a crucial point to understand the performance of the NL-means algorithm. The following mixing definition is a rather technical condition. In the case of images, it amounts to say that regions become more independent as it increases the distance between them. This is intuitively true for natural images. Definition 5.2.4 Let Z be a stochastic and stationary process {Zt | t = 1, 2, · · · , n}, and, for m < n, let Fnm be the σ − F ield induced in Ω by the r.v.’s Zj , m ≤ j ≤ n. Then the sequence Z is said to be β − mixing if for every A ∈ Fk1 and every B ∈ F∞ k+n |P (A ∩ B) − P (A)P (B)| ≤ β(n) with β(n) → 0, as n → ∞. The following theorem establishes the convergence of rn to r, see Roussas [31]. The theorem is established under the stationary and mixing hypothesis of {Yi , Xi }∞ i=1 and asymptotic conditions on the decay of φ, β(n) and K. This set of conditions will be denoted by H and it is more carefully detailed in the Appendix. Theorem 5.2.5 (Conditional expectation theorem) Let Zj = {Xj , Yj } for j = 1, 2, · · · be a strictly stationary and mixing process. For i ∈ I let X and Y be distributed as Xi and Yi . Let J be a compact subset J ⊂ Rp such that inf{fX (x); x ∈ J} > 0. Then, under hypothesis H, sup[ψn |rn (x) − r(x)|; x ∈ J] → 0 a.s, where ψn are positive norming factors. Let v be the observed noisy image and let i be a pixel. Taking φ as the identity we see that rn (v(Ñi )) converges to E[V (i) | V (Ñi ) = v(Ñi )] under stationary and mixing conditions of the sequence {V (i), V (Ñi )}∞ i=1 . In the case where an additive or multiplicative white noise model is assumed, the next result shows that this conditional expectation is in fact the function of V (Ñi ) that minimizes the mean square error with the original field U . Theorem 5.2.6 Let V, U, N1 , N2 be random fields on I such that V = U + N1 + g(U )N2 , where N1 and N2 are independent white noises. Let N be a neighborhood system on I . Then, the following statements are hold. (i) E[V (i) | V (Ñi ) = x] = E[U (i) | V (Ñi ) = x] for all i ∈ I and x ∈ Rp . (ii) The real value E[U (i) | V (Ñi ) = x] minimizes the following mean square error, min E[(U (i) − g ∗ )2 | V (Ñi ) = x] g ∗ ∈R (18) for all i ∈ I and x ∈ Rp . (iii) The expected random variable E[U (i) | V (Ñi )] is the function of V (Ñi ) that minimizes the mean square error min E[U (i) − g(V (Ñi ))]2 (19) g Remark 5.2.7 Given a noisy image observation v(i) = u(i) + n1 (i) + g(u(i))n2 (i), i ∈ I, where g is a real function and n1 and n2 are white noise realizations, then the NL-means algorithm is the function of v(Ñi ) that minimizes the mean square error with the original image u(i). 22 5.3 Experiments with NL-means The NL-means algorithm chooses for each pixel a different average configuration adapted to the image. As we explained in the previous sections, for a given pixel i, we take into account the similarity between the neighborhood configuration of i and all the pixels of the image. The similarity between pixels is measured as a decreasing function of the Euclidean distance of the similarity windows. Due to the fast decay of the exponential kernel, large Euclidean distances lead to nearly zero weights, acting as an automatic threshold. The decay of the exponential function and therefore the decay of the weights is controlled by the parameter h. In order to ensure that only pixels with a similar neighborhood have a large weight, h must be taken of the order of σ. (This is just an empirical observation, since, as indicated by the results in the previous section, the algorithm is consistent for every positive value of h). Empirical experimentation shows however that one can take a similarity window of size 7 × 7 or 9 × 9 for grey level images and 5 × 5 or even 3 × 3 in color images with little noise. These window sizes have shown to be large enough to be robust to noise and at the same time able to take care of the details and fine structure. Smaller windows are not robust enough to noise. Notice that in the limit case, one can take the window reduced to a single pixel i and get therefore back to the Yaroslavsky neighborhood filter. Figure 7: NL-means denoising experiment with a nearly periodic image. Left: Noisy image with standard deviation 30. Right: NL-means restored image. Due to the nature of the algorithm, the most favorable case for the NL-means is the periodic case. In this situation, for every pixel i of the image one can find a large set of samples with a very similar configuration, leading to a noise reduction and a preservation of the original image, see Figure 7 for an example. Another case which is ideally suitable for the application of the NL-means algorithm is the textural case. Texture images have a large redundancy. For a fixed configuration many similar samples can be found in the image. In Figure 8 one can see an example with a Brodatz texture. The Fourier Transform of the noisy and restored images show the ability of the algorithm to preserve the main features even in the case of high frequencies. The NL-means is not only able to restore periodic or texture images. Natural images also have enough redundancy to be restored by NL-means. For example in a flat zone, one can find many pixels lying in the same region and similar configurations. In a straight or curved edge a complete line of pixels with a similar configuration is found. In addition, redundancy of natural images allows us to find many similar configurations in far away pixels. Figures 9 and 10 show two examples on two well known standard processing images. The same algorithm applies to the restoration of color images and films, see Figure 11. 23 Figure 8: NL-means denoising experiment with a Brodatz texture image. Left: Noisy image with standard deviation 30. Right: NL-means restored image. The Fourier Transform of the noisy and restored images show how main features are preserved even at high frequencies. 5.4 Testing stationarity : a soft threshold optimal correction In this section, we describe a simple and useful statistical improvement of NL-means, with a technique similar to the wavelet thresholding. The stationarity assumption of Theorem 5.2.5 is not true everywhere, as each image may contain exceptional, non repeated structures. Such structures can be blurred out by the algorithm. The NL-means algorithm, and actually every local averaging algorithm, must involve a detection phase and special treatment of non stationary points. The principle of such a correction is quite simple and directly derived from other thresholding methods, like the SWT method. Let us estimate the original value at a pixel i, u(i), as the mean of the noisy grey levels v(j) for j ∈ J ⊂ I. One expects averaged pixels to have a non noisy grey level value, u(j), similar to u(i) in order to reduce the noise and restore the original value. Assuming this fact, û(i) = 1 X 1 X v(j) ≃ u(i) + n(j) → u(i) as |J| → ∞, |J| |J| j∈J j∈J because the average of noise values tends to zero. In addition, 1 X 1 X (v(j) − û(i))2 ≃ n(j)2 → σ 2 as |J| → ∞. |J| |J| j∈J j∈J If the averaged pixels have a non noisy grey level value close to u(i), as expected, then the variance of the average should be close to σ 2 . If it is a posteriori observed that this variance is much larger than σ 2 , this fact can hardly be caused by the noise only. This means that NL is averaging pixels whose original grey level values were very different in the original. At those pixels, a more conservative 24 Figure 9: NL-means denoising experiment with a natural image. Left: Noisy image with standard deviation 20. Right: Restored image. Figure 10: NL-means denoising experiment with a natural image. Left: Noisy image with standard deviation 35. Right: Restored image. estimate is required, and therefore the estimated value should be averaged with the noisy one. The next result tells us how to compute this average. Theorem 5.4.1 Let X and Y be two real random variables. Then, the linear estimate Ŷ , Ŷ = EY + Cov(X, Y ) (X − EY ) V arX minimizes the square error min E[(Y − (a + bX))2 ] a,b∈R In our case, X = Y + N where N is independent of Y, with zero mean and variance σ 2 . Thus, Ŷ = EY + which is equal to V arY (X − EY ), V arY + σ 2 µ ¶ V arX − σ 2 Ŷ = EX + max 0, (X − EX). V arX This strategy can be applied to correct any local smoothing filter. However, a good estimate of the mean and the variance at every pixel is needed. That is not the case for the local smoothing filters 25 Figure 11: NL-means denoising experiment with a color image. Left: Noisy image with standard deviation 15 in every color component. Right: Restored image. of Section 2. This strategy can instead be satisfactorily applied to the NL-means algorithm. As we have shown in the previous section, the NL-means algorithm converges to the conditional mean. The conditional variance can be also computed by the NL-means, by taking φ(x) = x2 in Theorem 5.2.5, and then computing the variance as EX 2 − (EX)2 . In Figure 12 one can see an application of this correction. Figure 12: Optimal correction experience. Left: Noisy image. Middle: NL-means solution. Right: NL-means corrected solution. The average with the noisy image makes the solution to be noisier, but details and fine structure are better preserved. 6 6.1 Discussion and Comparison NL-means as an extension of previous methods. As was said before, the gaussian convolution only preserves flat zones while contours and fine structure are removed or blurred. Anisotropic filters instead preserve straight edges but flat zones present many artifacts. One could think of combining these two methods to improve both results. A gaussian convolution could be applied in flat zones while an anisotropic filter could be applied on straight edges. Still, other types of filters should be designed to specifically restore corners or curved edges and 26 texture. The NL-means algorithm seems to provide a feasible and rational method to automatically take the best of each mentioned algorithm, reducing for every possible geometric configuration the image method noise. Although we have not computed explicitly the image method noise, Figure 13 illustrates how the NL-means algorithm chooses in each case a weight configuration corresponding to one of the previously analyzed filters. In particular, according to this set of experiments, we can consider that the consistency Theorems 2.1.1, 2.2.1 and 2.3.1 are all valid for this algorithm. (a) (b) (c) (d) (e) (f) Figure 13: On the right-hand side of each pair, we display the weight distribution used to estimate the central pixel of the right image by the NL-means algorithm. Image a: In flat zones, the weights are distributed as a convolution filter (as a gaussian convolution). Image b: In straight edges, the weights are distributed in the direction of the level line (as the mean curvature motion). Image c: On curved edges, the weights favor pixels belonging to the same contour or level line, which is a strong improvement w.r. to the mean curvature motion. Image d: In a flat neighborhood, the weights are distributed in a grey level neighborhood (as the neighborhood filters). In the case of Images e and f, the weights are distributed across the more similar configurations, even though are far away from the observed pixel. This shows a behavior similar to a nonlocal neighborhood filter or to an ideal Wiener filter. 27 6.2 Comparison In this section we shall compare the different algorithms based on three well defined criteria: the method noise, the mean square error and the visual quality of the restored images. Note that every criterion measures a different aspect of the denoising method. It is easy to show that only one criterion is not enough to judge the restored image, and so one expects a good solution to have a high performance under the three criteria. 6.2.1 Method noise comparison In previous sections we have defined the method noise and computed it for the different algorithms. Remember that the denoising algorithm is applied on the original (slightly noisy) image. A filtering parameter, depending mainly on the standard deviation of the noise, must be fixed for the most part of algorithms. Let us fix σ = 2.5: we can suppose that any digital image is affected by this amount of noise since it is not visually noticeable. The method noise tells us which geometrical features or details are preserved by the denoising process and which are eliminated. In order to preserve as much features as possible of the original image the method noise should look as much as possible like white noise. Figures 14-17 display the method noise of the different methods for a set of standard natural images. Let us comment them briefly. • The gaussian filter method noise highlights all important features of the image like texture, contours and details. All these features have a large Laplacian and are therefore modified by the application of the algorithm, see Theorem 2.1.1. • As announced in Theorem 2.2.1, the anisotropic filter method noise displays the corners and high frequency features. The straight edges are instead not to be seen : they have a low curvature. • The Total Variation method modifies most structures and details of the image. Even straight edges are not well preserved. • The iterated Total Variation refinement improves the total variation method noise. The straight edges and details are much better preserved. • The neighborhood filter preserves flat objects and contrasted edges, while edges with a low contrast are not kept. In any case, the contours, texture and details seem to be well preserved. • The TIHWT method noise is concentrated on the edges and high frequency features. These structures lead to coefficients of large enough value but lower than the threshold. They are removed by the algorithm. The average of the application to all translated versions reduces the method noise, and structures are hardly noticeable. • The TISWT method noise presents much more structure than the hard thresholding. Indeed, the method noise is not only based on the small coefficients but also on an attenuation of the large ones, leading to a high alteration of the original image. • It is difficult to find noticeable structure in the DCT empirical Wiener filter method noise. Only some contours are noticeable. In general, this filter seems to perform much better than all local smoothing filters and other frequency domain filters. Its results are similar to those of a hard stationary wavelet thresholding. • The NL-means method noise looks the more like white noise. 28 6.2.2 Visual quality comparison As commented before, the visual quality of the restored image is another important criterion to judge the performance of a denoising algorithm. Let us present some experiences on a set of standard natural images. The objective is to compare the visual quality of the restored images, the non presence of artifacts and the correct reconstruction of edges, texture and fine structure. Figures 18-21 present these experiences comparing the visual quality of previous methods. Figure 18 illustrates the fact that a non local algorithm is needed for the correct reconstruction of periodic images. Local smoothing filters and local frequency filters are not able to reconstruct the wall pattern. Only the NL-means algorithm and the global Fourier Wiener filter reconstruct the original texture. The Fourier Wiener filter is based on a global Fourier transform which is able to capture the periodic structure of the image in a few coefficients. Now, in practice, this is an ideal filter because the Fourier Transform of the original image is used. Figure 13 e) shows how the NL-means method chooses the correct weight configuration and explains the correct reconstruction of the wall pattern. Figure 19 illustrates the difficulty of local smoothing filters for recovering stochastic patterns. The high degree of noise present in the image makes the local comparisons of the neighborhood filter much dependent of noise. As a consequence, noise and texture are not well differentiated. The regularity assumption involved in the bounded variation makes it unsuitable for the restoration of textures which are filtered as noise. Figure 20 shows that the frequency domain filters are well adapted to the recovery of oscillatory patterns. Although some artifacts are noticeable in both solutions, the stripes are well reconstructed. The DCT transform seems to be more adapted to this type of texture and stripes are a little better reconstructed. The NL-means also performs well on this type of texture, due to its high degree of redundancy. Finally, Figure 21 shows how both the wavelets thresholdings and the iterated Total Variation are able to remove the noise in presence of geometry and texture. However, both methods can create the artifacts or over smoothing. 6.2.3 Mean square error comparison The mean square error is the square of the Euclidean distance between the original image and its estimate. This numerical quality measurement is the more objective one, since it does not rely on any visual interpretation. Table 1 shows the mean square error of the different denoising methods with the images presented in this paper. This error table seems to corroborate the observations made for the other criteria. One sees for example how the iterated Total Variation improves the original minimization, how the frequency domain filters have a lower mean square error than the local smoothing filters. One also sees that in presence of periodic or textural structures the Empirical Wiener Filter based on a DCT transform performs better than the wavelet thresholding, see also Figures 18 and 20. Note that, in presence of periodic or stochastic patterns, the NL-means mean square error is significantly more precise than the other algorithms. Of course, the errors presented in this table cannot be computed in a real denoising problem. Let us remark that a small error does not guarantee a good visual quality of the restored image. The mean square error by itself would not be meaningful and all previous quality criteria are also necessary to evaluate the performance of denoising methods. 7 Appendix The set H of assumptions necessary for the statement of Theorem 5.2.5 are: 29 Image Boat Lena Barbara Baboon Wall σ 8 20 25 35 35 GF 53 120 220 507 580 AF 38 114 216 418 660 TV 39 110 186 365 721 TV Iter 29 82 189 364 715 YNF 39 129 176 381 598 EWF 33 105 111 396 325 TIHWT 28 81 135 365 712 NL-means 23 68 72 292 59 Table 1: Mean square error table. A smaller mean square error indicates that the estimate is closer to the original image. The numbers have to be compared on each row. The square of the number on the right column gives the real variance of the noise. By comparing this square to the values on the same row, it is quickly checked that all studied algorithms indeed perform some denoising. A sanity check! In general, the comparison performance corroborates the previously mentioned quality criteria. p (H1) The sequence of random vectors Zi = {Yi , Xi }∞ i=1 where Yi is real valued and Xi is R valued form a strictly stationary sequence. (H2) The sequence {Zi }P is β − mixing and the sequence β(n) satisfies the following summability ∞ requirement: β ∗ = j=1 β(n) < ∞. (H3) Let α = α(n) be a positive integer and let µ = µ(n) be the largest positive integer for which 2αµ ≤ n. Then 1 lim sup[1 + 6e 2 β 1/(µ+1) (α)]µ < ∞. (H4) kxkp K(x) → 0, as x → ∞, where the norm kxk of x = (x1 , . . . , xp ) is defined by kxk = max(|x1 |, . . . , |xp |). (H5) i) φ is a real valued Borel function defined on R such that E|φ(Y )|s < ∞ for some s > 1. ii) sup[ Z R (H6) ′ |φ(y)|s fXY (x, y)dy; x ∈ Rp ] = C < ∞. p i) For any point x and x in R and for some positive constant C(independent of these points): |K(x) − K(x′ )| ≤ Ckx − x′ k ii) R kxkK(x)dx < ∞. (H7) For any point x in Rp , there are positive constants C(x) such that, for all x′ ∈ Rp and with J being as in (A7): i) ||fX (x) − fX (x′ )|| ≤ C(x)||x − x′ ||, sup[C(x); x ∈ J] < ∞. ii) ||ψ(x) − ψ(x′ )|| ≤ C(x)||x − x′ ||, sup[C(x); x ∈ J] < ∞, where r(x) = E[φ(Y ) | X = x] and ψ(x) = r(x)fX (x). (A8) There exists a compact subset J of Rp such that inf[fX (x); x ∈ J] > 0. 30 References [1] L. Alvarez and P-L. Lions and J-M. Morel, Image selective smoothing and edge detection by nonlinear diffusion (II), SIAM Journal of numerical analysis 29, pp. 845-866, 1992. [2] J.F. Aujol, G. Aubert, L. Blanc-Féraud, and A. Chambolle, Decomposing an image: Application to textured images and SAR images, 2003, INRIA Research Report 4704, to appear in JMIV. [3] B. Merriman, J. Bence and S. Osher, Diffusion generated motion by mean curvature, Computational Crystal Growers Workshop, Proceedings of the Geometry Center Workshop held in Minneapolis, Minnesota, February 22-28, 1992. [4] A. Chambolle and P. L. Lions, Image recovery via total variation minimization and related problems, Numer. Math. 76, pp. 167-188, 1997. [5] A. Chambolle, R. De Vore, N. Lee and B. Lucier, Nonlinear wavelet image processing: variational problems, compression, and noise removal through wavelet shrinkage, IEEE Tran. Image Proc., Vol 7, No 3, pp. 319-333, 1998. [6] T.F. Chan and H.M. Zhou, Total variation improved wavelet thresholding in image compression. Proceedings of ICIP’2000, Volume 2, pp. 391-394, 2000. [7] R.R. Coifman and D. Donoho, Translation-invariant de-noising, Wavelets and Statistics, pp. 125150, Springer Verlag, 1995. [8] R. R. Coifman and A. Sowa, Combining the calculus of variations and wavelets for image enhancement, Applied and Computational harmonic analysis, 9(1), pp. 1-18, 2000. [9] F. Catté, F. Dibos and G. Koepfler, A morphological scheme for mean curvature motion and applications to anisotropic diffusion and motion of level sets, SIAM J. Numerical Analysis 32 (6), pp. 1845-1909, 1995. [10] D. Donoho, De-noising by soft-thresholding, IEEE Transactions on Information Theory, 41, pp. 613-627, 1995. [11] D. Donoho and I. Johnstone, Ideal spatial adaptation via wavelet shrinkage, Biometrika, 81 pp. 425-455, 1994. [12] S. Durand and J. Froment, Reconstruction of wavelet coefficients using total variation minimization, SIAM Journal on Scientific computing, 24(5), pp. 1754-1767, 2003. [13] A. Efros and T. Leung, Texture synthesis by non parametric sampling, Proc. Int. Conf. Computer Vision (ICCV 99), Vol. 2, pp. 1033-1038, 1999. [14] S. Geman and D. Geman, Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images, IEEE Pat. Anal. Mach. Intell., 6, pp. 721-741, 1984. [15] F. Guichard, J.M Morel and R. Ryan, Image Analysis and P.D.E.’s. Preprint. [16] J. Heeger and J. R. Bergen, Pyramid based texture analysis synthesis, Computer Graphics, pp. 229-238. ACM SIGGRAPH, 1995. [17] E. Le Pennec and S. Mallat, Geometrical image compression with bandelets, Proceedings of the SPIE 2003, 5150, pp. 1273-1286, 2003. [18] E. Levina, Statistical Issues in Texture Analysis. Ph dissertation 2002. 31 [19] M. Lindenbaum and M. Fischer and A. M. Bruckstein, On Gabor Contribution To ImageEnhancement, Pattern Recognition 27, pp. 1-8, 1994. [20] S. Lintner and F. Malgouyres, Solving a variational image restoration model which involves L∞ constraints. University Paris 13, preprint no 2003-18, July 2003. [21] F. Malgouyres, A noise selection approach of image restoration, Applications in signal and image processing IX, volume 4478, pp. 34-41, 2001. [22] F. Malgouyres, Minimizing the total variation under a general convex constraint for image restoration, IEEE Trans. on Image Processing, 11(12), pp. 1450-1456, 2002. [23] S. Mallat, A wavelet tour of signal processing, Academic Press, 1997. [24] Y. Meyer, Oscillating Patterns in Image Processing and Nonlinear Evolution Equations, University Lecture Series Volume 22, AMS 2002. [25] E. A. Nadaraya, On estimating regression. Theory of Probability and its Applications 10: 186-190, 1964. [26] E. Ordentlich, G. Seroussi, S. Verdú, M. Weinberger and T. Weissman, A discrete universal denoiser and its application to binary images. Proc. IEEE ICIP, Vol 1, pp. 117-120, 2003. [27] S. Osher, M. Burger, D. Goldfarb, J. Xu and W. Yin, Using Geometry and iterated refinement for inverse problems (1): Total variation based image restoration, CAM-Report 04-13 UCLA, 2004. To appear in Multiscale Modeling and Simulation, SIAM. [28] S. Osher and L. Rudin. Feature oriented image enchancement using shock filters SIAM J. Numerichal Analysis, 27, pp. 919-940. 1990. [29] S.J. Osher, A. Sole, and L.A. Vese, Image decomposition and restoration using total variation minimization and the H −1 norm, Multiscale Modeling and Simulation: A SIAM Interdisciplinary Journal, 1(3), pp.349-370, 2003. [30] P. Perona and J. Malik, Scale space and edge detection using anisotropic diffusion, IEEE Trans. Patt. Anal. Mach. Intell., 12, pp. 629-639, 1990. [31] G.G Roussas, Nonparametric regression estimation under mixing conditions, Stochastic processes and their applications, Vol 36, pp. 107-116. 1990. [32] L. Rudin and S. Osher, Total Variation based image restoration with free local constraints, Proc. IEEE ICIP, Vol 1, pp. 31-35, 1994. [33] L. Rudin and S. Osher, Nonlinear total variation based noise removal algorithms, Physica D, 60, pp. 259-268, 1992. [34] C. Shannon. and W. Weaver, The mathematical theory of communication, University of Illinois Press, 1998. [35] S.M. Smith and J.M. Brady, Susan - a new approach to low level image processing, International Journal of Computer Vision, Volume 23 (1), pp. 45-78, 1997. [36] J. Starck, E. Candès, D. Donoho, The curvelet transform for image denoising IEEE Transactions on image processing, 11, pp. 670-684, 2000. 32 [37] J. Starck, E. Candès, D. Donoho, Very high quality image restoration by combining wavelets and curvelets Wavelet and applications in signal and image processing IX, volume 4478, pp. 9-19, 2001. [38] L. Vese and S. Osher, Modeling textures with total variation minimization and oscillating patterns in image processing Journal of Scientific Computing, 19 (1-3), pp. 553-572, 2003. [39] G. S. Watson, Smooth regression analysis. Sankhya - The Indian Journal of Statistics 26, pp. 359-372, 1964. [40] L.P. Yaroslavsky, Digital Picture Processing - An Introduction, Springer Verlag, 1985. [41] L. Yaroslavsky, M. Eden Fundamentals of Digital Optics. Birkhauser, Boston, 1996. [42] L. Yaroslavsky, Local Adaptive Image Restoration and enhancement with the use of DFT and DCT in a running window. Proceedings, Wavelet applications in signal and image processing IV, SPIE Proc. Series, v. 2825, pp. 1-13, 1996. 33 Figure 14: Image method noise. From left to right and from top to bottom: Gaussian convolution, Mean curvature motion, Total Variation, Iterated Total Variation, Neighborhood filter, Hard TIWT, Soft TIWT, DCT empirical Wiener filter and NL-means. 34 Figure 15: Zoom on image method noise of Figure 14. Top left: original image. Top right: Hard TIWT. Bottom left: DCT empirical Wiener filter. Bottom right: NL-means. 35 Figure 16: Image method noise. Top left: original image. Top right: Hard TIWT. Bottom left: DCT empirical Wiener filter. Bottom right: NL-means. 36 Figure 17: Image method noise. Top left: original image. Top right: Iterated Total variation. Bottom left: Neighborhood filter. Bottom right: NL-means. 37 Figure 18: Denoising experience on a periodic image. From left to right and from top to bottom: noisy image (standard deviation 35), Gauss filtering, Total variation, Neighborhood filter, Wiener filter (ideal filter), Hard TIWT, DCT empirical Wiener filtering, NL-means algorithm. Figure 19: Denoising experience on a natural image. From left to right and from top to bottom: noisy image (standard deviation 35), Neighborhood filter, Total variation, NL-means algorithm. 38 Figure 20: Denoising experience on a natural image. From left to right and from top to bottom: noisy image (standard deviation 25), DCT empirical Wiener filter, Hard TIWT, NL-means algorithm. 39 Figure 21: Denoising experience on a natural image. From left to right and from top to bottom: noisy image (standard deviation 25), Iterated Total variation, Hard TIWT, NL-means algorithm. 40