Denoising 3D Microscopy Images of Cell Nuclei
using Shape Priors on an Anisotropic Grid
Mathieu Bouyrie
1,2
, Cristina Manfredotti
1
, Nadine Peyri
´
eras
2
and Antoine Cornu
´
ejols
1
1
AgroParisTech, MMIP Department / INRA UMR518, 16 rue Claude Bernard, 75231 Paris Cedex 05, France
2
BioEmergences Lab, CNRS USR3695, Ave. de la Terrasse, 91198 Gif-sur-Yvette Cedex, France
Keywords:
3D Image Denoising, Biomedical Imaging, Mixed Poisson-Gaussian Noise, Isotropic Undecimated Wavelet
Transform, Stabilizing Transform, Hypothesis Testing, Optimization Algorithm, Cell Detection.
Abstract:
This paper presents a new multiscale method to denoise three-dimensional images of cell nuclei. The speci-
ficity of this method is its awareness of the noise distribution and object shapes. It combines a multiscale
representation called Isotropic Undecimated Wavelet Transform (IUWT) with a nonlinear transform, a statis-
tical test and a variational method, to retrieve spherical shapes in the image. Beyond extending an existing
2D approach to a 3D problem, our algorithm takes the sampling grid dimensions into account. We compare
our method to the two algorithms from which it is derived on a representative image analysis task, and show
that it is superior to both of them. It brings a slight improvement in the signal-to-noise ratio and a significant
improvement in cell detection.
1 INTRODUCTION
Three-dimensional (3D) confocal microscopy is a
technique that can image live specimens (in vivo)
and capture their entire volume (in toto). Imaging a
whole embryo at different time steps leads to 3D+time
movies displaying either cell nuclei or membranes, or
both. In this article we deal with images of cell nuclei.
From these datasets, numerical reconstruction algo-
rithms can extract, or “reconstruct”, the embryo’s ge-
ometric information by shape segmentation (Zanella
et al., 2010), and the cell trajectories and lineage by
cell tracking (Olivier et al., 2010). Then, based on
these numerical reconstructions, biologists can per-
form a quantitative and qualitative analysis to explain
the specimen’s developmental behavior.
Several improvements have been made to the orig-
inal concept of confocal microscopy, leading to two-
photon microscopy (Denk et al., 1990). This tech-
nique relies on modified RNA injected in the embryo
to make cell nuclei produce certain fluorescent pro-
teins. At the convergence of two infrared beams, these
proteins emit photons in the visible spectrum. The
light thus generated is collected by a camera. Moving
the focal point of the two beams on a Cartesian grid,
the entire volume of the embryo can be recorded as a
3D grayscale image. Figure 1 shows a 2D slice at a
fixed depth. As shown this figure, despite the advan-
Figure 1: Example of 2D slice from a 3D image of cell
nuclei in Ascidia (source: A. Rausch, CNRS). This image
presents higly deteriorated signal with low resolution.
tages of two-photon microscopy, the general process
still suffers from deterioration, which can become a
real problem when applying segmentation and track-
ing algorithms because image quality can have, in
fact, a great impact on the precision of these tasks.
Image deterioration takes many forms, and this paper
focuses on three of them: (1) a random distribution
whose expectation is a clean signal; (2) ectopic mark-
ing, i.e. adverse signal that can pollute other regions
than the nuclei of interest; (3) loss of signal caused by
moving deeper inside the embryo.
In this paper, we present a 3D generalization of a
Bouyrie, M., Manfredotti, C., Peyriéras, N. and Cornuéjols, A.
Denoising 3D Microscopy Images of Cell Nuclei using Shape Priors on an Anisotropic Grid.
DOI: 10.5220/0005706002910298
In Proceedings of the 5th International Conference on Pattern Recognition Applications and Methods (ICPRAM 2016), pages 291-298
ISBN: 978-989-758-173-1
Copyright
c
 2016 by SCITEPRESS – Science and Technology Publications, Lda. All rights reserved
291

