A non-local algorithm for image denoising
Antoni Buades, Bartomeu Coll
Dpt. Matemàtiques i Informàtica, UIB
Ctra. Valldemossa Km. 7.5,
07122 Palma de Mallorca, Spain
vdmiabc4@uib.es, tomeu.coll@uib.es

Abstract
We propose a new measure, the method noise, to evaluate and compare the performance of digital image denoising methods. We ﬁrst compute and analyze this method
noise for a wide class of denoising algorithms, namely the
local smoothing ﬁlters. Second, we propose a new algorithm, the non local means (NL-means), based on a non local averaging of all pixels in the image. Finally, we present
some experiments comparing the NL-means algorithm andthe local smoothing ﬁlters.

1. Introduction
The goal of image denoising methods is to recover theoriginal image from a noisy measurement,v(i) = u(i) + n(i),(1)

where v(i) is the observed value, u(i) is the “true” value andn(i) is the noise perturbation at a pixel i. The best simple
way to model the effect of noise on a digital image is to add
a gaussian white noise. In that case, n(i) are i.i.d. gaussian
values with zero mean and variance σ 2 .
Several methods have been proposed to remove the noiseand recover the true image u. Even though they may be very
different in tools it must be emphasized that a wide classshare the same basic remark : denoising is achieved by averaging. This averaging may be performed locally: the Gaussian smoothing model (Gabor [7]), the anisotropic ﬁltering(Perona-Malik [11], Alvarez et al. [1]) and the neighborhood ﬁltering (Yaroslavsky [16], Smith et al. [14], Tomasi
et al. [15]), by the calculus of variations: the Total Variation minimization (Rudin-Osher-Fatemi [13]), or in the frequency domain: the empirical Wiener ﬁlters (Yaroslavsky[16]) and wavelet thresholding methods (Coiffman-Donoho[5, 4]).

Jean-Michel Morel
CMLA, ENS Cachan
61, Av du Président Wilson
94235 Cachan, France
morel@cmla.ens-cachan.fr

Formally we deﬁne a denoising method Dh as a decompositionv = Dh v + n(Dh , v),
where v is the noisy image and h is a ﬁltering parameter which usually depends on the standard deviation of thenoise. Ideally, Dh v is smoother than v and n(Dh , v) looks
like the realization of a white noise. The decompositionof an image between a smooth part and a non smooth or
oscillatory part is a current subject of research (for example
Osher et al. [10]). In [8], Y. Meyer studied the suitable functional spaces for this decomposition. The primary scope of
this latter study is not denoising since the oscillatory partcontains both noise and texture.
The denoising methods should not alter the original image u. Now, most denoising methods degrade or remove theﬁne details and texture of u. In order to better understandthis removal, we shall introduce and analyze the method
noise. The method noise is deﬁned as the difference between the original (always slightly noisy) image u and its
denoised version.
We also propose and analyze the NL-means algorithm,
which is deﬁned by the simple formula(Ga ∗|u(x+.)−u(y+.)|2 )(0)1h2e−u(y) dy,
N L[u](x) =C(x) Ω

(Ga ∗|u(x+.)−u(z+.)|2 )(0)h2where x ∈ Ω, C(x) = Ω e−dz is anormalizing constant, Ga is a Gaussian kernel and h acts
as a ﬁltering parameter. This formula amounts to say thatthe denoised value at x is a mean of the values of all points
whose gaussian neighborhood looks like the neighborhoodof x. The main difference of the NL-means algorithm withrespect to local ﬁlters or frequency domain ﬁlters is the systematic use of all possible self-predictions the image canprovide, in the spirit of [6]. For a more detailed analysis on
the NL-means algorithm and a more complete comparison,see [2].
Section 2 introduces the method noise and computes its
mathematical formulation for the mentioned local smooth-ing ﬁlters. Section 3 gives a discrete deﬁnition of the NLmeans algorithm. In section 4 we give a theoretical resulton the consistency of the method. Finally, in section 5 wecompare the performance of the NL-means algorithm andthe local smoothing ﬁlters.

