<figure><img src="image_1.png" /></figure>

HAL
open science

# Image denoising in fluorescence microscopy using feature based gradient reconstruction

Suman Kumar Maji, Hussein Yahia

► To cite this version:

Suman Kumar Maji, Hussein Yahia. Image denoising in fluorescence microscopy using feature based gradient reconstruction. Journal of Medical Imaging, 2023, 10 (6), 10.1117/1.JMI.10.6.064004. hal-04339260

HAL Id: hal-04339260

https://inria.hal.science/hal-04339260v1

Submitted on 18 Dec 2023

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Image Denoising in Fluorescence Microscopy using Feature based Gradient Reconstruction

Suman Kumar Majiª, Hussein Yahiaᵇ

**Abstract.**

**Purpose:** The utility of fluorescence microscopy imaging comes with the challenge of low resolution acquisitions which severely limit information extraction and quantitative analysis. Image denoising is a technique that aims to remove noise from microscopy acquisitions by taking into account prior statistics of the corrupting noise. In this paper, we propose an image denoising technique for fluorescence microscopy imaging.

**Approach:** The proposed technique is based on the principle of multifractal feature extraction from a noisy sample followed by a reconstruction technique from these features. It is observed that by following a proper hierarchical classification procedure, meaningful features can be extracted from a noisy image. A denoised image is then estimated from this sparse feature set through proper formulation of an optimization problem.

**Results:** Experiments are performed on both synthetic image databases as well as on real fluorescence microscopy data. Superior denoising results, in comparison to multiple comparing techniques, validate the potential of the proposed approach.

**Conclusion:** The proposed method gives superior denoising results for low resolution fluorescence microscopy image acquisitions and can be used for post processing of data by biologists.

**Keywords:** Image Denoising, Fluorescence Microscopy, Poisson-Gaussian Noise, Singularity exponents.

*Suman Kumar Maji, smaji@iitp.ac.in*

## 1 Introduction

Over the decade, fluorescence microscopy has become a biologists tool of choice for *in vivo* analysis of cellular structures. However, due to various physical constraints like **short exposure time**, photo-toxicity, photo-bleaching, device imperfections, etc. the resulting image is a low resolution noisy acquisition. This seriously hampers further analysis like cell segmentation,¹ recognition,² reconstruction,³ etc of the acquired samples. Noise in fluorescent microscopy acquired samples is characterized by a **mixed Poisson-Gaussian (MPG) statistics**.⁴⁻⁸ Image denoising is a technique that specifically addresses noise removal in a sample through proper knowledge of the noise source and its distribution.

<footnote>ªDepartment of Computer Science & Engineering, Indian Institute of Technology Patna, 801106 Patna, India</footnote>

<footnote>ᵇTeam Geostat, INRIA Bordeaux Sud-Ouest, 200 rue de la Vieille Tour, 33405 Talence Cedex, France</footnote>
Multiple denoising strategies have been proposed in the literature for improving the quality of the low resolution fluorescent microscopy samples. These strategies range from Bayesian based formulations,⁹,¹⁰ wavelet based techniques,¹¹,¹² non-local patch based techniques⁴,¹³ to variance stabilization based approaches,¹⁴,¹⁵ to name a few. Bayesian based techniques rely on deriving an appropriate noise restoration model using Bayes principle, taking into account the statistics of the nature of the noise. In this regard, the denoising model in⁹ is particularly remarkable as it derives a convex optimization model using Bayes formulation, for fluorescent microscopy samples, based on the Poisson-Gaussian likelihood. Wavelet based denoising schemes operate on manipulating the wavelet coefficients within the framework of forward and inverse wavelet transform. The SURE-LET¹² approach is one such technique. It is an orthonormal wavelet thresholding based denoising scheme that takes into account the interscale dependencies between the wavelet coefficients for removing noise. This scheme has been extended to a PURE-LET¹¹ scheme, where in the authors have generalized the concept of Haar wavelet transform to fit into a series of linear expansion of thresholds. Patch based non-local denoising techniques, like,⁴ work on the principle of exploiting patch-based redundancy within a given sample. The NDSAFIR technique⁴ is one such method, where the authors take into account the spatio-temporal image patches of observed fluorescent samples for minimizing an objective non-local energy functional. Stein's unbiased risk estimator (SURE) is a similarity criteria that is often used in image denoising for predicting optimal solutions. In PG-URE⁵ the authors extend the concept of SURE to address the case of Poisson-Gaussian noise removal in fluorescent samples. Another popular category of image denoising technique are the variance stabilization (VST) based techniques.¹⁴,¹⁵ VST based methods make use of appropriate transformations in order to convert the signal dependent noisy component into signal independent noise. VST+BM3D¹⁵ and GAT+BM3D¹⁴ are two such popular denoising schemes
for fluorescent microscopy samples that make use of the Anscombe transform as a stabilizer for
the signal-dependent Poisson noise.

Deep learning based denoising techniques have also become popular in recent years.<sup>16,17</sup> A
content-aware image restoration based on deep learning is proposed in<sup>18</sup> for denoising fluorescence
microscopy samples. A three dimensional deep neural network based residual channel attention de-
noising network for fluorescence microscopy is proposed in.<sup>19</sup> In addition to these, self-supervised
denoising networks that work on the principle of learning only from the single noisy image is also
widely used for fluorescence microscopy denoising. The Self2Self<sup>20</sup> model is one such network
that uses Bernoulli sampling as a dropout strategy for denoising images. The Neighbor2Neighbor
(Nb2Nb)<sup>21</sup> is another approach that uses a random neighbor <u>up-sampler</u> to generate training data
for effective denoising.