2D algorithm for the denoising and restoration of cell
nuclei image called “mutiscale variance stabilizing
transforms” (MS-VST), a combination of a ”variance
stabilizing transforms” (VST) with the filter banks of
wavelets (Zhang et al., 2007). Differently from other
methods, it incorporates a prior bias over the shapes
of the objects of interest and over the noise statistic
assumed to follow a Mixed Poisson-Gaussian (MPG)
distribution. In this case, imaged objects are expected
to be nearly isotropic, which means that they can be
locally approximated by spheres. This method has al-
ready proven its efficacy in dealing with the above
three deterioration types on 2D astronomy images,
thus we expect equivalent results on 3D cell nuclei
since their shapes are similar to stars.
Another contribution of this paper is to take
into account the image grid dimensions, which are
anisotropic in the case of two-photon microscopy,
i.e. the height of the voxel (a 3D pixel) in z is different
from its width and length in (x,y). An isotropic shape
in real space is no longer isotropic when sampled on
this type of grid. One solution to this problem is re-
sampling the image to get an isotropic grid. This may
induce some interpolations between voxels, which
would create statistic distribution changes, making
the original method irrelevant, and approximation er-
rors before actually denoising the image. Therefore,
instead of re-sampling the image, our method adapts
2D MS-VST to the 3D anistropic grid.
The method presented here can be decomposed in
three main steps (Figure 2):
1. The image is decomposed in several scales, from
high-level (coarse-grained) to low-level (detailed
voxels). The differences between two scales are
called “detail images” (Figure 3). The transform
from the whole image to the set of detail images is
called Isotropic Undecimated Wavelet Transform
(IUWT).
2. “Significant voxels” are selected from the detail
images by a statistical significance test.
3. The denoised image is retrieved from the signifi-
cant voxels of the detail images.
To assess the improvements we brought to 2D
MS-VST introducing the third dimension and the
anisotropic grid, we compare the performances of
our method with one that does not consider grid
anisotropy and another that analyses 2D slices only.
2 ISOTROPIC UNDECIMATED
WAVELET TRANSFORM
Figure 3a shows how the Isotropic Undecimated
Wavelet Transform (IUWT) is constructed (Zhang
et al., 2007). This transform provides a sparse rep-
resentation of the images containing mainly isotropic
shapes. It means that such images produce an out-
put with many null voxels . IUWT distributes the im-
age across several scales, from coarse-grained (image
a
4
) to detailed voxels (image a
0
). An IUWT takes as
input J, the number of scales on which the image is
scattered. It constructs the sequence of scales {a
j
}
j
from an image u as following
1
(Figure 3a):
a
0
= u
a
j+1
= a
j
∗h
( j)
(1)
Where the convolution operator ∗, is the filtering pro-
cess. It is defined in section 2.1, and the sequence of
filters h
( j)
is defined in section 2.2. The ”detail” im-
ages {d
j
}
j
are the differences between two consecu-
tive scales. Given two consecutive scales a
j
and a
j+1
,
the ”detail” image d
j+1
is obtained as (Figure 3a):
d
j+1
= a
j
−a
j+1
∗h
( j)
. (2)
The minus sign represents the subtraction voxel by
voxel of the two images. The output of the trans-
form is the set of all the ”detail” images {d
j
}
j
and
the coarsest scale a
J
. The reconstruction of the image
from the ”detail” images d
j
and the coarsest scale a
J
(Figure 3b) can be directly deduced from Eq. 2:
a
j
= d
j+1
+ a
j+1
∗h
( j)
. (3)
The following section makes explicit all the steps nec-
essary to understand the whole process of the IUWT.
Our contribution is reported in section 2.3: we have
extended the IUWT to a non isotropic grid, i.e. a grid
with voxel dimensions in x, y, or z be equal.
2.1 Filtering 3D Grayscale Images
Grayscale images can be seen as 3D arrays of real
numbers. The indices designing the location along the
x, y, and z axes are denoted i
x
, i
y
, and i
z
. We denote
the set of all the three indices as i
= (i
x
,i
y
,i
z
). With
the description above, a voxel is denoted as u(i
x
,i
y
,i
z
)
or u(i).
A 3D filter is a function defined by (k
x
,k
y
,k
z
) 7→
h(k
x
,k
y
,k
z
), k
x
,k
y
and k
z
in Z. Filters can be applied
on images by the convolution product defined by:
(w∗h)(i) =
∑
k
x
,k
y
,k
z
w
o
(i
x
−k
x
,i
y
−k
y
,i
z
−k
z
)h(k
x
,k
y
,k
z
)
(4)
1
Throughout the paper, u is the noisy image.
ICPRAM 2016 - International Conference on Pattern Recognition Applications and Methods
292

