arXiv:1512.00717v1 [] 2 Dec 2015

MMSE Estimation for Poisson Noise Removal in
Images
Stanislav Pyatykh* and Jürgen Hesser

∗

December 3, 2015

Abstract
Poisson noise suppression is an important preprocessing step in several
applications, such as medical imaging, microscopy, and astronomical imaging. In this work, we propose a novel patch-wise Poisson noise removal
strategy, in which the MMSE estimator is utilized in order to produce
the denoising result for each image patch. Fast and accurate computation
of the MMSE estimator is carried out using k-d tree search followed by
search in the K-nearest neighbor graph. Our experiments show that the
proposed method is the preferable choice for low signal-to-noise ratios.

1

Introduction

Images corrupted with Poisson noise appear in many applications, such as medical imaging, fluorescence microscopy, and astronomical imaging. There, the
noise removal quality plays an essential role for the further image processing
steps. For Poisson noise, the noise variance equals the true image intensity,
which means that noise is signal-dependent and the signal-to-noise ratio equals
the square root of the image intensity. For this reason, the methods developed
for additive white Gaussian noise are not applicable for Poisson noise directly.
Many of currently existing approaches reduce the Poisson noise removal to
the Gaussian noise removal. This is an attractive strategy since the state of the
art algorithms for Gaussian noise suppression are very efficient: it was shown in
[9, 10] that they have almost achieved theoretical performance boundary. The
reduction of the denoising problem for Poisson noise to that for Gaussian noise
is based on the variance stabilization [2, 11, 15]. First, a pixel-wise variancestabilizing transformation (such as the Anscombe root transformation [1]) is
applied in order to transform the input image with Poisson noise into the image
with signal-independent Gaussian noise. Second, a method for Gaussian noise
removal is utilized (e.g. [4, 5, 3]). Finally, an inverse transformation is applied.
∗ The
authors are with the University Medical Center Mannheim, Heidelberg University, Theodor-Kutzer-Ufer 1-3, D-68167 Mannheim, Germany (e-mail:
stanislav.pyatykh@medma.uni-heidelberg.de; juergen.hesser@medma.uni-heidelberg.de).

1

In [11] and [6], it was pointed out that application of the algebraic inverse
introduces a bias and the exact unbiased inverse was developed.
Nevertheless, the variance stabilization becomes inaccurate as the image
intensity decreases, which limits the applicability of the above mentioned approaches. This fact led to development of denoising methods, which work with
Poisson-distributed data directly. Among them, [14] and [7] are the most efficient. First, these algorithms group similar image patches into clusters. Next,
patches of each cluster are represented using a small number of dictionary elements. Due to the projection onto the subspace defined by these dictionary
elements, the noise components orthogonal to this subspace are canceled out.
Therefore, the dimensionality of this subspace, i.e. the number of the dictionary
elements used to represent each patch, should be made as small as possible. In
[14], the dictionary elements are computed using Poisson principal component
analysis of each patch cluster, which is based on the Bregman divergence for
Poisson distribution. In [7], the dictionary elements are selected using a greedy
algorithm, which searches for the maximum of the likelihood of the denoised
patch given the input noisy patch.
In this paper, we propose a different technique for Poisson noise removal.
Instead of projecting image patches onto a subspace defined by a dictionary,
the denoised patch is computed as the minimum mean square error (MMSE)
estimator of the true patch given the noisy patch. If all possible noise-free image
patches were available for the algorithm as prior information, such approach
would provide the optimal estimate of the true patch. However, we show that
it outperforms existing methods also with computationally tractable amount of
prior information.
The MMSE estimation was considered in the image denoising context in
[9, 10], where it was used to compute the theoretical performance bounds of
denoising algorithms. Since the direct computation of the MMSE estimator is
very time-consuming, practical algorithms based on the MMSE estimation were
not developed in these works. An optimized denoising algorithm based on the
MMSE estimation was presented in [8] for the case of Gaussian noise, however,
its speed drops as signal-to-noise ratio decreases. In this paper, we propose a
fast and efficient method to compute the MMSE estimator for an image patch.
This leads to a practical image denoising algorithm, which is able to process
images in a reasonable time and outperforms existing approaches in terms of
the output image quality.