2. Method noiseDeﬁnition 1 (Method noise) Let u be an image and Dh adenoising operator depending on a ﬁltering parameter h.
Then, we deﬁne the method noise as the image differenceu − Dh u.
The application of a denoising algorithm should not alter the non noisy images. So the method noise should bevery small when some kind of regularity for the image isassumed. If a denoising method performs well, the method
noise must look like a noise even with non noisy images
and should contain as little structure as possible. Sinceeven good quality images have some noise, it makes senseto evaluate any denoising method in that way, without thetraditional “add noise and then remove it” trick. We shalllist formulas permitting to compute and analyze the method
noise for several classical local smoothing ﬁlters: the Gaussian ﬁltering [7], the anisotropic ﬁltering [1, 11], the Total
Variation minimization [13] and the neighborhood ﬁltering[16]. The formal analysis of the method noise for the frequency domain ﬁlters fall out of the scope of this paper.
These method noises can also be computed but their interpretation depends on the particular choice of the wavelet
basis.

2.1. The Gaussian ﬁltering
The image isotropic linear ﬁltering boils down to the
convolution of the image by a linear symmetric kernel. The
paradigm of such kernels is of course the gaussian kernel
|x|2

− 4h2
1
. In that case, Gh has standard
x → Gh (x) = (4πh2) edeviation h and it is easily seen that

Theorem 1 (Gabor 1960) The image method noise of the
convolution with a gaussian kernel Gh is
u − Gh ∗ u = −h2 ∆u + o(h2 ),for h small enough.
The gaussian method noise is zero in harmonic parts of
the image and very large near edges or texture, where theLaplacian cannot be small. As a consequence, the Gaussianconvolution is optimal in ﬂat parts of the image but edges
and texture are blurred.

