# SPITFIR(e): A supermaneuverable algorithm for fast denoising and deconvolution of 3D fluorescence microscopy images and videos

Sylvain Prigent, Hoai-Nam Nguyen, Ludovic Leconte, Cesar Augusto Valades-Cruz, Bassam Hajj, Jean Salamero, Charles Kervrann

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

## To cite this version:

Sylvain Prigent, Hoai-Nam Nguyen, Ludovic Leconte, Cesar Augusto Valades-Cruz, Bassam Hajj, et al. SPITFIR(e): A supermaneuverable algorithm for fast denoising and deconvolution of 3D fluorescence microscopy images and videos. Scientific Reports, 2023, 13, pp.1-34. 10.1038/s41598-022-26178-y. http://dx.doi.org/10.1038/s41598-022-26178-y

HAL Id: hal-03518618

https://inria.hal.science/hal-03518618v2

Submitted on 14 Jan 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 Attribution 4.0 International License
# SPITFIR(e): A supermaneuverable algorithm for fast denoising and deconvolution of 3D fluorescence microscopy images and videos

Sylvain Prigent<sup>1,2,+</sup>, Hoai-Nam Nguyen<sup>1,2,+</sup>, Ludovic Leconte<sup>1,2</sup>, Cesar Augusto Valades-Cruz<sup>1,2</sup>, Bassam Hajj<sup>3</sup>, Jean Salamero<sup>1,2</sup>, and Charles Kervrann<sup>*</sup>,<sup>1,2</sup>

## ABSTRACT

Modern fluorescent microscopy imaging is still limited by the optical aberrations and the photon budget available in the specimen. A direct consequence is the necessity to develop flexible and "off-road" algorithms in order to recover structural details and improve spatial resolution, which is critical when restraining the illumination to low levels in order to limit photo-damages. Here, we report SPITFIR(e) a flexible method designed to accurately and quickly restore 2D-3D fluorescence microscopy images and videos (4D images). We designed a generic sparse-promoting regularizer to subtract undesirable out-of-focus background and we developed a primal-dual algorithm for fast optimization. SPITFIR(e) is a "swiss-knife" method for practitioners as it adapts to any microscopy techniques, to various sources of signal degradation (noise, blur), to variable image contents, as well as to low signal-to-noise ratios. Our method outperforms existing state-of-the-art algorithms, and is more flexible than supervised deep-learning methods requiring ground truth datasets. The performance, the flexibility, and the ability to push the spatiotemporal resolution limit of sub-diffracted fluorescence microscopy techniques are demonstrated on experimental datasets acquired with various microscopy techniques from 3D spinning-disk confocal up to lattice light sheet microscopy.

## Introduction

Fluorescence microscopy provides a very powerful framework to biologists for observing, analyzing, and studying specific fluorescently-tagged structures and biological processes at very high spatial and temporal resolutions. Despite number of advantages, there are two major limitations of fluorescence microscopy. The first limitation is the presence of photon (shot) noise in acquired images. Shot noise is mainly due to the quantum nature of light, implying that the arrival of a photon on a sensor is a random event and thus the number of incident photons over a period of time is a random variable depending on the brightness of the light source. Moreover, in biological samples the signal-to-noise-ratio (SNR) is usually very low because low dose of illumination light is required to avoid photo-bleaching (i.e., progressive fading of the emission intensity) of fluorescent molecules and to preserve specimen integrity (photo-toxicity)<sup>1,2</sup>. Additionally, the quality of acquired fluorescence images is worsened by the blurring effects induced by several factors such as excitation wavelength, immersion medium refraction, specimen thickness, and the limited aperture of the objective which results in light diffraction through the optical system. The diffraction phenomenon implies that light emitted by an infinitely small point source appears wider at the focal plan and spreads into a specific pattern called "point spread function" (PSF).

Noise and blur not only degrade the image quality in terms of overall visualization but also have a negative influence on specimen analysis, including detection and segmentation. To improve the quality of images acquired by fluorescence microscopes, restoration (deconvolution and/or denoising) is then frequently applied as pre-processing step before quantitative analysis, and amounts to recovering the original unknown image *u* from the observed noisy and blurred image *f* represented as follows: $f = \mathcal{T}(h*u)$, where *h* denotes the 2D or 3D spatial response of the device representing the blur related to the optical system (e.g., PSF) assumed to be linear shift-invariant, * denotes the convolution operator, and $\mathcal{T}$ is a degradation operator modeling the measurement noise. If the image is not blurred but mainly corrupted by noise, *h* is represented by the Dirac (or impulse) function and the restoration problem translates into a denoising problem. Finally, the degradation operator (shot noise

<footnote n="1"><sup>1</sup>SERPICO Project-Team, Inria Centre Rennes-Bretagne Atlantique, 35042, Rennes Cedex, France</footnote>

<footnote n="2"><sup>2</sup>SERPICO/STED Team, UMR144 CNRS Institut Curie, PSL Research University, Sorbonne Universités, 75005, Paris, France</footnote>

<footnote n="3"><sup>3</sup>Laboratoire Physico-Chimie, Institut Curie, PSL Research University, Sorbonne Universités, CNRS UMR168, 75005, Paris, France</footnote>

<footnote>*charles.kervrann@inria.fr</footnote>

<footnote>+these authors contributed equally to this work</footnote>
and readout noise) is usually modeled as a mixed Poisson-Gaussian distribution which is non-stationary and signal-dependent.

In the last twenty years, many restoration methods have been investigated in order to "inverse" the model in 2D-3D fluorescence microscopy. The most popular restoration approaches are linear methods (e.g., Wiener filtering). Despite the simplicity and low-computation-requirement, they usually do not restore fine image details at frequency components that are beyond the bandwidth of the PSF (i.e., the support of its Fourier transform). Another popular technique is the Richardson-Lucy (RL) algorithm³⁻⁵. While well-dedicated to Poisson noise removal, RL deconvolution amplifies and creates structural noise and artifacts after a small number of iterations, which constitutes a major problem in fluorescence imaging⁶. To reduce these drawbacks and iteration number, several modern and fast RL implementations has been recently proposed⁷,⁸. Actually, the most flexible and robust deconvolution methods consist in minimizing an energy functional composed of a data fidelity term and a regularization term that encompasses prior knowledge (positivity, smoothness, sparsity,...) about the solution. The seminal Tikhonov-Miller (TM) approach⁹⁻¹¹ may be considered as the starting point to the development of regularization methods in fluorescence microscopy. While it was successful in many applications, the TM regularizer tends to produce blurred images since low gradient magnitudes are encouraged in the entire image, including at image discontinuities. To avoid over-smoothing caused by quadratic functionals, non-quadratic regularizers have been first studied, especially the Total Variation (TV) regularizer that penalizes the $L_1$ norm of the gradient. Nevertheless, TV creates "stair-casing" effects, which are particularly undesirable in fluorescence imaging¹²,¹³. To overcome this disadvantage, Lefkimmiatis et al.¹⁴ recently proposed a family of convex regularizers built upon the matrix norm of the Hessian (Schatten norm), computed at each point of the image. This piece-wise smoothness promoting regularizer combined with a non-quadratic data term, was recently evaluated in fluorescence microscopy¹⁵. Nevertheless, this regularizer does not impose sparsity prior on the fluorescent signals which is also an important feature in fluorescence imaging. This idea was investigated by Arigovindan et al.¹⁶ who reported good results 3D wide-field fluorescence images by using a non-convex sparsity-promoting regularizer; this regularizer (named here Log Sparse Hessian Variation (LHSV)) has been implemented into the entropy-regularized deconvolution algorithm (ER-Decon). Concurrently, other sparsity-promoting methods, focused on first-order differentiation-based regularizers have been investigated in medical imaging such as GraphNet (GN)¹⁷⁻¹⁹ and Sparse Variation (SV)²⁰, and in image processing (e.g., TV-$L_1$²¹,²²). Nevertheless, the sparsity-promoting regularizers reported in Table 2 have limitations. They were designed for 2D medical imaging¹⁷⁻²⁰ or 3D microscopy¹⁶,²³ but never for both. Moreover, they are based on the first-order¹⁷⁻²⁰ or second derivatives¹⁶,²³, and are formulated as non-convex¹⁶ or convex²³ optimization problems. While more tricky, non-convex energies are now efficiently minimized by performant algorithms such as ADMM²⁴,²⁵ and FISTA²⁶.

To overcome the main limitations of previous methods, we introduce a fast, supermaneuverable method (SPITFIR(e)) to deconvolve and/or denoise 2D or 3D multispectral fluorescence microscopy images and videos. Our method outperforms existing algorithms in three ways. (1) SPITFIR(e) includes a generic sparsity promoting regularizer to discard out-of-focus background, and adapts to most if not all fluorescence microscopy techniques even if the 3D PSF is not well known. (2) It reaches high restoration performance for variable levels and sources of noise and blur. It outperforms state-of-the-art deconvolution methods and is competitive when compared to supervised deep-learning methods. (3) The SPITFIR(e) parameters are automatically calibrated from the image contents and the practitioner needs only to specify if the desired image is "highly", "moderately" or "weakly" sparse. A "weakly" sparse image contains clutter and complex contents. We validated the overall performance of our approach on several challenging images acquired with lattice light sheet (LLSM), stimulated emission depletion (STED), multifocus microscopy (MFM), and spinning-disk confocal microscopy (CM). For the evaluation on different imaging techniques, we have chosen to apply SPITFIR(e) on a limited number of structures relatively difficult to optically resolve or for which the photon budget is critical, i.e., mitochondria, cytoskeleton elements and biomolecules associated to moving vesicles in living cells. We also developed CPU and GPU implementations which can be downloaded as open-source from the GitHub website (https://github.com/sylvainprigent/spitfire) under a AGPL-3.0 license. Finally SPITFIR(e) is interfaced with napari for ease of use (https://www.napari-hub.org/plugins/napari-sdeconv).

# Results

## Overview of SPITFIR(e)

SPITFIR(e) is an algorithm well-grounded in the regularization theory for inverse problems, that can robustly recover information from general multidimensional noisy and blurry fluorescence images. Unlike the traditional methods requiring the tedious manual adjustment of the parameter that balances the data fidelity term and the prior term, SPITFIR(e) includes a self-tuning and scale-invariant technique for automatic adaptation (Methods). To our knowledge, there are few methods able to select such parameters in the literature to process complex 3D microscopy images. Here, we mention an alternative approach²⁷ that exploits a single image and subsampling operations to compute Fourier Ring Correlation and further to determine restoration algorithm parameters and assess the most appropriate spatial resolution. In our approach, the practitioner can modify the default parameter setting and encourage more sparsity and/or smoothness from a preview of a 3 × 3 matrix depicting the restoration results on a 2D ROI with three levels – "weak", "moderate", "high" – of sparsity and smoothness.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** SPITFIR(e) overview. The flow chart describes the sequence of restoration steps applied to 3D, 2D+Time, or 3D+Time images based on responses of end-users.</figcaption></figure>

The minimization of the sum of the two convex terms is originally implemented in a computationally efficient way with a recent version of proximal-splitting optimization algorithm<sup>28,29</sup> (see Methods and Appendix). In what follows, we show that SPITFIR(e) can restore 3D fluorescence images under low excitation light excitation conditions by chaining denoising and deconvolution with the same algorithm with very small computing times compared to other competitive methods<sup>15,23</sup>.

**SPITFIR(e) is supermaneuverable and adapts to every fluorescence microscopy modality**

SPITFIR(e) is said "supermaneuverable" as it includes several strategies to adapt to microscopy specificities, and to particular spatial and temporal acquisition conditions (see flowchart in Fig. 1). For instance, the practitioner can apply the conventional 3D deconvolution strategy if the PSF is perfectly known. If this is not the case, a "plane by plane" (PbyP) strategy which consists in separately deconvolving each XY section of the 3D stack given a 2D Gaussian PSF model, can be applied. Surprisingly, the PbyP strategy could produce better visual results and is faster than direct 3D deconvolution with an imperfect or approximately measured 3D PSF. Furthermore, in the case of low-photon regimes or low-exposure times, the results may be significantly improved by applying the "denoise-before-deconvolve" approach. For instance, a 3D+Time (4D) volume can be denoised as a whole, while the planes of each 3D stack are independently deconvolved with a 2D Gaussian PSF model. An exemplary demonstration of the flexibility of SPITFIR(e) is given by the experiments performed with a lattice light sheet microscope<sup>30</sup> (LLSM) which allows to acquire 3D images for several hours with limited photobleaching. In LLSM imaging, an illumination angle makes the raw image stack skewed, post-processing is required to create a 3D "deskewed" stack<sup>30</sup>. Consequently, deconvolving 3D LLSM images is not straightforward as the PSF is not well characterized and the deskew operation affects the noise statistics. We then design the following SPITFIR(e) deconvolution pipeline according to the flowchart illustrated in Fig. 1. We hypothesized here that chaining the denoising and deconvolution steps could provide better restorations results. So, first we denoise the 3D+Time image sequence (4D denoising) with SPITFIR(e) (the PSF function *h* is a Dirac function). Then we deconvolve each 3D stack separately with SPITFIR(e) by using an anisotropic 3D Gaussian PSF model. The settings of the PSF are chosen given the image resolution along the XY and Z axes.

We first illustrate the efficiency of the "denoise-before-deconvolve" approach using microtubules and mitochondria internal membrane double labeling (Fig. 2, Live-cell LLSM). Live imaging mitochondria and their interactions with other intracellular components at high spatial resolution in 3D and high temporal frequency, is known to be challenging. Beyond the usual optical limits, labeled mitochondria are very sensitive to light<sup>31</sup>, which could induce artefactual fluorescence signals such as exaltation, probes dissociation and finally swelling and fragmentation of mitochondria and cell death. This may impair accurate 3D+Time observation of their behaviors in normal and stress conditions. For these reasons, imaging mitochondria dynamics and structural features will serve as one of the biological thread in the following parts of the manuscript.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** SPITFIR(e) on LLSM images.</figcaption></figure>

**Live-cell LLSM.** 60 planes 3D volumes of live RPE1 cells double stained with PKMR for Mitochondria (a and with Tubulin Tracker Deep Red for Microtubules (b) were acquired within 2.2s per stack using LLSM. Maximum intensity projection (MIP) of representative raw images labeled for Mitochondria and Microtubules, respectively, before (a, b; upper part) and after (middle part) SPITFIR(e) 4D denoising + 3D deconvolution (3D Gaussian PSF, $\sigma_{xy} = 1.5$ pixels and $\sigma_z = 1.0$ pixel) and after (lower part) 3D deconvolution (Richardson-Lucy (RL) algorithm with a 3D Gaussian PSF, $\sigma_{xy} = 1.5$ pixels and $\sigma_z = 1.0$; 10 iterations). Insets are zoomed area illustrating SPITFIR(e) improvement in spatial resolution and SNR as compared to RL processing. 3D angular views and MIP of composite images (red for Microtubules; green for Mitochondria) before (c, d; upper part) and after SPITFIR(e) (middle part) or RL (lower part). Comparative illustrations of 4D denoising SPITFIR(e), 4D denoising+3D deconvolution SPITFIR(e) and RL 3D deconvolution image improvement shown on a single X-Y plane (e) zoomed from the ring insets indicated in the three images in (d). Raw image is indicated as a blue lined circle (e; upper part). 4D denoised image is similarly indicated in magenta, while the 4D denoised + 3D deconvolved image is in yellow, and 3D RL deconvolution in gray. Intensity line profiles (e; lower part) were measured as indicated in the 4 circles (e; upper part) and plotted (e; lower part) for each processing step and both Mitochondria and Microtubules staining. Scale bars are indicated in the bottom right corners (10 µm). 3D+time series related to Raw and 4D denoised + 3D deconvolved SPITFIR(e) data are shown as S1 Video and S2 Video.
**Figure 2 (previous page)**

**Fixed-cell LLSM.** RPE1 cells labeled with SiR-Actin were fixed before immunostaining for α-tubulin (f, g, h). 60 planes 3D volumes were acquired within 2 s per stack using LLSM. Representative merged MIP images of actin filaments (red) and microtubules (green) staining for raw images, before (f), after Richardson-Lucy (RL) deconvolution (g, 3D Gaussian PSF, σ<sub>xy</sub> = 1.5 pixels and σ<sub>z</sub> = 1.0 pixel and 10 iterations) and after SPITFIR(e) 3D denoising and 3D deconvolution (h, 3D Gaussian PSF, σ<sub>xy</sub> = 1.5 pixels and σ<sub>z</sub> = 1.0 pixel). Zoomed area of a single X-Y plane (f-h, middle and right images; scale bars = 5 µm) was extracted from the insets (as indicated in f, g, h; left images) for qualitative comparison between SPITFIR(e) and RL. Fourier Image Resolution (FIRE) was estimated using Fourier Ring Correlation Plugin<sup>32</sup> for multiple single X-Y planes before (i), after RL 3D deconvolution (j) and after SPITFIR(e) 3D denoising + 3D deconvolution (k). SPITFIR(e) improvement was estimated separately for 33 X-Y single planes of microtubules from 6 RPE1 cells acquired using LLSM. Statistics analysis is represented in (l). It indicates the improvement of average FIRE (unpaired two-sample t-test).

We imaged live RPE1 cells in full volume labeled for mitochondria (Fig. 2(a)); Raw images, upper part) and microtubules (Fig. 2(b)); Raw images, upper part). We here linked 3D LLSM at low illumination regimes to efficiently preserve fluorescence from photobleaching and consequently from phototoxic effects, with 4D denoising and 3D image deconvolution by using SPITFIR(e) (Fig. 2(a, b)); middle images) and a comparison with 3D deconvolution using Richardson-Lucy (RL-Deconvolution) approach (Fig. 2(a, b)); bottom images). Results are shown as single labeling (Fig. 2(a, b)) and as merged images (Fig. 2(c, d)), with different angles of view.