Figure 2: Global denoising scheme in 2D. Step 1: Multiscale transform of the image, from coarse-grained (image a
4
in the
back) to low-level (image a
0
, in the front). Step 2: Finding relevant coefficients for the multiscale transform by hypothesis
test. Step 3: Retrieving a denoised image from the relevant coefficients.
Where w is a 3D image. F A filter can be seen as
a local mixture around each voxel of the image. In
practical cases, the filters are designed so that if k
x
,k
y
,
or k
z
is far enough from 0, then h(k
x
,k
y
,k
z
) is null.
This implies that this sum is computed on a bounded
window around the voxel i.
Note that some simple conditions, which we do
not detail here, must be imposed on the borders of the
image. Equation (4) can be used to make a new filter
from the convolution of two filters.
The filters that we use are separable and can be
written as h(k
x
,k
y
,k
z
) = h
x
(k
x
)h
y
(k
y
)h
z
(k
z
), where h
x
,
h
y
, and h
z
are some functions from Z to R.
2.2 Construction of the h
( j)
Sequence
The sequence of filters h
( j)
are separable and are
made from an original filter h, which is separable too
(h(k
x
,k
y
,k
z
) = h
x
(k
x
)h
y
(k
y
)h
z
(k
z
)). If δ is the direc-
tion x,y or z, the construction of h
( j)
δ
follows the prin-
ciple:
h
( j)
δ
(i) =

h
δ
(
i
2
j
) if
i
2
j
is integer.
0 otherwise
Hence we have
h
( j)
(k
x
,k
y
,k
z
) = h
( j)
x
(k
x
)h
( j)
y
(k
y
)h
( j)
z
(k
z
)
In two dimensions, (Zhang et al., 2007) designed
a filter h
1D
such that h
1D
(−2) = 1/16, h
1D
(−1) =
4/16, h
1D
(0) = 6/16, h
1D
(1) = 4/16, h
1D
(2) = 1/16,
h
1D
(k) = 0, if k is not between -2 and 2. Then they
set h
x
= h
y
= h
1D
. Extending the 2D filter in 3D is
not readily possible because of the non isotropy of the
grid: isotropic shapes in real space are not isotropic
on the grid of voxels.
2.3 Design of the h Filter to Deal with
an Anisotropic Grid
To preserve the ability to represent isotropic forms
from two-photon microscopy images, the filter h has
to be redesigned. As in 1D h
1D
=
1
16
[1,4,6,4,1] and
in 2D, h
x
= h
y
= h
1D
, the 3D definition would be
h
z
= h
1D
. But this simple combination of the same
h
1D
along the 3 directions does not take into account
the non-isotropy of the grid and will not properly rep-
resent isotropic shapes.
In (Cai, 1989), a parallel is made between a fil-
tering process using h
1D
and a diffusion process. In
particular, they discretized the diffusion process in
space and time and showed, with some pre-defined
time-step and spatial-step, that a diffusion process is
equivalent to a filtering process using h
1D
. We use
this result to redefine our filter fixing the spatial step.
Let the diffusion equation be defined as follows:
∂ f
∂t
=
∂
2
f
∂x
2
. (5)
In our case, the time t is artificial. To switch from
continuous problems to discrete ones, there exist nu-
merical schemes to solve Partial Differential Equa-
tions (PDE) like (5). We use the Finite Difference
Method (Smith, 1985). Let 4t be the time step and an
adjustable variable, and and 4x, the size of the voxel
along one direction. We aim at finding as a solution
of a PDE, a good approximation
˜
f (t
n
,i) where t
n
is
a discrete version of the time and i of the space. We
denote
˜
f
n
(i) =
˜
f (t
n
,i). With the condition
4x
2
4t
> 2,
the Finite Difference Method leads to:
˜
f
n+1
=
˜
f
n
∗
4t
4x
2

1,
4x
2
4t
−2, 1

(6)
The condition
4x
2
4t
> 2, is actually fulfilled with-
out problem: as already said, time is an artificial vari-
able, so we can fix the time step for the inequality to
be always verified. Moreover we set r
x
=
4x
2
4t
:
˜
f
2
=
˜
f
1
∗
1
r
x
[1,r
x
−2, 1]
=
˜
f
0
∗
1
r
x
[1,r
x
−2, 1] ∗
1
r
x
[1,r
x
−2, 1]
=
˜
f
0
∗
1
r
2
x
[1,2r
x
−4, r
2
x
−4r
x
+ 6, 2r
x
−4, 1] (7)
Denoising 3D Microscopy Images of Cell Nuclei using Shape Priors on an Anisotropic Grid
293