2

MMSE estimation for Image Patches

√
√
Consider an unknown noise-free image patch of size d × d rearranged into
vector x ∈ Rd and the corresponding patch of the input noisy image represented
as vector y ∈ (N ∪ 0)d . We model patch x using random vector X ∈ Rd with
probability density function fX (·), whereas patch y is modeled by random vector
Y ∈ (N ∪ 0)d with probability mass function fY (·). Assume that we know the
likelihood P (Yi = y|Xi = x), where Xi and Yi are the i-th elements of X and
2

Y respectively. In the case of Poisson noise,

y −x

x e /y! x > 0
P (Yi = y|Xi = x) = 0
x = 0, y > 0


1
x = 0, y = 0

(1)

Let us denote the peak intensity of the noise-free image by Imax . The problem to denoise patch y is typically formulated in terms of the peak signal-to-noise
ratio (PSNR): we look for estimate x̂ ∈ Rd , which maximizes
PSNR = 10 log10

2
Imax
MSE

(2)

or, equivalently, which minimizes the mean square error (MSE):
MSE = E(kx̂ − Xk2 |Y = y)
Z
= kx̂ − uk2 fX (u|Y = y)du.
Such x̂ is the MMSE estimator [9], which is computed as
Z
x̂ = E(X|Y = y) = ufX (u|Y = y)du.

(3)

(4)

Using Bayes’ theorem and the properties of marginal distributions, the latter
expression can be rewritten as
Z
1
x̂ =
ufY (y|X = u)fX (u)du
fY (y)
R
ufY (y|X = u)fX (u)du
= R
(5)
fY (y|X = u)fX (u)du
where fY (y|X = u) is the probability of y given u, which can be computed
using (1):
fY (y|X = u) =

d
Y

P (Yi = yi |Xi = ui ).

(6)

i=1