A reconstituted comparison of the successive steps of the SPITFIR(e) workflow, 4D denoised followed by 3D deconvolution (see Supplementary Fig. S1) and 3D RL deconvolution is shown in Fig. 2(e). For the sake of visual clarity, only one single plane of the zoomed area of a stack (as indicated in Fig. 2(d; upper, middle and bottom images) is used, and shown as zoomed insets in Fig. 2(e; upper part) from raw image to full SPITFIR(e) processing and compared to 3D RL deconvolution). More visual representations can be appreciated in supplementary Videos S1 and S2, in 3D+Time along both XY and Z axes. Interestingly, effect of chaining denoising and deconvolution with SPITFIR(e) is particularly visible on the microtubule staining which efficiency is low, while improvement in resolution is more obvious on mitochondria (compare insets of zoomed area between Fig. 2(b) upper and middle images and between Fig. 2(a) upper and middle images). As expected, SPITFIR(e) does not outperform RL improvement in resolution for high SNR images (green for mitochondria), while clearly improving both image quality and resolution for low SNR (red, microtubules). Line intensity plots through the center of the fluorescent filamentous structures (Fig. 2(e); upper part) demonstrate the recovery of the signal from the noise (Fig. 2(e); bottom part). For quantitative measurements of the performance of SPITFIR(e) vs RL using LLSM images we double labeled fixed cells for two filamentous intracellular structures, the microtubules (green) and the actin filaments (red). Merged images and zoomed extracts are shown on Fig. 2(f-h, Fixed-cell LLSM). We computed the Fourier Image Resolution (FIRE) using Fourier Ring Correlation Plugin<sup>32</sup> for multiple single X-Y planes before (Fig. 2(i)), after RL 3D deconvolution (Fig. 2(j)) and after SPITFIR(e) 3D denoising + 3D deconvolution (Fig. 2(k)). For both labeled structures SPITFIR(e) using 3D denoising before 3D deconvolution improves the resolution as compared to RL deconvolution, again with a more pronounced effect in the lowest SNR fluorescent channel (compare green and red curves in Fig. 2(i-k)). This clear improvement in the lowest SNR fluorescent channel is statistically relevant (33 X-Y planes, from 6 different cells) as confirmed by the results shown in Fig. 2(l).

With SPITFIR(e), we have also improved the resolution limits of current STED microscopy. STED microscopy provides subdiffraction resolution while preserving useful aspects of fluorescence microscopy, such as optical sectioning, and molecular specificity and sensitivity. The PSF shape depends on several factors such as optics, laser power, fluorescence response of the specimen, and thermal drift which can add deformation by sheering. STED nanoscopy is one of the very few light microscopy techniques<sup>33,34</sup> allowing to solve mitochondria cristae organization. STED microscopy provides subdiffraction resolution while preserving useful aspects of fluorescence microscopy, such as optical sectioning, and molecular specificity and sensitivity. The PSF shape depends on several factors such as optics, laser power, fluorescence response of the specimen, and thermal drift which can add deformation by sheering. In Fig. 3, we deconvolved with SPITFIR(e) two 2D images depicting live mitochondria and their alignment on microtubules, by applying SPITFIR(e) 2D with a 2D Gaussian PSF model and two different bandwidths σ<sub>xy</sub> of size 1.5 pixels and 1.0 pixel, respectively (Fig. 3(a)). Importantly, we tested and obtained very similar results with a Lorentzian 2D PSF with such bandwidth values (not shown as the differences are nearly indistinguishable), mainly because the PSF sizes are small in these experiments corresponding to few pixels when implemented in the discrete setting. We selected two regions of interest (ROI) depicting mitochondria cristae to better appreciate the restoration results in (Fig. 3(a, b)). Normalized fluorescence intensity line profiles, positioned at a distinct cristae region of mitochondria (Fig. 3(b) and Supplementary Fig. S2) show the improvement of SPITFIR(e) in terms of SNR and image quality, with no visual artifact and noise amplification. Two simulated PSFs sizes were proposed to evaluate the capacity of SPITFIR(e) to improve the image resolution. The first PSF (σ<sub>xy</sub> = 1.0 pixel) was similar than the expected experimental PSF of our microscope and the second PSF (σ<sub>xy</sub> = 1.5 pixels) corresponds to the situation where there is an issue of defocusing. Clearly in both cases we observe an image improvement of the two individual cristae line profiles. The fitting of a Gaussian model to intensity profiles (n = 75,
<figure><img src="image_5.png" /></figure>

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

<table><thead><tr><th rowspan="2"></th><th colspan="2">Mean [µm]</th><th rowspan="2">t test<br>p-val</th><th colspan="2">STD [µm]</th><th rowspan="2">F test<br>p-val</th></tr><tr><th>STED 2D</th><th>σ = 1 pixel</th><th>STED 2D</th><th>σ = 1 pixel</th></tr></thead><tbody><tr><td>All</td><td>0.0869</td><td>0.0725</td><td>0.0072 [**]</td><td>0.0402</td><td>0.022</td><td>&lt;0.0001 [****]</td></tr><tr><td>FWHM<sub>raw</sub> &gt; 60nm</td><td>0.1019</td><td>0.08</td><td>0.0003 [***]</td><td>0.0374</td><td>0.0203</td><td>&lt;0.0001 [****]</td></tr></tbody></table>

**Figure 3.** SPITFIR(e) increases the signal-to-noise of low-light images in STED imaging. Live RPE1 cells double stained with PKMO (mitochondria) and Tubulin Tracker Deep Red (microtubules) were imaged by STED nanoscopy before and after being processed by SPITFIR(e) (a). Partial overview of 2D STED raw data of RPE1 cells for mitochondria, microtubules and their superimposed images are shown from left to right (a, upper panel). 2D STED images were improved by SPITFIR(e) using two different 2D Gaussian model $\sigma_{xy} = 1.0$ pixel (a, middle panel) and $\sigma_{xy} = 1.5$ pixels (a; lower panel). Insets (1 and 2) show zoomed area where SPITFIR(e) image quality improvement is illustrated by comparison of two processing conditions (middle and lower panel) with the same area in the raw STED image. Inset(b) shows a magnified composite image where the 2D STED raw part is indicated as a blue lined rectangle and 2D denoised + 2D deconvolved parts are indicated in red ($\sigma_{xy} = 1.0$ pixel) and green ($\sigma_{xy} = 1.5$ pixels). (b) Normalized intensity line profiles were measured in the insets (b) for cristae regions. Yellow line indicated in the inset (b) serve to identify fluorescence profiles. Pixel size is equal to 25 nm. Scale bars: 1 µm in (a) and (b). Line profiles of individual cristae (n=75, from 8 different acquisitions) were fitted using a Gaussian model. Full width half maximum (FWHM) was estimated on raw STED and 2D denoised + 2D deconvolved with $\sigma_{xy} = 1.0$ pixel (c). SPITFIR(e) improvement is shown separately for all the line profiles analyzed and the line profiles from bigger cristae (FWHM<sub>raw</sub> = 60 nm), indication of lower SNR. Statistics analysis is including in (d). It indicates the improvement of average resolution (unpaired two-sample t-test) and variance of cristae resolution (F-test).
from 8 different acquisitions) allowed us to estimate the resolution at FWHM of raw and restored images (Fig. 3(b, c)). The obtained resolution is statistically improved from direct STED imaging when SPITFIR(e) is applied (with a smaller variance), as shown in (Fig. 3(d)), and correspond to twice the reported results reported in transmitted electron microscopy³⁵. Finally, we compared the resolution achieved by SPITFIR(e) and RL as shown in Supplementary Fig. S2. As examined by intensity profiles, resolutions reached by both algorithms are very similar. However, deconvolution artifacts are strongly reduced using SPITFIR(e) (Supplementary Fig. S2(c)).

## SPITFIR(e) reliably restores weak fluorescent signals under extreme low-light conditions

Additional demonstrations of the ability of SPITFIR(e) to retrieve the high spatial-frequency information under extreme light conditions in 3D multifocus microscopy and spinning-disk confocal microscopy, are presented below. We focus on the performance limits of 3D deconvolution knowing that it is possible to improve the results by denoising the images beforehand as described above.

First, we tested SPITFIR(e) to visualize the 3D dynamics of mitochondria in living cells under extreme low-light conditions with 3D multifocus microscopy (MFM; see Supplementary Fig. S3) for different exposure values by reducing the excitation light dose (see Fig. 4(a)). In this experiment, each 3D stack was restored with a 3D Gaussian PSF model. The results on 2D plane (6th plane in Fig. 4(a)) are shown in Fig. 4 for three different illumination values using the same sample. We considered three different levels of sparsity: "highly", "moderately", and "weakly" sparse. In Fig. 4(b), we display the increasing of SNR values for 100 time points computed from two ROIs selected in the foreground (red rectangle) and background (green rectangle) in Fig. 4(a, top left image). The results computed for the three different levels of sparsity demonstrate the ability of SPITFIR(e) to reduce noise while enhancing small structures and removing background. Overall, the most satisfying results with SPITFIR(e) are obtained with "moderate" sparsity, (Fig. 4(a) third column, and Supplementary Fig. S3(c) versus (d)).

We further evaluated the flexibility of SPITFIR(e) on spinning-disk confocal microscopy images. Multi-pinholes spinning-disk microscopes are compatible with acquisition regimes allowing high speed live cell imaging. Nevertheless, notably because of photobleaching effects, time exposure must be reduced in order to capture very fast intracellular events, such as membrane trafficking. We acquired images with different excitation light doses in order to increase the frame rate and capture fast movements while reducing photobleaching and photo-toxicity effects. To assess the overall sensitivity of SPITFIR(e) cells expressing Calcium Dependent Mannose 6 phosphate Receptor-eGFP (CD-M6PR-eGFP) were imaged at variable excitation light levels, reducing light intensity by a factor of 2 and 4, as shown in Fig. 4(d). The three time-lapse temporal series display bleaching curves which correlates with the excitation light intensity levels and SNR values (calculated from Fig. 4(c), and following the equation displayed in figure legend). As expected, when we reduce by half or one fourth the laser power (raw data denoted 1/2 and 1/4 in Fig. 4(f)) and orange and green curves in Fig. 4(d, e), the bleaching effect is reduced and the SNR is lower, as compared to data with illumination at 1. We applied the proposed SPITFIR(e) restoration confocal pipeline on the three 3D+Time series. As shown in Fig. 4(e), SPITFIR(e) restore the SNR for low (1/2 orange curve and 1/4 green curve) illumination conditions, while removing the time dependent photobleaching effect for high (1, blue curve) or intermediate (1/2 orange curve) laser power illumination. Beyond these quantitative measurements, we also explore the SPITFIR(e) image quality rendering in these experiments. Figure 4(f), display 3D+Time images 1, 1/2, 1/4 for three time points of the corresponding temporal series illustrating broadly the SPITFIR(e) restoration pipeline performance. While image quality is reduced for the lowest illumination intensity (1/4), the restored images are good enough for spot detection. More importantly, SPITFIR(e) preserves SNR and thus restore image quality even at high frame rate and high laser excitation (compare 1 and Spitfir(e)[1] at t<sub>84.1s</sub> in Fig. 4(f)). As other experiments complementarily depicting different intracellular structural features, Supplementary Fig. S4 shows the maximum intensity projection along the Z-axis of a single stack depicting Rab5 proteins coupled to Green Fluorescence Protein (Rab5-eGFP), labelling early endosomes, for which the original images are noisy and contains diffraction limited spots against a dark background, while Supplementary Fig. S5 shows the maximum intensity projection along the Z-axis of a stack depicting the thin filamentous structures of polymerized actin (mCherry-LifeAct). The latter image is apparently not as sparse as the Rab5-image. These two examples are challenging because the PSF is known approximately and the image contents are very different. The 3D PSF is assumed to be approximated by a 3D Gaussian function. In each figure, we reported the results obtained with four SPITFIR(e) strategies (3D deconvolution, PbyP deconvolution, 3D denoising + 3D deconvolution, 3D denoising + PbyP deconvolution) and different levels of sparsity. Applying denoising beforehand enables here to better enhance weak fluorescent signals while removing blur and background. From Supplementary Figs. S4-S5, it visually best results on the Rab5-eGFP and LifeActin-mCherry images are obtained when the sparsity is "high" and "moderate", respectively.

We conclude that even at low photon counts and low SNR, SPITFIR(e) helps to preserve sample integrity, while allowing to extract exploitable data for quantitative analysis.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** SPITFIR(e) image reconstruction applied to wide-field multifocus and spinning-disk confocal microscopy with different amounts of sparsity, exposure and illumination values.</figcaption></figure>
**Figure 4** (previous page).

**Multifocus microscopy.** Data are temporal series (100 time points) of 3D stacks composed of nine planes each depicting mitochondria in U2OS cells transfected with TOM20 (translocase of outer mitochondrial membrane) fused to GFP (GFP-TOM20). Exposure time: 50 ms. The sample is imaged with three different doses of illumination light. MFM allows simultaneous acquisition of a 3D stack of nine images equally spaced with (dz = 330 nm) focal 2D images (pixel size = 120 nm, see also Supplementary Fig. S3). (a) The image grid display the deconvolution results on the 6th plane for successive increasing laser power values (EXP 1, 2 and 3). The 3D stacks have been deconvolved with 3D Gaussian PSF model ($\sigma_{xy} = 1.5$ pixels and $\sigma_z = 1.5$ pixels) and with three different levels of sparsity ("weak", "moderate", and "high"). (b) The plots are SNR values calculated for the 6th plane at each time-point as follows: SNR = Power(signal) / Power(background) = $|\Omega(x_S)|^{-1} \sum_{x \in \Omega(x_S)} I(x)^2 / |\Omega(x_B)|^{-1} \sum_{x \in \Omega(x_B)} I(x)^2$, where $\Omega(x_S)$ (with size $|\Omega(x_S)|$) and $\Omega(x_B)$ (with size $|\Omega(x_B)|$) are two ROIs centered at pixels $x_S$ (red ROI) and $x_B$ (green ROI), respectively. Scale bar: 10 $\mu$m.

**Spinning-disk microscopy.** Data are temporal series (30 time points) of 3D stacks composed 14 planes each depicting CD-M6PR-eGFP Hela cells. Exposure time: 20 ms, $\Delta t$ between two stacks: 2.9 s. The sample is imaged 3 times with different doses of illumination light. (c) The figure depicts a MIP (Maximum Intensity Projection) of the 3D stack at $t_{0s}$, acquired with a relatively high laser power dose (intensity $I$) (Scale bar = 10 $\mu$m). The plots represent SNR values measured over time on the single 6th plane (time in seconds) for three time series acquired with $I/4$, $I/2$, $I$ of laser illumination, successively before (d) and after (e) applying SPITFIR(e). The SNR metric is defined as above with two ROIs centered at pixels $x_S$ (red ROI in c) and $x_B$ (green ROI in c), respectively. The figure (f) shows a zoomed area (yellow in c) of the 6th plane at time $t_{0s}$, $t_{43.5s}$ and $t_{84.1s}$ extracted from the full movies acquired with high dose, half and one fourth dose of illumination light before (I, $I/2$ $I/4$) and after applying 4D denoising + 3D deconvolution SPITFIR(e) (Spitfir(e) [I], Spitfir(e) [I/2]) and Spitfir(e)[I/4]) to the three 4D volumes, with a Gibson-Lanni PSF model truncated to 3 planes around the center of the PSF. Intensities "min" and "max" are indicated as grey level bars, under each thumbnail images. Image contrast was adapted for visualization purposes. Scale bar: 4.5 $\mu$m.

**SPITFIR(e) quantitatively outperforms state-of-the-art algorithms**

In fluorescence microscopy, restoration algorithms were dedicated to remove Poisson-Gaussian noise originating from the low-photons regimes (Poisson noise) and dark current induced by the electronic imaging detectors (Gaussian noise). In Methods, we formally demonstrate that a quadratic fidelity term combined with any regularization term is optimal to handle Poisson-Gaussian noise. This new result allows us to apply SPITFIR(e) on images acquired with low light doses (Fig. 4) and to quantitatively compare state-of-the-art deconvolution algorithms mainly developed for Gaussian noise removal, and applied here to fluorescence microscopy images.

First, we fairly benchmarked the performance of a large collection of competitive deconvolution algorithms given in Table 2, including the popular deconvolution methods such as RL, iterative constrained Tikhonov-Miller (ICTM), Gold-Meinel (GM)$^{36}$. As the implementation in 3D of several deconvolution algorithms is not always possible, we conducted experiments on four 2D images shown in Table 1. The four ground-truth images were normalized in the range [0, 1] and further blurred by considering a Gaussian point spread function (PSF) with different standard deviation values $\sigma_{xy}$. A Gaussian noise with zero mean and standard deviation $\tau$ was also added to these images in order to generate observed noisy and blurry data. Note that, in the case of RL algorithm, which is originally designed to deal with Poisson noise, the degraded images are re-scaled to the original dynamics of the underlying reference images. The RL deconvolution results are then re-normalized for a fair comparison with those obtained by the other methods. Moreover, before applying GM, the noisy images are smoothed with a Gaussian filter as pre-processing step because this algorithm assumes that the noise is negligible. In these experiments, we consider three different PSF sizes which correspond to $\sigma_{xy} \in \{1, 1.25, 1.5\}$ and three distinct noise levels $\tau \in \{0.01, 0.02, 0.04\}$ (intensities are in the range [0, 1]). In Table 1, we reported the best possible PSNR values for each method (with optimal parameter adjustment) and the true PSF size $\sigma_{xy}$ and noise standard deviation $\tau$. In this experiment, RL and GM algorithm were often ranked at the end of the benchmark, while SPITFIR(e) usually provided the best PSNR values.

