Impact Statement

This research paper describes a super-resolution method improving the spatial resolution of images acquired by common fluorescence microscopes and conventional blinking/fluctuating fluorophores. The problem is formulated in terms of a sparse and convex/nonconvex optimization problem in the covariance domain for which a well-detailed algorithmic and numerical description are provided. It is addressed to an audience working at the interface between applied mathematics and biological image analysis. The proposed approach is validated on several synthetic datasets and shows promising results also when applied to real data, thus paving the way for new future research directions.

2. Mathematical Modeling

For real scalars $ T,M>1 $ and $ t\hskip0.2em \in \left\{1,2,\dots, T\right\} $, let $ {\mathbf{Y}}_t\hskip0.2em \in {\unicode{x211D}}^{M\times M} $ be the blurred, noisy, and down-sampled image frame acquired at time $ t $. We look for a high-resolution image $ \mathbf{X}\hskip0.2em \in {\unicode{x211D}}^{L\times L} $ being defined as $ \mathbf{X}=\frac{1}{T}{\sum}_{t=1}^T{\mathbf{X}}_t $ with $ L= qM $ and defined on a $ q $-times finer grid, with $ q\hskip0.2em \in \unicode{x2115} $. Note that in the following applications, we typically set $ q=4 $. The image formation model describing the acquisition process at each $ t $ can be written as:

(1)$$ {\mathbf{Y}}_t={\mathcal{M}}_q\left(\mathcal{H}\left({\mathbf{X}}_t\right)\right)+\mathbf{B}+{\mathbf{N}}_{\mathbf{t}}, $$

where $ {\mathcal{M}}_q:{\unicode{x211D}}^{L\times L}\to {\unicode{x211D}}^{M\times M} $ is a down-sampling operator summing every $ q $ consecutive pixels in both dimensions, $ \mathcal{H}:{\unicode{x211D}}^{L\times L}\to {\unicode{x211D}}^{L\times L} $ is a convolution operator defined by the PSF of the optical imaging system and $ \mathbf{B}\hskip0.2em \in {\unicode{x211D}}^{M\times M} $ models the background, which collects the contributions of the out-of-focus (and the ambient) fluorescent molecules. Motivated by experimental observations showing that the blinking/fluctuating behavior of the out-of-focus molecules is not visible after convolution with wide de-focused PSFs, we assume that the background is temporally constant ($ \mathbf{B} $ does not depend on $ t $), while we allow it to smoothly vary in space. Finally, $ {\mathbf{N}}_t\hskip0.2em \in {\unicode{x211D}}^{M\times M} $ describes the presence of noise modeled here as a matrix of independent and identically distributed (i.i.d.) Gaussian random variables with zero mean and variance $ s\hskip0.2em \in {\unicode{x211D}}^{+} $ taking into account both the underlying electronic noise and the noise bias induced by $ \mathbf{B} $ (see Remark 1 for more details on the approximation considered). We assume that the molecules are located at the center of each pixel and that there is no displacement of the specimen during the imaging period, which is a reasonable assumption whenever short time acquisitions are considered.

Remark 1. A more appropriate model taking also into account the presence of signal-dependent Poisson noise in the data would be the following:

(2)$$ {\mathbf{Y}}_t=P\left({\mathcal{M}}_q\left(\mathcal{H}\left({\mathbf{X}}_t\right)\right)+\mathbf{B}\right)+{\mathbf{E}}_{\mathbf{t}}=P\left({\mathcal{M}}_q\left(\mathcal{H}\left({\mathbf{X}}_t\right)\right)\right)+P\left(\mathbf{B}\right)+{\mathbf{E}}_t,\hskip2em \forall t=1,2,\dots, T, $$

where, for $ \mathbf{W}\hskip0.2em \in {\unicode{x211D}}^{M\times M} $, $ P\left(\mathbf{W}\right) $ represents the realization of a multivariate Poisson variable of parameter $ \mathbf{W} $ and $ {\mathbf{E}}_{\mathbf{t}}\hskip0.2em \in {\unicode{x211D}}^{M\times M} $ models electronic noise with a matrix of i.i.d. Gaussian entries of zero mean and constant variance $ {\sigma}^2\hskip0.2em \in {\unicode{x211D}}_{+} $. Note that the second equality in (Reference Sage, Pham and Babcock2) holds due to the independence between $ {\mathcal{M}}_q\left(\mathcal{H}\left({\mathbf{X}}_t\right)\right) $ and $ \mathbf{B} $. Model (Reference Sage, Pham and Babcock2) is indeed the one we used for the generation of the simulated data, see Section 6.1. However, to simplify the reconstruction process, we simplified (Reference Sage, Pham and Babcock2) by assuming that $ \mathbf{B} $ has sufficiently large entries, so that $ P\left(\mathbf{B}\right) $ can be approximated as $ P\left(\mathbf{B}\right)\approx \hat{\mathbf{B}} $ with $ {\hat{\mathbf{B}}}_{i,j}\sim \mathcal{N}\left({\mathbf{B}}_{i,j},{\mathbf{B}}_{i,j}\right) $, where $ \left(i,j\right)\hskip0.2em \in {\left\{1,\dots, M\right\}}^2 $, thus considering:

(3)$$ {\mathbf{Y}}_t=P\left({\mathcal{M}}_q\left(\mathcal{H}\left({\mathbf{X}}_t\right)\right)\right)+\hat{\mathbf{B}}+{\mathbf{E}}_t,\hskip2em \forall t=1,2,\dots, T. $$

By now further approximating the variance of $ \hat{\mathbf{B}} $ with a constant $ b\hskip0.2em \in {\unicode{x211D}}_{+} $ to be interpreted as the average of $ \mathbf{B} $, we have that by simple manipulations:

$$ \hat{\mathbf{B}}+{\mathbf{E}}_t=\mathbf{B}+{\mathbf{N}}_t, $$

where the independence between $ \hat{\mathbf{B}} $ and $ {\mathbf{E}}_t $ has been exploited. We can thus retrieve (Reference Sage, Kirshner and Pengo1) from (Reference Betzig, Patterson and Sougrat3) by neglecting the Poisson noise dependence in $ P\left({\mathcal{M}}_q\left(\mathcal{H}\left({\mathbf{X}}_t\right)\right)\right) $ and that the variance of every entry of the random term $ {\mathbf{N}}_t $ is $ s={\sigma}^2+b $. A more detailed and less approximated modeling taking into account the signal-dependent nature of the noise in the data could represent a very interesting area of future research.

In vectorized form, model (Reference Sage, Kirshner and Pengo1) reads:

(4)$$ {\mathbf{y}}_t=\boldsymbol{\Psi} {\mathbf{x}}_t+\mathbf{b}+{\mathbf{n}}_t, $$

where $ \boldsymbol{\Psi} \hskip0.2em \in {\unicode{x211D}}^{M^2\times {L}^2} $ is the matrix representing the composition $ {\mathcal{M}}_q\circ \mathcal{H} $, while $ {\mathbf{y}}_t\hskip0.2em \in {\unicode{x211D}}^{M^2} $, $ {\mathbf{x}}_t\hskip0.2em \in {\unicode{x211D}}^{L^2} $, $ \mathbf{b}\hskip0.2em \in {\unicode{x211D}}^{M^2} $, and $ {\mathbf{n}}_t\hskip0.2em \in {\unicode{x211D}}^{M^2} $ are the column-wise vectorizations of $ {\mathbf{Y}}_t $, $ {\mathbf{X}}_t $, $ \mathbf{B} $, and $ {\mathbf{N}}_t $ in (Reference Sage, Kirshner and Pengo1), respectively.

For all $ t $ and given $ \boldsymbol{\Psi} $ and $ {\mathbf{y}}_t $, the problem can thus be formulated as