The denoising technique that we propose in this paper is based from ideas encountered in the
analysis of multifractal systems.<sup>22-25</sup> These are scale-free complex systems whose characteristic
quantities display a power-law behaviour and are characterized by values called singularity expo-
nents (SE's).<sup>26</sup> Natural scale-free complex systems exhibit multifractal traits<sup>27</sup> implying thereby
that signals acquired from such systems can be characterized by multifractals, an approach to
complex systems which has been validated over a large dataset.<sup>22</sup> The denoising scheme that
we propose is a two-step concept where in the first step we extract a sparse feature map through
multifractal analysis on the acquired noisy image. In the second step we build an optimization
problem to estimate a denoised image from the gradients confined to the extracted feature map.

The contribution of this paper are as follows:

* We formulate a novel denoising algorithm that uses multifractal analysis to extract meaning-
ful features from a noisy image, and then reconstruct a denoised image from these features.

Under the hypothesis of a multifractal system, the integrated signal's gradient of the original image has a power-law behaviour whose exponent at each point $x$ of the signal's domain, $h(x)$, known as a singularity exponent (SE), can be computed using a local correlation measure.

* We recall the computation of the local correlation measure and of the set of lowest exponents, $\mathcal{L}_{\infty}$, which is of fundamental significance in signal processing as it allows the reconstruction of a new signal from the SEs.

* We design a reconstruction algorithm to estimate a denoised image from the feature set $\mathcal{L}_{\infty}$.

* Experimental results over three image databases, under different levels of noise, show that the proposed method outperforms both classical denoising techniques as well as recent deep learning techniques in the field.

The rest of the paper is organised as follows: In section 2.1 we mathematically formulate the image acquisition problem and objective and briefly summarize the multifractal analysis of acquired samples. Here in subsection 2.3 we propose a technique for computing the quantity $h(x)$ required for estimating the feature set $\mathcal{L}_{\infty}$ and in subsection 2.4 we propose a method for estimating a denoised image from $\mathcal{L}_{\infty}$. In section 3 we present the experiment results and discuss them and finally conclude in section 4.
## 2 Proposed Method

### 2.1 Preliminaries & Objective

If $f: \Omega \to \mathbb{R}$ is the acquired noisy image defined on the image domain $\Omega \subset \mathbb{R}$ and $u: \Omega \to \mathbb{R}$ its noise free version, the image formation process for fluorescence microscopy can be written as:⁵,¹¹

$$f = \tau(u) + w \qquad (1)$$

where $\tau(\cdot) \sim P(\cdot|\zeta)$ represents the signal dependent noisy component with $\zeta \ge 0$ being the gain of the overall image acquisition system. $w \sim N(b, \sigma^2)$ represents the signal independent additive noise component with mean $b$ and variance $\sigma^2$. $\sigma, \zeta \ge 0$ determines the intensity of Poisson-Gaussian noise in $f$.⁵ ($\sigma, \zeta$) were empirically determined keeping into account the fact that higher the $\sigma$ value, higher the Gaussian noise. Similarly higher the $\zeta$ value, higher is the Poisson noise. Noise is generated with different sets of values of $(\sigma, \zeta)$, over the ground truth, following the above equation. We intend to recover an image $\hat{u} \approx u$ from the noisy fluorescent microscopy sample $f$.

### 2.2 Multifractal Analysis and Feature Extraction

As discussed earlier, natural scale-free systems exhibit self-similarity and multifractal behavior implying therefore that signals (or images) acquired through such systems will exhibit similar traits. If $f(x)$ (where $x$ is a vector representing the location of a pixel), which in our case is the acquired noisy image, is a signal with multifractal properties then it can be used to define a measure with multi-scale properties,²²⁻²⁴,²⁸ $\mu(x, r)$ ($x$: pixel location, $r$ scale) such that we have
the limiting behaviour:

$$h(\mathbf{x}) = \lim_{r \to 0} \frac{\log \mu(\mathbf{x}, r)}{\log r} \qquad (2)$$

which in turn implies that

$$\mu(\mathbf{x}, r) \approx \alpha(\mathbf{x}) r^{h(\mathbf{x})} \qquad (3)$$

where $\alpha(\mathbf{x})$ is a signal dependent constant. $h(\mathbf{x})$ is known as the singularity exponents or SE's. The transitions within a signal are well-recorded in $h(\mathbf{x})$,²⁶ which are in turn independent of the measure $\mu$, as long as it is multiscale. There are different ways to compute the $h(\mathbf{x})$, once $\mu$ is defined properly. The most evident algorithm consists in noting that equation (3) implies that, by taking the logarithms on both sides:

$$\log(\mu(\mathbf{x}, r)) = \log(\alpha(\mathbf{x})) + h(\mathbf{x}) \log(r) \quad (r \to 0) \qquad (4)$$

and then perform log-log regression for many different scales $r$. However, the accuracy of such method depends on the number of scales $r$ chosen for the regression; it has been shown²⁸ (and the references herein) that a most advanced methodology consists in evaluating $h(\mathbf{x})$ by use of a local correlation measure. This is justified by the existence of a local reconstruction of the whole signal at the points where $h(\mathbf{x})$ at the lowest values.²⁸ We postpone this computation to section 2.3 below. For the moment, supposing that $h(\mathbf{x})$ has been properly calculated, they are used to define a multifractal hierarchy of multiple fractals $\mathcal{L}_h$ as level sets: $\mathcal{L}_h = \{\mathbf{x} | h(\mathbf{x}) = h\}$,²⁵ from where it is possible to extract a particular set $\mathcal{L}_\infty$ consisting of components associated to the smallest
**Algorithm 1** Input: $h(\mathbf{x})$, Output: $\mathcal{L}_{\infty}$

**Step 1:** Compute the multifractal exponents $h(\mathbf{x})$ from Algorithm 2.

**Step 2:** Define density

$$ \delta_{\mathcal{F}_{\infty}}(\mathbf{x}) = \begin{cases} 1 & \text{if } h_{\infty} - \Delta h \le h(\mathbf{x}) \le h_{\infty} + \Delta h \\ 0 & \text{otherwise} \end{cases} $$

**Step 3:** $\Delta h$ is conventionally fixed (0.2 is usually a good choice).

**Step 4:** $\delta_{\mathcal{L}_{\infty}}(\mathbf{x})$ is the binary map that locates the $\mathcal{L}_{\infty}$ pixels.

possible value $h_{\infty}$ defined as:

$$ \mathcal{L}_{\infty} = \{\mathbf{x} | h(\mathbf{x}) = h_{\infty} = \min(h(\mathbf{x}))\} \qquad (5) $$

The $h_{\infty}$'s only encode the pixels where it observes sharp transitions in the image i.e., the boundaries of objects inside the image, and maintains its stability under noise. This can be visualised in Fig. 1 (middle row). The procedure for computing the set $\mathcal{L}_{\infty}$ has been listed in algorithm 1.

## 2.3 Local correlation measure for the computation of $h(\mathbf{x})$

As explained in the last subsection, a local correlation measure can be used to evaluate the SEs $h(\mathbf{x})$, instead of log-regression across the scales (which necessitates to define various multiscale versions of the original signal). The key observation lies in the existence of the set $\mathcal{L}_{\infty}$, which allows a local reconstruction of the signal from the most singular exponents (see section 2.4 and$^{28}$).

In this study, we propose a new local correlation measure, well adapted the noise reduction problem we have in sight.

In this paper we propose a method to compute $h(\mathbf{x})$ at every pixel location $\mathbf{x} = (i, j)$ in the image taking into consideration their 4 neighbouring pixels: denoting $f(i, j)$ the gray-level value at pixel $\mathbf{x} = (i, j)$, the 4 neighbouring pixel values form the set $\{f(i - 1, j), f(i + 1, j), f(i, j - 1), f(i, j + 1)\}$. In an image of size $N \times M$ ($N, M$: numbers of pixels in the $x$ and $y$ di-
<figure><img src="image_3.png" /><figcaption>Fig 1: **Top:** Noisy realizations. **Middle:** $\mathcal{L}_{\infty}$ computed over them and **bottom:** Proposed denoising results.</figcaption></figure>

rections respectively), and assuming the whole image of surface 1, the minimum resolution is
$r_0 = \frac{1}{\sqrt{N \times M}}$. The complete procedure to compute $h(\mathbf{x})$ at resolution $r_0$ is given in Algorithm 2 (where $*$ denotes convolution and $T$ stands for transposition). **Note that at minimum fixed resolution $r_0$, $\mu(\mathbf{x}, r_0)$ depends only on $\mathbf{x}$.**

**Algorithm 2** Input: Image $f$, Output: Singularity exponent $h(\mathbf{x})$ at $\mathbf{x} = (i, j)$

**Step 1:** $\epsilon_x = f \star [1 \ 0 \ -1])$; $\epsilon_y = f \star -[1 \ 0 \ -1]^T$, $\epsilon_x, \epsilon_y$ correspond to unit pixel displacements at maximum scale $r_0$.