In terms of visual assessment, we display the noise-free versus artificially blurred and noisy images (Image #1) depicting small fluorescent vesicles acquired with classical laser scanning confocal microscopy in Fig. 5(a) and Fig. 5(b)), respectively. In this example, the original image was corrupted by a Gaussian PSF ($\sigma_{xy} = 1$ pixel) and an additive white Gaussian noise with standard deviation $\tau = 0.01$ (corresponding to a value of PSNR = 35.35 dB). In Fig. 5(c) and Fig. 5(d), we can notice noise amplification with RL and GM, respectively, resulting in the apparition of undesired high-intensity pixels (i.e., the “night sky” effect) and unrealistic reconstructed structures since the high-frequency components were not correctly restored. ICTM produced deconvolution results without noise amplification artifacts but tends to over-smooth image details due to the quadratic nature of Tikhonov-Miller penalty (see Fig. 5(e)), The GN regularizer, which also incorporates Tikhonov penalty, provides indeed a similar deconvolved image with blurred details (see Fig. 5(h)) while comparing to the ICTM solution. In contrast to these over-smoothed results, TV (see Fig. 5(f)), TV-$L_1$ (see Fig. 5(g)) and SV (see Fig. 5(i)) regularizers generate artificial sharper edges (stair-casing effects), respectively. HV generates visually more pleasant deconvolution results (see Fig. 5(m)). Regarding smooth-approximation based regularizations, LSHV (see Fig. 5(j)) provides results with restored details which are slightly sharper than those obtained with ICTM (see Fig. 5(e)), but noise is not sufficiently removed. In the comparison
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Deconvolution of a confocal microscopy image with different methods. (a) The original 512 × 256 pixels image (Image #1) (source: CNRS UMR144 Institut Curie) is a single 2D plane extracted from a 77 time frame series depicting Rab11a-mCherry protein in M10 cells, was previously acquired in line scanning confocal microscopy. Scale bar: 10 µm. (b) Artificially corrupted image by a Gaussian PSF $\sigma_{xy} = 1.0$ pixel and Gaussian white noise with standard deviation $\tau = 0.01$ (b). (c) Image deconvolved with the Richardson-Lucy (RL)³ algorithm. (d) Image deconvolved with the Gold-Meinel (GM)³⁹ algorithm. (e) Image deconvolved with the iterative constrained Tikhonov-Miller (ICTM)⁹ algorithm. (f) Image deconvolved with the Total Variation (TV)⁴⁰ regularizer. (g) Image deconvolved with the Total-Variation-L1 (TV-L1)²¹. (h) Image deconvolved with the GraphNet (GN)¹⁷ regularizer. (i) Image deconvolved with the Total Variation-Huber (TVH)⁴¹ regularizer. (j) Image deconvolved with the Log Sparse Hessian Variation (LSHV)¹⁶ regularizer. (k) Image deconvolved with the Good's roughness (GR)⁴² regularizer. (l) Image deconvolved with the Sparse Variation (SV)²⁰ regularizer. (m) Image deconvolved with the Hessian Variation (HV)¹⁴, ¹⁵ regularizer. (n) Image deconvolved with the Sparse Hessian Variation (SHV) (our SPITFIR(e) method) regularizer (black box). Zoom-in views are displayed in yellow windows for comparison in details.</figcaption></figure>

with TVH (Fig. 5(i)) and LSHV (Fig. 5(j)), the smooth version of GR (Fig. 5(k)) yields visual similar deconvolution results. However, GR suffer from severe “white pixel” artifacts in the background region. These "white pixels" were discarded beforehand to fairly compute the PSNR values. Finally, SPITFIR(e) (SHV) avoids over-smoothing and over-sharpening effect while preserving details and removing noise in the background (see Fig. 5(n)). These quantitative results were confirmed on 3D microscopy images by considering additional performance criteria.

We next conducted experiments on 3D confocal images and compared the results of SPITFIR(e) to those produced by CARE⁴³, a supervised deep-learning-based deconvolution method (Supplementary Fig. S6). As CARE is learned from a training set of dedicated biological structures, we evaluate the algorithms on the stack provided by the authors⁴³ and depicting the envelopes of nuclei stained with GFP-LAP2b (developing eye of zebrafish (Daniorerio) embryos. For the sake of visual assessment, Supplementary Fig. S6 only shows the restoration of the 5th plane of an anisotropic 3D stack composed 18 images⁴³. As expected, the reconstruction quality is visibly superior with CARE. Nevertheless SPITFIR(e) with 3D Gaussian PSF ($\sigma_{xy} = 2.0$ pixels, $\sigma_z = 0.5$ pixel) and Gibson-Lanni model⁴⁴ (truncated to 3 planes around the center of the PSF) provided satisfactory results. The main differences are visible on small filaments better enhanced with CARE. The results can be visually improved with SPITFIR(e) if the volume is denoised beforehand and deconvolved with the PbyP strategy (Supplementary Fig. S6(a; right)).
<figure><img src="image_9.png" /><figcaption>Image #1 (512 × 256)
(CM)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Image #2 (227 × 238)
(TIRFM)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Image #3 (256 × 283)
(TIRFM)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Image #4 (1054 × 1028)
(SIM)</figcaption></figure>

<table><caption>**Table 1.** PSNR scores of deconvolution methods applied to four 2D images (confocal (CM), TIRFM, SIM). The best scores (±0.05 dB) are in bold style (including when the error is lower than 0.05 dB). Images #1-#3 depict fluorescently tagged proteins (Rab11-mCherry, Langerin-YFP, TfR-pHluorin (Transferin-receptor)) corresponding to small bright spots over a dark background, observed in spinning-disk confocal microscopy (CM) (Image #1, courtesy of PICT-IBiSA imaging platform) and total internal reflection fluorescence microscopy (TIRFM) (Images #2 and #3) respectively, previously acquired³⁷,³⁸. Image #4 depicts microtubules observed in Structured Illumination Microscopy (SIM) with high resolution (up to 100 nm). This image (CIL 36147) is taken from the Cell Image Library (CIL) (cellimagelibrary.org). The images were deconvolved with the following methods: Richardson-Lucy (RL)³ algorithm, Gold-Meinel (GM)³⁹ algorithm, iterative constrained Tikhonov-Miller (ICTM)⁹ algorithm, Total Variation (TV)⁴⁰ regularizer, Total-Variation-L1 (TV-L1)²¹, GraphNet (GN)¹⁷ regularizer, Total Variation-Huber (TVH)⁴¹ regularizer, Log Sparse Hessian Variation (LSHV)¹⁶ regularizer, Good's roughness (GR)⁴² regularizer, Sparse Variation (SV)²⁰ regularizer, Hessian Variation (HV)¹⁴,¹⁵ regularizer, and Sparse Hessian Variation (SHV) (our SPITFIR(e) method) regularizer. The PSNR values correspond to optimal settings of parameters for all competing methods, including SPITFIR(e). We deconvolved the four images several times with different parameters to objectively supply the best PSNR value as possible, for all methods. The PSNR values are averaged values conventionally computed from 3 realizations (standard deviation less than 0.05 dB).</caption><thead><tr><th></th><th></th><th>RL</th><th>GM</th><th>ICTM</th><th>TV</th><th>TVL1</th><th>GN</th><th>TVH</th><th>LSHV</th><th>GR</th><th>SV</th><th>HV</th><th>SHV (SPITFIR(e))</th></tr></thead><tbody><tr><td rowspan="12">Image #1 (CM)</td><td rowspan="3">σ<sub>xy</sub> = 1.0<br>τ = 0.01</td><td>30.01</td><td>23.73</td><td>39.62</td><td>40.21</td><td>40.23</td><td>39.71</td><td>36.66</td><td>35.48</td><td>36.01</td><td><b>42.13</b></td><td>41.78</td><td><b>42.15</b></td></tr><tr><td>25.12</td><td>19.89</td><td>37.04</td><td>38.04</td><td>38.06</td><td>37.10</td><td>32.98</td><td>32.75</td><td>32.93</td><td>39.45</td><td>39.34</td><td><b>39.55</b></td></tr><tr><td>20.52</td><td>15.89</td><td>34.40</td><td>35.76</td><td>35.76</td><td>34.43</td><td>28.56</td><td>28.52</td><td>28.47</td><td><b>36.68</b></td><td>36.43</td><td>36.61</td></tr><tr><td rowspan="3">σ<sub>xy</sub> = 1.25<br>τ = 0.01</td><td>31.02</td><td>28.92</td><td>38.52</td><td>38.88</td><td>38.86</td><td>38.61</td><td>35.88</td><td>33.96</td><td>34.55</td><td><b>40.84</b></td><td>40.42</td><td><b>40.87</b></td></tr><tr><td>25.95</td><td>24.11</td><td>36.35</td><td>36.81</td><td>36.80</td><td>36.41</td><td>32.56</td><td>31.88</td><td>32.20</td><td><b>38.65</b></td><td>38.28</td><td>38.57</td></tr><tr><td>21.02</td><td>19.36</td><td>33.90</td><td>34.87</td><td>34.90</td><td>33.99</td><td>28.42</td><td>28.17</td><td>28.27</td><td><b>36.09</b></td><td>35.80</td><td>35.96</td></tr><tr><td rowspan="3">σ<sub>xy</sub> = 1.5<br>τ = 0.01</td><td>31.85</td><td>32.89</td><td>37.50</td><td>37.95</td><td>37.89</td><td>37.61</td><td>35.20</td><td>32.70</td><td>33.28</td><td><b>39.87</b></td><td>39.12</td><td>39.61</td></tr><tr><td>26.86</td><td>27.92</td><td>35.51</td><td>36.04</td><td>36.03</td><td>35.57</td><td>32.15</td><td>31.05</td><td>31.44</td><td><b>37.78</b></td><td>37.19</td><td>37.69</td></tr><tr><td>21.67</td><td>22.65</td><td>33.25</td><td>33.86</td><td>33.86</td><td>33.34</td><td>28.23</td><td>27.73</td><td>27.84</td><td><b>35.34</b></td><td>34.83</td><td>35.08</td></tr><tr><td rowspan="12">Image #2 (TIRFM)</td><td rowspan="3">σ<sub>xy</sub> = 1.0<br>τ = 0.01</td><td>27.61</td><td>10.66</td><td>31.08</td><td>30.92</td><td><b>31.29</b></td><td>31.00</td><td>30.73</td><td>29.79</td><td>30.73</td><td>30.92</td><td>31.01</td><td>31.01</td></tr><tr><td>22.81</td><td>29.64</td><td>29.79</td><td>29.77</td><td><b>30.01</b></td><td>29.69</td><td>28.72</td><td>28.56</td><td>28.65</td><td>29.83</td><td>29.89</td><td><b>30.01</b></td></tr><tr><td>17.93</td><td>27.50</td><td>28.27</td><td>28.40</td><td>28.43</td><td>28.21</td><td>25.75</td><td>25.75</td><td>25.76</td><td>28.60</td><td>28.63</td><td><b>28.78</b></td></tr><tr><td rowspan="3">σ<sub>xy</sub> = 1.25<br>τ = 0.01</td><td>28.26</td><td>17.07</td><td>29.96</td><td>29.78</td><td><b>30.13</b></td><td>30.00</td><td>29.73</td><td>28.42</td><td>29.75</td><td>30.01</td><td>29.87</td><td>30.01</td></tr><tr><td>23.94</td><td>15.37</td><td>28.99</td><td>28.94</td><td><b>29.23</b></td><td>28.95</td><td>27.94</td><td>27.53</td><td>27.92</td><td>29.17</td><td>29.02</td><td>29.16</td></tr><tr><td>18.98</td><td>12.13</td><td>27.64</td><td>27.68</td><td>27.77</td><td>27.58</td><td>25.22</td><td>25.21</td><td>25.22</td><td>27.97</td><td>27.85</td><td><b>28.09</b></td></tr><tr><td rowspan="3">σ<sub>xy</sub> = 1.5<br>τ = 0.01</td><td>28.46</td><td>23.29</td><td>28.97</td><td>28.80</td><td>28.91</td><td>27.38</td><td>28.93</td><td><b>29.09</b></td><td>29.07</td><td><b>29.09</b></td><td>28.90</td><td><b>29.08</b></td></tr><tr><td>24.76</td><td>20.61</td><td>28.17</td><td>28.10</td><td>28.38</td><td>28.19</td><td>27.35</td><td>26.72</td><td>27.29</td><td><b>28.44</b></td><td>28.04</td><td>28.29</td></tr><tr><td>19.94</td><td>16.54</td><td>26.98</td><td>26.97</td><td>27.10</td><td>26.91</td><td>24.81</td><td>24.74</td><td>24.78</td><td><b>27.34</b></td><td>27.08</td><td><b>27.35</b></td></tr><tr><td rowspan="12">Image #3 (TIRFM)</td><td rowspan="3">σ<sub>xy</sub> = 1.0<br>τ = 0.01</td><td>28.19</td><td>13.56</td><td>33.30</td><td>33.05</td><td><b>33.40</b></td><td>33.22</td><td>32.64</td><td>32.16</td><td>32.68</td><td>33.07</td><td>33.27</td><td>33.27</td></tr><tr><td>23.21</td><td>11.75</td><td><b>31.96</b></td><td>31.50</td><td>31.74</td><td>31.89</td><td>30.27</td><td>30.26</td><td>30.26</td><td>31.49</td><td>31.83</td><td>31.82</td></tr><tr><td>18.26</td><td>28.65</td><td>30.40</td><td>30.20</td><td>30.23</td><td>30.34</td><td>26.86</td><td>26.85</td><td>26.87</td><td>30.20</td><td><b>30.50</b></td><td><b>30.51</b></td></tr><tr><td rowspan="3">σ<sub>xy</sub> = 1.25<br>τ = 0.01</td><td>29.11</td><td>19.47</td><td>32.42</td><td>32.14</td><td><b>32.49</b></td><td>32.28</td><td>31.86</td><td>30.99</td><td>31.83</td><td>32.22</td><td>32.34</td><td>32.34</td></tr><tr><td>24.31</td><td>16.88</td><td><b>31.24</b></td><td>30.82</td><td>31.19</td><td>31.15</td><td>29.67</td><td>29.51</td><td>29.63</td><td>30.81</td><td>31.11</td><td>31.14</td></tr><tr><td>19.28</td><td>12.99</td><td><b>29.76</b></td><td>29.42</td><td>29.51</td><td>29.68</td><td>26.47</td><td>26.46</td><td>26.46</td><td>29.60</td><td>29.68</td><td><b>29.76</b></td></tr><tr><td rowspan="3">σ<sub>xy</sub> = 1.5<br>τ = 0.01</td><td>29.69</td><td>25.03</td><td><b>31.51</b></td><td>31.36</td><td>31.62</td><td>31.39</td><td>31.13</td><td>29.98</td><td>31.12</td><td>31.36</td><td><b>31.48</b></td><td><b>31.48</b></td></tr><tr><td>25.31</td><td>21.46</td><td><b>30.60</b></td><td>30.18</td><td>30.55</td><td>30.49</td><td>29.15</td><td>28.80</td><td>29.14</td><td>30.33</td><td>30.42</td><td>30.40</td></tr><tr><td>20.23</td><td>16.86</td><td><b>29.26</b></td><td>28.86</td><td>28.98</td><td>29.19</td><td>26.17</td><td>26.15</td><td>26.08</td><td><b>29.21</b></td><td>29.05</td><td><b>29.24</b></td></tr><tr><td rowspan="12">Image #4 (SIM)</td><td rowspan="3">σ<sub>xy</sub> = 1.0<br>τ = 0.01</td><td>31.24</td><td>21.49</td><td>36.57</td><td>35.87</td><td>35.93</td><td>37.06</td><td>36.19</td><td>33.74</td><td>36.14</td><td>37.89</td><td>37.62</td><td><b>38.21</b></td></tr><tr><td>25.89</td><td>19.17</td><td>33.95</td><td>33.76</td><td>33.77</td><td>34.42</td><td>32.60</td><td>31.89</td><td>32.60</td><td>35.19</td><td>34.99</td><td><b>35.54</b></td></tr><tr><td>20.80</td><td>15.66</td><td>31.33</td><td>31.41</td><td>31.44</td><td>31.79</td><td>28.41</td><td>28.31</td><td>27.66</td><td>32.48</td><td>32.31</td><td><b>32.70</b></td></tr><tr><td rowspan="3">σ<sub>xy</sub> = 1.25<br>τ = 0.01</td><td>32.13</td><td>28.63</td><td>34.79</td><td>34.55</td><td>34.56</td><td>35.35</td><td>34.58</td><td>31.46</td><td>34.53</td><td>36.08</td><td>35.66</td><td><b>36.29</b></td></tr><tr><td>26.82</td><td>24.75</td><td>32.55</td><td>32.36</td><td>32.40</td><td>33.09</td><td>31.67</td><td>30.33</td><td>31.64</td><td>33.79</td><td>33.49</td><td><b>34.07</b></td></tr><tr><td>21.45</td><td>19.94</td><td>30.29</td><td>30.29</td><td>30.35</td><td>30.67</td><td>28.03</td><td>27.64</td><td>27.35</td><td>31.41</td><td>31.13</td><td><b>31.54</b></td></tr><tr><td rowspan="3">σ<sub>xy</sub> = 1.5<br>τ = 0.01</td><td>32.34</td><td>33.04</td><td>33.14</td><td>33.36</td><td>33.40</td><td>33.70</td><td>33.08</td><td>29.82</td><td>33.09</td><td>34.47</td><td>34.00</td><td><b>34.66</b></td></tr><tr><td>27.55</td><td>28.78</td><td>31.27</td><td>31.19</td><td>31.25</td><td>31.80</td><td>30.67</td><td>29.03</td><td>30.66</td><td>32.48</td><td>32.08</td><td><b>32.64</b></td></tr><tr><td>22.15</td><td>23.57</td><td>29.36</td><td>29.37</td><td>29.46</td><td>29.88</td><td>27.60</td><td>26.94</td><td>27.04</td><td>30.44</td><td>30.08</td><td><b>30.54</b></td></tr></tbody></table>
a

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

b

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

**Figure 6.** Selection of regularization and sparsity parameters. (a) Display of square ROIs with variable sizes superimposed on single X-Y slide of 3D LLSM image depicting σ subunit-eGFP of the AP2 complex in Hela cells (left). Scale bar: 10 µm. Illustrations of energies $E_{\lambda'}(f, u)$ with respect to $\lambda' \in [0, 1[$ for different ROIs, and three levels ("weak", "moderate", "high") of sparsity. The minimax value $\lambda''$ is in the interval $[0.2, 0.4]$ according to the sparsity level in all ROIs, suggesting that the maximum is invariant to scale (right). (b) Display of the $3 \times 3$ matrix associated to the yellow ROI (see (a), left) whose the rows and columns represent smoothness and sparsity. The default value $\lambda' = 1/2$ corresponds to matrix cells marked in red. The minimax value $\lambda''$ corresponds matrix cells marked in orange. We also display the "symmetric" solution corresponding to $\lambda' = 1 - \lambda''$ (matrix cells marked in green).
Discussion

An important challenge in live fluorescence imaging is to limit photo-damages on the observed sample by reducing the dose of excitation light, while at the same time collecting enough photons to produce informative images for quantitative analysis. To overcome the difficulties of acquisition of 2D+Time and 3D+Time images at low excitation dose and boost the signal corrupted by blur and different sources of noise (shot noise, readout noise), we developed SPITFIR(e), a supermaneuverable and easy-to-use solution able to generate spatially consistent and high-quality fluorescence images. The approach is based on the fast minimization of a convex energy composed of a novel sparse-promoting regularization term and a Poisson-Gaussian-aware quadratic fidelity term. Unlike previous methods which were often tested on a limited number of microscopy modalities¹⁵,¹⁶,²³ and need a careful adjustment of regularization parameters to obtain high-quality restoration results, SPITFIR(e) is a nearly automatized method which adapts to a broad range of microscopy techniques (here up to 6), techniques, whatever the detector technology (CCD, sCMOS, PMT, Avalanche photodiodes), and image dimensionality. With this unique algorithm, the practitioner can switch from 4D image denoising to 3D deconvolution, or can chain the two processing tasks to obtain high-quality images. It is only needed to supply the measured PSF or the parameters of the 2D or 3D theoretical PSF (e.g., Gaussian) model, and to specify the level of desired sparsity ("weak", "moderate", "high"). A preview of restoration results obtained on a 2D ROI, as illustrated in Fig. 6, is displayed to help the practitioner selecting the amount of sparsity and smoothness. We demonstrated on experimental data that SPITFIR(e) extracted high frequency information and achieved significant improvement in the output quality even with a poor calibrated PSF. Nevertheless, SPITFIR(e) could miss some details in noisy images taken at low excitation light intensity or short exposure time. As the lower light intensities approach the limits of ability of SPITFIR(e) to recover information reliably, a successful strategy consists in chaining denoising (3D or 4D) and deconvolution with a 2D or 3D PSF model. This idea was already presented by the authors². In live experiments, fast acquisition with a high frame rate is generally required to capture fast events¹, inducing low SNR values. In presence of Poisson shot noise and readout noise, our results demonstrated that the proposed strategy is robust and able to boost signals reliably, and confirms that the quadratic fidelity term is suitable for Poisson-Gaussian noise as theoretically proved (see Methods).

We have showed that SPITFIR(e) substantially reduces noise and provides better high-quality images than other currently available methods over a wide variety of experimental conditions. In the particular case of LLSM, which is designed already for long lasting fast imaging of 3D large and close to isotropic volumes, SPITFIR(e) shows all its power, since it not only improves resolution and reduces noise, but still allows an increase in the imaging rates by permitting to reduce the illumination dose. SPITFIR(e) also achieved high resolution in sub-diffracted and superresolution microscopy under low-light conditions, without perturbing the signal amplitudes. For two-dimensional live-cell STED imaging, line intensity profiles along mitochondria cristae showed that SPITFIR(e) reached an improved lateral resolution of a factor of 1.3 (see Fig. 3(c, d)). Little is known about the actual dynamics of mitochondrial cristae in multiple physiological processes because of the visualization challenge of cristae structure in living cells. Hessian-SIM modality³³ or STED nanoscopy³⁴,⁴⁵ are the only approaches so far, able to achieve this goal. While only unique 2D STED images of live mitochondria are shown here, SPITFIR(e) applied to STED images might be of high interest in this context since it enables to restore these fine structures even when STED depletion is incomplete and scanning imaging at very low illumination, making live STED relatively long-term imaging accessible.

More generally, unlike smoothness regularizers, our sparse-promoting regularizer is able to reveal small details in the raw data, while preserving signal peaks and removing background, with no discernible artifact in the final restored images. We also demonstrated that on 2D/3D benchmarks that SPITFIR(e) outperforms existing deconvolution algorithms. Performance has been quantitatively studied in terms of PSNR values on artificially degraded 2D images and in terms of FRC on 3D experimental data using Fourier Image Resolution (FIRE³²) (Fig. 2).

Finally, we showed that SPITFIR(e) is able to produce comparable results to supervised machine-learning techniques⁴³. Unlike content-dependent methods based on deep-learning algorithms requiring many training images, SPITFIR(e) is a method which utilizes sparsity and smoothness as a priori knowledge and adapts to any image contents.

In conclusion, our results on experimental data endorse that SPITFIR(e) may be considered as a "swiss-knife", able to handle different samples imaged with various fluorescence microscopes, adapt to various sources of signal degradation, image contents, and variable signal-to-noise ratios. This algorithm does not require a well calibrated PSF and is easily controlled by choosing three levels of desired sparsity, that is "high", "moderate", and "weak". SPITFIR(e) resolves 3D structures and dynamics of biological components observed at low excitation light intensity and thus holds the promise to push the spatiotemporal resolution limit of sub-diffracted fluorescence microscopy technique.
# Methods

## Variational model for image restoration

An observed 2D/3D image $f: \Omega \subset \mathbb{R}^d \to \mathbb{R}$, $d = 2,3$ is a blurred and noisy version of the underlying true image $u: \Omega \to \mathbb{R}$ (i.e., $u \in \mathbb{R}^\Omega$). A general approach for restoring an image consists in finding the minimizer of an energy functional, i.e.,

$$ \hat{u} = \arg \min_{u \in \mathbb{R}^\Omega} F(f, u; h, \mathcal{T}) + \lambda R_\rho(u) \quad (1) $$

where $F(f, u; h, \mathcal{T})$ and $R_\rho(u)$ are the data fidelity term and the regularization term, respectively, and $\lambda > 0$ and $\rho \in [0, 1]$ are the parameters controlling the amount of smoothness and sparsity in the image $u$. The fidelity data term is the distance between the restored image $u$ and the observed image $f$ where $h$ is the PSF function and $\mathcal{T}$ is the degradation operator that encompasses Poisson and Gaussian noise. In what follows, we give the explicit expressions of the two terms.

**Data fidelity term.** A fidelity data term is generally derived from the general formation model, for instance dedicated to low photon counts and low-light regimes in fluorescence microscopy (e.g.,$^{15}$). We formally demonstrate below that a conventional quadratic fidelity term, which is optimal when the images are corrupted by additive white Gaussian noise, is also appropriate for mixed Poisson-Gaussian noise.

Let us consider the following observation model$^{46}$:

$$ f(\mathbf{x}) = v(\mathbf{x}) + \eta(v(\mathbf{x})) + \varepsilon(\mathbf{x}), \quad (2) $$

with $v(\mathbf{x}) = (h \ast u)(\mathbf{x})$, where $\mathbf{x} \in \Omega$ is the pixel position in the domain $\Omega$, $\eta$ and $\varepsilon$ represent the signal-dependent Poisson noise component and the zero-mean Gaussian noise component respectively, such that:

$$ \alpha(v(\mathbf{x}) + \eta(v(\mathbf{x}))) \sim \mathcal{P}(\alpha v(\mathbf{x})), \quad \varepsilon(\mathbf{x}) \sim \mathcal{N}(0, \tau^2), \quad (3) $$

$\alpha > 0$ is the quantization factor of the photodetector, and $\tau^2 > 0$ represents the Gaussian noise variance. According to the properties of Poisson distribution ($\mathbb{E}[\cdot]$ and $\mathrm{Var}[\cdot]$ denote the mathematical expectation and the variance respectively), we have

$$ \begin{cases} \mathbb{E}[\alpha(v(\mathbf{x}) + \eta(v(\mathbf{x})))]] = \alpha v(\mathbf{x}) = \mathrm{Var}[\alpha(v(\mathbf{x}) + \eta(v(\mathbf{x})))] = \alpha^2 \mathrm{Var}[\eta(v(\mathbf{x}))] \\ \mathbb{E}[\alpha(v(\mathbf{x}) + \eta(v(\mathbf{x})))]] = \alpha v(\mathbf{x}) + \alpha \mathbb{E}[\eta(v(\mathbf{x}))], \end{cases} \quad (4) $$

which yields to $\mathbb{E}[\eta(v(\mathbf{x}))] = 0$ and $\mathrm{Var}[\eta(v(\mathbf{x}))] = v(\mathbf{x})/\alpha$. Hence,

$$ \mathbb{E}[f(\mathbf{x})] = v(\mathbf{x}) + \underbrace{\mathbb{E}[\eta(v(\mathbf{x}))]}_{= 0} + \underbrace{\mathbb{E}[\varepsilon(\mathbf{x})]}_{= 0}, \quad (5) $$

and the overall variance of $f(\mathbf{x})$ is the sum of $\mathrm{Var}[\eta(v(\mathbf{x}))]$ and $\mathrm{Var}[\varepsilon(\mathbf{x})] = \tau^2$:

$$ \mathrm{Var}[f(\mathbf{x})] = \mathbb{E}[(f(\mathbf{x}) - \mathbb{E}[f(\mathbf{x})])^2] = \mathbb{E}[(f(\mathbf{x}) - v(\mathbf{x}))^2] = v(\mathbf{x})/\alpha + \tau^2. \quad (6) $$

We conclude that $v(\mathbf{x})$ is defined as the expected value of the noisy observations $f(\mathbf{x})$ and $\mathrm{Var}[f(\mathbf{x})] = \alpha^{-1}v(\mathbf{x}) + \tau^2$ is the noise variance$^{46-48}$.

Assume an image $v = h \ast u$ corrupted by a mixed Poisson-Gaussian and denote $f$ the noisy image. As explained above, the local noise variance can be represented at a given location $\mathbf{x} \in \Omega$ as:

$$ V_{\mathrm{noise}}(\mathbf{x}) \triangleq \mathbb{E}[(f(\mathbf{x}) - (h \ast u)(\mathbf{x}))^2] = \alpha^{-1}v(\mathbf{x}) + \tau^2. \quad (7) $$

Let us assume that the average intensity is preserved, that is

$$ \int_{\Omega} v(\mathbf{x}) \, d\mathbf{x} = \int_{\Omega} f(\mathbf{x}) \, d\mathbf{x}. \quad (8) $$

It follows that, by pre-multiplying the intensities on both sides of the previous equation by a factor $\alpha^{-1}$ and adding the constant value $\tau^2$, we have:

$$ \int_{\Omega} \alpha^{-1} v(\mathbf{x}) + \tau^2 \, d\mathbf{x} = \int_{\Omega} \alpha^{-1} f(\mathbf{x}) + \tau^2 \, d\mathbf{x}. \quad (9) $$

Meanwhile, by integrating (7) over the image domain $\Omega$, we get

$$ \begin{aligned} \int_{\Omega} \mathbb{E}[(f(\mathbf{x}) - (h \ast u)(\mathbf{x}))^2] d\mathbf{x} &= \int_{\Omega} \alpha^{-1} v(\mathbf{x}) + \tau^2 d\mathbf{x}, \\ \int_{\Omega} \int_{\Omega} (f(\mathbf{x}) - (h \ast u')(\mathbf{x}))^2 p(u'(\mathbf{x})) d\mathbf{x} du'(\mathbf{x}) &= \int_{\Omega} \alpha^{-1} v(\mathbf{x}) + \tau^2 d\mathbf{x}, \end{aligned} \quad (10) $$
<table><caption>**Table 2.** Properties of sparse regularizers.</caption><thead><tr><th>Name</th><th>Regularizer<br>φ(Du(x))</th><th>Differentiation</th><th>Convexity<br>order</th><th>Smoothness</th></tr></thead><tbody><tr><td>TV-L<sub>1,ρ</sub><sup>21</sup></td><td>ρ||∇u(x)||<sub>2</sub> + (1 − ρ)|u(x)| d x</td><td>1</td><td>Convex</td><td>Non smooth</td></tr><tr><td>GraphNet (GN<sub>ρ</sub>)<sup>17</sup></td><td>ρ||∇u(x)||<sub>2</sub><sup>2</sup> + (1 − ρ)|u(x)|</td><td>1</td><td>Convex</td><td>Non smooth</td></tr><tr><td>Log Sparse Hessian Variation (LSHV<sub>ε</sub>)<sup>16</sup></td><td>log(ε + u<sup>2</sup>(x) + ||H u(x)||<sub>F</sub><sup>2</sup>)</td><td>2</td><td>Non convex</td><td>Smooth</td></tr><tr><td>Sparse Variation (SV<sub>ρ</sub>)<sup>20</sup></td><td>√ρ<sup>2</sup>||∇u(x)||<sub>2</sub><sup>2</sup> + (1 − ρ)<sup>2</sup>u(x)<sup>2</sup></td><td>1</td><td>Convex</td><td>Non smooth</td></tr></tbody></table>

where $p(u'(x))$ represents the probability distribution of $u$. Assume the following prior $p(u'(x)) = \mathbf{1}[f(x) = ((h \ast u)(x))]$ where $\mathbf{1}[\cdot]$ is the indicator function. Hence, we have by using (9)

$$ 
\begin{aligned}
\int_{\Omega} (f(\mathbf{x}) - (h \ast u)(\mathbf{x}))^2 d\mathbf{x} &= \int_{\Omega} \alpha^{-1} v(\mathbf{x}) + \tau^2 d\mathbf{x}, \\
&= \int_{\Omega} \alpha^{-1} f(\mathbf{x}) + \tau^2 d\mathbf{x}, \\
\|f - (h \ast u)\|_{2}^{2} &= \alpha^{-1} \|f\|_{1} + |\Omega| \tau^2.
\end{aligned}
\quad (11) $$

Now, starting from the seminal paper<sup>40</sup>, the restored image is found by solving the following optimization problem (Poisson-Gaussian noise):

$$ \min_{u} R(u) \quad \text{subject to} \quad \|f - (h \ast u)\|_{2}^{2} = \alpha^{-1} \|f\|_{1} + |\Omega| \tau^{2}. \quad (12) $$

It turns out that (12) is a constrained formulation with an equality constraint which is not convex. The corresponding formulation with inequality constraint is

$$ \min_{u} R(u) \quad \text{subject to} \quad \|f - (h \ast u)\|_{2}^{2} \le \alpha^{-1} \|f\|_{1} + |\Omega| \tau^{2}. \quad (13) $$

It has been established that, under additional assumptions, that (12) and (13) are equivalent. A Lagrange formulation can be then derived since the right-hand side of the inequality does not depend on the unknown image $u$:

$$ \min_{u} R(u) + \lambda' \|f - (h \ast u)\|_{2}^{2}, \quad (14) $$

where the parameter $\lambda' > 0$ balances the two energy terms and is an unknown function of the Poisson-Gaussian noise variance. The Karush-Tucker conditions guarantees that (15) and (13) are equivalent for a particular choice of $\lambda'$. Finally, if we set $\lambda' = \lambda^{-1} \ne 0$, we can equivalently reformulate the minimization problem as:

$$ \min_{u} \|f - (h \ast u)\|_{2}^{2} + \lambda R(u). \quad (15) $$

In conclusion, a quadratic fidelity term $F(\cdot)$ is appropriate for Gaussian-noise removal provided it is combined with regularizer $R(\cdot)$.

**Regularization term.** Commonly-used regularizers in image restoration have the following form:

$$ R(u) = \int_{\Omega} \phi(Du(\mathbf{x})) d\mathbf{x}, \quad (16) $$

where $D$ is a linear operator (called “regularization operator”) used to control the spatial distribution of $u$ and $\phi(\cdot)$ is a positive potential function usually related to a norm distance. A typical example is the Tikhonov-Miller regularization<sup>10</sup> defined as the squared norm of the gradient of $u$: $\phi(Du(\mathbf{x})) = \|\nabla u(\mathbf{x})\|_{2}^{2}$.

We propose a convex sparsity-promoting regularizer based on second-order derivatives to avoid the emergence of undesirable stair-casing effects. The Sparse Hessian Variation (SHV) regularizer is defined as follows

$$ \mathrm{SHV}_{\rho}(u) = \int_{\Omega} \underbrace{\sqrt{\rho^2 \|\mathcal{H}u(\mathbf{x})\|_{F}^{2} + (1-\rho)^2 u(\mathbf{x})^2}}_{\|D_{2,\rho} u(\mathbf{x})\|_{2}} d\mathbf{x}, \quad (17) $$
where $\| \cdot \|_F$ denotes the matrix Frobenius norm and $\rho \in [0, 1]$ is the weighting parameter. Moreover, since the Frobenius norm of the Hessian matrix is equal to the Euclidean norm of its vectorized version, the operator $D_{2,\rho} \in \mathbb{R}^{10}$ is defined (for $d=3$) as:

$$D_{2,\rho} u(\mathbf{x}) \triangleq \left( (1-\rho)u(\mathbf{x}), \rho \frac{\partial^2 u(\mathbf{x})}{\partial x_x}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial y_y}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial z_z}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial x_x}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial y_y}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial z_z}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial x_x}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial y_y}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial z_z}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial z_y} \right)^T. \quad (18)$$

