# Denoising 3D microscopy images of cell nuclei using shape priors on an anisotropic grid

Mathieu Bouyrie, Cristina Manfredotti, Nadine Peyriéras, Antoine Cornuéjols

► To cite this version:

Mathieu Bouyrie, Cristina Manfredotti, Nadine Peyriéras, Antoine Cornuéjols. Denoising 3D microscopy images of cell nuclei using shape priors on an anisotropic grid. ICPRAM 2016 International Conference on Pattern Recognition Applications and Method, Feb 2016, Rome, Italy. hal-01558022

HAL Id: hal-01558022

https://hal.science/hal-01558022

Submitted on 3 Jun 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Denoising 3D microscopy images of cell nuclei
using shape priors on an anisotropic grid

Mathieu Bouyrie<fnref n="1" />,<fnref n="2" />, Cristina Manfredotti<fnref n="1" />, Nadine Peyriéras<fnref n="2" /> and Antoine Cornuéjols<fnref n="1" />

**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 specificity 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 statistical 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 algorithms can extract, or “reconstruct”, the embryo’s geometric 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 perform a quantitative and qualitative analysis to explain the specimen’s developmental behavior.

Several improvements have been made to the original concept of confocal microscopy, leading to two-photon microscopy (Denk et al., 1990). This technique relies on modified RNA injected in the embryo to make cell nuclei produce certain fluorescent proteins. 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-

tages of two-photon microscopy, the general process still suffers from deterioration, which can become a real problem when applying segmentation and tracking 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 marking, 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 2D algorithm for the denoising and restoration of cell

<figure><img src="image_2.png" /><figcaption>Figure 1: Example of 2D slice from a 3D image of cell nuclei in *Ascidia* (source: A. Rausch, CNRS). This image presents highly deteriorated signal with low resolution.</figcaption></figure>

<footnote n="1">¹AgroParisTech, MMIP Department / INRA UMR518, 16 rue Claude Bernard, 75231 Paris Cedex 05, France</footnote>

<footnote n="2">²BioEmergences Lab, CNRS USR3695, Ave. de la Terrasse, 91198 Gif-sur-Yvette Cedex, France
{mathieu.bouyrie, cristina.manfredotti, antoine.cornuejols}@agroparistech.fr, nadine.peyrieras@inaf.cnrs-gif.fr</footnote>
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 already 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 errors before actually denoising the image. Therefore, instead of re-sampling the image, our method adapts 2D MS-VST to the 3D anisotropic 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 significant 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 representation of the images containing mainly isotropic shapes. It means that such images produce an output with many null voxels. IUWT distributes the image 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<fnref n="1" /> (Figure 3a):

$$a_0 = u \\ a_{j+1} = a_j * h^{(j)} \qquad (1)$$

Where the convolution operator *, is the filtering process. It is defined in section 2.1, and the sequence of filters $h^{(j)}$ is defined in section 2.2. The ”detail” images $\{d_j\}_j$ are the differences between two consecutive 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)}. \qquad (2)$$

The minus sign represents the subtraction voxel by voxel of the two images. The output of the transform 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)}. \qquad (3)$$