2.2. The anisotropic ﬁltering
The anisotropic ﬁlter (AF ) attempts to avoid the blurring
effect of the Gaussian by convolving the image u at x onlyin the direction orthogonal to Du(x). The idea of such ﬁlter
goes back to Perona and Malik [11]. It is deﬁned byDu(x)⊥)dt,
AFh u(x) = Gh (t)u(x + t|Du(x)|
for x such that Du(x) = 0 and where (x, y)⊥ = (−y, x)
and Gh is the one-dimensional 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 The image method noise of an anisotropic ﬁlter
AFh is1u(x) − AFh u(x) = − h2 |Du|curv(u)(x) + o(h2 ),2where the relation holds when Du(x) = 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. This method noise is zero wherever u behaves locally likea straight line and large in curved edges or texture (where
the curvature and gradient operators take high values). As aconsequence, the straight edges are well restored while ﬂatand textured regions are degraded.

2.3. The Total Variation minimization
The Total Variation minimization was introduced by
Rudin, Osher and Fatemi [13]. Given a noisy image v(x),these authors proposed to recover the original image u(x)as the solution of the minimization problem
TVFλ (v) = arg min T V (u) + λ |v(x) − u(x)|2 dx
uwhere T V (u) denotes the total variation of u and λ is agiven Lagrange multiplier. The minimum of the above minimization problem exists and is unique. The parameter λis related to the noise statistics and controls the degree ofﬁltering of the obtained solution.
Theorem 3 The method noise of the Total Variation minimization isu(x) − TVFλ (u)(x) = −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.

2.4. The neighborhood ﬁltering
We call neighborhood ﬁlter any ﬁlter which restores apixel by taking an average of the values of neighboring pixels with a similar grey level value. Yaroslavsky (1985) [16]averages pixels with a similar grey level value and belonging to the spatial neighborhood Bρ (x),|u(y)−u(x)|21h2
u(y)e−dy, (2)
YNFh,ρ u(x) =C(x) Bρ (x)|u(y)−u(x)|2h2
dy is the normalwhere x ∈ Ω, C(x) = Bρ (x) e−
ization factor and h is a ﬁltering parameter.
The Yaroslavsky ﬁlter is less known than more recent
versions, namely the SUSAN ﬁlter (1995) [14] and the Bilateral ﬁlter (1998) [15]. Both algorithms, instead of considering a ﬁxed spatial neighborhood Bρ (x), weigh the distanceto the reference pixel x,|y−x|2|u(y)−u(x)|21−h2
SNFh,ρ u(x) =u(y)e ρ2 e−dy,
C(x) Ω(3)2
 − |y−x||u(y)−u(x)|2−2h2
dy is the normalizawhere C(x) = Ω e ρ etion factor and ρ is now a spatial ﬁltering parameter. In practice, there is no difference between YNFh,ρ and SNFh,ρ . If
the grey level difference between two regions is larger thanh, both algorithms compute averages of pixels belonging to
the same region as the reference pixel. Thus, the algorithm
does not blur the edges, which is its main scope. In theexperimentation section we only compare the Yaroslavsky
neighborhood ﬁlter.
The problem with these ﬁlters is that comparing only
grey level values in a single pixel is not so robust whenthese values are noisy. Neighborhood ﬁlters also create artiﬁcial shocks which can be justiﬁed by the computation ofits method noise, see [3].

3. NL-means algorithm
Given a discrete noisy image v = {v(i) | i ∈ I}, theestimated value N L[v](i), for a pixel i, is computed as aweighted average of all the pixels in the image,
N L[v](i) =w(i, j)v(j),j∈Iwhere the family of weights {w(i, j)}j depend on the similarity between the pixels i and j, and satisfy the usual
conditions 0 ≤ w(i, j) ≤ 1 and j w(i, j) = 1.
The similarity between two pixels i and j depends on
the similarity of the intensity gray level vectors v(Ni ) andv(Nj ), where Nk denotes a square neighborhood of ﬁxed
size and centered at a pixel k. This similarity is measured

Figure 1. Scheme of NL-means strategy. Similar
pixel neighborhoods give a large weight, w(p,q1)
and w(p,q2), while much different neighborhoods
give a small weight w(p,q3).

as a decreasing function of the weighted Euclidean distance,v(Ni ) − v(Nj )22,a , where a > 0 is the standard deviationof the Gaussian kernel. The application of the Euclidean
distance to the noisy neighborhoods raises the following
equality
E||v(Ni ) − v(Nj )||22,a = ||u(Ni ) − u(Nj )||22,a + 2σ 2 .
This equality shows the robustness of the algorithm since in
expectation the Euclidean distance conserves the order ofsimilarity between pixels.
The pixels with a similar grey level neighborhood tov(Ni ) have larger weights in the average, see Figure 1.
These weights are deﬁned as,2

w(i, j) =j )||2,a
1 − ||v(Ni )−v(N
h2
e
,
Z(i)

where Z(i) is the normalizing constantZ(i) =

 − ||v(Ni )−v(Ni )||22,ah2
ej

and the parameter h acts as a degree of ﬁltering. It controls
the decay of the exponential function and therefore the decay of the weights as a function of the Euclidean distances.
The NL-means not only compares the grey level in a single point but the the geometrical conﬁguration in a whole
neighborhood. This fact allows a more robust comparison
than neighborhood ﬁlters. Figure 1 illustrates this fact, thepixel q3 has the same grey level value of pixel p, but theneighborhoods are much different and therefore the weightw(p, q3) is nearly zero.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2. Display of the NL-means weight distribution used to estimate the central pixel of every image. The
weights go from 1(white) to zero(black).

(i) E[V (i) | Xi = x] = E[U (i) | Xi = x] for all i ∈ I
and x ∈ Rp .

4. NL-means consistency
Under stationarity assumptions, for a pixel i, the NLmeans algorithm converges to the conditional expectationof i once observed a neighborhood of it. In this case, thestationarity conditions amount to say that as the size of theimage grows we can ﬁnd many similar patches for all thedetails of the image.
Let V be a random ﬁeld and suppose that the noisy image v is a realization of V . Let Z denote the sequence ofrandom variables Zi = {Yi , Xi } where Yi = V (i) is realvalued and Xi = V (Ni \{i}) is Rp valued. The NL-means
is an estimator of the conditional expectation r(i) = E[Yi |Xi = v(Ni \{i})].
Theorem 4 (Conditional expectation theorem) Let Z ={V (i), V (Ni \{i})} for i = 1, 2, . . . be a strictly stationary and mixing process.
Let N Ln denote the
NL-means algorithm applied to the sequence Zn ={V (i), V (Ni \{i})}ni=1 . Then,|N Ln (j) − r(j)| → 0

a.s

for j ∈ {1, . . . , n}.
The full statement of the hypothesis of the theorem and its
proof can be found in a more general framework in [12].
This theorem tells us that the NL-means algorithm corrects
the noisy image rather than trying to separate the noise (oscillatory) from the true image (smooth).
In the case that an additive white noise model is assumed,
the next result shows that the conditional expectation is thefunction of V (Ni \{i}) that minimizes the mean square error with the true image u.
Theorem 5 Let V, U, N be random ﬁelds on I such that
V = U + N , where N is a signal independent white noise.
Then, the following statements are hold.

(ii) The expected random variable E[U (i) | V (Ni \{i})]is the function of V (Ni \{i}) that minimizes the meansquare error
min E[U (i) − g(V (Ni \{i}))]2g

Similar optimality theoretical results have been obtained in[9] and presented for the denoising of binary images. Theoretical links between the two algorithms will be explored ina future work.

5. Discussion and experimentation
In this section we compare the local smoothing ﬁlters
and the NL-means algorithm under three well deﬁned criteria: the method noise, the visual quality of the restored
image and the mean square error, that is, the Euclidean difference between the restored and true images.
For computational purposes of the NL-means algorithm,
we can restrict the search of similar windows in a larger”search window” of size S × S pixels. In all the experimentation we have ﬁxed a search window of 21×21 pixels and asimilarity square neighborhood Ni of 7 × 7 pixels. If N 2 isthe number of pixels of the image, then the ﬁnal complexityof the algorithm is about 49 × 441 × N 2 .
The 7 × 7 similarity window has shown to be large
enough to be robust to noise and small enough to take careof details and ﬁne structure. The ﬁltering parameter h has
been ﬁxed to 10 ∗ σ when a noise of standard deviation σ
is added. Due to the fast decay of the exponential kernel,
large Euclidean distances lead to nearly zero weights acting
as an automatic threshold, see Fig. 2.

Figure 3. Denoising experience on a natural texture. From left to right: noisy image (standard deviation 35),
Gauss ﬁltering, anisotropic ﬁltering, Total variation, Neighborhood ﬁltering and NL-means algorithm.

Figure 4. Method noise experience on a natural image. Displaying of the image difference u−Dh (u). From left toright and from top to bottom: original image, Gauss ﬁltering, anisotropic ﬁltering, Total variation minimization,
Neighborhood ﬁltering and NL-means algorithm. The visual experiments corroborate the formulas of section 2.

In section 2 we have computed explicitly the method
noise of the local smoothing ﬁlters. These formulas are corroborated by the visual experiments of Figure 4. This ﬁgure displays the method noise for the standard image Lena,
that is, the difference u − Dh (u), where the parameter his been ﬁxed in order to remove a noise of standard deviation 2.5. The method noise helps us to understand the performance and limitations of the denoising algorithms, since
removed details or texture have a large method noise. We
see in Figure 4 that the NL-means method noise does notpresent any noticeable geometrical structures. Figure 2 explains this property since it shows how the NL-means algorithm chooses a weighting conﬁguration adapted to the
local and non local geometry of the image.
The human eye is the only one able to decide if thequality of the image has been improved by the denoising
method. We display some denoising experiences comparing the NL-means algorithm with local smoothing ﬁlters.
All experiments have been simulated by adding a gaussian
white noise of standard deviation σ to the true image. The
objective is to compare the visual quality of the restored
images, the non presence of artifacts and the correct recon-

struction of edges, texture and details.
Due to the nature of the algorithm, the most favorable
case for the NL-means is the textured or periodic case. In
this situation, for every pixel i, we can ﬁnd a large set ofsamples with a very similar conﬁguration. See Figure 2 e)
for an example of the weight distribution of the NL-means
algorithm for a periodic image. Figure 3 compares the performance of the NL-means and local smoothing ﬁlters for anatural texture.
Natural images also have enough redundancy to be restored by NL-means. Flat zones present a huge number ofsimilar conﬁgurations lying inside the same object, see Figure 2 (a). Straight or curved edges have a complete lineof pixels with similar conﬁgurations, see Figure 2 (b) and(c). In addition, natural images allow us to ﬁnd many similar conﬁgurations in far away pixels, as Figure 2 (f) shows.
Figure 5 shows an experiment on a natural image. This experience must be compared with Figure 4, where we display
the method noise of the original image. The blurred or degraded structures of the restored images coincide with thenoticeable structures of its method noise.
Finally Table 1 displays the mean square error for theFigure 5. Denoising experience on a natural image. From left to right and from top to bottom: noisy image(standard deviation 20), Gauss ﬁltering, anisotropic ﬁltering, Total variation, Neighborhood ﬁltering and NLmeans algorithm. The removed details must be compared with the method noise experience, Figure 4.