The idea behind this combination is to sparsify jointly the spatial distribution of image intensities and the second order derivatives of the image to encourage smooth variations between spatially-contiguous non-zero regions of the underlying image. The resulting image contain bright objects against a large smooth dark background, and no spurious edge. Unlike LSHV$^{16}$, SHV is convex and more flexible as it involves a parameter $\rho \in [0, 1]$ which is helpful to eliminate or preserve background. In SPITFIR(e), we consider only three values to control the amount of sparsity: $\rho = 0.1$ ("sparse"), $\rho = 0.6$ ("moderately" sparse), and $\rho = 0.9$ ("weakly" sparse). It is worth noting that, unlike conventional energy modelings in fluorescence microscopy, the background is not preliminarily computed with an ad-hoc method beforehand but is implicitly obtained when $u$ is estimated.

**Discrete formulation and optimization.** We present here the discretization setting required to derive the optimization algorithm. Let us consider a sampling grid of the following form in 3D ($d=3$) (expressions in 2D are similar).

$$\Lambda = \mathbb{Z}^3 \cap \Omega = \{1, 2, \dots, N_x\} \times \{1, 2, \dots, N_y\} \times \{1, 2, \dots, N_z\}.$$

The observed noisy and blurry image $f$ is represented by its digitized (discrete) version as follows: $f = \mathcal{T}(Hu)$, where $f, u \in \mathbb{R}_+^N$ with $N = N_x \times N_y \times N_z$, $H \in \mathbb{R}^{N \times N}$ is a matrix that models the point spread function of the microscope in the discrete setting, and $\mathcal{T}$ is degradation operator. Both $u$ and $f$ are assumed to be non negative. In the discrete setting, the blurring operator $H$ corresponds to a discrete convolution which can be efficiently computed by using fast Fourier transform (FFT)$^{49-52}$. To discretize $D_{2,\rho}$, we use standard finite differences for the Hessian operators with Neumann conditions on image boundaries. The 2D-3D deconvolution problem is then defined in the discrete setting as the minimizer of the following energy:

$$E(u) = \frac{1}{2} \|Hu - f\|_2^2 + \lambda \|D_{2,\rho}u\|_2 + \iota_{\mathcal{C}}(u), \quad (19)$$

where $\lambda > 0$ is the regularization parameter and $\iota_{\mathcal{C}}$ is the characteristic function of the convex subset $\mathcal{C}$ defined as: $\iota_{\mathcal{C}} = 0$ if $u \in \mathcal{C}$ and $+\infty$ otherwise, is very helpful to impose positivity constraint on the solution.

The optimization problem (19) is convex since the underlying energy functional is defined as the sum of convex terms, but it is non-smooth. To minimize the energy (19), several algorithms have been investigated, including FISTA$^{53}$, ADMM$^{24,25}$ and MM$^{54}$ algorithms. Instead of applying the aforementioned optimization methods, we investigated a first-order method to minimize the sum of convex functions, based on the proximal splitting approaches$^{28,29,55-60}$. The main idea consists in splitting the original problem into several simple sub-problems in the way that each single function of the sum can be processed separately (see details in Appendix). These operators are well-suited for large-scale problems arising in signal and image processing, because they only exploit first-order information of the function and thus enable fast and efficient computation. To solve the minimization problem (19), we adapted the full splitting approach described in$^{28,29}$. The key idea is to evaluate the gradient, proximity and linear operators individually in order to avoid implicit operations such as inner loops or inverse of linear operators.

**Selection of the regularization parameter.** The adjustment of the regularization parameter $\lambda$ in (19) is generally crucial in most image restoration algorithms, and may be time consuming. A non optimal choice of this parameter may over-smooth object borders, suppress structural details, generate artifacts or weakly reduce noise. The practical issue is then to automatically tune this parameter on a case-by-case basis to get the best performance as possible, from the input noisy image.

Instead of applying the usual principle of cross-validation, we investigated another strategy which is based on the comparison of several competing restored images (or regions of interest (ROIs)) over a range of parameters $\lambda$. All the intensity values are preliminary normalized with respect to the maximum intensity to get intensity values in the range $[0, 1]$. Henceforth, the energy (1) is rewritten as follows:

$$E_{\lambda'}(f, u) = \lambda' F(f, u; h, \mathcal{T}) + (1 - \lambda') R_{\rho}(u) \quad (20)$$

where $\lambda' \in [0, 1[$ is a positive constant such that $\lambda' = (1 + \lambda)^{-1}$. The two energy terms are expected be equally balanced if $\lambda' = 1/2$, which is considered as the default value in the algorithm. Nevertheless, we applied the minimax principle$^{61}$ to find $u^*$ that best minimizes, over all $\lambda'$ values, the maximum of $E_{\lambda'}(f, u)$:

$$(u^*, \lambda'^*) = \min_u \max_{\lambda'} E_{\lambda'}(f, u) = \max_{\lambda'} \min_u E_{\lambda'}(f, u) \triangleq \max_{\lambda'} E_{\lambda'}(f, u^*).$$
<table><caption>**Table 3.** Computing times of SPITFIR(e) (200 iterations for energy minimization).</caption><thead><tr><th>SPITFIR(e)<br>200 iterations</th><th>Image size</th><th>CPU<br>Intel Xeon<br>(2.8 GHz, 4 threads)</th><th>GPU<br>Nvidia Quadro M2000</th><th>Nvidia Tesla K80</th></tr></thead><tbody><tr><td>2D Denoising</td><td>512 × 512</td><td>349 ms</td><td>345 ms</td><td>278 ms</td></tr><tr><td>3D Denoising</td><td>512 × 512 × 22</td><td>8 s 74 ms</td><td>3 s 78 ms</td><td>1s 277ms</td></tr><tr><td>4D Denoising</td><td>512 × 512 × 14 × 32</td><td>8 min 23 s 789 ms</td><td>1 min 40 s 98 ms</td><td>45 s 744 ms</td></tr><tr><td>2D Deconvolution</td><td>512 × 512</td><td>767 ms</td><td>550 ms</td><td>428 ms</td></tr><tr><td>3D Deconvolution</td><td>512 × 512 × 22</td><td>23 s 574 ms</td><td>5 s 786 ms</td><td>2 s 713 ms</td></tr></tbody></table>

The optimization problem may be solved by examining the finite set of solutions $E_{\lambda'}(f, u^*)$ and selecting the $\lambda'$ valued associated to the highest energy $E_{\lambda'}(f, u^*)$. The good news is that Gennert and Yuille$^{61}$ demonstrated that there exists a unique value $\lambda'$ that maximizes $E_{\lambda'}(f, u^*)$ (convexity of $E_{\lambda'}(f, u^*)$). Hence, by initializing $\lambda' = 1/2$ and considering a finite set $\{\lambda'_m\}_{m=1}^M$ of increasing values such that $\lambda'_{[M/2]} = 1/2$, a fast search technique can be developed which consists in restoring the image with $\lambda' = \lambda'_{[M/2]+1}$ and $\lambda'_{[M/2]-1}$ and keeping the $\lambda'$ corresponding to the highest energy $E_{\lambda'}(f, u^*)$. The raw image is then restored with a highest (or lowest value) of $\lambda'$, and the SPITFIR(e) procedure is repeatedly applied on raw data (with $\lambda = (1 - \lambda')/\lambda'$) until the energy $E_{\lambda'}(f, u^*)$ decreases (or increases).

The optimal value solution $(u^*, \lambda'^*)$ is very helpful as it corresponds to the smoothest image that is acceptable (before oversmoothing) to the user. As illustrated in Fig. 6(a), we found that the optimal parameter $\lambda'^*$ varies little in ROIs $(\lambda'^* \in [0.2, 0.4]$ in Fig. 6(a)), including if the ROI is the entire image. This means that the ROI can be the entire image or a small patch containing meaningful fluorescent signals and manually selected by the practitioner. The aforementioned procedure is then run on a ROI (2D or 3D patch) by considering three levels of sparsity ("high", "moderate", "weak") given a finite set of values of $\lambda'$. This preview study is very low time consuming and computationally demanding. The results are displayed as a $3 \times 3$ matrix whose the rows and columns represent smoothness and sparsity (see Fig. 6(b)).

The practitioner then selects a given element in this $3 \times 3$ matrix to modify or to confirm the default parameter settings. In general, the 3x3 matrix summarizes the possible solutions, quantified on 3 levels for each parameter. This visualization is optional but is helpful to typically confirm the amount of sparsity and smoothness in the reconstructed images before processing large volumes. In a more automatic mode, the expert only specifies if the expected image is "highly", "moderately" or "weakly" sparse, without using the visualization module. The smoothing parameter is then automatically estimated following the minimax principle. To our knowledge, there is few practical rules to automatically set the two parameters to process 3D images and videos; an alternative approach to estimate algorithm parameters consists in computing FRC curves from a single image and subsampling operations, after each iteration of the algorithm (e.g, RL algorithm)$^{27}$.

## Implementation details and code availability

The input image is the raw (noisy and blurred) image. In the case of Lattice Light Sheet Microscopy, a necessary "deskew" operation is applied to raw data before to be supplied to SPITFIR(e).

Images have been processed on classical workstation for CPU calculation (CPU Intel Xeon 2.8GHz 4 threads). We use this setup since it is the most common setup in biology labs. For high speed processing we developed a GPU version of SPITFIR(e). We run the GPU on two types of GPU, a Nvidia Quadro M2000 that can be found in most workstations and a Nvidia Tesla K80 for computing grid. Typical computing times (CPU and GPU) for processing 2D-3D-4D images are reported in the Table 3. These computing times are small and take into account the full data processing steps (including data normalization and data copy to the GPU device). The code includes several CPU and GPU implementations and can be downloaded for free from the GitHub website (https://github.com/sylvainprigent/spitfire) along with accompanying documentation.

## Cell culture fluorescence labeling

The hTERT-immortalized RPE1 cells (Human Retinal Pigment Epithelial Cell) were purchased from ATCC (CRL-4000). RPE1 cells within 30 passage number were used and grown in DMEM-F12 medium without phenol red (Life Technologies) supplemented with 10% (vol/vol) fetal bovine serum (FBS) at 37°C in humidified atmosphere with 5% CO2. For LLSM (Lattice Light Sheet Microscopy) imaging, cells were seeded 5 to 6 hours before image acquisition at $150 \times 10^3$ per well of a 6 well plate, containing 3 to 4 coverslips with a 5mm diameter, treated as previously described in Chen *et al.*$^{30}$. Live cells were stained for microtubule and mitochondria using the Tubulin Tracker DR (Tub-Tracker Deep Red, ThermoFisher Scientific) probe (Ex. 652nm/Em. 669) at 10 ng/ml and 100ng/ml of PKMR (PK Mito Red), a cyclooctatetraene-conjugated-cyanine-3 dye (Ex. 549/Em. 569) kindly provided by Dr. Z. Chen from Peking University$^{31}$, respectively. Cell labeling was performed in two steps. Cells were first incubated for 15min with Tubulin Tracker DR and washed twice with LLSM medium (DMEM/F-12, without
phenol red, with 1% BSA, 1mM pyruvate, and 20mM HEPES). Coverslips were then transferred to the lattice light-sheet microscope (LLSM) sample holder and inserted into the imaging chamber containing 3.5 ml of the later medium and 0.1 mg/ml of PKMR.

Previously described Hela cells expressing Rab5-eGFP, CRISP genome edited SUM 159 breast carcinoma cell line, using eGFP-tagged σ unit of the AP2 complex (Adaptor Protein complex) were kindly provided by Dr. Tomas Kirchhausen (HMS, Boston, MS, USA) via Dr Ludger Johannes (Institut Curie, Paris, France), and CD-M6PR-EGFP expressing Hela cell line, by Dr Bernard Hoflack (TU-Dresden, Dresden, Germany). Human retinal pigment epithelial cells (RPE1) expressing mCherry-LifeAct were kindly provided by Dr Manuel Thery⁶². All these cells were cultured and prepared for imaging with the LLSM as described above.

For fixed samples, RPE1 cells were cultured and seeded as above, incubated for two additional hours with 1µl/ml of SiR-Actin (from a stock solution at 1mM; Ex.652/Em.674; Spirochorme AG, SWZ) and 1µl/ml of verapamil (from a stock solution at 1 mM) before washing in PBS and fixation with 4% paraformaldehyde freshly prepared in PBS from a 16% stock solution (Electron Microscopy Science, PA, USA). Cells were then permeabilized with 0.2% TX100 for 10 minutes, before incubation with anti α-tubulin DM1 IgG1-mAb for 45 minutes, diluted at 1/500 in PBS plus 0.1% FBS. After washing, goat anti-mouse IgGs secondary antibodies coupled to Alexa 488nm (diluted at 1/200 in PBS plus 0.1% FBS; Ex.496/Em.519; ThermoFisher Scientific, Bordeaux, France) were added for 30 minutes incubation. Coverslips were preserved at 4°C in PBS before transfer to the LLSM chamber containing 3.5 ml of the PBS at the room temperature.

For STED imaging RPE1 cells were plated on glass bottom Petri dishes (μ-dish 35mm-High, 1.5H, Ibidi, GmbH, Grafelfing, Germany). Microtubule labeling was performed as for LLSM imaging, while mitochondria were labeled by diluting 16000 × a 1mM stock solution of PKMITO Orange [gifted by Dr. Z. Chen from Peking University and Genvivo Biotech, Beijing, PR of China] in DMEM/F12 medium for a final 62.5 nM, for 15 minutes at 37°C before image acquisition.

For MFM imaging the U2OS cells, originating from human bone osteosarcoma with an epithelial morphology (obtained from ATCC (HTB-96)), were grown in DMEM culture medium (11880, by Thermo Fisher Scientific) without phenol red supplemented with 1% Glutamax, 1% Penicillin-Streptomycin and 10% fetal bovine serum (26140079, by Thermo Fisher Scientific). Cells are kept at 37°C in the presence of 5% CO2 in a humid environment. Cells were seeded on 25mm round coverglass for 12h in the same culture conditions. Cells were then transfected with TOM20-GFP DNA plasmid using X-tremeGENE HP DNA Transfection Reagent (Sigma). After 24 hours, cells are mounted and imaged. During the experiments cells were kept at 37°C with 5% CO2 in a humid environment using a Tokai Hit heating system (INUBG2E-PPZI).

Previously acquired datasets of M10 cells, stably expressing Langerin tagged with eYFP, or Rab11a tagged with mCherry³⁸ or transiently transfected with PhLuorin Transferrin-Receptor (TfR)³⁷ were also exploited in this study.

## Microscopy

**Lattice light-sheet microscopy (LLSM).** Lattice light-sheet microscopy was done on a commercialized version of a previously described setup³⁰ from 3i (Denver, USA). Cells were scanned incrementally through a 20 µm long light sheet in 600 nm steps using a fast piezoelectric flexure stage equivalent to ≈ 325 nm with respect to the detection objective and were imaged using a sCMOS camera (Orca-Flash 4.0; Hamamatsu, Bridgewater, NJ). Excitation was achieved with 488 nm (Sapphire Coherent), 560 nm (MPB Communications) or 642 nm (MPB Communications) diode lasers at 10-20% acousto-optic tunable filter transmittance with 50-200 mW (initial box power) through an excitation objective (Special Optics 28.6× 0.7 NA 3.74-mm water-dipping lens) and detected via a Nikon CFI Apo LWD 25× 1.1 NA water-dipping objective with a 2.5× tube lens with a final pixel size of 104 nm. Lattice light-sheet imaging was performed using an excitation pattern of outer NA equal to 0.55 and inner NA equal to 0.493. Composite volumetric datasets were obtained using ≈10 ms exposure/slice/channel at a time resolution of 2 seconds per total cell volume (about 60 planes). Fifty-six time points were acquired within 3 to 4 minutes (total raw data, only 1:45 min were analyzed for the movie) for live cell acquisition. Ten time-points were acquired within 15 to 20 seconds for fixed cell acquisition. Acquired data were deskewed, a necessary step to realigned image frames, using LLSpy, a python library to facilitate lattice light sheet data processing (copyright to T. Lambert, Harvard Medical School, Boston, USA; https://github.com/tlambert03/LLSpy). Deskewed images are then considered as Raw images. Napari, a multi-dimensional image viewer for Python (https://github.com/napari/napari), was used for 3D rendering. Maximum intensity projections were generated using ImageJ/Fiji 1.53c⁶³. Intensity profiles were plotted using Matlab2019b. Related supplementary Videos S1 and S2 were generated using napari and napari-animation plugin (https://github.com/napari/napari-animation). 3D Richardson-Lucy deconvolution was performed using deconvlucy from MATLAB imaging processing toolbox. Fourier Image Resolution (FIRE) was estimated using Fourier Ring Correlation³².

**Stimulated emission depletion microscopy (STED).** Image acquisition was performed with a STEDYCON module (Abbe-
rior Instruments, Göttingen, Germany) mounted at the camera port TCS SP8 STED microscope (Leica, Mannheim, Germany)
with a HC PL APO C2S ×100 oil objective (1.4 NA) used in 2D mode. Depletion was performed with a 775 nm pulsed laser
1-7 ns at 80% (413 mW at the objective lens position). Labeled mitochondria and microtubules were imaged with excitation at a wavelength of 594 nm and 640 nm, respectively. Nominal laser powers were adjusted at 7% for the 594 nm laser (3.15 µW at the objective lens position) and at 10% for the 640 nm (26 mW at the objective lens position). The time-gated fluorescence detection was done on a detector single photon counting avalanche photodiode between 578-627 nm and 650-700 nm with a pinhole settled at 1.1 AU. Imaging was executed with the 9 lines accumulation. Pixel dwell time was 10 µs and a pixel size of 25 nm × 25 nm. Images were generated using ImageJ/Fiji 1.53c3. Line profiles plots were measured with ImageJ. Then line profiles were fitted to a Gaussian model using Matlab2019b. Finally, Full width at half maximum (FWHM) is estimated from the fitting results. Unpaired two-sample *t*-test and variance *F*-test were performed using GraphPad Prism 9. 2D Richardson-Lucy deconvolution was performed using deconvlucy from MATLAB imaging processing toolbox.

**Multifocus microscopy (MFM).** Imaging was performed on a custom made MFM system with a Nikon 100x oil immersion (1.4 NA) detection objective. We recall briefly the main characteristics of the system⁶⁴. The emission path of a wide-field microscope (Nikon Ti) was modified by introducing several optical elements. At the output of the microscope a lens conjugated the back focal plane of the microscope to a multifocus grating (MFG). The motif of the phase grating is designed such that it splits with high efficiency the emission into 9 diffraction orders with equal intensities. The periodicity of the grating is distorted in order to add a different defocusing powers in the different diffraction orders. A chromatic correction module (CCM) based on the combination of a multifaceted prism and gratings corrects the chromatic aberration introduced by the MFG. A final lens is placed before the imaging camera (iXon Ultra897, Andor) in order to form the final images. Each diffraction order is focused on a different part of the camera, and yields an image of a different axial plane of the sample. As such, a single camera exposure results in a simultaneous acquisition of 9 focal planes. For live cell imaging, a 488nm laser was used to excite the sample. The fluorescence was filtered by a fluorescence filter (FF01-525/45-25, Semrock). The exposure time of the camera was set to 50ms and the laser power was adjusted to yield different signal to noise ratio.

**Other imaging modalities.** Fast confocal imaging in Fig. 4 and Supplementary Fig. S4-5 was performed with a spinning-disk confocal based on a Ti-Eclipse inverted microscope (Nikon, Tokyo, JP) equipped with X1 confocal head (Yokogawa, Tokyo, JP), Prime95B, an sCMOS camera (Teledyne Photometrics, Tucson, AZ). Thirty z-stack cell volumes were obtained, scanning incrementally in 300 nm z-step at 20 ms/frame (14 frames/volume), using a 100X, Plan Apo NA 1.4 oil objective (Nikon, Tokyo, JP) with a final pixel size of 110 nm. Excitation was achieved at 488-nm (100 mW) with a diode laser (GATACA systems, Massy, FR), at different illuminations, adjusting the acousto-optic tunable filter transmittance at 80% (I), 40%(I/2) and 20%(I/4). The acquisition is driven by the Metamorph 8.6 software (Molecular Devices, San Jose, CA). The original single 2D image in Fig. 5 was extracted from a previously acquired 3D+Time stack (512 × 256 × 1 plane × 74 frames) by laser line scanning microscopy (A1R, Nikon, Tokyo, JP) using a ×100, NA 1.35, APO lambda S silicone objective and excitation with a 561 nm laser (by courtesy of the PICT-Curie Imaging platform) and further denoised using nDSafir² [5 iterations and patch size of 11 × 11 pixels], in order to build a reference "ground truth". Spinning-disk confocal microscopy (CM) and Total Internal Reflection microscopy (TIRFM) of Table 1 and Fig. 5 were either used as described before in Gidon et al.³⁸, Pécot et al.⁶⁵ or previously acquired datasets³⁷,⁶⁶, reused in this study.

**Data availability**

Data generated during this study are available in figShare: https://doi.org/10.6084/m9.figshare.21378327.

**References**

1. Carlton, P. et al. Fast live simultaneous multi-wavelength 4-dimensional optical microscopy. Proc. Natl. Acad. Sci. USA **107**, 16016–16022 (2010).

2. Boulanger, J. et al. Patch-based nonlocal functional for denoising fluorescence microscopy image sequences. IEEE Trans. Medical Imaging **29**, 442–454 (2010).

3. Richardson, W. H. Bayesian-based iterative method of image restoration. J. Optical Society of America **62**, 55–59 (1972). URL http://www.osapublishing.org/abstract.cfm?URI=josa-62-1-55.

4. Lucy, L. B. An iterative technique for the rectification of observed distributions. The Astronomical J. **79**, 745–754 (1974).

5. Shepp, L. A. & Vardi, Y. Maximum likelihood reconstruction for emission tomography. IEEE Trans. Medical Imaging **1**, 113–122 (1982).

6. Sibarita, J.-B. Deconvolution Microscopy, 201–243 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2005). URL https://doi.org/10.1007/b102215.

7. Bruce, M. & Butte, M. Real-time gpu-based 3d deconvolution. Opt. Express **21**, 4766–4773 (2013).

8. Guo, M. & et al. Rapid image deconvolution and multiview fusion for optical microscopy. Nat. Biotechnol. **38**, 1337–1346 (2020).
9. van der Voort, H. T. M. & Strasters, K. C. Restoration of confocal images for quantitative image analysis. *J. Microscopy* **178**, 165-181 (1995). URL http://dx.doi.org/10.1111/j.1365-2818.1995.tb03593.x.

10. van Kempen, G. M. P., van der Voort, H. T. M., Bauman, J. G. J. & Strasters, K. C. Comparing maximum likelihood estimation and constrained Tikhonov-Miller restoration. *IEEE Engineering in Medicine and Biology Magazine* **15**, 76-83 (1996).

11. van Kempen, G. M. P., van Vliet, L. J., Verveer, P. J. & van der Voort, H. T. M. A quantitative comparison of image restoration methods for confocal microscopy. *J. Microscopy* **185**, 354-365 (1997). URL http://dx.doi.org/10.1046/j.1365-2818.1997.001-629.x.

12. Dey, N. et al. Richardson-Lucy algorithm with total variation regularization for 3d confocal microscope deconvolution. *Microsc. Res. Tech.* **69**, 260-266 (2006).

13. Hom, E. et al. Aida: An adaptive image deconvolution algorithm with application to multi-frame and three-dimensional data. *J. Opt. Soc. Am. A* **24**, 1580-1600 (2007).

14. Lefkimmiatis, S., Ward, J. P. & Unser, M. Hessian Schatten-norm regularization for linear inverse problems. *IEEE Trans. Image Processing* **22**, 1873-1888 (2013).

15. Ikoma, H., Broxton, M., Kudo, T. & Wetzstein, G. A convex 3d deconvolution algorithm for low photon count fluorescence imaging. *Scientific Reports* **8**, 111489 (2018).

16. Arigovindan, M. et al. High-resolution restoration of 3D structures from widefield images with extreme low signal-to-noise-ratio. Proc. Natl. Acad. Sci. USA **110**, 17344-17349 (2013). URL http://www.pnas.org/content/110/43/17344full.pdf.

17. Ng, B., Vahdat, A., Hamarneh, G. & Abugharbieh, R. Generalized sparse classifiers for decoding cognitive states in fMRI. In Wang, F., Yan, P., Suzuki, K. & Shen, D. (eds.) *Proceedings Int. MICCAI Workshop on Machine Learning in Medical Imaging (MLMI)*, 108-115 (Beijing, China, 2010).

18. Kandel, B. M., Wolk, D. A., Gee, J. C. & Avants, B. Predicting cognitive data from medical images using sparse linear regression. In Gee, J. C., Joshi, S., Pohl, K. M., Wells, W. M. & Zöllei, L. (eds.) *Int. Conf. Information Processing in Medical Imaging (IPMI 2013)*, 86-97 (Asilomar, CA, USA, 2013).

19. Grosenick, L., Klingenberg, B., Katovich, K., Knutson, B. & Taylor, J. E. Interpretable whole-brain prediction analysis with GraphNet. *NeuroImage* **72**, 304-321 (2013). URL http://www.sciencedirect.com/science/article/pii/S1053811912012487.

20. Eickenberg, M., Dohmatob, E., Thirion, B. & Varoquaux, G. Grouping total variation and sparsity: Statistical learning with segmenting penalties. In *Medical Image Computing and Computer-Assisted Intervention (MICCAI)*, 685-693 (Munich, Germany, 2015).

21. Chan, T. F. & Esedoglu, S. Aspects of total variation regularized L1 function approximation. *SIAM J. Applied Mathematics* **65**, 1817-1837 (2005). URL https://doi.org/10.1137/040604297.

22. Michel, V., Gramfort, A., Varoquaux, G., Eger, E. & Thirion, B. Total variation regularization for fMRI-based prediction of behavior. *IEEE Trans. Medical Imaging* **30**, 1328-1340 (2011).

23. Zhao, W. et al. Sparse deconvolution improves the resolution of live-cell super-resolution fluorescence microscopy. *Nat. Biotechnol.* (2021).

24. Afonso, M., Bioucas-Dias, J. & Figueiredo, M. An augmented la-grangian approach to the constrained optimization formulation of imaging inverse problems. *IEEE Trans. Image Processing* **20**, 681-695 (2011).

25. Boyd, S., Parikh, N., Chu, E., Peleato, B. & Eckstein, J. Distributed optimization and statistical learning via the alternating direction method of multipliers. *Found. Trends Mach. Learn.* **3**, 1-122 (2007).

26. Beck, A. & Teboulle, M. Fast iterative shrinkage-thresholding algorithm for linear inverse problems. *SIAM J. Imaging Sciences* **2**, 183-202 (2009).

27. Koho, S. et al. Fourier ring correlation simplifies image restoration in fluorescence microscopy. *Nat. Comm.* **10**, 3103 (2019).

28. Condat, L. A primal-dual splitting method for convex optimization involving lipschitzian, proximable and linear composite terms. *J. Optimization Theory and Applications* **158**, 460-479 (2013). URL http://dx.doi.org/10.1007/s10957-012-0245-9.

29. Condat, L. A generic proximal algorithm for convex optimization – application to total variation minimization. *IEEE Signal Processing Letters* **21**, 985-989 (2014).

30. Chen, B.-C. & et al. Lattice light-sheet microscopy: imaging molecules to embryos at high spatiotemporal resolution. *Science* **346**, 1257998 (2014). URL https://www.science.org/doi/10.1126/science.1257998.

31. Yang, Z. et al. Cyclooctatetraene-conjugated cyanine mitochondrial probes minimize phototoxicity in fluorescence and nanoscopic imaging. *Chem. Sci.* **11**, 8506-8516 (2020). URL http://dx.doi.org/10.1039/D0SC02837A.

32. Nieuwenhuizen, K., R. and Lidke et al. Measuring image resolution in optical nanoscopy. *Nat. Methods* **10**, 557-562 (2013).
33. Huang, X. et al. Fast, long-term, super-resolution imaging with hessian structured illumination microscopy. Nat. Biotechnol. 36, 451–459 (2018).

34. Till, S., Roesch, A., Riedel, D. & Jakobs, S. Live-cell sted nanoscopy of mitochondrial cristae. Sci. Rep. 9 (2019).

35. Lea, P. & Hollenberg, M. Mitochondrial structure revealed by high-resolution scanning electron microscopy. Am J Anat. 184 (1989).

36. Sage, D. et al. DeconvolutionLab2: an open-source software for deconvolution microscopy. Methods 115, 28–41 (2017).

37. Boulanger, J. et al. Fast high-resolution 3d total internal reflection fluorescence microscopy by incidence angle scanning and azimuthal averaging. Proc. Natl. Acad. Sci. USA 111, 17164–17169 (2014).

38. Gidon, A. et al. Rab11A/MyosinVb/Rab11-FIP2 complex frames two late recycling steps of langerin from erc to plasma membrane. Traffic 13 (2012).

39. Meinel, E. S. Origins of linear and nonlinear recursive restoration algorithms. J. Optical Society of America A 3, 787–799 (1986). URL http://josaa.osa.org/abstract.cfm?URI=josaa-3-6-787.

40. Rudin, L. I., Osher, S. & Fatemi, E. Nonlinear total variation based noise removal algorithms. Physica D: Nonlinear Phenomena 60, 259 – 268 (1992). URL http://www.sciencedirect.com/science/article/pii/016727899290242F.

41. Nikolova, M. Weakly constrained minimization: Application to the estimation of images and signals involving constant regions. J. Math. Imaging and Vision 21, 155–175 (2004). URL https://doi.org/10.1023/B:JMIV.0000035180.40477.bd.

42. Good, I. J. & Gaskins, R. A. Nonparametric roughness penalties for probability densities. Biometrika 58, 255–277 (1971).

43. Weigert, M., Schmidt, U. & et al., T. B. Content-aware image restoration: pushing the limits of fluorescence microscopy. Nat. Methods 15, 1090–1097 (2018).

44. Kirshner, H., Aguet, F., Sage, D. & Unser, M. 3-D PSF fitting for fluorescence microscopy: Implementation and localization application. Journal of Microscopy 249, 13–25 (2013).

45. Yang, X. et al. Mitochondrial dynamics quantitatively revealed by sted nanoscopy with an enhanced squaraine variant probe. Nat. Commun 11 (2020).

46. Foi, A., Trimeche, M., Katkovnik, V. & Egiazarian, K. Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data. IEEE Transactions on Image Processing 17, 1737–1754 (2008).

47. Mugnier, L., Fusco, T. & Conan, J.-M. Mistral: a myopic edge-preserving image restoration method, with application to astronomical adaptive-optics-corrected long-exposure images. J. Opt. Soc. Am. A 21, 1841–1854 (2004).

48. Soulez, F., Denis, L., Tourneur, Y. & Thiébaut, É. Blind deconvolution of 3D data in wide field fluorescence microscopy. In Int. Symp. Biomedical Imaging (ISBI) (Barcelone, Spain, 2012).

49. Heideman, M., Johnson, D. & Burrus, C. Gauss and the history of the fast Fourier transform. IEEE ASSP Magazine 1, 14–21 (1984).

50. Van Loan, C. Computational Frameworks for the Fast Fourier Transform (Society for Industrial and Applied Mathematics, 1992). URL http://epubs.siam.org/doi/abs/10.1137/1.9781611970999. http://epubs.siam.org/doi/pdf/ 10.1137/1.9781611970999.

51. Frigo, M. & Johnson, S. G. FFTW: an adaptive software architecture for the FFT. In IEEE Int. Conf. Acoustics, Speech and Signal Processing (ICASSP), vol. 3, 1381–1384 (1998).

52. Frigo, M. & Johnson, S. G. The design and implementation of FFTW3. Proceedings of the IEEE 93, 216–231 (2005).

53. Beck, A. & Teboulle, M. Fast gradient-based algorithm for constrained total variation image denoising and deblurring problems. IEEE Trans. Image Processing 18, 2419–2434 (2009).

54. Figueiredo, M., Bioucas-Dis, J. & Nowack, R. Majorization-minimization algorithms for wavelet-based image restoration. IEEE Trans. Image Processing 16, 2980–2991 (2007).

55. Mercier, B. Lectures on Topics in Finite Element Solution of Elliptic Problems (Springer-Verlag Berlin Heidelberg, 1979).

56. Eckstein, J. & Bertsekas, D. P. On the Douglas–Rachford splitting method and the proximal point algorithm for maximal monotone operators. Mathematical Programming 55, 293–318 (1992). URL http://dx.doi.org/10.1007/BF01581204.

57. Combettes, P. L. Solving monotone inclusions via compositions of nonexpansive averaged operators. Optimization 53, 475–504 (2004).

58. Combettes, P. L. & Wajs, V. R. Signal recovery by proximal forward-backward splitting. SIAM J. Multiscale Modeling and Simulation 4, 1168–1200 (2005). URL https://doi.org/10.1137/050626090.

59. Combettes, P. L. & Pesquet, J.-C. Proximal Splitting Methods in Signal Processing, 185–212 (Springer New York, New York, NY, 2011). URL http://dx.doi.org/10.1007/978-1-4419-9569-8_10.

60. Chambolle, A. & Pock, T. A first-order primal-dual algorithm for convex problems with applications to imaging. J. Math. Imaging and Vision 40, 120–145 (2011). URL http://dx.doi.org/10.1007/s10851-010-0251-1.

61. Gennert, M. & Yuille, A. Determining the optimal weights in multiple objective function optimization. In Proc. IEEE Int. Conf. Comp. Vision (ICCV), 87–89 (1988).
62. Vignaud, T. et al. Reprogramming cell shape with laser nano-patterning. *Cell Sci* **125**, 2134–2140 (2012).

63. Schindelin, J. & et al. Fiji: an open-source platform for biological-image analysis. *Nat. Methods* **9**, 676–682 (2012).

64. Hajj, B. et al. Whole-cell, multicolor superresolution imaging using volumetric multifocus microscopy. *Proc. Natl. Acad. Sci. USA* **111**, 17480–17485 (2014).

65. Pécot, T. et al. A quantitative approach for analyzing the spatio-temporal distribution of 3d intracellular events in fluorescence microscopy. *eLife* **7** (2018).

66. Basset, A., Boulanger, J., Salamero, J., Bouthemy, P. & Kervrann, C. Adaptive spot detection with optimal scale selection in fluorescence microscopy images. *IEEE Transactions on Image Processing* **24**, 4512–4527 (2015).

67. Combettes, P., Dung, D. & Vu, B. Proximity for sums of composite functions. *J. Math. Analysis and Applications* **380**, 680 – 688 (2011). URL http://www.sciencedirect.com/science/article/pii/S0022247X11002137.

68. Combettes, P., Condat, L., Pesquet, J. & Vū, B. A forward-backward view of some primal-dual optimization methods in image recovery. In *IEEE Int. Conf. Image Processing (ICIP)*, 4141–4145 (2014).

69. Moreau, J. J. Proximité et dualité dans un espace hilbertien. *Bulletin de la Société Mathématique de France* **93**, 273–299 (1965). URL http://www.numdam.org/item?id=BSMF_1965__93__273_0.

## Acknowledgements

This work was jointly supported by the French National Research Agency (France-BioImaging ANR-10-INBS-04-07 and ANR-10-INBS-04-01; DALLISH-ANR-16-CE23-0005; LabEx Cell(n)Scale (ANR-11-LABX-0038) as part of the Idex PSL ANR-10-IDEX-0001-02) and Innopsys company. This work was also supported by ITMO Cancer (18CQ091).

We would like to thank M. Maurin from Inserm U932 for his help with STED image acquisition, and the PICT imaging platform of Institut Curie, member of the France-BioImaging infrastructure (ANR-10-INBS-04-01) for maintaining spinning-disk and line-scanning confocal microscopes used in this study. We also thank Dr. Z. Chen from Peking University for providing us with mitochondrial fluorescent probes.

Finally, we warmly thank Dr. J. Boulanger from MRC Laboratory of Molecular Biology Cambridge, Dr. Z. Gan, Dr. P. Roudot, and Prof. Dr. G. Danuser from UT Southwestern Medical Center Dallas, Prof. Dr. J. Sedat from University of California San Francisco, Dr. M. Arigovindan from Indian Institute of Science, Prof. Dr. R. Heintzman from Friedrich Schiller University Jena, C. Gueudry and G. Petrazzuoli from Gataca Systems company, and V. Paveau and C. Cauchois from Innopsys company for fruitful discussions in early stages of this project.

## Author contributions statement

C.K. devised the project and the main conceptual ideas, supervised the project and was in charge of overall direction and planning. C.K., H.-N.N. and S.P. designed and implemented the presented SPITFIR(e) method, in discussions with J.S., C.A. V.-C. and L.L. J.S., C.A. V.-C. and L.L. provided the real datasets (LLSM, STED, CM), conceived experiments, and implemented the analysis workflows applied to real images. B.H. provided the real MFM datasets. H.-N.N. performed experiments on artificial datasets. S.P., C.A. V.-C., and L.L. performed experiments on real images (MFM, LLSM, STED, CM). C.K., J.S., C.A. V.-C., H.-N.N., and S.P. co-wrote the manuscript. All authors provided critical feedback and helped shape the research, analysis and manuscript.

## Additional information

### Competing interests

The authors declare no competing interests.

### Supplementary Information

APPENDIX (Technical details of the discrete formulation and optimization) and SUPPLEMENTARY FIGURES.

**Supplementary Video S1.** Composite movie corresponding to the time series of insets in Fig. 2 a and b. Upper part from left to right, microtubules, mitochondria and superimposed MIP data from LLSM after deskew treatment. Lower part, similar disposition after SPITFIR(e) 4D denoising and 3D deconvolution. Only 1:45 minutes of a 4 minutes series is shown. Scale bar: 5 µm.

**Supplementary Video S2.** Multi-Angle view of the 3D of superimposed image stacks from the same time series as in Supplementary Video 1. Upper part data from LLSM after deskew treatment, lower part after SPITFIR(e) 4D denoising and 3D deconvolution. Scale Bar = 5 µm.
SUPPLEMENTARY INFORMATION

SPITFIR(e): A supermaneuverable algorithm for fast
denoising and deconvolution of 3D fluorescence
microscopy images and videos

Sylvain Prigent<sup>1,2,+</sup>, Hoai-Nam Nguyen<sup>1,2,+</sup>, Ludovic Leconte<sup>1,2</sup>, Cesar Augusto

Valades-Cruz<sup>1,2</sup>, Bassam Hajj<sup>3</sup>, Jean Salamero<sup>1,2</sup>, Charles Kervrann<sup>1,2,*</sup>

<footnote n="1">¹ SERPICO Project-Team, Inria Centre Rennes-Bretagne Atlantique, 35042, Rennes Cedex, France</footnote>

<footnote n="2">² SERPICO/STED Team, UMR144 CNRS Institut Curie, PSL Research University, Sorbonne Universités, 75005, Paris, France</footnote>

<footnote n="3">³ Laboratoire Physico-Chimie, Institut Curie, PSL Research University, Sorbonne Universités, CNRS UMR168, 75005, Paris, France</footnote>

<footnote>* charles.kervrann@inria.fr</footnote>

<footnote>+ these authors contributed equally to this work.</footnote>
# APPENDIX

## Discrete formulation

The observed noisy and blurry image $f$ is represented by its digitized (discrete) version and not by its continuously defined counterpart. The continuous model is not appropriate for discrete images even though the estimation of the continuous image $u$ from discrete samples of $f$ is in principle possible. Consider a discrete formulation by assuming that the images $u$ and $f$ are non negative and sampled according to the sampling grid

$$ \Lambda = \mathbb{Z}^2 \cap \Omega = \{1, 2, \dots, N_x\} \times \{1, 2, \dots, N_y\} \times \{1, 2, \dots, N_z\}. \qquad (21) $$

The observed noisy and blurry image $f$ is represented by its digitized (discrete) version as follows:

$$ f = \mathcal{T}(h * u) $$

where $f, u \in \mathbb{R}_+^N$ with $N = N_x \times N_y \times N_z$, $H \in \mathbb{R}^{N \times N}$ is a matrix that models the point spread function of the microscope in the discrete setting, and $\mathcal{T}$ is the degradation operator. For a coordinate $(i, j, k) \in \Lambda$, we denote by $u_{i,j,k}$ (resp. $f_{i,j,k}$) the value of $u$ (resp. $f$) at position $(i, j, k) \in \Lambda$. A discrete version of these images are therefore given by $\{u_{i,j,k}\}_{1 \le i \le N_x, 1 \le j \le N_y, 1 \le k \le N_z}$ and $\{f_{i,j,k}\}_{1 \le i \le N_x, 1 \le j \le N_y, 1 \le k \le N_z}$. Denote $\mathcal{X} = \mathbb{R}^N$ with $N = N_x \times N_y \times N_z$, a finite dimensional vector space equipped with a standard inner (scalar) product

$$ \langle w, w' \rangle_{\mathcal{X}} = \sum_{i=1}^{N_x} \sum_{j=1}^{N_y} \sum_{k=1}^{N_z} w_{i,j,k} w'_{i,j,k}. \qquad (22) $$

The induced norm by the defined inner product is given by

$$ \|w\|_{\mathcal{X}_2} = \sqrt{\langle w, w \rangle_{\mathcal{X}}} = \left( \sum_{i=1}^{N_x} \sum_{j=1}^{N_y} \sum_{k=1}^{N_z} w_{i,j,k}^2 \right)^{\frac{1}{2}}. \qquad (23) $$

In the discrete setting, the blurring operator $H$ corresponds to a discrete convolution which can be efficiently computed by using fast Fourier transform (FFT)49-52. To discretize $D_{2,p}$, we use standard finite differences to approximate the second derivatives along the three dimensions of the volume at voxel $(i, j, k)$, with Neumann conditions on image boundaries:

$$ (D_{2,p} u)_{i,j,k} = \begin{pmatrix} (1-\rho)u_{i,j,k} \\ \rho(\Delta_{1,1}^2 u)_{i,j,k} \\ \rho(\Delta_{2,2}^2 u)_{i,j,k} \\ \rho(\Delta_{3,3}^2 u)_{i,j,k} \\ \rho(\Delta_{1,2}^2 u)_{i,j,k} \\ \rho(\Delta_{1,3}^2 u)_{i,j,k} \\ \rho(\Delta_{2,3}^2 u)_{i,j,k} \\ \rho(\Delta_{1,1}^2 u)_{i,j,k} \\ \rho(\Delta_{3,2}^2 u)_{i,j,k} \end{pmatrix} \in \mathbb{R}^{10}, \qquad (24) $$

where

$$ (\Delta_{1,1}^2 u)_{i,j,k} = \begin{cases} u_{i+1,j,k} - 2u_{i,j,k} + u_{i-1,j,k} & \text{if } 1 < i < N_x, \\ 0 & \text{otherwise,} \end{cases} $$

$$ (\Delta_{2,2}^2 u)_{i,j,k} = \begin{cases} u_{i,j+1,k} - 2u_{i,j,k} + u_{i,j-1,k} & \text{if } 1 < j < N_y, \\ 0 & \text{otherwise,} \end{cases} $$

$$ (\Delta_{3,3}^2 u)_{i,j,k} = \begin{cases} \delta^2(u_{i,j,k+1} - 2u_{i,j,k} + u_{i,j,k-1}) & \text{if } 1 < k < N_z, \\ 0 & \text{otherwise,} \end{cases} $$

$$ (\Delta_{1,2}^2 u)_{i,j,k} = (\Delta_{2,1}^2 u)_{i,j,k} = \begin{cases} (u_{i+1,j+1,k} - u_{i+1,j,k} - u_{i,j+1,k} + u_{i,j,k}) & \text{if } i < N_x, j < N_y \text{ and } k < N_z, \\ 0 & \text{otherwise.} \end{cases} $$

$$ (\Delta_{1,3}^2 u)_{i,j,k} = (\Delta_{3,1}^2 u)_{i,j,k} = \begin{cases} \delta(u_{i+1,j,k+1} - u_{i+1,j,k} - u_{i,j,k+1} + u_{i,j,k}) & \text{if } i < N_x, j < N_y \text{ and } k < N_z, \\ 0 & \text{otherwise.} \end{cases} $$

$$ (\Delta_{2,3}^2 u)_{i,j,k} = (\Delta_{3,2}^2 u)_{i,j,k} = \begin{cases} \delta(u_{i,j+1,k+1} - u_{i,j+1,k} - u_{i,j,k+1} + u_{i,j,k}) & \text{if } i < N_x, j < N_y \text{ and } k < N_z, \\ 0 & \text{otherwise.} \end{cases} $$
and $\delta$ is the ratio of the lateral-to-axial step sizes.

The discrete operators can be used to define the discrete SHV regularizer as:

$$ \|D_{2,\rho}u\|_2 = \sum_{i=1}^{N_x} \sum_{j=1}^{N_y} \sum_{k=1}^{N_z} \|(D_{2,\rho}u)_{i,j}\|_2, \qquad (26) $$

where the $L_1$-norm acts now on the discrete domain $\Lambda$. The 3D deconvolution problem is defined in the discrete setting as the minimizer of the following energy:

$$ E(u) = \frac{1}{2} \|Hu - f\|_2^2 + \lambda \|D_{2,\rho}u\|_2 + \iota_{\mathcal{C}}, \qquad (27) $$

where $\lambda > 0$ is the regularization parameter and $\iota_{\mathcal{C}}$ is the indicator of a convex set $\mathcal{C}$ such as:

$$ \mathcal{C} = \{u: \underline{u}_{i,j,k} \le u_{i,j,k} \le u_{\max}, 1 \le i \le N_x, 1 \le j \le N_y, 1 \le k \le N_z\}, \qquad (28) $$

where the upper bound $u_{\max} > 0$ is the maximal intensity value allowed and $\underline{u}_{i,j,k} \ge 0$ is an estimated lower bound of the pixel intensity which is spatially varying and then adapted to each pixel location. The spatially varying constraint on the lower bound of pixel intensity not only guarantees positivity of the solution but also helps to avoid over-sparsifying effect. In our experiment, for the sake of simplicity, we use the following lower bound $\underline{u}_{i,j,k} = \max(0, \tilde{f}_{i,j,k} - c\tau)$ where $c > 0$ and $\tilde{f}$ is a smoothed version of the observed noisy image $f$ by a low-pass (Gaussian) filter.

## Energy minimization and splitting algorithms

We notice that the objective function

$$ E(u) = \frac{1}{2} \|Hu - f\|_2^2 + \lambda \|D_{2,\rho}u\|_2 + \iota_{\mathcal{C}}, \qquad (29) $$

is a sum of linear composite functions as $u \mapsto \sum_{n=1}^3 \mathcal{F}_m(L_m u)$, where each $\mathcal{F}_m$ is a convex function and each $L_m$ is a linear bounded operator. Formally, we can write $\mathcal{F}_1 = \iota_{\mathcal{C}}$, $L_1 = \text{Id}$, $\mathcal{F}_2 = \lambda \|\cdot\|_2$, $L_2 = D_{2,\rho}$ and $\mathcal{F}_3 = \frac{1}{2} \|\cdot\|^2$, $L_3 = H(\cdot) - f$. Generic primal-dual proximal approaches can be used to minimize this linear combination of convex functions as proposed in$^{67,68}$, but it is not optimal since the smoothness of the quadratic terms $\frac{1}{2}\|Hu - f\|_2$ is not exploited. In order to solve the problem (29), the design of an appropriate algorithm requires therefore to take into account the specific form of the corresponding energy, i.e., the sum of a simple convex function $\mathcal{F} = \iota_{\mathcal{C}}$, a more sophisticated composite function $\mathcal{G} \circ L = \lambda \|D_{2,\rho}(\cdot)\|_2$ (here, $\mathcal{G} = \lambda \|\cdot\|_2$ and $L = D_{2,\rho}$) and a smooth function $\mathcal{H} = \frac{1}{2} \|H(\cdot) - f\|_2$.

In what follows, we present a first-order method to minimize the sum of convex functions based on the proximal splitting approaches$^{28,29,55-60}$. It consists in decomposing (splitting) the original problem into several simple sub-problems in the way that each single function of the sum can be processed separately. Indeed, smooth function involves its gradient operator, while non-smooth function implies its Moreau proximity operator$^{69}$. These operators are well-suited for large-scale problems arising in signal and image processing, because they only exploit first-order information of the function and thus enable fast and efficient computation.

Let us recall first that the proximity operator of a convex function $\mathcal{J}: \mathbb{R}^N \to \mathbb{R}$ is defined as

$$ \operatorname{prox}_{\tau \mathcal{J}}(u) = \arg \min_{w \in \mathbb{R}^N} \mathcal{J}(w) + \frac{1}{2\tau} \|u - w\|_2^2, \qquad (30) $$

where $\tau > 0$ is a control parameter. From this definition, it easy to verify that the proximity operator of the function $\mathcal{F}(u) = \iota_{\mathcal{C}}(u)$ is nothing else than the projection onto the convex subset $\mathcal{C}$ as the following

$$ \begin{aligned} \operatorname{prox}_{\tau \mathcal{F}}(u) &= \arg \min_{w \in \mathbb{R}^N} \iota_{\mathcal{C}}(w) + \frac{1}{2\tau} \|u - w\|_2^2 \\ &= \arg \min_{w \in \mathcal{C}} \|u - w\|_2^2. \end{aligned} \qquad (31) $$

If we denote $\operatorname{proj}_{\mathcal{C}}$ the projection operator on $\mathcal{C}$, its closed-form expression is given by

$$ (\operatorname{proj}_{\mathcal{C}}(u))_{i,j,k} = \max(\underline{u}_{i,j,k}, \min(u_{\max}, u_{i,j,k})). \qquad (32) $$

Moreover, the quadratic function $\mathcal{H}(u) = \frac{1}{2} \|Hu - f\|_2^2$ possesses an analytic form for its associated proximity operator

$$ \begin{aligned} \operatorname{prox}_{\tau \mathcal{H}}(u) &= \arg \min_{w \in \mathbb{R}^N} \frac{1}{2} \|Hw - f\|_2^2 + \frac{1}{2\tau} \|u - w\|_2^2 \\ &= (\tau H^\top H + \text{Id})^{-1} (\tau H^\top H f + u), \end{aligned} \qquad (33) $$

where the symbol $\top$ denotes the adjoint of a linear operator and $H^\top: \mathbb{R}^N \to \mathbb{R}^N$ satisfies $\langle Hw, w' \rangle_{\mathbb{R}^N} = \langle w, H^\top w' \rangle_{\mathbb{R}^N}$. The evaluation of $\operatorname{prox}_{\tau \mathcal{H}}(u)$ corresponds to the inverse of a linear system that is not always possible in practice due to the high dimensionality of the problem.
For this reason, the optimization methods which involve the gradient of $\mathcal{H}$ are more appropriate since they do not require any inverse operator. In the comparison with $\mathcal{F}$ and $\mathcal{H}$, the calculation of the proximity operator in the case of the composite function $\mathcal{G} \circ L(u) = \lambda \|D_{2,p} u\|_2$ theoretically possible but is challenging because of the presence of $D_{2,p}$ which is not diagonal.

To solve the minimization problem (29), we adopt the full splitting approach described in$^{28,29}$. The key idea of this approach is to evaluate the gradient, proximity and linear operators individually in order to avoid implicit operations such as inner loops or inverse of linear operators. Accordingly, only “simple” computations are considered such as the gradient $\nabla \mathcal{H}$, the proximity operator of $\mathcal{F}$ and $\mathcal{G}$, the linear mapping $L$ and its adjoint operators $L^\top$. The corresponding proximal algorithm for the problem is written under the following general form at iteration $\ell$:

$$u^{(\ell+1)} = \operatorname{prox}_{\gamma \mathcal{F}} \{u^{(\ell)} - \gamma (L^{\top} z^{(\ell)} + \nabla \mathcal{H}(u^{(\ell)}))\}, \quad (34a)$$

$$z^{(\ell+1)} = \operatorname{prox}_{\zeta \mathcal{G}^*} \{z^{(\ell)} + \zeta L (2u^{(\ell+1)} - u^{(\ell)})\}, \quad (34b)$$

where $\gamma, \zeta > 0$ are proximal parameters, $\mathcal{G}^*$ denotes the Legendre-Fenchel conjugate of $\mathcal{G}$ and its proximity operator $\operatorname{prox}_{\zeta \mathcal{G}^*}$ can be directly computed from $\operatorname{prox}_{\mathcal{G}^*} \zeta \mathcal{G}^*$ by using the Moreau's identity $v = \zeta \operatorname{prox}_{\mathcal{G}^*} (\frac{v}{\zeta}) + \operatorname{prox}_{\zeta \mathcal{G}^*} (v)$. Following$^{28,29}$, to guarantee the convergence of the proposed algorithm, the parameters $\gamma$ and $\zeta$ must fulfill the condition

$$\gamma \left( \frac{1}{2} + \zeta \|L^\top L\| \right) < 1, \quad (35)$$

where $\|\cdot\|$ denotes the operator norm. The proofs of the convergence can be found in$^{28}$. We also note that the proposed algorithm belongs to the class of primal-dual algorithms which provide not only the solution of the primal problem (a.k.a. the original minimization problem) but also the solution of its dual problem.

Since the closed-form of $\operatorname{prox}_{\gamma \mathcal{F}}$ is already given, it remains to define the analytic expression of other terms in (34a) and (34b). We start with the gradient of the quadratic function $\mathcal{H}$ which is straightforwardly obtained by

$$\nabla \mathcal{H}(u) = H^\top (Hu - f). \quad (36)$$

Next, we notice that the regularization operator $L = D_{2,p}$ is a linear mapping, then its adjoint operator $L^\top = D_{2,p}^\top$ is defined using the equation $\langle u, D_{2,p}^\top v \rangle = \langle D_{2,p} u, v \rangle$, which implies ($d=3$)

$$(D_{2,p}^\top v)_{i,j,k} = (1-\rho)v_{i,j,k}^1 + \rho[(\Delta_{1,1}^{(2)})^\top v^2]_{i,j,k} + \rho[(\Delta_{2,2}^{(2)})^\top v^3]_{i,j,k} + \rho[(\Delta_{3,3}^{(2)})^\top v^4]_{i,j,k} \\ + \rho[(\Delta_{1,2}^{(2)})^\top v^5]_{i,j,k} + \rho[(\Delta_{1,3}^{(2)})^\top v^6]_{i,j,k} + \rho[(\Delta_{2,3}^{(2)})^\top v^7]_{i,j,k}, \quad (37)$$

where the involving adjoint operators are given below:

$$[(\Delta_{1,1}^{(2)})^\top w]_{i,j,k} = \begin{cases} w_{i+1,j,k} - 2w_{i,j,k} + w_{i-1,j,k} & \text{if } 1 < i < N_x \\ 0 & \text{otherwise} \end{cases}$$

$$[(\Delta_{2,2}^{(2)})^\top w]_{i,j,k} = \begin{cases} w_{i,j+1,k} - 2w_{i,j,k} + w_{i,j-1,k} & \text{if } 1 < j < N_y \\ 0 & \text{otherwise} \end{cases}$$

$$[(\Delta_{3,3}^{(2)})^\top w]_{i,j,k} = \begin{cases} w_{i,j,k+1} - 2w_{i,j,k} + w_{i,j,k-1} & \text{if } 1 < k < N_z \\ 0 & \text{otherwise} \end{cases}$$

$$[(\Delta_{1,2}^{(2)})^\top w]_{i,j,k} = [(\Delta_{2,1}^{(2)})^\top w]_{i,j,k} = \begin{cases} w_{i,j,k} & \text{if } i=1, j=1, 1 \le k \le N_z \\ w_{i,j,k} - w_{i,j-1,k} & \text{if } i=1, 1 < j < N_y, 1 \le k \le N_z; \\ -w_{i,j-1,k} & \text{if } i=1, j=N_y, 1 \le k \le N_z; \end{cases} \\ w_{i,j,k} - w_{i-1,j,k} - w_{i,j-1,k} + w_{i-1,j-1,k} & \text{if } 1 < i < N_x, j=1, 1 \le k \le N_z; \\ -w_{i,j-1,k} + w_{i-1,j-1,k} & \text{if } 1 < i < N_x, j=N_y, 1 \le k \le N_z; \\ -w_{i-1,j,k} & \text{if } i=N_x, j=1, 1 \le k \le N_z; \\ -w_{i-1,j,k} + w_{i-1,j-1,k} & \text{if } i=N_x, 1 < j < N_y, 1 \le k \le N_z; \\ w_{i-1,j-1,k} & \text{if } i=N_x, j=N_y, 1 \le k \le N_z. \end{cases}$$
$$
[\Delta_{1,3}^{(2)}^\top w]_{i,j,k} = [\Delta_{3,1}^{(2)}^\top w]_{i,j,k} =
\begin{cases}
w_{i,j,k} & \text{if } i=1, k=1, 1 \le j \le N_y \\
w_{i,j,k} - w_{i,j,k-1} & \text{if } i=1, 1 < k < N_z, 1 \le j \le N_y; \\
-w_{i,j,k-1} & \text{if } i=1, k=N_z, 1 \le j \le N_y;
\end{cases}
\\
\begin{cases}
w_{i,j,k} - w_{l-1,j,k} & \text{if } 1 < i < N_x, k=1, 1 \le j \le N_y; \\
w_{i,j,k} - w_{i-1,j,k} - w_{i,j,k-1} + w_{i-1,j,k-1} & \text{if } 1 < i < N_x, 1 < k < N_z, 1 \le j \le N_y; \\
-w_{i,j,k-1} + w_{i-1,j,k-1} & \text{if } 1 < i < N_x, k=N_z, 1 \le j \le N_y;
\end{cases}
\\
\begin{cases}
-w_{l-1,j,k} & \text{if } i=N_x, k=1, 1 \le j \le N_y; \\
-w_{l-1,j,k} + w_{l-1,j,k-1} & \text{if } i=N_x, 1 < k < N_z, 1 \le j \le N_y; \\
w_{l-1,j,k-1} & \text{if } i=N_x, k=N_z, 1 \le j \le N_y.
\end{cases}
$$

$$
[\Delta_{2,3}^{(2)}^\top w]_{i,j,k} = [\Delta_{3,2}^{(2)}^\top w]_{i,j,k} =
\begin{cases}
w_{i,j,k} & \text{if } j=1, k=1, 1 \le i \le N_x \\
w_{i,j,k} - w_{i,j,k-1} & \text{if } j=1, 1 < k < N_z, 1 \le i \le N_x; \\
-w_{i,j,k-1} & \text{if } j=1, k=N_z, 1 \le i \le N_x;
\end{cases}
\\
\begin{cases}
w_{i,j,k} - w_{i,j-1,k} & \text{if } 1 < j < N_y, k=1, 1 \le i \le N_x; \\
w_{i,j,k} - w_{i,j-1,k} - w_{i,j,k-1} + w_{i,j-1,k-1} & \text{if } 1 < j < N_y, 1 < k < N_z, 1 \le i \le N_x; \\
-w_{i,j,k-1} + w_{i,j-1,k-1} & \text{if } 1 < j < N_y, k=N_z, 1 \le i \le N_x;
\end{cases}
\\
\begin{cases}
-w_{i,j-1,k} & \text{if } j=N_y, k=1, 1 \le i \le N_x; \\
-w_{i,j-1,k} + w_{i,j-1,k-1} & \text{if } j=N_y, 1 < k < N_z, 1 \le i \le N_x; \\
w_{i,j-1,k-1} & \text{if } j=N_y, k=N_z, 1 \le i \le N_x.
\end{cases}
$$

From equations (37), one can deduce the following upper bound:

$$
\|D_{2,\rho}^{\top} D_{2,\rho}\| \le (1-\rho)^2 + 64\rho^2, \quad (38a)
$$

that are used for choosing the proximal parameters $\gamma$ and $\zeta$ according to (35).

The last term we want to deal with is the proximity operator prox<sub>ζ,q</sub>*. We also note that the proposed primal-dual algorithm does not necessitate evaluating the proximity operator of the composite function *G* ◦ *L* as in the case of generic proximal algorithms, but only prox<sub>G</sub> is required. Since *G* is related to the mixed norm whose the proximity operator is defined as:

$$
(\operatorname{prox}_{\tau\|\cdot\|_2}(v))_{i,j,k}^{\kappa} = \max(0, \|L_{i,j,k}v\|_2 - \tau) \frac{v_{i,j,k}^{\kappa}}{\|L_{i,j,k}v\|_2}, \quad (39)
$$

where $L_{i,j,k} : v \in \mathcal{X}^{10} \mapsto (v_{i,j,k}^{\kappa})_{1 \le k \le 10} \in \mathbb{R}^{10}$ is a linear operator ($d=3$). By using the Moreau's identity, we obtain the closed-form expression of $\text{prox}_{\zeta,q^*}$ as the following:

$$
(\operatorname{prox}_{\zeta q^*}(v))_{i,j,k}^{\kappa} = \frac{v_{i,j,k}^{\kappa}}{\max\left(1, \frac{\|L_{i,j,k}v\|_2}{\lambda}\right)}, \quad (40)
$$

which shows that prox<sub>ζq*</sub> is independent from ζ and moreover it is an pointwise operator. These properties allow therefore fast and efficient computation by exploiting the intrinsic parallelism of multicore processors.

Energy model for 4D denoising

A temporal series of 3D noisy images $f: \Omega \subset \mathbb{R}^4 \rightarrow \mathbb{R}$ is a noisy version of the underlying true 3D image sequence $u: \Omega \rightarrow \mathbb{R}$ modeled as follows: $f = \mathcal{T}(u)$. It follows that the denoising problem is formulated as the minimization of an energy functional defined as

$$
\int_{\Omega} (u(\mathbf{x}) - f(\mathbf{x}))^2 + \lambda \|D_{2,\rho} u(\mathbf{x})\|_2 d\mathbf{x} + \iota_{\mathcal{C}}(u) \quad (41)
$$
where

$$ D_{2,p}u(\mathbf{x}) := \left( (1-\rho)u(\mathbf{x}), \rho \frac{\partial^2 u(\mathbf{x})}{\partial x_x}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial y_y}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial z_z}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial t_t}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial x_x}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial y_y}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial z_z}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial t_t}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial x_x}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial y_y}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial z_z}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial t_t}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial x_x}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial y_y}, \rho \frac{\partial^2 u(\mathbf{x})}{\partial z_z} \right) \in \mathbb{R}^{14}. \quad (42) $$

