POISSON NL MEANS: UNSUPERVISED NON LOCAL MEANS FOR POISSON NOISE
Charles-Alban Deledalle, Florence Tupin
Institut Telecom, Telecom ParisTech
CNRS LTCI
Paris, France
Lo
¨
ıc Denis
Observatoire de Lyon, CNRS CRAL
UCBL, ENS de Lyon, Universit´e de Lyon
Lyon, France
ABSTRACT
An extension of the non local (NL) means is proposed for
images damaged by Poisson noise. The proposed method is
guided by the noisy image and a pre-filtered image and is
adapted to the statistics of Poisson noise. The influence of
both images can be tuned using two filtering parameters. We
propose an automatic setting to select these parameters based
on the minimization of the estimated risk (mean square error).
This selection uses an estimator of the MSE for NL means
with Poisson noise and Newton’s method to find the optimal
parameters in few iterations.
Index Terms— Non local means, Poisson noise, mean
square error, SURE, PURE, Newton’s method
1. INTRODUCTION
Poisson noise appears in low-light conditions when the num-
ber of collected photons is small, such as in optical mi-
croscopy or astronomy. Poisson noise is signal-dependent
and, thus, needs to adapt the usual denoising approaches.
NL means have been proposed by Buades et al. in [2] to
denoise images damaged by additive white Gaussian noise.
While local filters lead to biases and resolution loss, NL tech-
niques are known to efficiently reduce noise and preserve
structures. Instead of combining neighboring pixels, the NL
means average similar pixels. Let k
s
be the observed noisy
value at site s and λ
s
its underlying noise-free value. NL
means define the estimate
ˆ
λ
s
as a weighted average:
ˆ
λ
s
=
t
w
s,t
k
t
t
w
s,t
(1)
where t is a pixel index and w
s,t
is a data-driven weight de-
pending on the similarity between pixels with indexes s and
t. In practice, the pixels t are located in a search window
centered on s. For robustness reasons, pixel similarity is
evaluated by comparing surrounding patches around s and
t. Patch-similarity is classically defined by the Euclidean
Thanks to Vincent Duval for its comments and for the reference [1].
distance, leading to the following weight expression:
w
s,t
=exp
−
b
(k
s+b
− k
t+b
)
2
α
(2)
where
s+b
and
t+b
denote the b-th pixels in the patches B
s
and
B
t
centered on s and t,andα is a filtering parameter.
In case of low signal-to-noise ratio images, it has been
shown that the performances of the NL means can be im-
proved by refining the weights using a pre-estimate
ˆ
θ of the
noise-free image [3, 4, 5, 6]. A general expression of refined
NL means is:
ˆ
λ
s
=
t
w
s,t
k
t
t
w
s,t
(3)
with w
s,t
=exp
−
F
s,t
α
−
G
s,t
β
,
F
s,t
=
b
f(k
s+b
,k
t+b
)
and G
s,t
=
b
g
ˆ
θ
s+b
,
ˆ
θ
t+b
where α and β are filtering parameters, and f and g are two
similarity criteria suitable respectively to compare noisy data
and pre-estimated data. A typical choice for f and g is the
squared difference: f (x, y)=g(x, y)=(x−y)
2
. The choice
of the filtering parameters for α and/or β is a critical task al-
ready explored in [3, 6] and [7]. According to our knowledge,
there is no known method to jointly set α and β.
2. PATCH-SIMILARITIES UNDER POISSON NOISE
Let k be a noisy observation following a Poisson distribution
with parameters described by a noise-free value λ:
p(k|λ)=
λ
k
e
−λ
k!
. (4)
We take inspiration from the probabilistic approach of [6, 8]
to extend the refined NL means in (3) to the Poisson noise
degradation model. The squared difference classically used
801978-1-4244-7993-1/10/$26.00 ©2010 IEEE ICIP 2010
Proceedings of 2010 IEEE 17th International Conference on Image Processing September 26-29, 2010, Hong Kong
for f can be replaced by:
f
L
(k
1
,k
2
)=− log
max
λ
p(k
1
|λ
1
= λ)p(k
2
|λ
2
= λ)
max
λ
p(k
1
|λ
1
= λ)max
λ
p(k
2
|λ
2
= λ)
= k
1
log k
1
+ k
2
log k
2
− (k
1
+ k
2
)log
k
1
+ k
2
2
. (5)
This similarity evaluates the likelihood ratio corresponding
to the hypothesis that k
1
and k
2
share an identical λ param-
eter against the hypothesis that their λ parameters are inde-
pendent. The squared difference generally used for g can be
replaced by the symmetric Kullback-Leibler divergence:
g
KL
ˆ
θ
1
,
ˆ
θ
2
= D
KL
ˆ
θ
1
ˆ
θ
2
=
ˆ
θ
1
−
ˆ
θ
2
log
ˆ
θ
1
ˆ
θ
2
. (6)
This criterion is a good candidate to define similarities be-
tween estimated values since it can be considered as a statis-
tical test of the hypothesis λ
1
= λ
2
.
The setting of the parameters α and β in the case of Pois-
son noise is maybe a more critical problem than in other de-
noising tasks. In [2, 3] and [6], the authors propose to define
the parameters according to the variance or the quantiles of
the similarity criteria. Unfortunately, in the case of Poisson
noise, these quantities depend on the unknown image λ since
the noise is signal-dependent. Van De Ville et al. propose
a risk minimization approach for Gaussian noise [7]. Their
method selects the parameters minimizing the risk (without
any specific assumption on the underlying image λ). This
kind of approach seems relevant in the case of Poisson noise.
3. AUTOMATIC SETTING OF PARAMETERS
B ASED ON RISK MINIMIZATION
The parameters of the denoising technique can be selected as
those that minimize the expected MSE:
E
1
N
λ −
ˆ
λ
2
=
1
N
s
λ
2
s
+ E
ˆ
λ
2
s
− 2E
λ
s
ˆ
λ
s
(7)
with N the image size and E[.] the expectation operator. The
MSE requires the knowledge of the noise-free image λ but
can still be estimated from the noisy image k alone. Since
the first term λ
2
s
in (7) is independent of
ˆ
λ
s
, it can be omit-
ted when minimizing the MSE with respect to the denoising
parameters. The Stein’s unbiased risk estimator (SURE) is an
estimator of the MSE under Gaussian noise [9]. It is based on
an estimator of E[λ
s
ˆ
λ
s
] which does not require λ.SUREhas
already been used successfully on images damaged by addi-
tive white Gaussian noise for wavelet filtering [10] and NL
means filtering [7]. The main result in [7] is that SURE for
NL means can be obtained in closed form. For Poisson noise,
we can use the result of Chen [1] to follow the same approach:
E [λ
s
h(k)
s
]=E
k
s
h(k)
s
(8)
with
k
t
=
k
t
if t = s
k
t
− 1 otherwise
where h(.)
s
denotes the estimated value on site s obtained by
the application of the estimator h on the given noisy image.
By injecting (8) in (7), we obtain the Poisson unbiased risk
estimator (PURE) as introduced in [11] for wavelet denoising.
Considering now that the estimator h is the refined NL means
given in eq. (3), PURE can be expressed as follows:
R(
ˆ
λ)=
1
N
s
λ
2
s
+
ˆ
λ
2
s
− 2k
s
λ
s
. (9)
The value
λ
s
refers to the denoised value obtained by the ap-
plication of the NL means on the noisy image
k,i.e:
λ
s
=
t
w
s,t
k
t
t
w
s,t
(10)
with
w
s,t
=exp
−
F
s,t
α
−
G
s,t
β
and
F
s,t
=
b
f(k
s+b
, k
t+b
).
According to (7), (8) and (9), it is straightforward to show that
R(
ˆ
λ) is unbiased. However, note that (10) holds by assuming
that G
s,t
(i.e the pre-estimate
ˆ
θ) does not depend on the noise
component of k. To satisfy this assumption, the noise vari-
ance in
ˆ
θ has to be reduced significantly.
In terms of time complexity, we note as in [7] that the
computation time is unchanged since the computation of
PURE can be incorporated within the core of the NL means.
Moreover, the scan of the patches of
k can be avoided thanks
to the following relation:
F
s,t
= F
s,t
+
⎧
⎪
⎪
⎨
⎪
⎪
⎩
f(
k
s
, k
s
) − f(k
s
,k
s
), if s = t,
f(
k
s
,k
t
) − f(k
s
,k
t
)
+f (k
2s−t
, k
s
) − f(k
2s−t
,k
s
), if s ∈B
t
,
f(
k
s
,k
t
) − f(k
s
,k
t
), otherwise.
Selecting parameters that minimize PURE gives param-
eters close to that minimizing the MSE. In the case of the
classical NL means, the authors of [7, 11] compute the op-
timal parameters by exhaustive search. Optimization tech-
niques can be applied to reach the optimal parameters in few
iterations. We suggest to optimize R(
ˆ
λ) for the refined NL
means by using Newton’s method on the joint filtering pa-
rameters α and β. Newton’s method iteratively refines α and
β with the update:
α
(n+1)
β
(n+1)
=
α
(n)
β
(n)
− H
−1
∇ (11)
with H
−1
∇ =
∂
2
R(
ˆ
λ
(n)
)
∂α
2
∂
2
R(
ˆ
λ
(n)
)
∂α∂β
∂
2
R(
ˆ
λ
(n)
)
∂β∂α
∂
2
R(
ˆ
λ
(n)
)
∂β
2
−1
∂R(
ˆ
λ
(n)
)
∂α
∂R(
ˆ
λ
(n)
)
∂β
.
where n is the current iteration index. To perform the opti-
mization procedure in (11), the first and second order differ-
802
20
25
30
35
40
45
Risk
MSE
PURE
−10
−8
−6
−4
−2
0
2
Variations of the risk
∂ MSE / ∂ α
∂ PURE / ∂ α
Zero
5 10 15 20
0
2
4
6
α
Slope variations of the risk
∂
2
MSE / ∂ α
2
∂
2
PURE / ∂ α
2
Zero
MSE
PURE
∂ MSE / ∂ β
∂ PURE / ∂ β
Zero
5 10 15 20
β
∂
2
MSE / ∂ β
2
∂
2
PURE / ∂ β
2
Zero
Fig. 1. The risk (MSE and PURE) and their two first order
variations (from top to bottom) with respect to the parameters
α (left) and β (right).
entials are required. Their expressions are given by substitut-
ing x and y by α or β in the following equations:
∂R(
ˆ
λ)
∂x
=
2
N
s
ˆ
λ
s
∂
ˆ
λ
s
∂x
−
2
N
s
k
s
∂λ
s
∂x
,
∂
2
R(
ˆ
λ)
∂x∂y
=
2
N
s
ˆ
λ
s
∂
2
ˆ
λ
s
∂x∂y
+
2
N
s
∂
ˆ
λ
s
∂x
∂
ˆ
λ
s
∂y
−
2
N
s
k
s
∂
2
λ
s
∂x∂y
,
∂
ˆ
λ
s
∂x
=
X
s,t
w
s,t
(k
t
−
ˆ
λ
s
)
x
2
w
s,t
,
∂
2
ˆ
λ
s
∂x
2
, =
X
2
s,t
w
s,t
(k
t
−
ˆ
λ
s
)
x
4
w
s,t
− 2
∂
ˆ
λ
s
∂x
(X
s,t
+ x)w
s,t
x
2
w
s,t
,
∂
2
ˆ
λ
s
∂x∂y
=
X
s,t
Y
s,t
w
s,t
(k
t
−
ˆ
λ
s
)
x
2
y
2
w
s,t
−
∂
ˆ
λ
s
∂x
Y
s,t
w
s,t
y
2
w
s,t
−
∂
ˆ
λ
s
∂y
X
s,t
w
s,t
x
2
w
s,t
where X = F (resp. Y = F )whenx = α (resp. y = α)
and X = G (resp. Y = G)whenx = β (resp. y = β). The
differentials for
λ are the same with respect to k, w and F .
Newton’s method finds in few iterations the best trade-
off between the information brought by the noisy image and
the pre-estimated image to define the weights. For instance,
Peppers (256 × 256)
Noisy 3.14 13.14 17.91 23.92
MA filter 19.20 20.93 21.11 21.16
PURE-LET [11]
19.33 24.29 27.27 30.79
NL means [2]
18.12 23.33 26.98 30.64
Refined NL means
19.84 24.48 27.56 30.93
Poisson NL means
19.90 25.32 28.07 31.06
α
opt
(209) (13.6) (10.05) (9.21)
β
opt
(0.72) (1.31) (2.76) (7.64)
#iterations
(13.5) (8.02) (7.03) (6.90)
Cameraman (256 × 256)
Noisy 3.28 13.27 18.03 24.05
MA filter 18.71 20.15 20.29 20.33
PURE-LET [11]
19.67 24.32 26.87 30.36
NL means [2]
18.17 23.53 26.77 29.39
Refined NL means
19.80 24.53 27.18 29.44
Poisson NL means
19.89 25.07 27.42 29.47
α
opt
(62.1) (9.48) (8.81) (7.34)
β
opt
(0.51) (1.19) (3.57) (16.19)
#iterations
(11.0) (6.80) (7.60) (11.3)
Tab le 1 . PSNR values averaged over ten realisations using
different methods on images damaged by Poisson noise with
different levels of degradation. The averaged optimal param-
eters and the averaged number of iterations of the proposed
Poisson NL means are given.
β will get a high value when the pre-estimated image has a
poor quality, resulting in weights determined only from the
noisy image. Reciprocally, α will get a high value when the
pre-estimated image has a high quality: the weights will be
determined only from the well pre-estimated image.
4. EXPERIMENTS AND RESULTS
The proposed extension of the NL means (Poisson NL means)
is applied with a search window of size 21 × 21 and patches
of size 7 × 7. Newton’s method is performed until PURE
does not change between two successive iterations. The pre-
estimated image is obtained by a moving average (MA) filter
with a 13×13 disk kernel. Using the optimization of [12], the
computational time is of about 10s per iteration on a 256×256
image and C implementation on an Intel Core 2 Duo 64-bit
CPU 3.00GHz.
Figure 1 shows the risk and its two first order differentials
with respect to α and β. These curves have been computed
by applying the proposed method on a 150 × 150 noisy im-
age for different values of the parameters. The MSE and its
differentials have been computed from the noise-free image
and finite differences. PURE and its differentials have been
evaluated using the expressions given in Section 3.
Table 1 gives the peak signal-to-noise ratio (PSNR) val-
ues obtained by different denoising methods averaged over
803
(a) (b) (c)
Fig. 2. (a) Original images damaged by Poisson noise, de-
noised images obtained by (b) PURE-LET [11] and (c) the
proposed Poisson NL means.
ten realizations of two 256 × 256 reference images damaged
by synthetic Poisson noise with different degradation levels.
The MA filter is applied with a 13 × 13 disk kernel. We
compare our method with the PURE based wavelet approach
(PURE-LET) proposed in [11]. Three versions of NL means
are applied. NL means denotes the classical one i.e f (x, y)=
(x − y)
2
and β = ∞. We call the refined NL means when
f(x, y)=g(x, y)=(x−y)
2
. Poisson NL means denotes our
proposed method with f and g defined as in Section 2. For all
NL means versions, the optimal parameters are obtained by
PURE minimization using Newton’s method. The table gives
the optimal parameters α
opt
and β
opt
and the number of iter-
ations for the Poisson NL means. The refined and the Poisson
NL means both use the pre-estimated image obtained by the
MA filter. Poisson NL means provides generally better per-
formance with on average 6 to 14 iterations. The parameters
behave as predicted with respect to the relative qualities of the
noisy image and the pre-estimated image.
Figure 2 presents visual results
1
obtained by the PURE-
LET approach
2
and Poisson NL means on two images. The
first image is an image degraded by synthetic Poisson noise
and the second one is an image
3
of a mitochondrion sensed
in low-light conditions by confocal fluorescence microscopy
[13]. Even if both filters seem to preserve well the resolu-
tion while reducing the noise, the Poisson NL means seem to
provide a more regular result with fewer processing artifacts.
1
More results available at http://perso.telecom-paristech.
fr/
˜
deledall/poisson_nlmeans.php
2
we are grateful to F. Luisier for providing the results of PURE-LET
3
image courtesy of Y. Tourneur
5. CONCLUSION
Inspired by the methodology of [6], an extension of the NL
means has been proposed for images damaged by Poisson
noise. It is based on probabilistic similarities to compare
noisy patches and patches of a pre-estimated image. A risk
estimator for NL means, based on the idea of [7], has been
derived for Poisson noise. This risk estimator is used in an
optimization method to automatically select the filtering pa-
rameters in few iterations. Numerical results as well as visual
results support the efficiency of the proposed method.
6. REFERENCES
[1] L.H.Y. Chen, “Poisson approximation for dependent trials,”
The Annals of Probability, vol. 3, no. 3, pp. 534–545, 1975.
[2] A. Buades, B. Coll, and J.M. Morel, “A Non-Local Algorithm
for Image Denoising,” in Pr oc. IEEE Computer Society Conf.
CVPR, 2005, vol. 2, pp. 60–65.
[3] C. Kervrann and J. Boulanger, “Optimal spatial adaptation for
patch-based image denoising,” IEEE Transactions on Image
Processing, vol. 15, no. 10, pp. 2866–2878, 2006.
[4] T. Brox, O. Kleinschmidt, and D. Cremers, “Efficient Nonlocal
Means for Denoising of Textural Patterns,” IEEE Transactions
on Image Processing, 2007.
[5] B. Goossens, H. Luong, A. Piˇzurica, and W. Philips, “An im-
proved non-local denoising algorithm,” in Proc. Int. Workshop
on LNLA in Image Processing, Lausanne, Switzerland, 2008.
[6] C.A. Deledalle, L. Denis, and F. Tupin, “Iterative Weighted
Maximum Likelihood Denoising with Probabilistic Patch-
Based Weights,” IEEE Transactions on Image Processing,vol.
18, no. 12, pp. 2661–2672, Dec. 2009.
[7] D. Van De Ville and M. Kocher, “SURE-Based Non-Local
Means,” IEEE Signal Processing Letters, vol. 16, no. 11, pp.
973–976, November 2009.
[8] F. Alter, Y. Matsushita, and X. Tang, “An intensity similarity
measure in low-light conditions,” Lecture Notes in Computer
Science, vol. 3954, pp. 267, 2006.
[9] C.M. Stein, “Estimation of the mean of a multivariate normal
distribution,” The Annals of Statistics, pp. 1135–1151, 1981.
[10] D.L. Donoho and I.M. Johnstone, “Adapting to Unknown
Smoothness Via Wavelet Shrinkage.,” Journal of the ameri-
can statistical association, vol. 90, no. 432, 1995.
[11] F. Luisier, C. Vonesch, T. Blu, and M. Unser, “Fast interscale
wavelet denoising of Poisson-corrupted images,” Signal Pro-
cessing, vol. 90, no. 2, pp. 415–427, 2010.
[12] J. Darbon, A. Cunha, T.F. Chan, S. Osher, and G.J. Jensen,
“Fast nonlocal filtering applied to electron cryomicroscopy,”
in Pr oc. IEEE Int. Symp. on Biomedical Imaging: From Nano
to Macro, 2008, pp. 1331–1334.
[13] S. Pelloux, J. Robillard, R. Ferrera, A. Bilbaut, C. Ojeda,
V. Saks, M. Ovize, and Y. Tourneur, “Non-beating HL-1 cells
for confocal microscopy: application to mitochondrial func-
tions during cardiac preconditioning,” Progress in biophysics
and molecular biology, vol. 90, no. 1-3, pp. 270–298, 2006.
804
![Fig. 2. (a) Original images damaged by Poisson noise, denoised images obtained by (b) PURE-LET [11] and (c) the proposed Poisson NL means.](/figures/fig-2-a-original-images-damaged-by-poisson-noise-denoised-yugdfmym.webp)