(a) construction of transform from the image (b) reconstruction of the image from the ”detail” images
and the coarsest scale
Figure 3: Illustration of the IUWT algorithm for J = 4: Direct Transform and Reconstruction. In green input of the algorithms,
in red output.
If r
x
= 4, then equation (7) becomes
˜
f
0
convoluted by
h
1D
filter. This shows how a filtering process can be
interpreted as a discrete diffusion process. h
1D
can
then be generalized with:
h
x
=
1
r
2
x
[1,2r
x
−4, r
2
x
−4r
x
+ 6, 2r
x
−4, 1]
Using the same method, h
x
, h
y
and h
z
can be de-
fined with r
y
=
4y
2
4t
and r
z
=
4z
2
4t
. 4x, 4y and 4z, are
the size of the voxel along x, y, and z, and are fixed.
The whole 3D filter is actually parameterized by only
one parameter: 4t which must verify:
4x
2
4t
> 2 and
4y
2
4t
> 2 and
4z
2
4t
> 2 (8)
which is equivalent to
4t <
1
2
MIN(4x,4y,4z)
2
(9)
We finally fix
2
4t =
1
4
MIN(4x,4y,4z)
2
and we
build the 3D filter.
3 MUTISCALE VARIANCE
STABILIZING TRANSFORMS
MS-VST is a nonlinear transform dealing with MPG
noise. In the biomedical imaging state of the art,
methods which consider explicitly the characteris-
tics of microscopy imaging noise are rare (Boulanger
et al., 2010) and (Zhang et al., 2007). As noise model
we consider the mixed Poisson-Gaussian Noise:
u(i) = P (v(i)) + ε((i)). (10)
2
In this way, (9) is verified, and we get the h
1D
filter in
at least 1 Direction
where v is the clean image and v(i) the expected
value of u(i), P (v(i)) is a poisson noise of parame-
ter v(i): the probability of u(i) taking the value m is
P(u(i) = m) = e
−v(i)
v(i)
m
m!
. ε(i) ∼ N (0,σ
2
), is a noise
distributed as a normal distribution of variance σ
2
, un-
known a-priori and expectation equal to 0.
As explained before, ectopic marking can be as
harmful as noise. We use the method presented in
(Zhang et al., 2007) to select the scale of the filtered
objects to filter the ectopic marking as well.
In the following, we recall how IUWT is used in
the literature to clean images deteriorated by Mixed
Poisson-Gaussian (MPG) Noise. First we realize an
hypothesis test to select significant voxels in the ”de-
tail” images using a technique called Variance Stabi-
lizing Transform (second step of Figure 2). Once we
get all significant voxels on the ”detail” images, the
third step of Figure 2 describes an iterative method to
retrieve the image from them.
3.1 Hypothesis Test on an IUWT
Let d
j
(i
) be a voxel of a ”detail” image. We want
to check if it is relevant or not. To do so, we de-
fine a test function on hypotheses: we want to test
if E(d
j
(i)) = 0. E(d
j
(i)) is the noisy image’s IUWT
and it corresponds to the clean image’s IUWT.
In the following paragraph we suppose that
E(d
j
(i)) = 0, and consequently (Zhang et al., 2007):
q
(a
0
∗h
[ j]
)(i) + c
j
−
q
(a
0
∗h
0[ j+1]
)(i) + c
0
j+1
D
−→
v(i)→∞
N

0,σ
2
j

(11)
ICPRAM 2016 - International Conference on Pattern Recognition Applications and Methods
294

where N is the normal distribution,
σ
2
j
=
1
4
τ
2
(h
[ j]
)
(τ
1
(h
[ j]
))
2
+
1
4
τ
2
(h
0[ j+1]
)
(τ
1
(h
0[ j+1]
))
2
−
hh
[ j]
,h
0[ j+1]
i
2τ
1
(h
[ j]
)τ
1
(h
0[ j+1]
)
and
hh
[ j]
,h
0[ j+1]
i =
∑
k
x
,k
y
,k
z
h
[ j]
(k
x
,k
y
,k
z
)h
0[ j+1]
(k
x
,k
y
,k
z
)
τ
p
(h
[ j]
) =
∑
k
x
,k
y
,k
z

h
[ j]
(k
x
,k
y
,k
z
)

p
τ
p
(h
0[ j+1]
) =
∑
k
x
,k
y
,k
z

h
0[ j+1]
(k
x
,k
y
,k
z
)