To solve the optimization problem, we split the original problem (41) into several simple sub-problems as explained above. The resulting individual functions involved in the sum are then minimized separately. The algorithm is expected to be faster since FFT is not required in the implementation for image denoising. Nevertheless, the amount of data is much more larger since the whole 4D image sequence is denoised at once. In practice, a very long image sequence can be segmented into sub-sequences with a small overlap, and is processed independently and in parallel. The overlapping images can be averaged to reduce possible artifacts at the end.
SUPPLEMENTARY FIGURES

<figure><img src="image_15.png" /><figcaption>Supplementary Figure S1. Workflow for image restoration extracted from the global flow chart in Fig. 1 applied to 4D-denoising, background subtraction, and 3D deconvolution of mitochondria and microtubules images shown in Fig. 2.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Supplementary Figure S2.** Comparison of 2D image improvement in STED imaging. Images of Mitochondria (PKMO) in Live RPE1 cells (used in Fig. 3) before (a), after being processed by Richardson-Lucy 2D deconvolution (b) and SPITFIR(e) 2D denoising + 2D deconvolution (c). (b) Richardson-Lucy(RL) 3D deconvolution using deconvlucy from MATLAB imaging processing toolbox (2D Gaussian PSF and $\sigma_{xy} = 1$ pixel; 10 iterations). (c) SPITFIR(e) 2D denoising + 2D deconvolution (2D Gaussian PSF and $\sigma_{xy} = 1$ pixel). Insets are zoomed area illustrating SPITFIR(e) improvement (b; center) in signal to noise ratio (SNR) and image quality as compared to RL treatment (c; center). Intensity line profiles were measured as indicated in the insets and plotted for 2D STED (a; right), RL (b; right) and SPITFIRE(e) (c; left). Scale bars = 1 µm.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Supplementary Figure S3.** 3D multifocus microscopy reconstruction using SPITFIR(e). Data are temporal series (100 time points) of 3D stacks composed of nine planes each depicting mitochondria in U2OS cells transfected with TOM20 (translocase of outer mitochondrial membrane) fused to GFP (GFP-TOM20). Exposure time: 50 ms. 3D MFM stack (a) is acquired at the same time point. (b) Z-planes are equally spaced with dz = 330 nm and pixel size = 120 nm. 3D angular views are shown before (c) and after (d) SPITFIR(e) using "moderate" sparsity level. 3D average rendering was performed in napari. Data corresponds to EXP 2 on Fig. 4 (a). Scale bar: 5 µm.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Supplementary Figure S4.** Deconvolution of 3D spinning-disk confocal microscopy volume (14 planes) depicting live Hela cells expressing Rab5-eGFP proteins (exposure time: 100 ms). (a) Display of the maximum intensity projection along the Z-axis of raw data (scale bar: 8 µm). (b) Results with SPITFIR(e) with different amounts of sparsity ("weak" ($\rho$ = 0.9), "moderate" ($\rho$ = 0.6), "high" ($\rho$ = 0.1), automatic selection of the regularization parameter, and different strategies with SPITFIR(e): 3D Deconvolution, PbyP Deconvolution, 3D Denoising + 3D Deconvolution, 3D Denoising + PbyP deconvolution ($\sigma_{xy}$ = 1.5 pixels and $\sigma_{z}$ = 0.5 pixel).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Supplementary Figure S5.** Deconvolution of 3D spinning-disk confocal microscopy volume (24 planes) depicting live RPE1 cells expressing mCherry-LifeAct (exposure time: 20 ms). (a) Display of the maximum intensity projection along the Z-axis of raw data (scale bar: 10 µm). (b) Results with SPITFIR(e) with different levels of sparsity ("weak" ($\rho$ = 0.9), "moderate" ($\rho$ = 0.6), "high" ($\rho$ = 0.1), automatic selection of the regularization parameter and different strategies with SPITFIR(e): 3D Deconvolution, PbyP Deconvolution, 3D Denoising + 3D Deconvolution, 3D Denoising + PbyP deconvolution ($\sigma_{xy}$ = 1.5 pixels and $\sigma_{z}$ = 0.5 pixel).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Supplementary Figure S6.** Deconvolution of a 3D spinning-disk confocal microscopy image depicting the envelopes of nuclei stained with GFP-LAP2b (developing eye of zebrafish (Daniorerio) embryos, 60X/1.3-NA objective, 488 nm, *z* step = 2 µm, **source**:⁴³). (a) Comparison of CARE and SPITFIR(e) (3D deconvolution with a Gibson-Lanni PSF model (generated from⁴³), 3D deconvolution with a 3D Gaussian PSF model (σ<sub>xy</sub> = 2.0 pixels and σ<sub>z</sub> = 0.5 pixel), 3D denoising + PbyP deconvolution (σ<sub>xy</sub> = 2.0 pixels)). We display the results obtained on the 5th plane (scale bar: 25 µm). (b) Zoom-in views (5th plane) of deconvolution results obtained with CARE, RL, and SPITFIR(e) 3D (Gibson-Lanni PSF model).</figcaption></figure>