The following section makes explicit all the steps necessary 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 $\tilde{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(\tilde{i})$.

A 3D filter is a function defined by $(k_x, k_y, k_z) \mapsto h(k_x, k_y, k_z)$, $k_x, k_y$ and $k_z$ in $\mathbb{Z}$. Filters can be applied on images by the convolution product defined by:

$$(w * h)(i) = \sum_{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) \qquad (4)$$

Where $w$ is a 3D image. 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$,

<footnote n="1">¹Throughout the paper, $u$ is the noisy image.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: Global denoising scheme in 2D. Step 1: Multiscale transform of the image, from coarse-grained (image a₄ in the back) to low-level (image a₀, 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.</figcaption></figure>

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 $\mathbb{Z}$ to $\mathbb{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 $\delta$ is the direction x,y or z, the construction of $h_{\delta}^{(j)}$ follows the principle:

$$h_{\delta}^{(j)}(i) = \begin{cases} h_{\delta}\left(\frac{i}{2j}\right) & \text{if } \frac{i}{2j} \text{ is integer.} \\ 0 & \text{otherwise} \end{cases}$$

Hence we have

$$h^{(j)}(k_x, k_y, k_z) = h_x^{(j)}(k_x) h_y^{(j)}(k_y) h_z^{(j)}(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} = \frac{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 represent isotropic shapes.

In (Cai, 1989), a parallel is made between a filtering 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:

$$\frac{\partial f}{\partial t} = \frac{\partial^2 f}{\partial x^2} \qquad (5)$$

In our case, the time $t$ is artificial. To switch from continuous problems to discrete ones, there exist numerical schemes to solve Partial Differential Equations (PDE) like (5). We use the Finite Difference Method (Smith, 1985). Let $\Delta t$ be the time step and an adjustable variable, and $\Delta x$, the size of the voxel along one direction. We aim at finding as a solution of a PDE, a good approximation $\tilde{f}(t_n, i)$ where $t_n$ is a discrete version of the time and $i$ of the space. We denote $\tilde{f}_n(i) = \tilde{f}(t_n, i)$. With the condition $\frac{\Delta x^2}{\Delta t} > 2$, the Finite Difference Method leads to:

$$\tilde{f}_{n+1} = \tilde{f}_n * \frac{\Delta t}{\Delta x^2} \left[ 1, \frac{\Delta x^2}{\Delta t} - 2, 1 \right] \quad (6)$$

The condition $\frac{\Delta x^2}{\Delta t} > 2$, is actually fulfilled without problem: as already said, time is an artificial variable, so we can fix the time step for the inequality to be always verified. Moreover we set $r_x = \frac{\Delta x^2}{\Delta t}$:

$$\begin{aligned} \tilde{f}_2 &= \tilde{f}_1 * \frac{1}{r_x} [1, r_x - 2, 1] \\ &= \tilde{f}_0 * \frac{1}{r_x} [1, r_x - 2, 1] * \frac{1}{r_x} [1, r_x - 2, 1] \\ &= \tilde{f}_0 * \frac{1}{r_x^2} [1, 2r_x - 4, r_x^2 - 4r_x + 6, 2r_x - 4, 1] \end{aligned} \quad (7)$$

If $r_x = 4$, then equation (7) becomes $\tilde{f}_0$ convoluted by $h_{1D}$ filter. This shows how a filtering process can be
<figure><img src="image_4.png" /><figcaption>Figure 3: Illustration of the IUWT algorithm for $J = 4$: Direct Transform and Reconstruction. In green input of the algorithms, in red output.</figcaption></figure>

interpreted as a discrete diffusion process. $h_{1D}$ can then be generalized with:

$$h_x = \frac{1}{r_x^2} [1, 2r_x - 4, r_x^2 - 4r_x + 6, 2r_x - 4, 1]$$

Using the same method, $h_x$, $h_y$ and $h_z$ can be defined with $r_y = \frac{\Delta y^2}{\Delta t}$ and $r_z = \frac{\Delta z^2}{\Delta t}$. $\Delta x$, $\Delta y$ and $\Delta z$ 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: $\Delta t$ which must verify:

$$\frac{\Delta x^2}{\Delta t} > 2 \text{ and } \frac{\Delta y^2}{\Delta t} > 2 \text{ and } \frac{\Delta z^2}{\Delta t} > 2 \quad (8)$$

which is equivalent to

$$\Delta t < \frac{1}{2} \text{MIN}(\Delta x, \Delta y, \Delta z)^2 \quad (9)$$

We finally fix<fnref n="2" /> $\Delta t = \frac{1}{4} \text{MIN}(\Delta x, \Delta y, \Delta z)^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 characteristics 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(\hat{i}) = \mathcal{P}(v(\hat{i})) + \varepsilon((\hat{i})). \quad (10)$$

where $v$ is the clean image and $v(\hat{i})$ the expected value of $u(\hat{i})$, $\mathcal{P}(v(\hat{i}))$ is a poisson noise of parameter $v(\hat{i})$: the probability of $u(\hat{i})$ taking the value $m$ is $\mathcal{P}(u(\hat{i}) = m) = e^{-v(\hat{i})} \frac{v(\hat{i})^m}{m!}$. $\varepsilon(\hat{i}) \sim \mathcal{N}(0, \sigma^2)$, is a noise distributed as a normal distribution of variance $\sigma^2$, unknown 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 "detail" images using a technique called Variance Stabilizing 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(\hat{i})$ be a voxel of a "detail" image. We want to check if it is relevant or not. To do so, we define a test function on hypotheses: we want to test if $\mathbb{E}(d_j(\hat{i})) = 0$. $\mathbb{E}(d_j(\hat{i}))$ is the noisy image's IUWT and it corresponds to the clean image's IUWT.

In the following paragraph we suppose that $\mathbb{E}(d_j(\hat{i})) = 0$, and consequently (Zhang et al., 2007):

$$\sqrt{(a_0 * h^{[j]})(\hat{i}) + c_j} - \sqrt{(a_0 * h^{[j+1]})(\hat{i}) + c'_{j+1}} \xrightarrow{\mathcal{D}} \mathcal{N}(0, \sigma_j^2) \quad (11)$$

<footnote n="2">²In this way, (9) is verified, and we get the $h_{1D}$ filter in at least 1 Direction</footnote>
where $\mathcal{N}$ is the normal distribution,

$$ \sigma_j^2 = \frac{1}{4} \frac{\tau_2(h^{[j]})}{(\tau_1(h^{[j]})^2)} + \frac{1}{4} \frac{\tau_2(h'^{[j+1]})}{(\tau_1(h'^{[j+1]})^2)} - \frac{\langle h^{[j]}, h'^{[j+1]} \rangle}{2\tau_1(h^{[j]})\tau_1(h'^{[j+1]})} $$

and

$$ \langle h^{[j]}, h'^{[j+1]} \rangle = \sum_{k_x, k_y, k_z} h^{[j]}(k_x, k_y, k_z) h'^{[j+1]}(k_x, k_y, k_z) $$

$$ \tau_p(h^{[j]}) = \sum_{k_x, k_y, k_z} (h^{[j]}(k_x, k_y, k_z))^p $$

$$ \tau_p(h'^{[j+1]}) = \sum_{k_x, k_y, k_z} (h'^{[j+1]}(k_x, k_y, k_z))^p $$

The symbol $\xrightarrow{\mathcal{D}}$ means the convergence in probability 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'^{[j+1]}$ are defined by ($\delta$ is either x, or y, or z):

$$ h_{\delta}^{[j]} = h_{\delta}^{(0)} * h_{\delta}^{(1)} * \dots * h_{\delta}^{(j)}, \quad h^{[j]} = h_x^{[j]} h_y^{[j]} h_z^{[j]} $$

and

$$ h'_{\delta}^{[j+1]} = h_{\delta}^{[j+1]} * h_{\delta}^{(j+1)}, \quad h'^{[j+1]} = h_x'^{[j+1]} h_y'^{[j+1]} h_z'^{[j+1]} $$

This function is called a Variance Stabilization Transform, and it is a method that decorrelates noise from signal. In our case, the MPG Noise becomes Gaussian. Equation (11) works with any $c_j$, $c'_{j+1}$. The authors determined the best $c_j$ and $c'_{j+1}$ so that the convergence is faster: $c_j = \frac{7}{8} \frac{\tau_2(h^{[j]})}{\tau_1(h^{[j]})^2} - \frac{1}{2} \frac{\tau_3(h^{[j]})}{\tau_2(h^{[j]})}$ and $c'_{j+1} = \frac{7}{8} \frac{\tau_2(h'^{[j+1]})}{\tau_1(h'^{[j+1]})^2} - \frac{1}{2} \frac{\tau_3(h'^{[j+1]})}{\tau_2(h'^{[j+1]})}$. They experimentally discovered that the convergence was reached for very low values of $v(i)$ (i.e. for $v(i) \ge 0.1$).

All these results are true if $\mathbb{E}(d_j(i)) = 0$. Thus we can now define a function for hypothesis testing: we want to test whether $\mathbb{E}(d_j(i))$ is null. If $\mathbb{E}(d_j(i))$ is null then (11) says that $l_{i,j} = \sqrt{(a_0 * h^{[j]})(i) + c_j - \sqrt{(a_0 * h'^{[j+1]})(i) + c'_{j+1}}}$ follows a normal distribution with mean 0 and known variance. We ensure then that, under the hypothesis $\mathbb{E}(d_j(i)) = 0$,

$$ \mathbb{P}(|X| < |l_{i,j}|) < FPR \quad (12) $$

FPR being the false positive rate, a parameter of the hypothesis test, and $\mathbb{P}(|X| < |l_{i,j}|) = 1 - \text{erf}(\frac{|l_{i,j}|}{\sqrt{2}\sigma_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 distribution. If (12) is verified then we accept the hypothesis, 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 $\mathcal{M}$ by

$$ \begin{aligned} \mathcal{M} &= \{ (i, j), d_j(i) \text{ not significantly null} \} \\ &= \left\{ (i, j), \text{erf} \left( \frac{|l_{i,j}|}{\sqrt{2}\sigma_j} \right) > 1 - FPR \right\} \end{aligned} $$

## 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 $\mathcal{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 $\mathcal{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. $\mathcal{M}$ can be then redefined as follows:

$$ \mathcal{M} = \left\{ (i, j), \text{erf} \left( \frac{|l_{i,j}|}{\sqrt{2}\sigma_j} \right) > 1 - FPR \text{ and } J_m \le j \le J_M \right\} \quad (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 Figure 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 parameter of the Poisson distribution) are non-negative. This problem suggests a more subtle reconstruction model.

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 reconstructed images are non-negative. To do so, they establish a cost which constrains the number of positive coefficients added to the set of significant coefficients, to keep it as small as possible. If $s$ is the IUWT representation of the denoised image, the cost is $\|s\|_1 = \sum_{i,j} |s_j(\hat{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(\hat{i}) = d_j(\hat{i})$ for all tuple $(\hat{i}, j)$ in $\mathcal{M}$.

The final model is the solution of the following optimization problem:

$$ \text{minimize } \|s\|_1 \text{ subject to } \begin{cases} s_j(\hat{i}) = d_j(\hat{i}) \text{ for all } \hat{i} \text{ in } \mathcal{M} \\ (Rs)(\hat{i}) \ge 0 \text{ for all } \hat{i} \end{cases} $$

where $\mathcal{R}$ is the reconstruction operator defined by (3) and given in Figure 3b.

The solution of this model is obtained by some advanced iterative optimization algorithm described in (Chambolle and Pock, 2011) which can be adapted to our problem as shown in Algorithm (1). This algorithm, 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 $\mathcal{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 operator is defined as follows:

$$ (\text{ST}_{\tau}(s))_j(\hat{i}) = \max(|s_j(\hat{i})| - \tau, 0)\text{sign}(s_j(\hat{i})) \quad (14) $$

The $\mathcal{P}_{\mathcal{M}}$ operator constrains its argument $s$ to be equal to $d$ for the tuples $(\hat{i}, j)$ in $\mathcal{M}$:

$$ (\mathcal{P}_{\mathcal{M}}(s))_j(\hat{i}) = \begin{cases} d_j(\hat{i}) & \text{if } (\hat{i}, j) \in \mathcal{M} \\ s_j(\hat{i}) & \text{otherwise.} \end{cases} $$

The $\mathcal{P}_{\mathcal{O}^{-}}$ is the projection operator on the set of arrays with non-positive elements:

$$ (\mathcal{P}_{\mathcal{O}^{-}}(u))(\hat{i}) = \begin{cases} u(\hat{i}) & \text{if } u(\hat{i}) \le 0 \\ 0 & \text{otherwise.} \end{cases} $$

Finally the $\mathcal{R}^*$ operator is given by the following formula:

$$ (R^*u)_0 = u \\ (R^*u)_{j+1} = (R^*u)_j * h^{(j)} $$

We chose this algorithm because of its very fast convergence (experimentally it converges within 50 iterations).

**Data:** Noisy Image $u$, $r_s$, $J_m$ and $J_M$, False Positive rate FPR

**Result:** Denoised Image $\tilde{v}$. Set $\mathcal{M}$ as defined in (13);

$$ \tau \leftarrow \sqrt{\frac{1}{r_s(J_M+2)^2}}; \\ \sigma \leftarrow \sqrt{\frac{r_s}{(J_M+2)^2}}; \\ u_0 \leftarrow u; \\ (s_0)_j(\hat{i}) \leftarrow 0 \text{ for all } \hat{i} \text{ and } j; \\ \bar{s}_0 \leftarrow s_0; $$

**repeat**

$$ \left| \begin{array}{l} u_{k+1} \leftarrow \mathcal{P}_{\mathcal{O}^{-}}(u_k + \sigma \mathcal{R} \bar{s}_k); \\ s_{k+1} \leftarrow \mathcal{P}_{\mathcal{M}}(\text{ST}_{\tau}(s_k - \tau \mathcal{R}^* u_{k+1})); \\ \bar{s}_{k+1} \leftarrow 2s_{k+1} - s_k; \\ k \leftarrow k+1; \end{array} \right. $$

**until convergence;**

$$ \tilde{v} \leftarrow \mathcal{R} s_{\infty}; $$

**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 image; 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 embryo” to have a common ground truth. We tested all three methods on two different experiments. In the first experiment, the methods were prevented from removing the background and we measured the resulting signal-to-noise ratio with respect to the noiseless ground-truth image containing the cells and the background. 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 distributed 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 certain radius. In each center, we installed a full ellipse with random size, random orientation, and a deformed 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
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 compared to the groundtruth image. The main parameter is the FPR parameter explained in Section 3.1, ranging 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$ defined 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 preset before applying the denoising method. To do so, one can directly filter the desired scale from the image 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 reconstruction operator $\mathcal{R}$, which is itself a function of the size of the image, the dimensions of its voxels, and the number 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 image 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 compare the different results. Figure 4 shows the signal-to-noise ratio (SNR) of the output image after applying each denoising method. SNR is defined by $-10 \log_{10} (\|\tilde{v} - v\|_2^2 / \|v\|_2^2)$ (where $v$ is the ground truth and $\tilde{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 explained 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 detector” software (Olivo-Marin, 2002) which detects nuclei in the denoised image. This software actually relies on IUWT to detect cells in 3D images. As it was run on artificial images, we already know the true positions and number of cells. We set the parameter of all the methods as follows<fnref n="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 default 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 pixels (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 better results compared to the other two.

## 5 CONCLUSION

In this paper, we introduced a method that tackles noise problems in cell images coming from two-photon microscopy. More than just denoising images 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 stabilization 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 results 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 resulting images. We expect better results with the fine-

<footnote n="3">³for the scale parameters, the outcome does not change much if we choose $J_m = 1$</footnote>
<figure><img src="image_5.png" /><figcaption>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 anisotropy, a 3D method without anisotropy, the 2D method. In blue/red/yellow: three different scale ranges.</figcaption></figure>

<table><caption>Table 1: Number of false detections on 2 example images after applying the 3 denoising methods</caption><thead><tr><th>Method</th><th>3D-Grid Anisotropy</th><th>3D-Grid Isotropy</th><th>2D</th></tr></thead><tbody><tr><td>Image 1 (SNR<sub>0</sub> = 2.68)</td><td>1</td><td>1</td><td>1</td></tr><tr><td>Image 2 (SNR<sub>0</sub> = 0.075)</td><td>7</td><td>14</td><td>17</td></tr></tbody></table>

tuning of the parameters of our method, which is the goal of future work.

This opens up more perspectives, such as comparing our method to other methods not using IUWT. Another interesting working direction would be to combine cell detection and denoising. Indeed the cell detection 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é Doursat (BioEmergences 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 microscopy image sequences. *Medical Imaging, IEEE Transactions on*, 29(2):442–454.

* Cai, L. (1989). Spline smoothing: A special case of diffusion smoothing. In *Proceedings of the Alvey Vision 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é, N., Pop, S., Provoost, T., Meas-Yedid, V., Pankajakshan, 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ébarre, 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 recognition*, 35(9):1989–1996.

* Smith, G. D. (1985). *Numerical solution of partial differential equations: finite difference methods*. Oxford university press.

* Zanella, C., Campana, M., Rizzi, B., Melani, C., Sanguinetti, G., Bourgine, P., Mikula, K., Peyrieras, N., and Sarti, A. (2010). Cells segmentation from 3-d confocal 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.