**Step 2:** $e(\mathbf{x}) = \epsilon_x^2 + \epsilon_y^2$.

**Step 3:** $\mu(\mathbf{x}, r_0) = \sqrt{e(\mathbf{x})}$.

**Step 4:** $h(\mathbf{x}) = \frac{\log \mu(\mathbf{x}, r_0)}{\log r_0}$.

## 2.4 Reconstructing the Denoised Image

In this section, we propose a novel formulation to reconstruct a denoised image $\hat{u} \approx u$ from the sparse feature set $\mathcal{L}_{\infty}$. Such a reconstruction comes from the possibility on integrating gradients over the set $\mathcal{L}_{\infty}$ (see$^{25,28}$). However in this section we will introduce a new reconstruction algorithm which makes use of a L2-L1 minimization with sparse gradients instead of reconstruction in
<figure><img src="image_4.png" /><figcaption>Fig 2: Result over synthetic data corrupted with MPG noise. Top row: $\sigma, \zeta = 10^{-2.5}, 10^{-3}$ Middle row: $\sigma, \zeta = 10^{-1}, 10^{-2}$ and Bottom row: $\sigma, \zeta = 10^{-2}, 10^{-1.5}$.</figcaption></figure>

<table><caption>Table 1: Quantitative evaluation using PSNR (refer equation (11)) and SSIM (refer equation (12)) metric over standard image dataset. Best cases are highlighted in bold.</caption><thead><tr><th colspan="2"></th><th colspan="2">Param</th><th colspan="7">Mean PSNR</th></tr><tr><th></th><th></th><th>σ</th><th>ζ</th><th>Input</th><th>PG-URE</th><th>SURE</th><th>VST+</th><th>GAT+</th><th>Self2Self</th><th>Nb2Nb</th><th>Proposed</th></tr><tr><th></th><th></th><th></th><th></th><th></th><th colspan="2">-LET</th><th>BM3D</th><th>BM3D</th><th></th><th></th><th></th></tr></thead><tbody><tr><td rowspan="4">Classic5</td><td rowspan="4">Set12</td><td>10<sup>-2.5</sup></td><td>10<sup>-3</sup></td><td>24.56</td><td>33.33</td><td>32.90</td><td>32.19</td><td>32.01</td><td>30.20</td><td>30.02</td><td><strong>34.91</strong></td></tr><tr><td>10<sup>-1</sup></td><td>10<sup>-2</sup></td><td>18.48</td><td>25.04</td><td>23.84</td><td>23.90</td><td>25.32</td><td>24.78</td><td>24.14</td><td><strong>26.84</strong></td></tr><tr><td>10<sup>-2</sup></td><td>10<sup>-1.5</sup></td><td>13.88</td><td>22.824</td><td>21.96</td><td>20.73</td><td>21.53</td><td>20.19</td><td>20.87</td><td><strong>24.90</strong></td></tr><tr><td>10<sup>-2.5</sup></td><td>10<sup>-3</sup></td><td>25.99</td><td>33.92</td><td>31.09</td><td>30.79</td><td>30.42</td><td>29.85</td><td>29.64</td><td><strong>35.46</strong></td></tr><tr><td rowspan="4">BSD68</td><td rowspan="4">Set12</td><td>10<sup>-1</sup></td><td>10<sup>-2</sup></td><td>18.80</td><td>25.16</td><td>24.14</td><td>25.18</td><td>25.38</td><td>24.64</td><td>24.48</td><td><strong>26.45</strong></td></tr><tr><td>10<sup>-2</sup></td><td>10<sup>-1.5</sup></td><td>13.45</td><td>22.41</td><td>20.48</td><td>22.82</td><td>22.57</td><td>21.48</td><td>21.98</td><td><strong>24.67</strong></td></tr><tr><td>10<sup>-2.5</sup></td><td>10<sup>-3</sup></td><td>25.37</td><td>31.23</td><td>28.25</td><td>30.22</td><td>29.77</td><td>28.35</td><td>28.85</td><td><strong>33.89</strong></td></tr><tr><td>10<sup>-1</sup></td><td>10<sup>-2</sup></td><td>17.39</td><td>23.27</td><td>22.61</td><td>23.64</td><td>23.78</td><td>22.37</td><td>22.13</td><td><strong>24.87</strong></td></tr><tr><td rowspan="4">Classic5</td><td rowspan="4">Set12</td><td>10<sup>-2</sup></td><td>10<sup>-1.5</sup></td><td>0.417</td><td>0.657</td><td>0.622</td><td>0.644</td><td>0.653</td><td>0.620</td><td>0.624</td><td><strong>0.688</strong></td></tr><tr><td>10<sup>-2.5</sup></td><td>10<sup>-3</sup></td><td>0.768</td><td>0.934</td><td>0.924</td><td>0.918</td><td>0.905</td><td>0.893</td><td>0.896</td><td><strong>0.954</strong></td></tr><tr><td>10<sup>-1</sup></td><td>10<sup>-2</sup></td><td>0.588</td><td>0.784</td><td>0.695</td><td>0.771</td><td>0.770</td><td>0.707</td><td>0.701</td><td><strong>0.811</strong></td></tr><tr><td>10<sup>-2</sup></td><td>10<sup>-1.5</sup></td><td>0.437</td><td>0.662</td><td>0.547</td><td>0.654</td><td>0.645</td><td>0.631</td><td>0.635</td><td><strong>0.694</strong></td></tr><tr><td rowspan="4">BSD68</td><td rowspan="4">Set12</td><td>10<sup>-2.5</sup></td><td>10<sup>-3</sup></td><td>0.748</td><td>0.929</td><td>0.891</td><td>0.905</td><td>0.911</td><td>0.898</td><td>0.904</td><td><strong>0.943</strong></td></tr><tr><td>10<sup>-1</sup></td><td>10<sup>-2</sup></td><td>0.568</td><td>0.765</td><td>0.700</td><td>0.755</td><td>0.748</td><td>0.702</td><td>0.701</td><td><strong>0.798</strong></td></tr><tr><td>10<sup>-2</sup></td><td>10<sup>-1.5</sup></td><td>0.444</td><td>0.654</td><td>0.588</td><td>0.648</td><td>0.632</td><td>0.611</td><td>0.621</td><td><strong>0.688</strong></td></tr></tbody></table>