p
The symbol
D
−→
v(i)→∞
means the convergence in proba-
bility distribution: for v(i) large enough, the left hand
side of the operator ”nearly follows” a distribution
given by the right hand side. Finally the filters h
[ j]
and h
0[ j+1]
are defined by (δ is either x, or y, or z):
h
[ j]
δ
= h
(0)
δ
∗h
(1)
δ
∗. . . ∗h
( j)
δ
, h
[ j]
= h
[ j]
x
h
[ j]
y
h
[ j]
z
and
h
0[ j+1]
δ
= h
[ j+1]
δ
∗h
( j+1)
δ
, h
0[ j+1]
= h
0[ j+1]
x
h
0[ j+1]
y
h
0[ j+1]
z
This function is called a Variance Stabilization Trans-
form, and it is a method that decorrelates noise from
signal. In our case, the MPG Noise becomes Gaus-
sian. Equation (11) works with any c
j
, c
0
j+1
. The
authors determined the best c
j
and c
0
j+1
so that the
convergence is faster: c
j
=
7
8
τ
2
(h
[ j]
)
τ
1
(h
[ j]
)
2
−
1
2
τ
3
(h
[ j]
)
τ
2
(h
[ j]
)
and
c
0
j+1
=
7
8
τ
2
(h
0[ j+1]
)
τ
1
(h
0[ j+1]
)
2
−
1
2
τ
3
(h
0[ j+1]
)
τ
2
(h
0[ j+1]
)
. They experimentally
discovered that the convergence was reached for very
low values of v(i) (i.e. for v(i) ≥ 0.1).
All these results are true if E(d
j
(i)) = 0. Thus we
can now define a function for hypothesis testing: we
want to test whether E(d
j
(i)) is null. If E(d
j
(i)) is
null then (11) says that l
i, j
=
q
(a
0
∗h
[ j]
)(i) + c
j
−
q
(a
0
∗h
0[ j+1]
)(i) + c
0
j+1
follows a normal distribu-
tion with mean 0 and known variance. We ensure then
that, under the hypothesis E(d
j
(i)) = 0,
P(|X|< |l
i, j
|) < FPR (12)
FPR being the false positive rate, a parameter of the
hypothesis test, and P(|X| < |l
i, j
|) = 1 −erf

|l
i, j
|
√
2σ
j

.
The larger the FPR, the more the noise is removed,
and so are the structures in the image. Here, erf is the
classical error function of the centered normal distri-
bution. If (12) is verified then we accept the hypothe-
sis, we reject it otherwise.
In this way, we get all the coefficients d
j
(i) that
are not significantly null and we define the set M by
M = {(i, j), d
j
(i) not significantly null}
=
(
(i, j),erf
|l
i, j
|
√
2σ
j
!
> 1 −FPR
)
3.2 Scale Selection to Improve
Denoising Efficiency and
Background Removal
The background, due to ectopic marking, is a signal
whose scale is larger than that of the cell nuclei. For
this reason, removing all tuples (i
, j) for j > J
M
from
M (Zhang et al., 2007), is a good way to remove the
background. Here J
M
is a scale defining the maximal
size an object can have in an image.
Another useful property of scale selection is that it
allows one to improve the set of relevant coefficients
by imposing that the objects present in the image have
a minimal size given by another scale index J
m
. We
can then remove all tuples (i, j) for j < J
m
from M .
Even though J
m
and J
M
are parameters of the
method, they can be set very easily with a pre-
visualization, before applying the denoising method.
To do so, we realize a IUWT of our noisy image,
we put d
j
(i) = 0 for all (i, j) such that j > J
M
and,
j < J
m
, and we reconstruct the image with the new
coefficients. M can be then redefined as follows:
M =
(
(i, j),erf
|l
i, j
|
√
2σ
j
!
> 1 −FPR and J
m
≤ j ≤ J
M
)
(13)
In practice, defining the J parameter of the IUWT
when it is larger than J
M
+ 1 is useless since all the
images with scale greater than J
M
are set to 0 anyway.
From now on in this paper, we will impose J = J
M
+1.
3.3 Reconstructing an Image from
Relevant Coefficients
We now investigate how to reconstruct a denoised
image from relevant coefficients (third step of Fig-
ure 2). Once all the significant coefficients have been
acquired, the naive method would be to set all other
coefficients to 0 and then reconstruct the image. This
solution is however imperfect: indeed, some of the
voxels in the reconstructed image may have negative
intensities. By definition of the mixed Poisson Noise,
the black and white image intensities (i.e. the parame-
ter of the Poisson distribution) are non-negative. This
problem suggests a more subtle reconstruction model.
Denoising 3D Microscopy Images of Cell Nuclei using Shape Priors on an Anisotropic Grid
295