Since the probability density function for image patches fX (·) cannot be
accurately represented in some analytic form, direct computation of the integrals
in (5) is not possible in practice. In order to approximate these integrals, one can
use a large number of natural image patches u(1) , . . . , u(NP ) , which are treated
as realizations of random vector X [9]:
PNP (k)
fY (y|X = u(k) )
k=1 u
x̂app = P
.
(7)
NP
(k) )
k=1 fY (y|X = u
This expression gives a computationally feasible approximation for the MMMSE
estimator, which can be utilized for denoising of patch y.
3

3

Efficient Patch Denoising

In order to approximate x̂ by x̂app accurately enough, one has to use about
NP = 108 natural image patches. Therefore, direct evaluation of the sums
in (7) is very computationally expensive and, therefore, inapplicable. In this
section, we present a fast and practical method to compute x̂app .

3.1

Offline step

Several data structures are constructed based on patches u(1) , . . . , u(NP ) at the
offline (pre-processing) step.
First, patches u(1) , . . . , u(NP ) are normalized since the intensity range of
images from which they are extracted does not correspond to the intensity range
of the input image. Specifically, each patch u(k) is divided by the average of all
intensities occurring in u(1) , . . . , u(NP ) , resulting in normalized patch u(k) :
I=

NP X
d
1 X
(k)
ui
dNP
i=1
k=1

(k)

u

(k)

=u

/I

k = 1, . . . , NP .

(8)

Next, we apply the k-means clustering to patch set u(1) , . . . , u(NP ) . For each
cluster j, it provides the cluster centroid c(j) and the number of patches assigned
to this cluster nj . Since each patch u(k) can be approximated by the centroid
of the cluster it is assigned to, (7) can be replaced with
PNC
x̂app =

(j)
(j)
j=1 c nj fY (y|X = c )
PNC
(j)
j=1 nj fY (y|X = c )

(9)

where NC is the number of clusters.
Then, centroids c(1) , . . . , c(NC ) are used to construct a set of NT randomized
k-d trees as described in [13]. For each node of a randomized k-d tree, the
splitting dimension is chosen randomly from several dimensions with the largest
variance, so that each k-d tree represents a different space partitioning.
Additionally, the K-nearest neighbor graph is constructed using centroids
c(1) , . . . , c(NC ) as nodes. For each centroid c(j) , we compute the neighbor list
containing the indexes of K centroids from c(1) , . . . , c(NC ) , which are the nearest
to c(j) in the Euclidean metric (excluding c(j) itself).

4

3.2

Online step

Algorithm 1
1: s ← d-dimensional zero vector
2: w ← 0
Pd
1
3: µy ← d
i=1 yi
4: Q ← empty priority queue
5: Mark all clusters as unprocessed
6: for each k-d tree T do
7:
{jm }L
y/µy
m=1 ← leaf of T containing

8:
ProcessClusters {jm }L
m=1
9: end for
10: while Q is not empty and w has not converged do
11:
Pop element j ∗ with the highest priority from Q
(j ∗ )
12:
{jm }K
m=1 ← neighbors of c

13:
ProcessClusters {jm }K
m=1
14: end while
15: x̂app ← s/w
16: return x̂app
Sub-Algorithm ProcessClusters {jm }):
1: for each unprocessed cluster jm do
2:
s ← s + njm fY (y|X = µy c(jm ) ) · µy c(jm )
3:
w ← w + njm fY (y|X = µy c(jm ) )
4:
Push jm into Q with priority fY (y|X = µy c(jm ) )
5:
Mark cluster jm as processed
6: end for
In order to compute estimator (9) efficiently, we utilize the fact that patches
c(j) similar to input patch y have the largest likelihood fY (y|X = c(j) ) and,
therefore, the largest contribution to the sums in (9). On the other hand,
likelihood fY (y|X = c(j) ) is almost zero for patches c(j) , which are significantly
different from y, so that these patches have almost no influence on the result.
Therefore, we start the calculation of (9) from patches c(j) , which are the closest
to y, and continue by processing patches farther away from y until convergence
of the result.
This procedure is presented in Algorithm 1, which takes noisy patch y as
the input and computes denoised patch x̂app . Vector s ∈ Rd represents the
numerator of (9), whereas w ∈ R stores the denominator. First, centroids in
the neighborhood of y are found by traversing each k-d tree (lines 6–9). The
traversal starts from the root and goes to either the left or the right child
depending on the location of y relative to the splitting plane until a leaf node
is reached. Centroids {c(jm ) }L
m=1 assigned to this node are then used to update
s and w and added to the priority queue. Next, centroids farther away from y
are found using the K-nearest neighbor graph (lines 10–14). At each iteration
of the loop, the centroid with the largest likelihood is removed from the priority
5

Figure 1: Test images. From left to right: Peppers, Bridge, Boat, Couple, Hill,
Mandril, Pirate.
queue and its previously unprocessed neighbors are examined. The loop repeats
until the convergence criterion is met (line 10), namely, until the relative change
of w during the last M iterations is less than 10−12 .
Note that centroids c(j) were constructed using normalized patches u(k) .
Therefore, noisy patch y should be normalized in order to be utilized for the
k-d tree search (line 7) and, vise versa, centroids c(j) should be multiplied
by the mean of y to compute s, w, and the likelihood (see the sub-algorithm
ProcessClusters).
In order to denoise the input image, we extract its overlapping patches in a
sliding window manner. Next, these patches are denoised independently from
each other using Algorithm 1 and returned to their original positions. Since
they are overlapped, there are several estimates for each pixel, which are then
averaged to produce the final denoised image.

4
4.1

Experiments
Setup