Fourier space as in.²⁵,²⁸ This choice is motivated by our primary goal of noise reduction, for which sparse methods have been shown to be particularly effective (²⁹⁻³¹).

Let $f_x, f_y$ be the gradients of the noisy image $f$ and $\delta_{\mathcal{L}_{\infty}}$ be the binary map of the set $\mathcal{L}_{\infty}$. We intend to recover $u$ from the sparse gradient set $v = \{\hat{f}_x, \hat{f}_y\}$, where $\hat{f}_x = f_x \cdot \delta_{\mathcal{L}_{\infty}}$ and $\hat{f}_y = f_y \cdot \delta_{\mathcal{L}_{\infty}}$. Since $\mathcal{L}_{\infty}$ encodes the most meaningful pixels of the image, the gradients corresponding to $\mathcal{L}_{\infty}$ will be the most informative. A reconstruction from the set $v$ will therefore lead to a proper estimation
of the whole image and noise propagation will be automatically minimal since a lot of noisy pixels are dropped. Keeping this in mind, we therefore propose to minimize the following optimization problem:

$$ \hat{u} = \underset{u}{\operatorname{argmin}} \frac{1}{2} \|\nabla u - \mathbf{v}\|_2^2 + \lambda \|\nabla u\|_1 \quad (6) $$

where $\lambda$ is a positive regularization parameter. With Half-Quadratic splitting, using an intermediate variable $z$ and an additional positive regularization term $\beta$, we get:

$$ \hat{u} = \underset{u,z}{\operatorname{argmin}} \frac{1}{2} \|\nabla u - \mathbf{v}\|_2^2 + \lambda \{\|z\|_1 + \frac{\beta}{2} \|\nabla u - z\|_2^2\} \quad (7) $$