The method proposed in (Zhang et al., 2007) is
based on constrained optimization. The idea is to add
few non-null intensities in the transform so that the re-
constructed images are non-negative. To do so, they
establish a cost which constrains the number of pos-
itive coefficients added to the set of significant coef-
ficients, to keep it as small as possible. If s is the
IUWT representation of the denoised image, the cost
is ksk
1
=
∑
i, j
|s
j
(i)|. We want to minimize the cost
while respecting two constraints:
• the reconstruction of the image from s has to be
non-negative
• s
j
(i) = d
j
(i) for all tuple (i, j) in M .
The final model is the solution of the following opti-
mization problem:
minimize ksk
1
subject to

s
j
(i) = d
j
(i) for all i in M
(Rs)(i) ≥ 0 for all i
where R is the reconstruction operator defined by (3)
and given in Figure 3b.
The solution of this model is obtained by some ad-
vanced iterative optimization algorithm described in
(Chambolle and Pock, 2011) which can be adapted to
our problem as shown in Algorithm (1). This algo-
rithm, can be viewed as the combination of a gradient
descent in the primal space (i.e. IUWT span), and a
gradient ascent in the dual space (i.e. image space). It
takes as argument, the noisy image, the set M , and
the ratio between the gradient descent step and the
gradient ascent step. It returns the denoised image.
In this algorithm, the + and the − operators define
the element-wise addition and difference for the two
operands. The multiplication by a real number is an
element-wise operation too. The soft-tresholding op-
erator is defined as follows:
(ST
τ
(s))
j
(i) = max(|s
j
(i)|−τ,0)sign(s
j
((i))) (14)
The P
M
operator constrains its argument s to be equal
to d for the tuples (i, j) in M :
(P
M
(s))
j
(i) =

d
j
(i) if (i, j) ∈ M
s
j
(i) otherwise.
The P
O
−
is the projection operator on the set of arrays
with non-positive elements:
(P
O
−
(u))(i) =

u(i) if u(i) ≤ 0
0 otherwise.
Finally the R
∗
operator is given by the following for-
mula:
(R
∗
u)
0
= u
(R
∗
u)
j+1
= (R
∗
u)
j
∗h
( j)
We chose this algorithm because of its very fast con-
vergence (experimentally it converges within 50 iter-
ations).
Data: Noisy Image u, r
s
, J
m
and J
M
, False
Positive rate FPR
Result: Denoised Image ˜v.
Set M as defined in (13);
τ ←
q
1
r
s
(J
M
+2)
2
;
σ ←
q
r
s
(J
M
+2)
2
;
u
0
← u;
(s
0
)
j
(i) ← 0 for all i and j;
s
0
← s
0
;
repeat
u
k+1
← P
O
−
(u
k
+ σR s
k
);
s
k+1
← P
M
(ST
τ
(s
k
−τR
∗
u
k+1
));
s
k+1
← 2s
k+1
−s
k
;
k ←k +1;
until convergence;
˜v ← R s
∞
;
Algorithm 1: Reconstructing Image Algorithm.
4 EXPERIMENTAL PROTOCOL
To assess the validity of our method, in particular the
importance of combining both 3D and anisotropy, we
compared it with two other partial methods: 1) the
2D version from which it was designed (Zhang et al.,
2007) applied to each fixed-depth slice of the 3D im-
age; 2) another version in 3D assuming an isotropic
grid, i.e. the same voxel size along every direction.
For a better comparison, we built an “artificial em-
bryo” to have a common ground truth. We tested all
three methods on two different experiments. In the
first experiment, the methods were prevented from re-
moving the background and we measured the result-
ing signal-to-noise ratio with respect to the noiseless
ground-truth image containing the cells and the back-
ground. In the second experiment, the three methods
were allowed to remove the background signal and
we applied the “spot detector” algorithm by (Olivo-
Marin, 2002). This software was run on the ICY
graphical interface (De Chaumont et al., 2012).
4.1 Building the Artificial Embryo
To simulate 3D images of an embryo, we first dis-
tributed a few centers over the whole 3D volume such
that neighboring cells were separated by a minimal
distance, while confined within a sphere of a cer-
tain radius. In each center, we installed a full el-
lipse with random size, random orientation, and a de-
formed contour to represent a nucleus. As for the
background, we applied a highly smoothing filter to
the output of the previous step. Then, we added the
background to the image of the cell nuclei voxel by
ICPRAM 2016 - International Conference on Pattern Recognition Applications and Methods
296