The standard test images of size 256 × 256 shown in Fig. 1 were utilized in the
experiments. Since the variance equals the expected value for Poisson noise, the
peak intensity of the original image defines the noise level in the noisy image.
On one hand, it was shown in [7] that the strategy referred to as binning is
preferred when the peak intensity is less than 1. According to this approach,
the input image is first down-sampled by aggregating neighboring pixels, which
results in a smaller image with a higher peak intensity. This image is then
processed by any Poisson denoising algorithm and upscaled to the original size
using interpolation. On the other hand, when the peak intensity is higher than
5, the methods based on the variance stabilization are preferable. Therefore, we
varied the peak intensity between 1 and 5 in our experiments.
We compared our approach with the SPDA [7] and the BM3D with the
exact unbiased inverse Anscombe transformation [11], which are currently the
best existing techniques. Their code is available online and was taken with
the default parameter settings. The C++ implementation of our algorithm,
which was utilized in the experiments, is available at http://physics.medma.uniheidelberg.de/cms/content/poisson-denoising. All overlapping patches from the
BSD500 image set [12] were taken as u(1) , . . . , u(NP ) , which resulted in a set of
NP ≈ 7 · 107 patches.
6

4.2

Choice of Parameters

√
√
Increase of the patch size d × d makes it possible to handle more complex
image structures but requires a larger number of patches NP to represent the
natural image patch distribution. In our experiments, which were conducted
with NP ≈ 7 · 107 , we used patches of size 14 × 14.
Regarding the number of clusters NC , the larger it is, the more accurately
estimator (9) approximates (7). Nevertheless, NC cannot be arbitrary large
since it should be computationally feasible to construct the k-d trees and the
K-nearest neighbor graph. Due to this reason, we set NC to 106 .
Increase of the neighbor list size K results in a larger number of clusters,
which should be processed at each iteration of the loop in lines 10–14 of Algorithm 1. On the other hand, it reduces the number of loop iterations, which are
required to achieve sufficient denoising accuracy, and makes it possible to set
smaller M in the convergence criterion. We utilized K = 2d and M = 10 in our
experiments.
The number of k-d trees NT and the k-d tree leaf size L do not significantly
affect the results and were set to NT = 64 and L = 32.

4.3

Results

The results are presented in Table 1. For small peak intensities, the presented
method outperforms the existing techniques by about 0.3 dB. As the peak intensity increases, the results of our algorithm become comparable to those of
[11].
It took, on average, 14 minutes for our method implemented in C++ to
process one 256 × 256 image on a PC with CPU Intel i7-5820K 3.3 GHz.

5

Discussion and Conclusion

In this work, we present a new noise suppression method and demonstrate its
applicability to Poisson noise removal. Compared with the approaches [7, 14,
11], our algorithm has a different assumption about the original image. Instead
of utilizing image self-similarity or sparse patch representation, patches of the
original image are assumed to be sampled from the distribution representing
natural image patches. As a result, denoised patches are computed using the
MMSE estimator, i.e. as the conditional expectation of the original patch given
the noisy patch (4). As pointed out in [10], such approach would provide optimal
denoising results if the number of patches NP tended to infinity and the patch
size d was increased to the image size. Nevertheless, it is useful also with fixed
NP and d. As the experiments showed, it is the preferable method for peak
intensities below 4, for which the variance stabilization is not efficient enough.
Since a large number of patches is necessary to represent the image patch
distribution, the choice of suitable data structures is of primary importance in
order to calculate the MMSE estimator with a reasonable computational effort
and make it practically applicable. In this work, it is demonstrated that the k-d
7

Table 1: The denoising results. Each value is the average of five noise realizations. The best result is selected with the bold font. The PSNR is given in
decibels.
Method
SPDA [7]
BM3D [11]
proposed
SPDA [7]
BM3D [11]
proposed
SPDA [7]
BM3D [11]
proposed
SPDA [7]
BM3D [11]
proposed
SPDA [7]
BM3D [11]
proposed

Peak

1

2

3

4

5

Peppers
19.96
19.86
20.38
21.20
21.94
22.26
21.91
23.35
23.37
22.06
24.03
23.92
22.31
24.67
24.40

Bridge
19.21
19.22
19.55
20.14
20.31
20.65
20.43
21.00
21.21
20.55
21.50
21.60
20.60
21.87
21.90

Boat
20.07
20.01
20.24
21.01
21.07
21.28
21.34
21.86
21.93
21.51
22.36
22.32
21.60
22.84
22.68