which is further organised into independent sub-problems $p_1$ and $p_2$:

$$ p_1: z^{(k+1)} \leftarrow \underset{z}{\operatorname{argmin}} \|z\|_1 + \frac{\beta}{2} \|\nabla u^{(k)} - z\|_2^2 $$

$$ p_2: u^{(k+1)} \leftarrow \underset{u}{\operatorname{argmin}} \|\nabla u - \mathbf{v}\|_2^2 + \lambda \|\nabla u - z^{(k+1)}\|_2^2 $$

Problem $p_1$ is in the proximal form $\underset{z}{\operatorname{argmin}} \|z\|_1 + \frac{\beta}{2}\|x - z\|_2^2$ whose solution correspond to a soft thresholding shrinkage operator:³⁰

$$ z = \operatorname{shrink}_{l_p}(x, \beta) = \max \left\{ 0, |x| - \frac{|x|^{p-1}}{\beta} \right\} \frac{|x|}{x} \quad (8) $$

Considering $z$ as a 2-component vector, the solution to sub-problem $p_1$ is given as:

$$ z^{(k+1)} = \begin{cases} z_x^{(k+1)} = \operatorname{shrink}_{l_1}(\nabla_x u^{(k)}, \beta) \\ z_y^{(k+1)} = \operatorname{shrink}_{l_1}(\nabla_y u^{(k)}, \beta) \end{cases} \quad (9) $$
<table><caption>Table 2: Quantitative evaluation using PSNR (refer equation (11)) and SSIM (refer equation (12)) metric on FMD datasets. Best cases are highlighted in bold.</caption><thead><tr><th rowspan="2">Sample</th><th colspan="8">Mean PSNR</th></tr><tr><th>Raw</th><th>PG-URE</th><th>SURE</th><th>VST+</th><th>GAT+</th><th>DnCNN</th><th>N2N</th><th>Proposed</th></tr><tr><th></th><th></th><th></th><th></th><th>-LET</th><th>BM3D</th><th>BM3D</th><th></th><th></th></tr></thead><tbody><tr><td>Confocal Zebrafish</td><td>22.76</td><td>30.63</td><td>27.75</td><td>31.46</td><td>30.86</td><td>32.24</td><td>33.04</td><td><strong>35.66</strong></td></tr><tr><td>Twophoton Mice Brain</td><td>22.37</td><td>33.46</td><td>28.38</td><td>34.25</td><td>33.22</td><td>34.58</td><td>34.91</td><td><strong>36.06</strong></td></tr><tr><td>Widefield BPAE</td><td>20.06</td><td>27.05</td><td>28.05</td><td>27.73</td><td>27.78</td><td>28.76</td><td><strong>29.69</strong></td><td>29.11</td></tr></tbody></table><table><thead><tr><th rowspan="2">Sample</th><th colspan="8">Mean SSIM</th></tr><tr><th>Raw</th><th>PG-URE</th><th>SURE</th><th>VST+</th><th>GAT+</th><th>DnCNN</th><th>N2N</th><th>Proposed</th></tr></thead><tbody><tr><td>Confocal Zebrafish</td><td>0.427</td><td>0.834</td><td>0.856</td><td>0.883</td><td>0.888</td><td>0.898</td><td>0.911</td><td><strong>0.947</strong></td></tr><tr><td>Twophoton Mice Brain</td><td>0.278</td><td>0.907</td><td>0.798</td><td>0.915</td><td>0.908</td><td>0.922</td><td>0.924</td><td><strong>0.949</strong></td></tr><tr><td>Widefield BPAE</td><td>0.759</td><td>0.885</td><td>0.927</td><td>0.892</td><td>0.911</td><td>0.914</td><td><strong>0.968</strong></td><td>0.944</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Fig 3: Result over FMD dataset: single-channel (gray) image of zebrafish embryo under confocal microscopy</figcaption></figure>

Problem $p_2$ is in the regular quadratic form and can be solved using Euler-Lagrange. The final solution using FFT will be:

$$u^{(k+1)} = \mathbf{F}^{-1} \left\{ \frac{\mathbf{F}(\operatorname{div}(\mathbf{v}) + \operatorname{div}(z^{(k+1)}))}{(1 + \lambda)\Delta} \right\} \quad (10)$$

where $\mathbf{F}$ denotes Fourier transform, $\Delta = \nabla^2$ the Laplace operator. The results of reconstruction, under different levels of noise, are shown in Fig. 1 (bottom row) and Fig. 2.
<figure><img src="image_6.png" /><figcaption>Fig 4: Result over FMD dataset: single-channel (gray) image of mice brain under twophoton microscopy</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig 5: Result over FMD dataset: multi-channel (color) image of BPAE cells under widefield microscopy</figcaption></figure>

3 Experiment & Results

We have used two types of data: Synthetic data and fluorescence microscopy acquisition data. For comparison purpose we have considered the four popular classical denoising techniques designed for removing MPG noise: SURE-LET,¹¹ PG-URE,⁵ VST+BM3D¹⁵ and GAT+BM3D.¹⁴ We have also used two recent deep learning based self-supervised denoising methods that learns only from the input noisy image, thereby making them ideal for comparison. They are: Self2Self²⁰ and Nb2Nb.²¹ We have also used DnCNN³² and Noise2Noise³³ models for comparison on Fluorescence Microscopy Denoising (FMD) dataset.⁸ The parameter settings for all these methods have been kept the same as suggested by their authors in the papers and in the codes. For the purpose of
quantitative evaluation, we have considered the peak signal-to-noise ratio (PSNR) and structural similarity index (SSIM) metrics, defined as follows:

$$ \text{PSNR} = 10 \log_{10} \frac{u_{\text{max}}^2}{(\|\hat{u} - u\|^2/(m \times n))} \qquad (11) $$

where $u_{\text{max}}$ is the maximum intensity of the noise-free image $u$ and $m \times n$ the size of the image $u$.

$$ \text{SSIM}(u, \hat{u}) = \frac{(2\mu_u\mu_{\hat{u}} + c_1)(2\sigma_{u\hat{u}} + c_2)}{(\mu_u^2 + \mu_{\hat{u}}^2 + c_1)(\sigma_u^2 + \sigma_{\hat{u}}^2 + c_2)} \qquad (12) $$

where $\mu_u, \mu_{\hat{u}}$ are the average of $u$ and $\hat{u}$; $\sigma_u^2, \sigma_{\hat{u}}^2$ are the variance of $u$ and $\hat{u}$; $\sigma_{u\hat{u}}$ the covariance of $u$ and $\hat{u}$ respectively; $c_1$ and $c_2$ are two normalizing parameters.

## 3.1 Synthetic Data

We have considered three image datasets: the Classic5, Set12 and BSD68 which consists of 5, 12 and 68 synthetic images respectively. We have generated multiple noisy versions of these datasets, by adding MPG noise of level $\{\sigma, \zeta\} = \{10^{-2.5}, 10^{-3}\}, \{10^{-1}, 10^{-2}\}$ and $\{10^{-2}, 10^{-1.5}\}$, over which we perform denoising. In Fig. 2, we present the visual comparison results of our technique with the competing methods. As can be seen from the excerpts (highlighted in green bounding boxes), SURE-LET is capable of retaining the details to a certain extent, but introduces artefacts in the denoised output. PG-URE, VST + BM3D and GAT+BM3D does better denoising than SURE-LET, but are mostly accompanied by over smoothing which leads to loss in details. For the cameraman image one can see that the background building behind the tripod, in the case of boat image the paint linings corresponding to the boat name and in the case of foreman image the black spots in the building behind are best reconstructed by the proposed method. The proposed
technique is therefore able to retain the fine details of the image while rendering better visual
quality than the others in terms of noise removal. This is also justified quantitatively by the results
of the mean PSNR and mean SSIM, calculated over all the images of the datasets, under different
levels of noise, as is shown in table. 1.

3.2 Fluorescence Microscopy Data

The dataset that we have used is the fluorescence microscopy denoising dataset (FMD).<sup>8</sup> The vi-
sual results of denoising are shown in fig. 3, fig. 4 and fig. 5 respectively. The ground truth images
shown in these figures are obtained (as described by the authors) image averaging over 50 realiza-
tions of the noisy image. In fig. 3 the denoising results of a single channel (gray) image of zebrafish
embryo captured under confocal microscopy is shown. The image is of size 512 × 512 pixels. It can
be observed from the zoomed section that the details of the embryo are best preserved in the output
of the proposed method. Noise2Noise and VST+BM3D also tries to retain the details to a certain
extent, but are not that accurate in comparison to the output of the proposed method. PG-URE and
GAT+BM3D gives a blurred impression whereas SURE-LET and DnCNN gives oversmoothed
results. Fig. 4 shows the denoising results of a single channel (gray) image of mice brain captured
under two photon microscopy. The image is of size 512 × 512 pixels. Careful observation shows
that the proposed method better reconstructs the details and sharpness of the white spots as well as
the blob like structures inside the cells, than the Noise2Noise model which is slightly smoothened
in its output. DnCNN, VST+BM3D and PG-URE also gives better denoising but has smoothened
outputs compared to the proposed model and Noise2Noise. The rest gives over smoothened outputs.
Quantitative evaluation of fig. 3 and fig. 4, as shown in table 3, also justifies the superior denoising
results of the proposed method.
Fig. 5 shows the denoising results over a multi channel (color) image of bovine pulmonary artery endothelial (BPAE) cells captured under widefield microscopy. The image is of size 512 x 512 x 3 pixels. It is clear from the results that the proposed method gives the best denoising while retaining the finer details. The zoomed section shows that the filament like structures are well retained as well as the circular background patch. The rest of the denoising method fails to retain such detail and sharpness. The color contrast is also improved in the denoising results. It is to be noted that the FMD dataset ground truth images are obtained by image averaging over 50 realizations of the noisy image. The ground truth is therefore also contaminated with noise, albeit less than the noisy acquisition. For a high noise acquisition like BPAE (which has multiple noisy channels) this is more pronounced than compared to single channel acquisitions of zebrafish and mice brain. PSNR and SSIM of Noise2Noise (N2N) is highest for BPAE because N2N has been trained on the FMD dataset (where it has considered the less noisy ground truth as the clean image during training) and has therefore recovered the low noise also in its reconstruction. The proposed method (as well as PG-URE, SURE-LET, VST+BM3D, GAT+BM3D and DnCNN) cleans the background noise completely. However its PSNR and SSIM values will be lower because it mis-matches with the low noisy ground truth. The ground truth of BPAE is therefore not suitable for quantitative analysis using PSNR and SSIM.

### 3.3 Run Time

In the table below, we provide the running time for all the comparing methods for real data all of whose image size are 512 x 512 pixels. It is to be noted that PG-URE, SURE-LET, VST+BM3D, GAT+BM3D and Proposed methods are Matlab based methods and utilizes the CPU. We tested it on a standard core i5 8th Gen laptop with 8GB RAM. DnCNN and N2N are deep learning models
<table><caption>Table 3: Running time comparison</caption><thead><tr><th rowspan="3">Sample</th><th colspan="6">Running Time (in seconds)</th><th colspan="2">GPU Running Time</th></tr><tr><th colspan="5">CPU Running Time</th><th rowspan="2">DnCNN</th><th rowspan="2">N2N</th></tr><tr><th>PG-URE</th><th>SURE-LET</th><th>VST+BM3D</th><th>GAT+BM3D</th><th>Proposed</th></tr></thead><tbody><tr><td>Widefield BPAE</td><td>16.56</td><td>4.54</td><td>8.23</td><td>4.15</td><td>3.9</td><td>0.19</td><td>0.08</td></tr></tbody></table>