voxel. Finally, the image was deteriorated with MPG
noise using a random generator.
4.2 Results
The three tested methods are parametrized by four
quantities. For each image, for each method, we tuned
manually these parameters to get the best result. The
best result is the one that gives the best measure com-
pared to the groundtruth image. The main parameter
is the FPR parameter explained in Section 3.1, rang-
ing in [0, 1]. It defines the filtering coefficient: the
higher its value, the better the structures present in
the image are detected, but so is the noise, too. In all
our experiments, it was set to values lower than 0.3.
The second and third parameters are J
m
and J
M
de-
fined in Section 3.2 controlling the scale at which the
structures of the image are filtered. These parameters
only depend on the size of the objects and can be pre-
set before applying the denoising method. To do so,
one can directly filter the desired scale from the im-
age and check the output. The final quantity, and the
most difficult to tune, is the convergence parameter
from Section 3.3. It only depends on the reconstruc-
tion operator R which is itself a function of the size of
the image, the dimensions of its voxels, and the num-
ber of scales on which it processes an image (here
J
M
+1). In our experiments, we designed it manually,
and found a good default value (5e-5) for a fixed im-
age size, which also worked at any scale (i.e. with any
J
M
).
4.2.1 A Comparison without Background
Subtraction
The images were then denoised and we could com-
pare the different results. Figure 4 shows the signal-
to-noise ratio (SNR) of the ouput image after ap-
plying each denoising method. SNR is defined by
−10log
10

k˜v −vk
2
2
/kvk
2
2