Image
Couple
20.05
19.93
20.26
20.97
21.00
21.22
21.22
21.65
21.81
21.34
22.26
22.26
21.45
22.65
22.55

Average
Hill
20.61
20.69
20.98
21.81
21.63
22.05
22.13
22.48
22.78
22.16
23.06
23.23
22.22
23.51
23.61

Mandril
18.40
18.14
18.43
18.92
18.79
18.98
19.05
19.28
19.35
19.11
19.58
19.56
19.16
19.82
19.72

tree search followed by the search in the K-nearest neighbor graph allows fast
and accurate calculation of the MMSE estimator.

References
[1] Francis J Anscombe. The transformation of poisson, binomial and negativebinomial data. Biometrika, pages 246–254, 1948.
[2] Jérôme Boulanger, Charles Kervrann, Patrick Bouthemy, Peter Elbau,
Jean-Baptiste Sibarita, and Jean Salamero. Patch-based nonlocal functional for denoising fluorescence microscopy image sequences. Medical
Imaging, IEEE Transactions on, 29(2):442–454, 2010.
[3] Harold Christopher Burger, Christian Schuler, and Stefan Harmeling.
Learning how to combine internal and external denoising methods. In
Pattern Recognition, pages 121–130. Springer, 2013.
[4] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian. Image denoising by
sparse 3-D transform-domain collaborative filtering. IEEE Transactions on
Image Processing, 16(8):2080–2095, 2007.

8

Pirate
20.06
20.12
20.49
21.22
21.28
21.60
21.60
22.04
22.29
21.86
22.61
22.79
22.04
23.03
23.09

19.77
19.71
20.05
20.75
20.86
21.15
21.10
21.67
21.82
21.23
22.20
22.24
21.34
22.63
22.56

[5] Michael Elad and Michal Aharon. Image denoising via sparse and redundant representations over learned dictionaries. Image Processing, IEEE
Transactions on, 15(12):3736–3745, 2006.
[6] A. Foi. Clipped noisy images: Heteroskedastic modeling and practical
denoising. Signal Processing, 89(12):2609–2629, 2009.
[7] Raja Giryes and Michael Elad. Sparsity-based poisson denoising with dictionary learning. Image Processing, IEEE Transactions on, 23(12):5057–
5069, 2014.
[8] Chul Lee, Chulwoo Lee, and Chang-Su Kim. An mmse approach to nonlocal
image denoising: Theory and practical implementation. Journal of Visual
Communication and Image Representation, 23(3):476–490, 2012.
[9] Anat Levin and Boaz Nadler. Natural image denoising: Optimality and
inherent bounds. In Computer Vision and Pattern Recognition (CVPR),
2011 IEEE Conference on, pages 2833–2840. IEEE, 2011.
[10] Anat Levin, Boaz Nadler, Fredo Durand, and William T Freeman. Patch
complexity, finite pixel correlations and optimal denoising. In Computer
Vision–ECCV 2012, pages 73–86. Springer, 2012.
[11] Markku Mäkitalo and Alessandro Foi. Optimal inversion of the anscombe
transformation in low-count poisson image denoising. Image Processing,
IEEE Transactions on, 20(1):99–109, 2011.
[12] D. Martin, C. Fowlkes, D. Tal, and J. Malik. A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics. In Proc. 8th Int’l Conf. Computer Vision, volume 2, pages 416–423, July 2001.
[13] Marius Muja and David G Lowe. Scalable nearest neighbor algorithms for
high dimensional data. Pattern Analysis and Machine Intelligence, IEEE
Transactions on, 36(11):2227–2240, 2014.
[14] Joseph Salmon, Zachary Harmany, Charles-Alban Deledalle, and Rebecca
Willett. Poisson noise reduction with non-local PCA. Journal of mathematical imaging and vision, 48(2):279–294, 2014.
[15] Bo Zhang, Jalal M Fadili, and Jean-Luc Starck. Wavelets, ridgelets, and
curvelets for poisson noise removal. Image Processing, IEEE Transactions
on, 17(7):1093–1108, 2008.

9