and utilize the GPU, for which they are faster but computationally expensive in terms of resource utilization.

Algorithm 1 always achieves real-time. Half quadratic splitting depends of the size of the input image. For a 10000 × 10000 pixel image, GPU implementation of Half quadratic splitting can be reduced down to 1 second cpu time on a high level desktop (we tested it on Nvidia Quadro NVS 290 with Cuda).

## 4 Conclusion

This paper presents a denoising technique for recovering low resolution fluorescence microscopy samples. The denoising technique is motivated from the concept of multifractal feature extraction from an observed noisy sample, and subsequent reconstruction from the gradient information limited to those features. The feature extraction method records only the most informative pixels and drops the remaining pixels. A denoised image is then reconstructed from the gradients of this feature set by solving a half-quadratic based optimization problem. Experiments were performed on both synthetic image databases and real fluorescence microscopy data. Results, both qualitative and quantitative, have proved the superior denoising potential of the proposed technique over multiple competing techniques in the field.
## Disclosures

The authors have no relevant financial interests in the manuscript and no other potential conflicts of interest to disclose.

## Code, Data, and Materials Availability

The datasets analysed during the current study are available in the CBSD68 repository, https://paperswithcode.com/dataset/cbsd68 and FMD repository.⁸

## References

1 A. Gharipour and A. W.-C. Liew, "Segmentation of cell nuclei in fluorescence microscopy images: An integrated framework using level set segmentation and touching-cell splitting," *Pattern Recognition* **58**, 1-11 (2016).

2 S. Almasi, A. Ben-Zvi, B. Lacoste, et al., "Joint volumetric extraction and enhancement of vasculature from low-snr 3-d fluorescence microscopy images," *Pattern Recognition* **63**, 710-718 (2017).

3 K. Kolev, N. Kirchgeßner, S. Houben, et al., "A variational approach to vesicle membrane reconstruction from fluorescence imaging," *Pattern Recognition* **44**(12), 2944-2958 (2011).

4 J. Boulanger, C. Kervrann, P. Bouthemy, et al., "Patch-based non-local functional for denoising fluorescence microscopy image sequences," *IEEE Trans. on Medical Imaging* **29**(2) (2010).

5 Y. L. Montagner, E. D. Angelini, and J. C. Olivo-Marin, "An unbiased risk estimator for image denoising in the presence of mixed poisson-gaussian noise," *IEEE Trans. Image Process.* **23**, 1255-68 (2014).
6 S. K. Maji, C. Dargemont, J. Salamero, et al., “Joint denoising-deconvolution approach for fluorescence microscopy,” *Proc. IEEE ISBI* (2016).

7 S. Maji and J. Boulanger, “A variational model for poisson gaussian joint denoising deconvolution,” *Proc IEEE Int Symp Biomed Imaging*, 1527–1530 (2021).

8 Y. Zhang, Y. Zhu, E. Nichols, et al., “A poisson-gaussian denoising dataset with real fluorescence microscopy images,” in *Proceedings of the IEEE conference on computer vision and pattern recognition*, (2019).

9 E. Chouzenoux, A. Jezierska, J.-C. Pesquet, et al., “A convex approach for image restoration with exact poisson–gaussian likelihood,” *SIAM Journal on Imaging Sciences* 8(4), 2662–2682 (2015).

10 Q. Gao, S. Eck, J. Matthias, et al., “Bayesian joint super-resolution, deconvolution, and denoising of images with poisson-gaussian noise,” *Proc. IEEE Int. Symp. on Biomed. Imag. (ISBI)*, 938–942 (2018).

11 F. Luisier, T. Blu, and M. Unser, “Image denoising in mixed poisson-gaussian noise,” *IEEE Trans. Image Process.* **20**, 696–708 (2011).

12 T. Blu and F. Luisier, “The sure-let approach to image denoising,” *IEEE Trans. Image Process.* **16**, 2778–2786 (2007).

13 A. Baudes, “A non-local algorithm for image denoising,” *Proceedings of the IEEE conference on computer vision and pattern recognition*, 60–65 (2005).

14 M. Makitalo and A. Foi, “Optimal inversion of the generalized anscombe transformation for poisson-gaussian noise,” *IEEE Trans. Image Process.* **22**, 91–103 (2013).
15 L. Azzari and A. Foi, "Variance stabilization for noisy+estimate combination in iterative poisson denoising," *IEEE Signal Process. Lett.* **23**, 1086–1090 (2016).

16 C. Belthangady and L. A. Royer, "Applications, promises, and pitfalls of deep learning for fluorescence image reconstruction," *Nature Methods* **16**, 1215–1225 (2019).

17 E. Moen, D. Bannon, T. Kudo, et al., "Deep learning for cellular image analysis," *Nature Methods* **16**, 1233–1246 (2019).

18 M. Weigert, U. Schmidt, T. Boothe, et al., "Content-aware image restoration: pushing the limits of fluorescence microscopy," *Nature Methods* **15**, 1090–1097 (2018).

19 J. Chen, H. Sasaki, and e. a. H. Lai, "Three-dimensional residual channel attention networks denoise and sharpen fluorescence microscopy image volumes," *Nature Methods* **18**, 678–687 (2021).