(where v is the ground
truth and ˜v the denoised image. The higher the SNR,
the less noisy the image. All three methods were
tested on three images with different noise levels
and three different original SNR parameters (called
SNR
0
), which produced three graphs. We also used
three different scale ranges, where J
m
= 0, 1, 2 and
J
M
= 8. Results show that our method is only a slight
improvement over the isotropic version, but is much
better than the 2D version. Comparing the two 3D
methods, we also noticed that the best scale selection
is not the same when changing filters. This can be ex-
plained by the fact that the size of a cell does not fit a
whole scale, but stands between two scales.
4.2.2 A Comparison with Background
Subtraction
To subtract the background, we used the “spot detec-
tor” software (Olivo-Marin, 2002) which detects nu-
clei in the denoised image. This software actually re-
lies on IUWT to detect cells in 3D images. As it was
run on artificial images, we already know the true po-
sitions and number of cells. We set the parameter of
all the methods as follows
3
: J
m
= 2, J
M
= 3, r
s
= 5e-5,
and for each method, we adjusted the FPR parameter
so that the number of detected cells in the software
were the number of actual cells present in the image
(our ground-truth image contains 235 cells). In this
way we only had to compare only one quantity since
the number of false positives is equal to the number
of false negatives. For each output, all the parameters
of the detection algorithm but one were set to the de-
fault values. The tuned parameter of detection was the
scale, i.e. size of the objects we were looking for, and
was set to 3. For each image and algorithm, the spatial
precision of the detection was no greater than 3.5 pix-
els (for a cell whose actual size is around 15 voxels).
We ran this experiment on two images whose original
signal-to-noise ratio were respectively 2.68 and 0.075
(Table 1). As one can see, the differences between all
the methods, while favorable to our method, are not
significant on the less noisy image. However, when
image quality worsens, our method shows much bet-
ter results compared to the other two.
5 CONCLUSION
In this paper, we introduced a method that tack-
les noise problems in cell images coming from two-
photon microscopy. More than just denoising im-
ages of cells, it is also able to remove the background
due to ectopic marking. Relying on a known multi-
scale transform called IUWT and a variance stabiliza-
tion nonlinear transform, our method is able to deal
with Mixed Poisson-Gaussian noise in two-photon
microscopy. This method comes from an extension of
a 2D original method and we were able to show that it
leads to a significant improvement over it. Moreover,
we were able to take into account the grid anisotropy,
which is of paramount importance to obtain better re-
sults in cell detection. Even if the difference is not
outstanding as measured by the signal-to-noise ratio,
the results are quite satisfying when viewing the re-
sulting images. We expect better results with the fine-
3
for the scale parameters, the outcome does not change
much if we choose J
m
= 1
Denoising 3D Microscopy Images of Cell Nuclei using Shape Priors on an Anisotropic Grid
297

Image 1 (SNR
0
= 2.68) Image 2 (SNR
0
= 3.49) Image 3 (SNR
0
= 4.16)
Figure 4: Signal-to-noise ratio of the output of the three denoising methods. From left to right: application to three different
images. Inside each bar chart, from left to right: our 3D method with grid anisotroy, a 3D method without anisotropy, the 2D
method. In blue/red/yellow: three different scale ranges.
Table 1: Number of false detections on 2 example images after applying the 3 denoising methods.
Method 3D-Grid Anisotropy 3D-Grid Isotropy 2D
Image 1 (SNR
0
= 2.68) 1 1 1
Image 2 (SNR
0
= 0.075) 7 14 17
tuning of the parameters of our method, which is the
goal of future work.
This opens up more perspectives, such as compar-
ing our method to other methods not using IUWT. An-
other interesting working direction would be to com-
bine cell detection and denoising. Indeed the cell de-
tection technique we used was also involving IUWT:
it should thus be possible to merge denoising and cell
detection into a single overall method.
ACKNOWLEDGEMENTS
The authors would like to thank Jalal Fadili (GR-
EYC, CNRS) for his useful advices to go on with
this project. Thanks to Ren
´
e Doursat (BioEmer-
gences Lab, CNRS) for proof-reading the manuscript.
Thanks the Paris Complex System Institute (IS-
CPIF, CNRS) for facilities and material. Thanks to
A. Rausch and D. Fabresges (BioEmergences Lab,
CNRS) for providing the microscopy images in this
paper. Thanks to J. Sublime (MMIP, INRA) and P.A.
Murena for their comments. Thanks to the Idex Paris
Saclay Foundation for funding.
REFERENCES
Boulanger, J., Kervrann, C., Bouthemy, P., Elbau, P.,
Sibarita, J.-B., and Salamero, J. (2010). Patch-based
nonlocal functional for denoising fluorescence mi-
croscopy image sequences. Medical Imaging, IEEE
Transactions on, 29(2):442–454.
Cai, L. (1989). Spline smoothing: A special case of dif-
fusion smoothing. In Proceedings of the Alvey Vi-
sion Conference, AVC 1989, Reading, UK, September,
1989, pages 1–4.
Chambolle, A. and Pock, T. (2011). A first-order primal-
dual algorithm for convex problems with applications
to imaging. J. Math. Imaging Vis., 40(1):120–145.
De Chaumont, F., Dallongeville, S., Chenouard, N., Herv
´
e,
N., Pop, S., Provoost, T., Meas-Yedid, V., Panka-
jakshan, P., Lecomte, T., Le Montagner, Y., et al.
(2012). Icy: an open bioimage informatics platform
for extended reproducible research. Nature methods,
9(7):690–696.
Denk, W., Strickler, J. H., and Webb, W. W. (1990). Two-
Photon Laser Scanning Fluorescence Microscopy.
Science, 248:73–76.
Olivier, N., Luengo-Oroz, M. A., Duloquin, L., Faure,
E., Savy, T., Veilleux, I., Solinas, X., D
´
ebarre, D.,
Bourgine, P., Santos, A., et al. (2010). Cell lineage
reconstruction of early zebrafish embryos using label-
free nonlinear microscopy. Science, 329(5994):967–
971.
Olivo-Marin, J.-C. (2002). Extraction of spots in biological
images using multiscale products. Pattern recogni-
tion, 35(9):1989–1996.
Smith, G. D. (1985). Numerical solution of partial differ-
ential equations: finite difference methods. Oxford
university press.
Zanella, C., Campana, M., Rizzi, B., Melani, C., San-
guinetti, G., Bourgine, P., Mikula, K., Peyrieras, N.,
and Sarti, A. (2010). Cells segmentation from 3-d con-
focal images of early zebrafish embryogenesis. Image
Processing, IEEE Transactions on, 19(3):770–781.
Zhang, B., Fadili, M., Starck, J.-L., and Olivo-Marin, J.-C.
(2007). Multiscale variance-stabilizing transform for
mixed-poisson-gaussian processes and its applications
in bioimaging. In Image Processing, 2007. ICIP 2007.
IEEE International Conference on, volume 6, pages
VI–233. IEEE.
ICPRAM 2016 - International Conference on Pattern Recognition Applications and Methods
298