$$ \mathrm{find}\hskip1em \mathbf{x}=\frac{1}{T}\sum \limits_{t=1}^T\;{\mathbf{x}}_t\hskip0.2em \in {\unicode{x211D}}^{L^2},\mathbf{b}\hskip0.2em \in {\unicode{x211D}}^{M^2}\mathrm{and}\hskip0.5em s>0\hskip1em \mathrm{s}.\mathrm{t}.\hskip1em {\mathbf{x}}_t\hskip1em \mathrm{s}\mathrm{olves}\;(4). $$

In order to exploit the statistical behavior of the fluorescent emitters, we reformulate the model in the covariance domain. This idea was previously exploited by the SOFI approach⁽Reference Dertinger, Colyer, Iyer, Weiss and Enderlein¹¹⁾ and was shown to significantly reduce the full-width-at-half-maximum (FWHM) of the PSF. In particular, the use of second-order statistics for a Gaussian PSF corresponds to a reduction factor of the FWHM of $ \sqrt{2} $.

To formulate the model, we consider the frames $ {\left({\mathbf{y}}_t\right)}_{t\hskip0.2em \in \left\{1,\dots, T\right\}} $ as $ T $ realizations of a random variable $ \mathbf{y} $ with covariance matrix defined by:

(5)$$ {\mathbf{R}}_{\mathbf{y}}={\unicode{x1D53C}}_{\mathbf{y}}\left\{\left(\mathbf{y}-{\unicode{x1D53C}}_{\mathbf{y}}\left\{\mathbf{y}\right\}\right){\left(\mathbf{y}-{\unicode{x1D53C}}_{\mathbf{y}}\left\{\mathbf{y}\right\}\right)}^{\intercal}\right\}, $$

where $ {\unicode{x1D53C}}_{\mathbf{y}}\left\{\cdot \right\} $ denotes the expected value computed w.r.t. to the unknown law of $ \mathbf{y} $. We estimate $ {\mathbf{R}}_{\mathbf{y}} $ by computing the empirical covariance matrix, that is,

$$ {\mathbf{R}}_{\mathbf{y}}\approx \frac{1}{T-1}\sum \limits_{t=1}^T\left({\mathbf{y}}_t-\overline{\mathbf{y}}\right){\left({\mathbf{y}}_t-\overline{\mathbf{y}}\right)}^{\intercal }, $$

where $ \overline{\mathbf{y}}=\frac{1}{T}{\sum}_{t=1}^T{\mathbf{y}}_t $ denotes the empirical temporal mean. From (Reference Rust, Bates and Zhuang4) and (Reference Li and Vaughan5), we thus deduce the relation:

(6)$$ {\mathbf{R}}_{\mathbf{y}}=\boldsymbol{\Psi} {\mathbf{R}}_{\mathbf{x}}{\boldsymbol{\Psi}}^{\intercal }+{\mathbf{R}}_{\mathbf{n}}, $$

where $ {\mathbf{R}}_{\mathbf{x}}\hskip0.2em \in {\unicode{x211D}}^{L^2\times {L}^2} $ and $ {\mathbf{R}}_{\mathbf{n}}\hskip0.2em \in {\unicode{x211D}}^{M^2\times {M}^2} $ are the covariance matrices of $ {\left({\mathbf{x}}_t\right)}_{t\hskip0.2em \in \left\{1,\dots, T\right\}} $ and $ {\left({\mathbf{n}}_t\right)}_{t\hskip0.2em \in \left\{1,\dots, T\right\}}, $ respectively. As the background is stationary by assumption, the covariance matrix of $ \mathbf{b} $ is zero. Recalling now that the emitters are uncorrelated by assumption, we deduce that $ {\mathbf{R}}_{\mathbf{x}} $ is diagonal. We thus set $ {\mathbf{r}}_{\mathbf{x}}:= \operatorname{diag}\left({\mathbf{R}}_{\mathbf{x}}\right)\hskip0.2em \in {\unicode{x211D}}^{L^2} $. Furthermore, by the i.i.d. assumption on $ {\mathbf{n}}_t $, we have that $ {\mathbf{R}}_{\mathbf{n}}=s{\mathbf{I}}_{{\mathbf{M}}^{\mathbf{2}}} $, where $ s\hskip0.2em \in {\unicode{x211D}}_{+} $ and $ {\mathbf{I}}_{{\mathbf{M}}^{\mathbf{2}}} $ is the identity matrix in $ {\unicode{x211D}}^{M^2\times {M}^2} $. Note that the model in equation (6) is similar to the SPARCOM one presented in Ref. (Reference Solomon, Eldar, Mutzafi and Segev18), with the difference that here we consider also noise contributions by including in the model the diagonal covariance matrix $ {\mathbf{R}}_{\mathbf{n}} $. Finally, the vectorized form of the model in the covariance domain can thus be written as:

$$ {\mathbf{r}}_{\mathbf{y}}=\left(\boldsymbol{\Psi} \odot \boldsymbol{\Psi} \right){\mathbf{r}}_{\mathbf{x}}+s{\mathbf{v}}_{\mathbf{I}}, $$

where $ \odot $ denotes the Khatri–Rao (column-wise Kronecker) product, $ {\mathbf{r}}_{\mathbf{y}}\hskip0.2em \in {\unicode{x211D}}^{M^4} $ is the column-wise vectorization of $ {\mathbf{R}}_{\mathbf{y}} $ and $ {\mathbf{v}}_{\mathbf{I}}=\mathrm{vec}\left({\mathbf{I}}_{{\mathbf{M}}^{\mathbf{2}}}\right) $.

3. COL0RME, Step I: Support Estimation for Precise Molecule Localization

Similarly to SPARCOM⁽Reference Solomon, Eldar, Mutzafi and Segev¹⁸⁾, our approach makes use of the fact that the solution $ {\mathbf{r}}_{\mathbf{x}} $ is sparse, while including further the estimation of $ s>0 $ for dealing with more challenging scenarios. In order to compare specific regularity a priori constraints on the solution, we make use of different regularization terms, whose importance is controlled by a regularization hyperparameter $ \lambda >0 $. By further introducing some non-negativity constraints for both variables $ {\mathbf{r}}_{\mathbf{x}} $ and $ s $, we thus aim to solve:

(7)$$ \underset{{\mathbf{r}}_{\mathbf{x}}\ge 0,\hskip0.3em s\ge 0}{\arg \hskip0.1em \min}\mathcal{F}\left({\mathbf{r}}_{\mathbf{x}},s\right)+\mathcal{R}\left({\mathbf{r}}_{\mathbf{x}};\lambda \right), $$

where the data fidelity term is defined by:

(8)$$ \mathcal{F}\left({\mathbf{r}}_{\mathbf{x}},s\right)=\frac{1}{2}\parallel {\mathbf{r}}_{\mathbf{y}}-\left(\boldsymbol{\Psi} \odot \boldsymbol{\Psi} \right){\mathbf{r}}_{\mathbf{x}}-s{\mathbf{v}}_{\mathbf{I}}{\parallel}_2^2, $$

and $ \mathcal{R}\left(\cdot; \lambda \right) $ is a sparsity-promoting penalty. Ideally, one would like to make use of the $ {\mathrm{\ell}}_0 $ norm to enforce sparsity. However, as it is well-known, solving the resulting noncontinuous, nonconvex, and combinatorial minimization problem is an NP-hard problem. A way to circumvent this difficulty consists in using the continuous exact relaxation of the $ {\mathrm{\ell}}_0 $ norm (CEL0) proposed by Soubies et al. in Ref. (Reference Soubies, Blanc-Féraud and Aubert22). The CEL0 regularization is continuous, nonconvex and preserves the global minima of the original $ {\mathrm{\ell}}_2-{\mathrm{\ell}}_0 $ problem while removing some local ones. It is defined as follows:

(9)$$ \mathcal{R}\left({\mathbf{r}}_{\mathbf{x}};\lambda \right)={\Phi}_{\mathrm{CEL}0}\left({\mathbf{r}}_{\mathbf{x}};\lambda \right)=\sum \limits_{i=\unicode{x1D7D9}}^{L^2}\lambda -\frac{\parallel {\mathbf{a}}_i{\parallel}^2}{2}{\left(|{\left({\mathbf{r}}_{\mathbf{x}}\right)}_i|-\frac{\sqrt{2\lambda }}{\parallel {\mathbf{a}}_i\parallel}\right)}^2{\unicode{x1D7D9}}_{\left\{|{\left({\mathbf{r}}_{\mathbf{x}}\right)}_i|\le \frac{\sqrt{2\lambda }}{\parallel {\mathbf{a}}_i\parallel}\right\}}, $$

where $ {\mathbf{a}}_i={\left(\boldsymbol{\Psi} \odot \boldsymbol{\Psi} \right)}_i $ denotes the $ i $th column of the operator $ \mathbf{A}:= \boldsymbol{\Psi} \odot \boldsymbol{\Psi} $.

A different, convex way of favoring sparsity consists in taking as regularizer the $ {\mathrm{\ell}}_1 $ norm, that is,

(10)$$ \mathcal{R}\left({\mathbf{r}}_{\mathbf{x}};\lambda \right)=\lambda \parallel {\mathbf{r}}_{\mathbf{x}}{\parallel}_1. $$

Besides convexity and as it is well-known, the key difference between using the $ {\mathrm{\ell}}_0 $ and the $ {\mathrm{\ell}}_1 $-norm is that the $ {\mathrm{\ell}}_0 $ provides a correct interpretation of sparsity by counting only the number of the nonzero coefficients, while the $ {\mathrm{\ell}}_1 $ depends also on the magnitude of the coefficients. However, its use as a sparsity-promoting regularizer is nowadays well-established (see, e.g., Ref. (Reference Candès, Wakin and Boyd23)) and also used effectively in other microscopy applications, such as SPARCOM⁽Reference Solomon, Eldar, Mutzafi and Segev¹⁸⁾.

Finally, in order to model situations where piece-wise constant structures are considered, we consider a different regularization term favoring gradient-sparsity using the total variation (TV) regularization defined in a discrete setting as follows:

(11)$$ \mathcal{R}\left({\mathbf{r}}_{\mathbf{x}};\lambda \right)=\lambda TV\left({\mathbf{r}}_{\mathbf{x}}\right)=\lambda \sum \limits_{i=1}^{L^2}{\left({\left|{\left({\mathbf{r}}_{\mathbf{x}}\right)}_i-{\left({\mathbf{r}}_{\mathbf{x}}\right)}_{n_{i,1}}\right|}^2+{\left|{\left({\mathbf{r}}_{\mathbf{x}}\right)}_i-{\left({\mathbf{r}}_{\mathbf{x}}\right)}_{n_{i,2}}\right|}^2\right)}^{\frac{1}{2}}, $$

where $ \left({n}_{i,1},{n}_{i,2}\right)\hskip0.2em \in {\left\{1,\dots, {L}^2\right\}}^2 $ indicate the locations of the horizontal and vertical nearest neighbor pixels of pixel $ i $, as shown in Figure 2. For the computation of the TV penalty, Neumann boundary conditions have been used.

Figure 2. The one-sided nearest horizontal and vertical neighbors of the pixel $ i $ used to compute the gradient discretization in (Reference Dertinger, Colyer, Iyer, Weiss and Enderlein11).

To solve (Reference Gustafsson7), we use the alternate minimization algorithm between $ s $ and $ {\mathbf{r}}_{\mathbf{x}} $⁽Reference Attouch, Bolte, Redont and Soubeyran²⁴⁾, see the pseudo-code reported in Algorithm 1. Note that, at each $ k\ge 1 $, the update for the variable $ s $ can be efficiently computed through the following explicit expression:

$$ {s}^{k+1}=\frac{1}{M^2}{{\mathbf{v}}_{\mathbf{I}}}^{\intercal}\left({\mathbf{r}}_{\mathbf{y}}-\left(\boldsymbol{\Psi} \odot \boldsymbol{\Psi} \right){{\mathbf{r}}_{\mathbf{x}}}^k\right). $$

Concerning the update of $ {\mathbf{r}}_{\mathbf{x}}\hskip-0.2em $, different algorithms were used depending on the choice of the regularization term in (Reference Denoyelle, Duval, Peyré and Soubies9)–(Reference Dertinger, Colyer, Iyer, Weiss and Enderlein11). For the CEL0 penalty (Reference Denoyelle, Duval, Peyré and Soubies9), we used the iteratively reweighted $ {\mathrm{\ell}}_1 $ algorithm (IRL1)⁽Reference Ochs, Dosovitskiy, Brox and Pock²⁵⁾, following Gazagnes et al. ⁽Reference Gazagnes, Soubies and Blanc-Féraud²⁶⁾ with fast iterative shrinkage-thresholding algorithm (FISTA)⁽Reference Beck and Teboulle²⁷⁾ as inner solver. If the $ {\mathrm{\ell}}_1 $ norm (Reference Holden, Uphoff and Kapanidis10) is chosen, FISTA is used. Finally, when the TV penalty (Reference Dertinger, Colyer, Iyer, Weiss and Enderlein11) is employed, the primal-dual splitting method in Ref. (Reference Condat28) was considered.

Following the description provided by Attouch et al. in Ref. (Reference Attouch, Bolte, Redont and Soubeyran24), convergence of Algorithm 1 can be guaranteed only if an additional quadratic term is introduced in the objective function of the second minimization subproblem. Nonetheless, empirical convergence was observed also without such additional terms.

To evaluate the performance of the first step of the method COL0RME using the different regularization penalties described above, we created two noisy simulated datasets, with low background (LB) and high background (HB), respectively and used them to apply COL0RME and estimate the desired sample support. More details on the two datasets are available in the following Section 6.1. The results obtained using the three different regularizers are reported in Figure 3. In this example, we chose the regularization parameter $ \lambda $ heuristically, while more details about the selection of the parameter are given in the Section 5.1.

Figure 3. (a) Noisy simulated dataset with low background (LB) and stack size: $ T=500 $ frames, (b) Noisy simulated high-background (HB) dataset, with $ T=500 $ frames. From left to right: Superimposed diffraction limited image (temporal mean of the stack) with 4× zoom on ground truth support (blue), CEL0 reconstruction, $ {\mathrm{\ell}}_1 $ reconstruction and total variation (TV) reconstruction.

Despite its continuous and smooth reconstruction, we observe that the reconstruction obtained by the TV regularizer does not provide precise localization results. For example, the separation of the two filaments on the top-right corner is not visible and while the junction of the other two filaments on the bottom-left should appear further down, we clearly see that those filaments are erroneously glued together. Nonetheless, the choice of an appropriate regularizer tailored to favor fine structures as the ones observed in the GT image constitutes a challenging problem that should be addressed in future research.

The Jaccard indices (JIs) of both the results obtained when using the CEL0 and $ {\mathrm{\ell}}_1 $ regularizer, that allow for more precise localization, have been computed. The JI is a quantity in the range $ \left[0,1\right] $ computed as the ratio between correct detections (CD) and the sum of correct detections, false positives (FP) and false negatives (FN), that is, $ \mathrm{JI}:= \mathrm{CD}/\left(\mathrm{CD}+\mathrm{FN}+\mathrm{FP}\right), $ up to a tolerance $ \delta >0 $, measure in nm. A correct detection occurs when one pixel at most $ \delta $ nm away from a ground truth pixel is added to the support. In order to match the pixels from the estimated support to the ones from the ground truth, we employ the standard Gale–Shapley algorithm⁽Reference Gale and Shapley²⁹⁾. Once the matching has been performed, we can simply count the number of ground truth pixels which have not been detected (false negatives) and also the number of pixels in the estimated support which have not been matched to any ground truth pixel (false positives).

Figure 4 reports the average JI computed from 20 different noise realizations, as well as, an error bar (vertical lines) that represent the standard deviation, for several stack sizes. According to the figure, a slightly better JI is obtained when the CEL0 regularizer is being used, while an increase in the number of frames, when both regularizers being used, leads to better JI, hence better localization. As the reader may notice, such quantitative assessment could look inconsistent with the visual results reported in Figure 3. By definition, the JI tends to assume higher values whenever more CD are found even in presence of more FP (as it happens for the CEL0 reconstruction), while it gets more penalized when FN happen, as they affect the computation “twice,” reducing the numerator and increasing the denominator.