Image
Lena
Baboon

GF
120
507

AF
114
418

TVF
110
365

YNF
129
381

NL
68
292

Table 1. Mean square error table. A smaller meansquare error indicates that the estimate is closer
to the original image.

denoising experiments given in the paper. This numericalmeasurement is the most objective one, since it does notrely on any visual interpretation. However, this error is notcomputable in a real problem and a small mean square error
does not assure a high visual quality. So all above discussed
criteria seem necessary to compare the performance of algorithms.

References[1] L. Alvarez, P.-L. Lions, and J.-M. Morel. Image selective
smoothing and edge detection by nonlinear diffusion (ii).
Journal of numerical analysis, 29:845–866, 1992.
[2] A. Buades, B. Coll, and J. Morel. On image denoising methods. Technical Report 2004-15, CMLA, 2004.
[3] A. Buades, B. Coll, and J. Morel. Neighborhood ﬁlters andpde’s. Technical Report 2005-04, CMLA, 2005.
[4] R. Coifman and D. Donoho. Wavelets and Statistics, chapter
Translation-invariant de-noising, pages 125–150. Springer
Verlag, 1995.
[5] D. Donoho. De-noising by soft-thresholding. IEEE Transactions on Information Theory, 41:613–627, 1995.

[6] A. Efros and T. Leung. Texture synthesis by non parametricsampling. In Proc. Int. Conf. Computer Vision, volume 2,
pages 1033–1038, 1999.
[7] M. Lindenbaum, M. Fischer, and A. Bruckstein. On gabor
contribution to image enhancement. Pattern Recognition,
27:1–8, 1994.
[8] Y. Meyer. Oscillating Patterns in Image Processing and
Nonlinear Evolution Equations, volume 22. AMS University Lecture Series, 2002.
[9] E. Ordentlich, G. Seroussi, M. W. S. Verdú, and T. Weissman. A discrete universal denoiser and its application to binary images. In Proc. IEEE ICIP, volume 1, pages 117–120,
2003.
[10] S. Osher, A. Sole, and L. Vese. Image decomposition andrestoration using total variation minimization and the h−1
norm. Multiscale Modeling and Simulation, 1(3):349–370,
2003.
[11] P. Perona and J. Malik. Scale space and edge detection using
anisotropic diffusion. IEEE Trans. Patt. Anal. Mach. Intell.,
12:629–639, 1990.
[12] G. Roussas. Nonparametric regression estimation under
mixing conditions. Stochastic processes and their applications, 36:107–116, 1990.
[13] L. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation
based noise removal algorithms. Physica D, 60:259–268,
1992.
[14] S. Smith and J. Brady. Susan - a new approach to low level
image processing. International Journal of Computer Vision, 23(1):45–78, 1997.
[15] C. Tomasi and R. Manduchi. Bilateral ﬁltering for gray andcolor images. In Proceedings of the Sixth Internatinal Conference on Computer Vision, pages 839–846, 1998.
[16] L. Yaroslavsky. Digital Picture Processing - An Introduction. Springer Verlag, 1985.