20 Y. Quan, M. Chen, T. Pang, et al., "Self2self with dropout: Learning self-supervised denoising from single image," in *Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition*, 1890–1898 (2020).

21 T. Huang, S. Li, X. Jia, et al., "Neighbor2neighbor: Self-supervised denoising from single noisy images," in *Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition*, 14781–14790 (2021).

22 A. Turiel and N. Parga, "The multifractal structure of contrast changes in natural images: From sharp edges to textures," *Neural Computation* **12**, 763–793 (2000).

23 A. Turiel, H. Yahia, and C. J. Pérez-Vicente, "Microcanonical multifractal formalism -a geometrical approach to multifractal systems: Part 1. singularity analysis," *J. Phys. A: Math. Theor.* **41**(015501) (2008).
24 A. Turiel, C. J. Pérz-Vicente, and J. Grazzini, "Numerical methods for the estimation of multifractal singularity spectra on sampled data: A comparative study," *Journal of Computational Physics* **216**, 362–390 (2006).

25 H. Badri, H. Yahia, and K. Daoudi, "Fast and accurate texture recognition with multilayer convolution and multifractal analysis," *European Conf. Comp. Vision* (2014).

26 S. K. Maji and H. M. Yahia, "Edges, transitions and criticality," *Pattern Recognition* **47**, 2104–2115 (2014).

27 K. Falconer, *Techniques in fractal geometry*, John Wiley (1997).

28 Yahia, H., Schneider, N., Bontemps, S., et al., "Description of turbulent dynamics in the interstellar medium: multifractal-microcanonical analysis - i. application to herschel observations of the musca filament," *A&A* **649**, A33 (2021).

29 R. Chartrand, "Exact reconstruction of sparse signals via nonconvex minimization," in *IEEE Signal Processing Letters*, **14**(10), 707–710 (2007).

30 R. Chartrand, "Fast algorithms for nonconvex compressive sensing: Mri reconstruction from very few data," *Proc. IEEE ISBI* (2009).

31 G. Gasso, A. Rakotomamonjy, and S. Canu, "Recovering sparse signals with a certain family of non-convex penalties and dc programing," in *IEEE Transactions on Signal Processing*, **57**(12), 486–498 (2009).

32 K. Zhang, W. Zuo, Y. Chen, et al., "Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising," *IEEE Transactions on Image Processing* **26**(7), 3142–3155 (2017).
33 V. Mannam, Y. Zhang, Y. Zhu, et al., “Instant image denoising plugin for imagej using convolutional neural networks,” *Biophotonics Congress: Biomedical Optics 2020 (Translational, Microscopy, OCT, OTS, BRAIN)*, MW2A.3, Optica Publishing Group (2020).

**Suman Kumar Maji** received his B.Tech. degree in Electronics & Communication Engineering from West Bengal University of Technology, India in 2006, Postgraduate degree in Telecommunication Networks from the Indian Institute of Technology Kharagpur, India in 2008 and Ph.D. degree in Computer Science from INRIA Bordeaux France in 2013.

From 2014 to 2015 he worked as a Research Engineer with the Institute of Hematology, University Paris 7 and INSERM. He is currently working as Assistant Professor with the Department of Computer Science and Engineering, Indian Institute of Technology Patna, India. His research interests are in the area of Medical Imaging, Bioinformatics, Machine Learning and Image processing.

Dr. Maji has authored several conferences and journal papers and is the recipient of various research fellowships and awards such as the European CORDIS Doctoral Fellowship (2010), Region Aquitaine OPTAD Research Fellowship (2010), FRM Research Fellowship (2014), SERB Early Career Research Award from DST, Govt of India (2017), etc.

**Hussein Yahia** (M'87) received the Doctorat de 3eme cycle degree from the University of Paris-Sud, Orsay, France, in 1987, and the Habilitation à Diriger des Recherches (HDR) degree from Paris 13 University, Villetaneuse, France, in 2003.

He is currently the Head of the Geostat (Geometry and Statistics in Acquisition Data) Research Team with the French National Public Research Institute in Computer Science and Applied Mathematics (INRIA), Bordeaux, France. He specializes in the analysis of complex signals and systems
using approaches from statistical physics; he is also studying sparse signal representations and fil-
tering with specific application to signals in astrophysics in collaboration with the Astrophysics
Laboratory of Bordeaux (analysis of turbulent properties in the interstellar medium). H.Yahia is
involved in many national, European, and international contracts and has been supervising over
15 PhD. students. He has authored or coauthored more than 100 publications in international
peer-reviewed top journals and conferences.

**List of Figures**

1 **Top:** Noisy realizations. **Middle:** $\mathcal{L}_{\infty}$ computed over them and **bottom:** Proposed denoising results.

2 **Result over synthetic data corrupted with MPG noise.** Top row: $\sigma, \zeta = 10^{-2.5}, 10^{-3}$
Middle row: $\sigma, \zeta = 10^{-1}, 10^{-2}$ and Bottom row: $\sigma, \zeta = 10^{-2}, 10^{-1.5}$.

3 **Result over FMD dataset: single-channel (gray) image of zebrafish embryo**
**under confocal microscopy**

4 **Result over FMD dataset: single-channel (gray) image of mice brain under**
**twophoton microscopy**

5 **Result over FMD dataset: multi-channel (color) image of BPAE cells under**
**widefield microscopy**

**List of Tables**

1 Quantitative evaluation using PSNR (refer equation (11)) and SSIM (refer equation (12)) metric over standard image dataset. Best cases are highlighted in bold.
2 Quantitative evaluation using PSNR (refer equation (11)) and SSIM (refer equation (12)) metric on FMD datasets. Best cases are highlighted in bold.

3 Running time comparison