Figure 4. Jaccard Index values with tolerance $ \delta =40\;\mathrm{nm} $ for the low-background (LB) and high-background (HB) datasets, for different stack sizes and regularization penalty choices. The tolerance, $ \delta =40 $ nm, is set so that we allow the correct detections, that needed to be counted for the computation of the Jaccard Index, to be found not only in the same pixel but also to any of the 8-neighboring pixels.

3.1. Accurate noise variance estimation

Along with the estimations of the emitter’s temporal sparse covariance matrix, the estimation of the noise variance in the joint model (Reference Gustafsson7) allows for much more precise results even in challenging acquisition conditions. In Figure 5, we show the relative error between the computed noise variance $ s $ and the constant variance of the electronic noise $ {\sigma}^2 $ used to produce simulated LB and HB data. The relative error is higher in the case of the HB dataset, something that is, expected, as in our noise variance estimation $ s $ there is a bias coming from the background (see Remark 1). In the case of the LB dataset, as the background is low, the bias is sufficiently small so that it is barely visible in the error graph. In our experiments, a Gaussian noise with a corresponding SNR of approximately 16 dB is being used, while the value of $ {\sigma}^2 $ is in average equal to $ 7.11\times {10}^5 $ for the LB dataset and $ 7.13\times {10}^5 $ for the HB dataset. Note that, in general, the estimation of the noise variance $ s $ obtained by COL0RME is very precise.

Figure 5. The relative error in noise variance estimation, defined as: Error = $ \frac{\mid s-{\sigma}^2\mid }{\mid {\sigma}^2\mid } $, where $ {\sigma}^2 $ is the constant variance of the electronic noise. The error is computed for 20 different noise realizations, presenting in the graph the mean and the standard deviation (error bars).

4. COL0RME, Step II: Intensity Estimation

From the previous step, we obtain a sparse estimation of $ {\mathbf{r}}_{\mathbf{x}}\hskip0.2em \in {\unicode{x211D}}^{L^2} $. Its support, that is, the location of nonzero variances, can thus be deduced. This is denoted in the following by $ \Omega := \left\{i:{\left({\mathbf{r}}_{\mathbf{x}}\right)}_i\ne 0\right\}\subset \left\{1,\dots, {L}^2\right\} $. Note that this set corresponds indeed to the support of the desired $ \mathbf{x} $, hence in the following we will use the same notation to denote both sets.

We are now interested in enriching COL0RME with an additional step where intensity information of the signal $ \mathbf{x} $ can be retrieved in correspondence with the estimated support $ \Omega $. To do so, we thus propose an intensity estimation procedure for $ \mathbf{x} $ restricted only to the pixels of interest. Under this modeling assumption, it is thus reasonable to consider a regularization term favoring smooth intensities on $ \Omega $, in agreement to the intensity typically found in real images.

In order to take into account the modeling of blurry and out-of-focus fluorescent molecules, we further include in our model (Reference Rust, Bates and Zhuang4) a regularization term for smooth background estimation. We can thus consider the following joint minimization problem:

(12)$$ \underset{\mathbf{x}\hskip0.2em \in {\unicode{x211D}}_{+}^{\mid \Omega \mid },\mathbf{b}\hskip0.2em \in {\unicode{x211D}}_{+}^{M^2}}{\arg \hskip0.1em \min}\frac{1}{2}\parallel {\boldsymbol{\Psi}}_{\boldsymbol{\Omega}}\mathbf{x}-\left(\overline{\mathbf{y}}-\mathbf{b}\right){\parallel}_2^2+\frac{\mu }{2}\parallel {\nabla}_{\Omega}\mathbf{x}{\parallel}_2^2+\frac{\beta }{2}\parallel \nabla \mathbf{b}{\parallel}_2^2, $$

where the data term models the presence of Gaussian noise, $ \mu, \beta >0 $ are regularization parameters and the operator $ {\boldsymbol{\Psi}}_{\boldsymbol{\Omega}}\hskip0.2em \in {\unicode{x211D}}^{M^2\times \mid \Omega \mid } $ is a matrix whose $ i $th column is extracted from $ \boldsymbol{\Psi} $ for all indices $ i\hskip0.2em \in \Omega $. Finally, the regularization term on $ \mathbf{x} $ is the squared norm of the discrete gradient restricted to $ \Omega $, that is,

$$ \parallel {\nabla}_{\Omega}\mathbf{x}{\parallel}_2^2:= \sum \limits_{i\hskip0.2em \in \Omega}\sum \limits_{j\hskip0.2em \in \mathcal{N}(i)\cap \Omega}{\left({x}_i-{x}_j\right)}^2, $$

where $ \mathcal{N}(i) $ denotes the 8-pixel neighborhood of $ i\hskip0.2em \in \Omega $. Note that, according to this definition, $ {\nabla}_{\Omega}\mathbf{x} $ denotes a (redundant) isotropic discretization of the gradient of $ \mathbf{x} $ evaluated for each pixel in the support $ \Omega $. Note that this definition coincides with the standard one for $ \nabla \mathbf{x} $ restricted to points in the support $ \Omega $.

The non-negativity constraints on $ \mathbf{x} $ and $ \mathbf{b} $ as well as the one restricting the estimation of $ \mathbf{x} $ on $ \Omega $ can be relaxed using suitable smooth penalty terms, so that, finally, the following optimization problem can be addressed:

(13)$$ \underset{\mathbf{x}\hskip0.2em \in {\unicode{x211D}}^{L^2},\mathbf{b}\hskip0.2em \in {\unicode{x211D}}^{M^2}}{\arg \hskip0.1em \min}\frac{1}{2}\parallel \boldsymbol{\Psi} \mathbf{x}-\left(\overline{\mathbf{y}}-\mathbf{b}\right){\parallel}_2^2+\frac{\mu }{2}\parallel \nabla \mathbf{x}{\parallel}_2^2+\frac{\beta }{2}\parallel \nabla \mathbf{b}{\parallel}_2^2+\frac{\alpha }{2}\left(\parallel {\mathbf{I}}_{\boldsymbol{\Omega}}\mathbf{x}{\parallel}_2^2+\sum \limits_{i=1}^{L^2}\;{\left[\phi \left({\mathbf{x}}_i\right)\right]}^2+\sum \limits_{i=1}^{M^2}\;{\left[\phi \left({\mathbf{b}}_i\right)\right]}^2\right), $$

where the parameter $ \alpha \gg 1 $ can be chosen arbitrarily high to enforce the constraints, $ {\mathbf{I}}_{\boldsymbol{\Omega}} $ is a diagonal matrix acting as characteristic function of $ \Omega $, that is, defined as:

and $ \phi :\unicode{x211D}\to \unicode{x211D} $ is used to penalize negative entries, being defined as:

(14)$$ \phi (z):= \left\{\begin{array}{cc}0& \mathrm{if}\;z\ge 0,\\ {}z& \mathrm{if}\;z<0,\end{array}\right.\hskip2em \forall z\hskip0.2em \in \unicode{x211D}. $$

We anticipate here that considering the unconstrained problem (Reference Gustafsson, Culley, Ashdown, Owen, Pereira and Henriques13) instead of the original, constrained, one (Reference Geissbuehler, Bocchio, Dellagiacoma, Berclaz, Leutenegger and Lasser12), will come in handy for the design of an automatic parameter selection strategy, as we further detail in Section 5.2.

To solve the joint-minimization problem (Reference Gustafsson, Culley, Ashdown, Owen, Pereira and Henriques13), we use the Alternate Minimization algorithm, see Algorithm 2. In the following subsections, we provide more details on the solution of the two minimization subproblems.

4.1. First subproblem: update of x

In order to find at each $ k\ge 1 $ the optimal solution $ {\mathbf{x}}^{k+1}\hskip0.2em \in {\unicode{x211D}}^{L^2} $ for the first subproblem, we need to solve a minimization problem of the form:

(15)$$ {\mathbf{x}}^{k+1}=\underset{\mathbf{x}\hskip0.2em \in {\unicode{x211D}}^{L^2}}{\arg \hskip0.1em \min }g\left(\mathbf{x};{\mathbf{b}}^k\right)+h\left(\mathbf{x}\right), $$

where, for $ {\mathbf{b}}^k\hskip0.2em \in {\unicode{x211D}}^{M^2} $ being fixed at each iteration $ k\ge 1 $, $ g\left(\cdot; {\mathbf{b}}^k\right):{\unicode{x211D}}^{M^2}\to {\unicode{x211D}}_{+} $ is a proper and convex function with Lipschitz gradient, defined as:

(16)$$ g\left(\mathbf{x};{\mathbf{b}}^k\right):= \frac{1}{2}\parallel \boldsymbol{\Psi} \mathbf{x}-\left(\overline{\mathbf{y}}-{\mathbf{b}}^k\right){\parallel}_2^2+\frac{\mu }{2}\parallel \nabla \mathbf{x}{\parallel}_2^2, $$

and where the function $ h:{\unicode{x211D}}^{L^2}\to \unicode{x211D} $ encodes the penalty terms:

(17)$$ h\left(\mathbf{x}\right)=\frac{\alpha }{2}\left(\parallel {\mathbf{I}}_{\boldsymbol{\Omega}}\mathbf{x}{\parallel}_2^2+\sum \limits_{i=1}^{L^2}\;{\left[\phi \left({\mathbf{x}}_i\right)\right]}^2\right). $$

Solution of (Reference Yahiatene, Hennig, Müller and Huser15) can be obtained iteratively, using, for instance, the proximal gradient descent algorithm, whose iteration can be defined as follows:

(18)$$ {\mathbf{x}}^{n+1}={\mathbf{prox}}_{h,\tau}\left({\mathbf{x}}^n-\tau \nabla g\left({\mathbf{x}}^n\right)\right),\hskip1em n=1,2,\dots, $$

where $ \nabla g\left(\cdot \right) $ denotes the gradient of $ g $, $ \tau \hskip0.2em \in \left(0,\frac{1}{L_g}\right] $ is the algorithmic step-size chosen inside a range depending on the Lipschitz constant of $ \nabla g $, here denoted by $ {L}_g $, to guarantee convergence. The proximal update in (Reference Solomon, Eldar, Mutzafi and Segev18) can be computed explicitly using the computations reported in Appendix A. One can show in fact that, for each $ \mathbf{w}\hskip0.2em \in {\unicode{x211D}}^{L^2} $ there holds element-wise:

(19)$$ {\left({\mathrm{prox}}_{h,\tau}\left(\mathbf{w}\right)\right)}_i={\mathrm{prox}}_{h,\tau}\left({\mathbf{w}}_i\right)=\left\{\begin{array}{cc}\frac{{\mathbf{w}}_i}{1+\alpha \tau {\mathbf{I}}_{\boldsymbol{\Omega}}\left(i,i\right)}& \hskip-0.5em \mathrm{if}\hskip0.5em {\mathbf{w}}_i\ge 0,\\ {}\frac{{\mathbf{w}}_i}{1+\alpha \tau \left({\mathbf{I}}_{\boldsymbol{\Omega}}\left(i,i\right)+1\right)}& \mathrm{if}\hskip0.5em {\mathbf{w}}_i<0.\end{array}\right. $$

Remark 2. As the reader may have noted, we consider the proximal gradient descent algorithm (Reference Solomon, Eldar, Mutzafi and Segev18) for solving (Reference Yahiatene, Hennig, Müller and Huser15), even though both functions $ g $ and $ h $ in (Reference Deng, Sun, Lin, Ma and Shaevitz16) and (Reference Solomon, Mutzafi, Segev and Eldar17) respectively, are smooth and convex, hence, in principle, (accelerated) gradient descent algorithms could be used. Note, however, that the presence of the large penalty parameter $ \alpha \gg 1 $ would significantly slow down convergence speed in such case as the step size $ \tau $ in this case would be constrained to the smaller range $ \left(0,\frac{1}{L_g+\alpha}\right] $. By considering the penalty contributions in terms of their proximal operators, this limitation does not affect the range of $ \tau $ and convergence is still guaranteed⁽ Reference Combettes and Wajs³⁰⁾ in a computationally fast way through the update (Reference Dardikman-Yoffe and Eldar19).

4.2. Second subproblem: update of b

As far as the estimation of the background is concerned, the minimization problem we aim to solve at each $ k\ge 1 $ takes the form:

(20)$$ {\mathbf{b}}^{k+1}=\underset{\mathbf{b}\hskip0.2em \in {\unicode{x211D}}^{M^2}}{\arg \hskip0.1em \min }r\left(\mathbf{b};{\mathbf{x}}^{k+1}\right)+q\left(\mathbf{b}\right), $$

where

$$ r\left(\mathbf{b};{\mathbf{x}}^{k+1}\right):= \frac{1}{2}\parallel \mathbf{b}-\left(\overline{\mathbf{y}}-\boldsymbol{\Psi} {\mathbf{x}}^{k+1}\right){\parallel}_2^2+\frac{\beta }{2}\parallel \nabla \mathbf{b}{\parallel}_2^2,\hskip2em q\left(\mathbf{b}\right):= \frac{\alpha }{2}\sum \limits_{i=1}^{M^2}\;{\left[\phi \left({\mathbf{b}}_i\right)\right]}^2. $$

Note that $ r\left(\cdot; {\mathbf{x}}^{k+1}\right):{\unicode{x211D}}^{M^2}\to {\unicode{x211D}}_{+} $ is a convex function with $ {L}_r $-Lipschitz gradient and $ q:{\unicode{x211D}}^{M^2}\to {\unicode{x211D}}_{+} $ encodes (large, depending on $ \alpha \gg 1 $) penalty contributions. Recalling Remark 2, we thus use again the proximal gradient descent algorithm for solving (Reference Goulart, Blanc-Féraud, Debreuve and Schaub20). The desired solution $ \hat{\mathbf{b}} $ at each $ k\ge 1 $ can thus be found by iterating:

(21)$$ {\mathbf{b}}^{n+1}={\mathbf{prox}}_{q,\delta}\left({\mathbf{b}}^n-\delta \nabla r\left({\mathbf{b}}^n\right)\right),\hskip1em n=1,2,\dots, $$

for $ \delta \hskip0.2em \in \left(0,\frac{1}{L_r}\right] $. The proximal operator $ {\mathbf{prox}}_{q,\delta}\left(\cdot \right) $, has an explicit expression and it is defined element-wise for $ i=1,\dots, {M}^2 $ as:

(22)$$ {\left({\mathbf{prox}}_{q,\delta}\left(\mathbf{d}\right)\right)}_i={\mathrm{prox}}_{q,\delta}\left({\mathbf{d}}_i\right)=\left\{\begin{array}{cc}\hskip-2.5em {\mathbf{d}}_i& \mathrm{if}\;{\mathbf{d}}_i\ge 0,\\ {}\frac{{\mathbf{d}}_i}{1+\alpha \delta}& \hskip0.18em \mathrm{if}\;{\mathbf{d}}_i<0.\end{array}\right. $$

4.3. Intensity and background estimation results

Intensity estimation results can be found in Figure 6 where (Reference Gustafsson, Culley, Ashdown, Owen, Pereira and Henriques13) is used for intensity/background estimation on the supports $ {\Omega}_{\mathcal{R}} $ estimated from the first step of COL0RME using $ \mathcal{R}= $CEL0, $ \mathcal{R}={\mathrm{\ell}}_1 $, and $ \mathcal{R}= $TV. We are referring to them as COL0RME-CEL0, COL0RME-$ {\mathrm{\ell}}_1 $, and COL0RME-TV, respectively. The colormap ranges are different for the coarse-grid and fine-grid representations, as explained in Section 6.1 The result on $ {\Omega}_{TV} $, even after the second step does not allow for the observation of a few significant details (e.g., the separation of the two filament on the bottom left corner) and that is, why it will not further discussed.

Figure 6. On top: Diffraction limited image $ \overline{\mathbf{y}}=\frac{1}{T}{\sum}_{t=1}^T{\mathbf{y}}_{\mathbf{t}} $, with T = 500 (4× zoom) for the low-background (LB) dataset and for the high-background (HB) dataset, ground truth (GT) intensity image. (a) Reconstructions for the noisy simulated dataset with LB. (b) Reconstruction for the noisy simulated dataset with HB. From left to right: intensity estimation result on estimated support using CEL0 regularization, $ {\mathrm{\ell}}_1 $ regularization and TV regularization. For all COL0RME intensity estimations, the same colorbar, presented at the bottom of the figure, has been used.

A quantitative assessment for the other two regularization penalty choices, $ {\Omega}_{\mathrm{CEL}0} $ and $ {\Omega}_{{\mathrm{\ell}}_1} $, is available in Figure 7. More precisely, we compute the peak-signal-to-noise-ratio (PSNR), given the following formula:

(23)$$ {\mathrm{PSNR}}_{\mathrm{dB}}=10{\log}_{10}\left(\frac{{\operatorname{MAX}}_{\mathbf{R}}^2}{\mathrm{MSE}}\right),\hskip2em \mathrm{MSE}=\frac{1}{L^2}\sum \limits_{i=1}^{L^2}{\left({\mathbf{R}}_i-{\mathbf{K}}_i\right)}^2, $$

where $ \mathbf{R}\hskip0.2em \in {\unicode{x211D}}^{L^2} $ is the reference image, $ \mathbf{K}\hskip0.2em \in {\unicode{x211D}}^{L^2} $ the image we want to evaluate using the PSNR metric and $ {\operatorname{MAX}}_{\mathbf{R}} $ the maximum value of the image $ \mathbf{R} $. In our case, the reference image is the ground truth intensity image: $ {\mathbf{x}}^{GT}\hskip0.2em \in {\unicode{x211D}}^{L^2} $. The higher the PSNR, the better the quality of the reconstructed image.

Figure 7. COL0RME peak-signal-to-noise-ratio (PSNR) values for two different datasets (low-background and high-background datasets), stack sizes and regularization penalty choices. The mean and the standard deviation of 20 different noise realizations are presented.

According to Figures 6 and 7, when only a few frames are considered (e.g., $ T=100 $ frames, high temporal resolution), the method performs better using the CEL0 penalty for the support estimation. However, when longer temporal sequences are available (e.g., $ T=500 $ or $ T=700 $ frames) the method performs better using the $ {\mathrm{\ell}}_1 $-norm instead. In addition to this, for both penalizations, PSNR improves as the number of temporal frames increases.

Background estimation results are available in Figure 8 where (Reference Gustafsson, Culley, Ashdown, Owen, Pereira and Henriques13) is used for intensity/background estimation on the supports $ {\Omega}_{\mathcal{R}} $, with $ \mathcal{R}= $ CEL0 and $ \mathcal{R}={\mathrm{\ell}}_1 $, that have been already estimated in the first step. In the figure, there is also the constant background generated by the SOFI Simulation Tool⁽Reference Girsault, Lukes and Sharipov³¹⁾, the software we used to generate our simulated data (more details in Section 6.1). Although the results look different due to the considered space-variant regularisation on $ \mathbf{b} $, the variations are very little. The estimated background is smooth, as expected, while higher values are estimated near the simulated filaments and values closer to the true background are found away from them.

Figure 8. (a) Low-background (LB) dataset: Diffraction limited image $ \overline{\mathbf{y}}=\frac{1}{T}{\sum}_{t=1}^T{\mathbf{y}}_{\mathbf{t}} $ with T = 500 (4× zoom), background estimation result on estimated support using CEL0 and $ {\mathrm{\ell}}_1 $ regularization, ground truth (GT) background image. (b) High-background (HB) dataset: Diffraction limited image $ \overline{\mathbf{y}}=\frac{1}{T}{\sum}_{t=1}^T{\mathbf{y}}_{\mathbf{t}} $ with T = 500 (4× zoom), Background estimation result on estimated support using CEL0 and $ {\mathrm{\ell}}_1 $ regularization, GT background image. Please note the different scales between the diffraction limited and background images for a better visualization of the results.

5. Automatic Selection of Regularization Parameters

We describe in this section two parameter selection strategies addressing the problem of estimating the regularization parameters $ \lambda $ and $ \mu $ appearing in the COL0RME support estimation problem (Reference Gustafsson7) and intensity estimation one (Reference Geissbuehler, Bocchio, Dellagiacoma, Berclaz, Leutenegger and Lasser12), respectively. The other two regularization parameters $ \beta $ and $ \alpha $ do not need fine tuning. They are both chosen arbitrary high, so as with large enough $ \beta $ to allow for a very smooth background and with very high $ \alpha $ to respect the required constraints (positivity for both intensity and background and restriction to the predefined support only for the intensity estimation).

5.1. Estimation of support regularization parameter λ

The selection of the regularization parameter value $ \lambda $ in (Reference Gustafsson7) is critical, as it determines the sparsity level of the support of the emitters. For its estimation, we start by computing a reference value $ {\lambda}_{max} $, defined as the smallest regularization parameter for which the identically zero solution is found. It is indeed possible to compute such a $ {\lambda}_{max} $ for both regularization terms CEL0 and $ {\mathrm{\ell}}_1 $ (see Refs. (Reference Soubies32) and (Reference Koulouri, Heins and Burger33)). Once such values are known, we thus need to find a fraction $ \gamma \hskip0.2em \in \left(0,1\right) $ of $ {\lambda}_{max} $ corresponding to the choice $ \lambda ={\gamma \lambda}_{max} $. For the CEL0 regularizer the expression for $ {\lambda}_{max} $ (see Proposition 10.9 in Ref. (Reference Soubies32)) is:

(24)$$ {\lambda}_{max}^{CEL0}:= \underset{1\le i\le {L}^2}{\mathit{\max}}\frac{{\left\langle {\mathbf{a}}_i,{\mathbf{r}}_{\mathbf{y}}\right\rangle}^2}{2\parallel {\mathbf{a}}_i{\parallel}^2}, $$

where $ {\mathbf{a}}_i={\left(\boldsymbol{\Psi} \odot \boldsymbol{\Psi} \right)}_i $ denotes the $ i $th column of the operator $ \mathbf{A}:= \boldsymbol{\Psi} \odot \boldsymbol{\Psi} $. Regarding the $ {\mathrm{\ell}}_1 $-norm regularization penalty, $ {\lambda}_{max} $ is given as follows:

(25)$$ {\lambda}_{max}^{{\mathrm{\ell}}_1}:= \parallel {\mathbf{A}}^{\intercal }{\mathbf{r}}_{\mathbf{y}}{\parallel}_{\infty }=\underset{1\le i\le {L}^2}{\max}\left\langle {\mathbf{a}}_i,{\mathbf{r}}_{\mathbf{y}}\right\rangle . $$

As far as $ {\mathrm{\ell}}_1 $ is used as regularization term in (Reference Gustafsson7), we report in Figure 9 a graph showing how the PSNR value of the final estimated intensity image (i.e., after the application of the second step of COL0RME) varies for the two datasets considered depending on $ \lambda $. It can be observed that for a large range of values $ \lambda $, the final PSNR remains almost the same. Although this may look a bit surprising at a first sight, we remark that such a robust result is due, essentially, to the second step of the algorithm where false localizations related to an underestimation of $ \lambda $ can be corrected through the intensity estimation step. Note, however, that in the case of an overestimation of $ \lambda $, points contained in the original support are definitively lost so no benefit is obtained from the intensity estimation step, hence the overall PSNR decreases.

Figure 9. The peak-signal-to-noise-ratio (PSNR) value of the final COL0RME image, using the $ {\mathrm{\ell}}_1 $-norm regularizer for support estimation, for different $ \gamma $ values, evaluating in both the low-background (LB) and high-background (HB) datasets. The mean and the standard deviation of 20 different noise realization are presented.

When the CEL0 penalty is used for support estimation, a heuristic parameter selection strategy can be used to improve the localization results, but also to avoid the fine parameter tuning. More specifically, the nonconvexity of the model can be used by considering an algorithmic restarting approach to improve the support reconstruction quality. In short, a value of $ \lambda $ can be fixed, typically $ \lambda ={\gamma \lambda}_{\mathrm{max}}^{CEL0} $ with $ \gamma \approx 5\times {10}^{-4} $, so as to achieve a very sparse reconstruction. Then, the support estimation algorithm can be run and iteratively repeated with a new initialization (i.e., restarted) several times. While keeping $ \lambda $ fixed along this procedure, a wise choice of the initialization depending, but not being equal to the previous output can be used to enrich the support, see Appendix C for more details. Nonconvexity is here exploited by changing, for a fixed $ \lambda $, the initialization at each algorithmic restart, so that new local minimizers (corresponding to possible support points) can be computed. The final support image can thus be computed as the superposition of the different solutions computed at each restarting. In such a way, a good result for a not-finely-tuned value of $ \lambda $ can be computed.

5.2. Estimation of intensity regularization parameter $ \mu $ by discrepancy principle

In this section, we provide some details on the estimation of the parameter $ \mu $ in (Reference Geissbuehler, Bocchio, Dellagiacoma, Berclaz, Leutenegger and Lasser12), which is crucial for an accurate intensity estimation. Recall that the problem we are looking at in this second step is

(26)$$ \mathrm{find}\hskip1em \mathbf{x}\hskip0.2em \in {\unicode{x211D}}^{L^2}\hskip1em \mathrm{s}.\mathrm{t}.\hskip1em \overline{\mathbf{y}}=\boldsymbol{\Psi} \mathbf{x}+\mathbf{b}+\overline{\mathbf{n}}, $$

where the quantities correspond to the temporal averages of the vectorized model in (Reference Rust, Bates and Zhuang4), so that $ \overline{\mathbf{n}}=\frac{1}{T}{\sum}_{t=1}^T{\mathbf{n}}_{\mathbf{t}} $. The temporal realizations $ {\mathbf{n}}_{\mathbf{t}} $ of the random vector $ \mathbf{n} $ follow a normal distribution with zero mean and covariance matrix $ s{\mathbf{I}}_{M^2} $, where $ s $ has been estimated in the first step of the algorithm, see Section 3.1. Consequently, the vector $ \overline{\mathbf{n}} $ follows also a normal distribution with zero mean and covariance matrix equal to $ \frac{s}{T}{\mathbf{I}}_{M^2} $. As both $ s $ and $ T $ are known, we can use the discrepancy principle, a well-known a posteriori parameter-choice strategy (see, e.g., Refs. (Reference Hansen34) and (Reference Gfrerer35)), to efficiently estimate the hyper-parameter $ \mu $. To detail how the procedure is applied to our problem, we write $ {\mathbf{x}}_{\mu } $ in the following to highlight the dependence of $ \mathbf{x} $ on $ \mu $. According to the discrepancy principle strategy, the regularization parameter $ \mu $ is chosen so that the residual norm of the regularized solution satisfies:

(27)$$ \parallel \overline{\mathbf{y}}-\Psi {\hat{\mathbf{x}}}_{\mu }-\hat{\mathbf{b}}{\parallel}_2^2={\nu}_{DP}^2\parallel \overline{\mathbf{n}}{\parallel}_2^2, $$

where $ {\hat{\mathbf{x}}}_{\mu}\hskip0.2em \in {\unicode{x211D}}^{L^2} $ and $ \hat{\mathbf{b}}\hskip0.2em \in {\unicode{x211D}}^{M^2} $ are the solutions of (Reference Geissbuehler, Bocchio, Dellagiacoma, Berclaz, Leutenegger and Lasser12). The expected value of $ \parallel \overline{\mathbf{n}}{\parallel}_2^2 $ is:

(28)$$ \unicode{x1D53C}\left\{\parallel \overline{\mathbf{n}}{\parallel}_2^2\right\}={M}^2\frac{s}{T}, $$

which can be used as an approximation of $ \parallel \overline{\mathbf{n}}{\parallel}_2^2 $ for $ {M}^2 $ big enough. The scalar value $ {\nu}_{DP}\approx 1 $ is a “safety factor” that plays an important role in the case when a good estimate of $ \parallel \overline{\mathbf{n}}{\parallel}_2 $ is not available. In such situations, a value $ {\nu}_{DP} $ closer to $ 2 $ is used. As detailed in Section 3.1, the estimation of $ s $ is rather precise in this case, hence we fix $ {\nu}_{DP}=1 $ in the following.

We can now define the function $ f\left(\mu \right):{\unicode{x211D}}_{+}\to \unicode{x211D} $ as:

(29)$$ f\left(\mu \right)=\frac{1}{2}\parallel \overline{\mathbf{y}}-\boldsymbol{\Psi} {\hat{\mathbf{x}}}_{\mu }-\hat{\mathbf{b}}{\parallel}_2^2-\frac{\nu_{DP}^2}{2}\parallel \overline{\mathbf{n}}{\parallel}_2^2. $$

We want to find the value $ \hat{\mu} $ such that $ f\left(\hat{\mu}\right)=0 $. This can be done iteratively, using the Newton’s method whose iterations read:

(30)$$ {\mu}_{n+1}={\mu}_n-\frac{f\left({\mu}_n\right)}{f^{\prime}\left({\mu}_n\right)},\hskip1em n=1,2,\dots . $$

In order to be able to compute easily the values $ f\left(\mu \right) $ and $ {f}^{\prime}\left(\mu \right) $, the values $ {\hat{\mathbf{x}}}_{\mu}\hskip0.2em \in {\unicode{x211D}}^{L^2} $, $ \hat{\mathbf{b}}\hskip0.2em \in {\unicode{x211D}}^{M^2} $, and $ {\hat{\mathbf{x}}}_{\mu}^{\prime }=\frac{\partial }{\partial \mu }{\hat{\mathbf{x}}}_{\mu}\hskip0.2em \in {\unicode{x211D}}^{L^2} $ need to be computed, as it can be easily noticed by writing the expression of $ {f}^{\prime}\left(\mu \right) $ which reads:

(31)$$ {f}^{\prime}\left(\mu \right)=\frac{\partial }{\partial \mu}\left\{\frac{1}{2}\parallel \overline{\mathbf{y}}-\boldsymbol{\Psi} {\hat{\mathbf{x}}}_{\mu }-\hat{\mathbf{b}}{\parallel}_2^2\right\}={\left({\hat{\mathbf{x}}}_{\mu}^{\prime}\right)}^{\intercal }{\boldsymbol{\Psi}}^{\intercal}\left(\overline{\mathbf{y}}-\boldsymbol{\Psi} {\hat{\mathbf{x}}}_{\mu }-\hat{\mathbf{b}}\right). $$

The values $ {\hat{\mathbf{x}}}_{\mu } $ and $ \hat{\mathbf{b}} $ can be found by solving the minimization problem (Reference Geissbuehler, Bocchio, Dellagiacoma, Berclaz, Leutenegger and Lasser12). As far as $ {\hat{\mathbf{x}}}_{\mu}^{\prime } $ is concerned, we report in Appendix B the steps necessary for its computation. We note here, however, that in order to compute such a quantity, the relaxation of the support/non-negativity constraints by means of the smooth quadratic terms discussed above is fundamental. One can show that $ {\hat{\mathbf{x}}}_{\mu}^{\prime } $ is the solution of the following minimization problem:

(32)$$ {\hat{\mathbf{x}}}_{\mu}^{\prime }=\underset{\mathbf{x}\hskip0.2em \in {\unicode{x211D}}^{L^2}}{\arg \hskip0.1em \min}\frac{1}{2}\parallel \boldsymbol{\Psi} \mathbf{x}{\parallel}_2^2+\frac{\mu }{2}\parallel \nabla \mathbf{x}+\mathbf{c}{\parallel}_2^2+\frac{\alpha }{2}\left(\parallel {\mathbf{I}}_{\boldsymbol{\Omega}}\mathbf{x}{\parallel}_2^2+\parallel {\mathbf{I}}_{{\hat{\mathbf{x}}}_{\mu }}\mathbf{x}{\parallel}_2^2\right), $$

where $ \mathbf{c} $ is a known quantity defined by $ \mathbf{c}=\frac{1}{\mu}\nabla {\hat{\mathbf{x}}}_{\mu } $, and the diagonal matrix $ {\mathbf{I}}_{{\hat{\mathbf{x}}}_{\mu }}\hskip0.2em \in {\unicode{x211D}}^{L^2\times {L}^2} $ identifies the support of $ {\hat{\mathbf{x}}}_{\mu } $ by:

$$ {\mathbf{I}}_{{\hat{\mathbf{x}}}_{\mu }}\left(i,i\right)=\left\{\begin{array}{cc}0& \mathrm{if}\;{\left({\hat{\mathbf{x}}}_{\mu}\right)}_{\mathrm{i}}\ge 0,\\ {}1& \hskip0.2em \mathrm{if}\;{\left({\hat{\mathbf{x}}}_{\mu}\right)}_{\mathrm{i}}<0.\end{array}\right. $$

We can find $ {\hat{\mathbf{x}}}_{\mu}^{\prime } $ by iterating

(33)$$ {{\mathbf{x}}^{\prime}}_{\mu}^{n+1}={\mathbf{prox}}_{\overline{h},\tau}\left({{\mathbf{x}}^{\prime}}_{\mu}^n-\tau \nabla \overline{g}\left({{\mathbf{x}}^{\prime}}_{\mu}^n\right)\right),\hskip1em n=1,2,\dots, $$

where

(34)$$ \overline{g}\left(\mathbf{x}\right):= \frac{1}{2}\parallel \boldsymbol{\Psi} \mathbf{x}{\parallel}_2^2+\frac{\mu }{2}\parallel \nabla \mathbf{x}+\mathbf{c}{\parallel}_2^2,\hskip2em \overline{h}\left(\mathbf{x}\right):= \frac{\alpha }{2}\left(\parallel {\mathbf{I}}_{\boldsymbol{\Omega}}\mathbf{x}{\parallel}_2^2+\parallel {\mathbf{I}}_{{\hat{\mathbf{x}}}_{\mu }}\mathbf{x}{\parallel}_2^2\right). $$

For $ \mathbf{z}\hskip0.2em \in {\unicode{x211D}}^{L^2} $, the proximal operator $ {\mathbf{prox}}_{\overline{h},\tau}\left(\mathbf{z}\right) $ can be obtained following the computations in Appendix A:

(35)$$ {\left({\mathbf{prox}}_{\overline{h},\tau}\left(\mathbf{z}\right)\right)}_i={\mathrm{prox}}_{\overline{h},\tau}\left({\mathbf{z}}_i\right)=\frac{{\mathbf{z}}_i}{1+\alpha \tau \left({\mathbf{I}}_{\boldsymbol{\Omega}}\left(i,i\right)+{\mathbf{I}}_{{\hat{\mathbf{x}}}_{\mu }}\Big(i,i\Big)\right)}, $$

while

(36)$$ \nabla \overline{g}\left({\mathbf{x}}^{\prime}\right)=\left({\boldsymbol{\Psi}}^{\intercal}\boldsymbol{\Psi} +\mu {\nabla}^{\intercal}\nabla \right){\mathbf{x}}^{\prime }+{\nabla}^{\intercal}\nabla {\hat{\mathbf{x}}}_{\mu }, $$

and the step $ \tau \hskip0.2em \in \left(0,\frac{1}{L_{\overline{g}}}\right] $, with $ {L}_{\overline{g}}=\parallel {\boldsymbol{\Psi}}^{\intercal}\boldsymbol{\Psi} +\mu {\nabla}^{\intercal}\nabla {\parallel}_2 $ the Lipschitz constant of $ \nabla \overline{g} $. A pseudo-code explaining the procedure we follow to find the optimal $ \hat{\mu} $ can be found in Algorithm 3. Finally, in Figure 10, a numerical example is available to show the good estimation of the parameter $ \hat{\mu} $.

Algorithm 3 Discrepancy Principle.

Require: $ \overline{\mathbf{y}}\hskip0.2em \in {\unicode{x211D}}^{M^2},{\mathbf{x}}^0\hskip0.2em \in {\unicode{x211D}}^{L^2},{\mathbf{b}}^0\hskip0.2em \in {\unicode{x211D}}^{M^2},{\mu}_0,\beta >0 $,$ \alpha \gg 1 $.

repeat.

Find $ {\hat{\mathbf{x}}}_{\mu_n},\hat{\mathbf{b}} $ using Algorithm 2

Find $ {\hat{\mathbf{x}}}_{\mu_n}^{\prime } $ solving (Reference Soubies32)

Compute $ f\left({\mu}_n\right),{f}^{\prime}\left({\mu}_n\right) $ from (Reference Gale and Shapley29) and (Reference Girsault, Lukes and Sharipov31)

$ {\mu}_{n+1}\leftarrow {\mu}_n-\frac{f\left({\mu}_n\right)}{f^{\prime}\left({\mu}_n\right)} $

until convergence.

return $ \hat{\mu} $.

Figure 10. The solid blue line shows the peak-signal-to-noise-ratio (PSNR) values computed by solving (Reference Gustafsson, Culley, Ashdown, Owen, Pereira and Henriques13) for several values of $ \mu $ within a specific range. Tha data used are the high-background (HB) dataset with $ T=500 $ frames (Figure 12c) and the $ {\mathrm{\ell}}_1 $-norm regularization penalty. The red cross shows the PSNR value $ \hat{\mu} $ obtained by applying the Discrepancy Principle. We note that such value is very close to one maximizing the PSNR metric.

7. Discussion and Conclusions

In this paper, we propose and discuss the model and the performance of COL0RME, a method for super-resolution microscopy imaging based on the sparse analysis of the stochastic fluctuations of molecules’ intensities. Similarly to other methods exploiting temporal fluctuations, COL0RME relaxes all the requirements for special equipment (microscope and fluorescent dyes) and allows for live-cell imaging, due to the good temporal resolution and the low power lasers employed. In comparison with competing methods, COL0RME achieves higher spatial resolution than other methods exploiting fluctuations while having a sufficient temporal resolution. COL0RME is based on two different steps: a former one where accurate molecule localization and noise estimation are achieved by solving nonsmooth convex/nonconvex optimization problems in the covariance domain and the latter where intensity information is retrieved in correspondence with the estimated support only. Our numerical results show that COL0RME outperforms competing approaches in terms of localization precision. To the best of our knowledge, COL0RME is the only super-resolution method exploiting temporal fluctuations which is capable of retrieving intensity-type information, signal and spatially varying background, which are of fundamental interest in biological data analysis. For both steps, automatic parameter selection strategies are detailed. Let us remark that such strategy of intensity estimation could be applied to the other competing super-resolution methods in the literature. Several results obtained on both simulated and real data are discussed, showing the superior performance of COL0RME in comparison with analogous methods such as SPARCOM, LSPARCOM and SRRF. Possible extensions of this work shall address the use of intensity information estimated by COL0RME for 3D reconstruction in, for example, MA-TIRF acquisitions. Furthermore, a systematic study to assess quantitatively the spatial resolution achieved by COL0RME under different scenarios (different background levels, different PSNRs, and number of frames) is envisaged.