# Optical diffraction tomography from single-molecule localization microscopy

Thanh-An Pham, Emmanuel Soubies, Ferréol Soulez, Michael Unser

► To cite this version:

Thanh-An Pham, Emmanuel Soubies, Ferréol Soulez, Michael Unser. Optical diffraction tomography from single-molecule localization microscopy. Optics Communications, 2021, 499, pp.127290. 10.1016/j.optcom.2021.127290. hal-03290069

HAL Id: hal-03290069

https://hal.science/hal-03290069

Submitted on 28 Sep 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Optical Diffraction Tomography from Single-Molecule Localization Microscopy

Thanh-an Pham<sup>a,*</sup>, Emmanuel Soubies<sup>b</sup>, Ferréol Soulez<sup>c</sup>, Michael Unser<sup>a</sup>

<sup>a</sup>Biomedical Imaging Group, École polytechnique fédérale de Lausanne, Lausanne, Switzerland.

<sup>b</sup>IRIT, Université de Toulouse, CNRS, Toulouse, France.

<sup>c</sup>Univ. Lyon, Univ. Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon, Saint-Genis-Laval, France.

## Abstract

Single-molecule localization microscopy (SMLM) is a powerful method for the imaging of cellular structures. This modality delivers nanoscale resolution by sequentially activating a subset of fluorescent molecules and by extracting their super-resolved positions from the microscope images. The emission patterns of individual molecules can be distorted by the refractive-index (RI) map of the sample, which reduces the accuracy of the molecule localization if not accounted for. In this work, we show that one can exploit those sample-induced aberrations to reveal the structural information of the specimen. Our work is related to the optical diffraction tomography in that we aim to recover the RI map. To that end, we propose an optimization framework in which we reconstruct the RI map and optimize the positions of the molecules in a joint fashion. The benefits of our method are twofold. On one side, we effectively recover the RI map of the sample. On the other side, we further improve the molecule localization—the primary purpose of SMLM. We validate our joint-optimization framework on simulated data. Our results lay the foundation of an exciting and novel extension of SMLM.

**Keywords:** computational imaging, joint optimization, refractive index reconstruction.

## 1. Introduction

Single-molecule localization microscopy (SMLM) is a method of choice for the observation of biological phenomena at nanoscale resolution [1, 2, 3]. It breaks the diffraction limit of conventional fluorescence microscopy by sequentially activating and localizing a subset of fluorescent molecules. As such, SMLM is a prime example of computational microscopy where suitable acquisitions and algorithmic reconstruction are combined so as to enhance the capabilities of traditional systems. Although SMLM acquisitions are two-dimensional, innovative point-spread functions (PSF) whose shapes vary with depth have been designed to encode the axial position of molecules. These include the popular astigmatism [4] or double-helix [5] PSFs. Therefore, in addition to efficient localization algorithms, well-calibrated models of these PSFs are essential to reach the promised nanoscale resolution [6].

The standard practice is to estimate these PSFs from acquisitions of sub-resolved objects (e.g., fluorescent microspheres) [7, 8]. However, this strategy ignores sample-induced distortions. Indeed, the heterogeneity of biological specimens—through variations in their refractive index (RI)—induces a scattering of the emitted light. This distorts the recorded emission patterns and compromises the accurate localization of molecules. To mitigate this effect, Xu et al. [9] proposed an algorithm to jointly localize fluorescent molecules and estimate an *in situ* PSF model that

has the ability to capture sample-induced aberrations and, hence, to improve localization accuracy.

If we could estimate both the RI and the position of molecules from the SMLM acquisition stack alone, then we would have a unique combination of structural (RI) and functional (fluorescence) information about the sample [10]. To our knowledge, such a reconstruction of both RI and fluorescence density from the same fluorescent dataset (i.e., without phase measurements) has been investigated only recently by Xue and Waller [11]. They consider two-layers samples where the bottom layer contains fluorescence-labeled objects and the top layer contains non-labelled objects. In this context, they demonstrated that the 3D RI map of the non-labelled objects can be reconstructed from defocused fluorescence images that are collected by sequentially stimulating small regions of the fluorescence-labeled layer of the sample. Moreover, they showed that the obtained RI map can be exploited to obtain the scattered PSF and improve the fluorescence signal through deconvolution. This setting differs from the exploitation of the individual emission of fluorophores in SMLM that we propose here in two respects. First, RI and fluorescence objects are mixed (i.e., no two-layers samples). Second, fluorescence measurements are recorded at two distinct focal planes (i.e., biplane SMLM modality).

In SMLM, the recovery of the RI has been addressed in [12, 13]. These two works exploit the fact that SMLM data can be seen as measurements of an optical-diffraction tomography (ODT) system with point-source illuminations

<footnote>*Corresponding author: thanh-an.pham@epfl.ch</footnote>
inside the sample. In [12], the authors assumed that the phase of the measurements was accessible, an assumption which is not met in practice. Moreover, their proposed approach relies on a linear model whose validity is limited to weakly scattering samples [14]. In our preliminary work [13], we introduced a refined RI-reconstruction approach that can handle phaseless measurements under the assumption that the positions of the fluorescent molecules are perfectly known. This method is based on the exact (nonlinear) Lippmann-Schwinger equation.

## 1.1. Contributions and Roadmap

In this work, we extend [13] in two ways. First, we consider a more realistic image-formation model (described in Section 2) that integrates background fluorescence as well as the shot noise inherent to fluorescence microscopy. Second, and more importantly, we consider that the positions of the molecules are known only approximately, as opposed to [13], and then take advantage of our model-based scheme to refine them.

To cope with this more challenging scenario, we propose a joint-optimization framework in Section 3. Our method simultaneously reconstructs the RI and refines the positions and amplitudes of the molecules. The benefits of our framework are twofold. On one side, we are accurately estimating the structural information (RI) from SMLM acquisitions. On the other side, we significantly improve the localization of the molecules—the primary objective of SMLM. We validate our framework on simulated data in Section 4.

## 1.2. Notations

Scalar and continuously defined functions are denoted by italic letter (e.g., $k_b \in \mathbb{R}$, $u \in L_2(\mathbb{R})$). The complex conjugate of $v \in \mathbb{C}$ is denoted by $v^*$. Vectors and matrices are denoted by bold lowercase and bold uppercase letters, respectively (e.g., $\mathbf{x} \in \mathbb{R}^N$, $\mathbf{G} \in \mathbb{C}^{N \times N}$). The $m$th element of a vector $\mathbf{v} \in \mathbb{R}^M$ is denoted as $v_m$ or $[\mathbf{v}]_m$. Similarly, the $n$th column of a matrix $\mathbf{X} \in \mathbb{R}^{M \times N}$ is denoted as $\mathbf{x}_n$ or $[\mathbf{X}]_n$. The notation $\mathbf{G}^H$ refers to the conjugate transpose of the matrix $\mathbf{G} \in \mathbb{C}^{N \times N}$. The matrix $\mathbf{I}_N \in \mathbb{R}^{N \times N}$ is the identity and $\mathbf{diag}(\mathbf{f}) \in \mathbb{R}^{N \times N}$ is a diagonal matrix formed out of the entries of $\mathbf{f} \in \mathbb{R}^N$. The notation $\mathbf{1}_M = (1, 1, \dots, 1) \in \mathbb{R}^M$ stands for an $M$-length vector of ones. Similarly, $\mathbf{0}_M$ denotes a vector of $M$ zeros. Finally, $\odot$ and $\otimes$ stand for the Hadamard product and the pointwise division, respectively.

# 2. Image-Formation Model

## 2.1. SMLM Meets ODT

The space-varying refractive index of the sample under consideration is represented by the function $\eta: \Omega \to \mathbb{R}$ with $\Omega \subseteq \mathbb{R}^3$. The sample is populated with $L$ fluorophores located at spatial position $\{\mathbf{x}_l \in \Omega\}_{l=1}^L$. Without loss of generality, we consider an SMLM acquisition stack where a

<figure><img src="image_2.png" /><figcaption>Figure 1: Biplane single-molecule localization microscopy. A fluorophore emits fluorescent light which scatters through the sample. Then, an optical system records the intensity of the total field at two different focal planes.</figcaption></figure>

single fluorophore is activated on each frame. Indeed, because fluorophores are incoherent sources, the image produced upon activation of multiple emitters is simply the sum of the individual contribution of each emitter [15].

When activated, the $l$th fluorophore at position $\mathbf{x}_l \in \Omega$ emits a spherical wave with intensity $a_l > 0$, which leads to

$$u_{\text{in}}(\mathbf{x}; \mathbf{x}_l, a_l) = a_l \frac{\exp(jk_b \|\mathbf{x} - \mathbf{x}_l\|_2)}{4\pi \|\mathbf{x} - \mathbf{x}_l\|_2}, \quad (1)$$

where $j$ is the imaginary unit and $k_b = \frac{2\pi\eta_b}{\lambda}$ is the wavenumber determined by the emission wavelength $\lambda$ and the RI $\eta_b > 1$ of the surrounding medium. The spherical wave acts as an “incident” field that illuminates from within the sample. As such, it scatters through the sample and produces a field $u_l: \mathbb{R}^3 \to \mathbb{C}$ that satisfies the Lippmann-Schwinger equation

$$u_l(\mathbf{x}) = u_{\text{in}}(\mathbf{x}; \mathbf{x}_l, a_l) + \int_{\Omega} g(\mathbf{x} - \mathbf{z}) f(\mathbf{z}) u_l(\mathbf{z}) \, d\mathbf{z}, \quad (2)$$

where $f(\mathbf{x}) = k_b^2 (\frac{\eta(\mathbf{x})^2}{\eta_b^2} - 1)$ is the scattering potential and $g: \mathbb{R}^3 \to \mathbb{C}$ is the Green function that corresponds to the centered spherical wave $u_{\text{in}}(\mathbf{x}; 0, 1)$ [16]. The intensity of the field $u_l$ at the camera plane $\Gamma$ is then recorded by an optical system to form the $l$th SMLM frame $\mathbf{y}_l \in \mathbb{R}^M$. Formally, we have, $\forall l \in \{1, \dots, L\}$, that

$$\mathbf{y}_l = \text{Pois} (|P u_l|_{\Gamma}^2 + \mathbf{b}_l), \quad (3)$$

where Pois denotes Poisson's distribution (shot noise), $u_l|_{\Gamma}$ denotes the restriction of $u_l$ to $\Gamma$, and $\mathbf{b}_l \in \mathbb{R}^M$ is a background signal that can originate from autofluorescence or spurious out-of-focus fluorophore emissions. Finally, $P: \mathbb{C}^2 \to \mathbb{R}^M$ is a linear operator that models both the effect of the optical system (i.e., pointwise multiplication with the pupil function in the Fourier domain) and the sampling on the $M$ camera pixels.
**Algorithm 1 Joint-Optimization Framework**

Require: $\mathbf{f}^0 \in \mathbb{R}_{\ge 0}^N$, $[\mathbf{x}_1^0 \cdots \mathbf{x}_L^0] \in \Omega^L$, $\mathbf{a}^0 \in \mathbb{R}_{> 0}^L$

1: $t = 0$

2: **while** (Not converged) **do**

3: Select a subset $\mathcal{L} \subset \{1, \dots, L\}$
   ▷ Update amplitudes and positions

4: **for** $l \in \mathcal{L}$ **do**

5: $(a_l^{t+1}, \mathbf{x}_l^{t+1}) = \text{Refine} (\mathcal{J}_l(\mathbf{f}^t, \cdot, \cdot); a_l^t, \mathbf{x}_l^t)$

6: **end for**

▷ Update the scattering potential:

7: $\mathbf{f}^{t+1} = \text{aFBS} (\sum_{l \in \mathcal{L}} \mathcal{J}_l(\cdot, \mathbf{x}_l^{t+1}, a_l^{t+1}); \mathbf{f}^t)$

8: $t \leftarrow t + 1$

9: **end while**

Return: $\mathbf{f}^t, \mathbf{X}^t, \mathbf{a}^t$

2.2. Discrete Forward Model

Let us rasterize $\Omega$ into $N$ voxels of length $h$. Following [17, 18, 19], we define the discrete forward model by

$$
\mathbf{H}: \mathbb{R}_{\ge 0}^N \times \Omega \times \mathbb{R}_{> 0} \rightarrow \mathbb{R}^M \\
(\mathbf{f}, \mathbf{x}_l, a_l) \mapsto \mathbf{B} |\mathbf{P}[\mathbf{A}(\mathbf{f}), \mathbf{I}_M] \mathbf{s}_{\text{in}}(\mathbf{x}_l, a_l)|^2 \quad (4)
$$

with

$$
\mathbf{A}(\mathbf{f}) = \tilde{\mathbf{G}} \operatorname{diag}(\mathbf{f}) (\mathbf{I}_N - \mathbf{G} \operatorname{diag}(\mathbf{f}))^{-1}, \quad (5)
$$

$$
\mathbf{s}_{\text{in}}(\mathbf{x}_l, a_l) = [(\mathbf{u}_{\text{in},l}^{\Omega})^T, (\mathbf{u}_{\text{in},l}^{\Gamma})^T]^T. \quad (6)
$$

Here, $\mathbf{f} \in \mathbb{R}^N$ is a sampled version of $f$ within $\Omega$. The vectors $\mathbf{u}_{\text{in},l}^{\Omega} \in \mathbb{C}^N$ and $\mathbf{u}_{\text{in},l}^{\Gamma} \in \mathbb{C}^M$ are the sampled versions of $u_{\text{in}}(\cdot; \mathbf{x}_l, a_l)$ within $\Omega$ and $\Gamma$, respectively. We denote by $\{\mathbf{x}_n^{\Omega}\}_{n=1}^N$ and $\{\mathbf{x}_m^{\Gamma}\}_{m=1}^M$ the sampling points within $\Omega$ and $\Gamma$. The matrix $\mathbf{G} \in \mathbb{C}^{N \times N}$ encodes the convolution with the Green function in (2). Similarly, $\tilde{\mathbf{G}} \in \mathbb{C}^{M \times N}$ is a matrix that, given the total field within $\Omega$, gives the scattered field at the measurement plane $\Gamma$. Next, $\mathbf{P} \in \mathbb{C}^{M \times M}$ is the discrete version of $P$ and $|\cdot|^2$ denotes the pointwise-squared magnitude. A full description of $\mathbf{G}$, $\tilde{\mathbf{G}}$, and $\mathbf{P}$ is provided in [19]. Finally, the matrix $\mathbf{B} \in \mathbb{R}^{M \times M}$ encodes a convolution with a Gaussian filter of length $\sigma_b = 0.7h$. It accounts for the mismatch between our physical model derived from the scalar diffraction theory and the vectorial nature of light [20, 21].

In this work, we adopt a biplane configuration [22] that
involves two pupil functions with separate focal planes. To
keep the notation simple, we shall use a single matrix **P**
to represent the effect of the two pupil functions (i.e., two
focal planes). Given the discrete forward model (4), the
image formation model (3) writes as, $\forall l \in \{1, \dots, L\}$,

$$
\mathbf{y}_l = \text{Pois}(\mathbf{H}(\mathbf{f}, \mathbf{x}_l, a_l) + \mathbf{b}_l). \quad (7)
$$

**Remark 1.** Although we consider a biplane modality in our experiments, the proposed joint optimization framework (Section 3) can be deployed with any number of focal planes. In this proof-of-concept work, we considered two

focal planes because i) it corresponds to a standard SMLM modality ii) it helps to compensate for the lack of phase measurements.

3. Joint Recovery of the Molecule Localization and Refractive Index

3.1. Joint-Optimization Framework

Our goal is to jointly recover the distribution of the re-
fractive index and the localization of fluorescent molecules.
To that end, we propose to solve the minimization problem

$$
(\mathbf{f}_*, \mathbf{X}_*, \mathbf{a}_*) \in \underset{\mathbf{x} \in \Omega^L, \mathbf{a} \in \mathbb{R}_{>0}^L}{\operatorname{arg\,min}} \sum_{l=1}^{L} \mathcal{J}_l(\mathbf{f}, \mathbf{x}_l, a_l) + \tau \mathcal{R}(\mathbf{f}), \quad (8)
$$

where, for all $\mathbf{f} \in \mathbb{R}_{\ge 0}^N$, $\mathbf{x} \in \Omega$, and $a > 0$,

$$
\mathcal{J}_l(\mathbf{f}, \mathbf{x}, a) = \mathcal{D}_{\text{KL}}(\mathbf{H}(\mathbf{f}, \mathbf{x}, a) + \mathbf{b}_l; \mathbf{y}_l). \quad (9)
$$

The matrix **X** = [$\mathbf{x}_1 \cdots \mathbf{x}_L$] ∈ Ω^L and the vector **a** = (a_1, ..., a_L) ∈ R_{>0}^L are the concatenation of positions and amplitudes of the fluorophores, respectively. The functional R: R^N → R_{≥0} is a regularization term that introduces prior knowledge on the RI distribution, and τ > 0 is a tradeoff parameter. In this work, we use the total-variation (TV) regularization [23], although alternatives such as the Hessian-Schatten norm [24] or learnt regularizers [25, 26, 27, 28] can be easily plugged into our framework. The data-fidelity term D_KL is the Kullback-Leibler divergence [29] defined as, ∀(z, y) ∈ R_{≥0}^M × R_{≥0}^M,

$$
\mathcal{D}_{\text{KL}}(\mathbf{z}; \mathbf{y}) = \mathbf{z}^T \mathbf{1}_M - \mathbf{y} \odot \log(\mathbf{z} + \beta), \quad (10)
$$

where $\beta > 0$ is a stabilizing parameter. Note that the
Kullback-Leibler divergence corresponds to the Poisson
negative log-likelihood up to some constant term.

To optimize (8), we alternate between an update of
the refractive index and an update of the amplitudes and
positions of the fluorophores (Algorithm 1), inspired by the
self-calibrating reconstruction techniques developed for oth-
er modalities [30, 31]. Updates are performed on a sub-
set of molecules (Line 3) in a stochastic fashion. In Al-
gorithm 1, aFBS ($\sum_{l \in \mathcal{L}} \mathcal{J}_l(\cdot, \mathbf{x}_l^{t+1}, a_l^{t+1}); \mathbf{f}^t$) refers to the
minimization of $\sum_{l \in \mathcal{L}} \mathcal{J}_l(\cdot, \mathbf{x}_l^{t+1}, a_l^{t+1})$ with the algorithm
aFBS initialized with $\mathbf{f}^t$. We use the same notation for the
refinement step at Line 5. Details on the algorithms de-
ployed for each sub-problem are provided in Sections 3.2
and 3.3. We implemented the method within the Global-
BioIm framework [32].

3.2. Update of Molecule Amplitudes and Positions

For the refinement procedure in Line 5 of Algorithm 1,
we again adopt an alternating scheme between an update
of the amplitude and the position, as summarized in Al-
gorithm 2. In the Sections 3.2.1 and 3.2.2, we describe
the Newton and gradient update steps used to refine the
amplitude and position, respectively.
**Algorithm 2** Refinement procedure for the lth molecule

Require: $\mathbf{x}_l^0 \in \Omega, a_l^0 > 0, T_{\text{mol}} \in \mathbb{N}$

1: $t = 0$

2: **while** (Not converged or $t < T_{\text{mol}}$) **do**

3: $a_l^{t+1} = \text{NewtonUpdate}(\mathcal{J}_l(\mathbf{f}, \mathbf{x}_l^t, \cdot); a_l^t)$

4: $\mathbf{x}_l^{t+1} = \text{GradientUpdate}(\mathcal{J}_l(\mathbf{f}, \cdot, a_l^{t+1}); \mathbf{x}_l^t)$

5: $t \leftarrow t + 1$

6: **end while**

Return: $\mathbf{x}_l^t, a_l^t$

### 3.2.1. Amplitudes

Let $\mathbf{f} \in \mathbb{R}_{\ge 0}^N$ and $\mathbf{X} \in \Omega^L$ be fixed. First of all, one can see from (4) that, for $a_l > 0$,

$$ \mathbf{H}(\mathbf{f}, \mathbf{x}_l, a_l) = a_l^2 \mathbf{H}(\mathbf{f}, \mathbf{x}_l, 1), \quad (11) $$

which is very helpful to reduce the computational cost of our joint-optimization procedure. Indeed, denoting $\mathbf{v}_l = \mathbf{H}(\mathbf{f}, \mathbf{x}_l, 1)$, we have that

$$ \begin{aligned} \mathcal{J}_l(\mathbf{f}, \mathbf{x}_l, a_l) &= \mathcal{D}_{\text{KL}}(a_l^2 \mathbf{v}_l + \mathbf{b}_l; \mathbf{y}_l) \\ &= (a_l^2 \mathbf{v}_l + \mathbf{b}_l)^T \mathbf{1}_M \\ &\quad - \mathbf{y}_l \odot \log(a_l^2 \mathbf{v}_l + \mathbf{b}_l + \beta). \end{aligned} \quad (13) $$

The function $\mathcal{J}_l$ is twice differentiable with respect to $a$. Its first derivative is given by

$$ \partial_a \mathcal{J}_l(\mathbf{f}, \mathbf{x}_l, a_l) = 2a_l \sum_{m=1}^{M} v_{lm} \left( 1 - \frac{y_{lm}}{a_l^2 v_{lm} + b_{lm} + \beta} \right). \quad (14) $$

Its second derivative reads as

$$ \begin{aligned} \partial_a^2 \mathcal{J}_l(\mathbf{f}, \mathbf{x}_l, a_l) = & 2 \sum_{m=1}^{M} v_{lm} \left( 1 - \frac{y_{lm}}{a_l^2 v_{lm} + b_{lm} + \beta} \right) \\ & + \sum_{m=1}^{M} \frac{(2a_l v_{lm})^2 y_{lm}}{(a_l^2 v_{lm} + b_{lm} + \beta)^2}. \end{aligned} \quad (15) $$

As such, we can perform a Newton update on $a_l$ as

$$ a_l^{t+1} = a_l^t - s \frac{\partial_a \mathcal{J}_l(\mathbf{f}, \mathbf{x}_l, a_l^t)}{\partial_a \mathcal{J}_l(\mathbf{f}, \mathbf{x}_l, a_l^t)}, \quad (16) $$

where $s$ is the length of a step computed via line-search so as to satisfy Wolfe's conditions [33].

### 3.2.2. Positions

Let $\mathbf{f} \in \mathbb{R}_{\ge 0}^N$ and $\mathbf{a} \in \mathbb{R}_{> 0}^L$ be fixed. We want to perform a gradient update on the position $\mathbf{x}_l$ of the lth molecule. However, one can see that the spherical wave in (1) is not differentiable whenever $\mathbf{x} = \mathbf{x}_l$. Consequently, we prefer to consider the smoothed version of the spherical wave

$$ u_{\text{in}}^{\text{smth}}(\mathbf{x}; \mathbf{x}_l, a_l) = a_l \frac{\exp(jk_b \|\mathbf{x} - \mathbf{x}_l\|_{2,\epsilon})}{4\pi \|\mathbf{x} - \mathbf{x}_l\|_{2,\epsilon}}, \quad (17) $$

**Algorithm 3** aFBS

Require: $\mathbf{f}^0 \in \mathbb{R}_{\ge 0}^N, T_{\text{RI}} \in \mathbb{N}, \gamma > 0, \alpha \in [0, 1]$

1: $t = 0, \mathbf{w}^0 = \mathbf{f}^0, v^0 = 1$

2: **while** (Not converged or $t < T_{\text{RI}}$) **do**

3: $\mathbf{g} = \sum_{l \in \mathcal{L}} \nabla \mathcal{J}_l(\cdot, \mathbf{x}_l, a_l)(\mathbf{w}^t)$

4: $\mathbf{f}^{t+1} = \text{prox}_{\gamma \tau \mathcal{R}}(\mathbf{w}^t - \gamma \mathbf{g})$

5: $v^{t+1} = \frac{1 + \sqrt{1 + 4(v^t)^2}}{2}$

6: $\mathbf{w}^{t+1} = \mathbf{f}^t + \alpha \frac{v^{t-1}}{v^{t+1}}(\mathbf{f}^t - \mathbf{f}^{t+1})$

7: $t \leftarrow t + 1$

8: **end while**

Return: $\mathbf{f}^t$

where $\|\cdot\|_{2,\epsilon} = \sqrt{\|\cdot\|_2^2 + \epsilon}$ with $0 < \epsilon \ll 1$. Then, the gradient of $\mathcal{J}_l$ with respect to $\mathbf{x}$, evaluated at $\mathbf{x}_l$, is given by

$$ \begin{aligned} \nabla_{\mathbf{x}} \mathcal{J}_l(\mathbf{f}, \mathbf{x}_l, a_l) = & 2 \mathbf{J}_{\text{sin},l}^H [\mathbf{A}(\mathbf{f}), \mathbf{I}_M]^H \mathbf{P}^H \mathbf{P} [\mathbf{A}(\mathbf{f}), \mathbf{I}_M] \mathbf{s}_{\text{in},l} \\ & \odot \mathbf{B}^T \nabla_{\mathbf{z}} \mathcal{D}_{\text{KL}} (\mathbf{H}(\mathbf{f}, \mathbf{x}_l, a_l) + \mathbf{b}_l), \end{aligned} \quad (18) $$

where $\mathbf{s}_{\text{in},l} = \mathbf{s}_{\text{in}}(\mathbf{x}_l, a_l) \in \mathbb{C}^{N+M}$. The gradient of $\mathcal{D}_{\text{KL}}$ in (10) with respect to the first variable $\mathbf{z}$ is given by

$$ \nabla_{\mathbf{z}} \mathcal{D}_{\text{KL}}(\mathbf{z}; \mathbf{y}) = \mathbf{1}_M - \mathbf{y}_l \otimes (\mathbf{z} + \beta). \quad (19) $$

Finally, it remains to provide the expression of the Hermitean transpose of the Jacobian matrix of $\mathbf{s}_{\text{in}}(\cdot, a_l)$, evaluated at $\mathbf{x}_l$, which we denote $\mathbf{J}_{\mathbf{s}_{\text{in},l}}^H \in \mathbb{C}^{3 \times (N+M)}$. Its qth column is given by

$$ [\mathbf{J}_{\mathbf{s}_{\text{in},l}}^H]_q = [\mathbf{s}_{\text{in},l}]_q^* \left( jk_b + \frac{1}{\|\mathbf{r}_q - \mathbf{x}_l\|_{2,\epsilon}} \right) \frac{(\mathbf{r}_q - \mathbf{x}_l)}{\|\mathbf{r}_q - \mathbf{x}_l\|_{2,\epsilon}}. \quad (20) $$

Let us emphasize that $\mathbf{r}_q = \mathbf{x}_q^\Omega$ ($\Omega$ sampling points) for $q \le N$ and $\mathbf{r}_q = \mathbf{x}_q^\Gamma$ ($\Gamma$ sampling points) for $N < q \le N + M$. Equipped with this closed-form gradient, we can deploy a projected-gradient update on $\mathbf{x}_l$ as

$$ \mathbf{x}_l^{t+1} = \mathcal{P}_\Omega (\mathbf{x}_l^t - s \nabla_{\mathbf{x}} \mathcal{J}_l(\mathbf{f}, \mathbf{x}_l^t, a_l)), \quad (21) $$

where $s$ is a step-size computed via a backtracking line-search [34]. The projector $\mathcal{P}_\Omega : \mathbb{R}^3 \to \Omega$ constrains the fluorophore positions to remain in $\Omega$.

### 3.3. Update of the Refractive Index

When the positions $\mathbf{X} \in \Omega^L$ and amplitudes $\mathbf{a} \in \mathbb{R}_{> 0}^L$ are fixed, the RI update consists in solving

$$ \mathbf{f}^* \in \arg \min_{\mathbf{f} \in \mathbb{R}_{\ge 0}^N} \sum_{l=1}^{L} \mathcal{J}_l(\mathbf{f}, \mathbf{x}_l, a_l) + \tau \mathcal{R}(\mathbf{f}). \quad (22) $$

It corresponds to an inverse-scattering problem from intensity measurements [35, 36]. To solve (22), we deploy a relaxed variant [37] of the accelerated forward-backward splitting (aFBS) algorithm [38, 39] (Algorithm 3). It requires the computation of two quantities.
1. The gradient of $J_l(\cdot, \mathbf{x}_l, a_l)$ which involves the Jacobian of $\mathbf{A}(\mathbf{f})$ in (4) whose expression is provided in [18].

2. The proximal operator of $\mathcal{R}$ which, for TV, can be efficiently evaluated by using the fast gradient-projection algorithm [40].

## 3.4. Initialization Strategies

### 3.4.1. Initialization of the Refractive Index

In phaseless diffraction tomography, the light-field refocusing method is a standard tool to obtain an initial guess of the RI distribution [35, 41]. However, this initialization requires coherent light sources with known geometry, which prevents its use on SMLM data. We therefore adopt an alternative approach that comprises two steps. We first replicate the widefield image (sum of the SMLM stack) along the axial direction and then blur the obtained volume with a Gaussian filter. The rationale behind this choice is that we can only expect to recover the RI where fluorophore emissions have propagated, that is, at the vicinity of fluorescent molecules. We then define $\mathbf{f}^0$ as a scaled version of this filtered volume so that its values belong to an admissible range of RI (see Figure 3).

### 3.4.2. Single-Molecule Localization

Any SMLM localization software can be used to compute the initial positions $\{\mathbf{x}_l^0\}_{l=1}^L$. However, we found that existing software packages for a biplane modality were not performing well on our simulated dataset. We believe that this is due to the high thickness of the sample together with the small number of acquisitions.

Therefore, we adopted a simple yet efficient method. We localize the position of the $l$th fluorophore based on cross-correlations between the measurements $\mathbf{y}_l$ and a set $\{\mathbf{k}_p\}_{p=1}^P$ of PSF models in $\in \mathbb{R}^M$. We define them as the output of the forward model with no scatterer, like in

$$ \mathbf{k}_p = \mathbf{H}(\mathbf{0}_N, \mathbf{x}_p^{\text{psf}}, 1), \quad (23) $$

where the positions $\mathbf{x}_p^{\text{psf}} = (0, 0, p\Delta z)$ for $p = \{-P, \dots, P\}$ vary along the axial direction. We then initialize

$$ \mathbf{x}_l^0 = (\hat{m}h, \hat{n}h, \hat{p}\Delta z) \quad (24) $$

where

$$ (\hat{m}, \hat{n}, \hat{p}) = \arg \max_{m,n,p} [\mathbf{y}_l \star \mathbf{k}_p^{\vee}]_{m,n}, \quad (25) $$

with $\star$ the 2D discrete convolution and $(\cdot)^\vee$ the reflection operator. Once localized, we initialize the amplitude as

$$ a_l^0 = \left( \|\mathbf{M}_{\mathbf{x}_l^0}(\mathbf{y}_l - \hat{\mathbf{b}}_l)\|_1 / \|\mathbf{k}_p\|_1 \right)^{\frac{1}{2}}, \quad (26) $$

where $\mathbf{M}_{\mathbf{x}_l^0} \in \mathbb{R}^{R \times M}$ crops a region-of-interest centered at $\mathbf{x}_l^0$ and $\hat{\mathbf{b}}_l$ denotes the estimated background (see Section 3.4.3).

### 3.4.3. Background Estimation

To estimate $\{\mathbf{b}_l\}_{l=1}^L$, we apply a simple algorithm suitable for a background that slowly varies in space and time. In SMLM, this is a common assumption [42]. Our procedure proceeds in two steps.

1. For each measurement $\mathbf{y}_l$, we mask an area around the estimated position $\mathbf{x}_l$ and inpaint it using the function `regionfill` of Matlab<fnref n="1" /> to obtain $\bar{\mathbf{y}}_l$.

2. We apply a spatio-temporal (3D) median filter along the stack of masked and inpainted measurements $\bar{\mathbf{Y}} = [\bar{\mathbf{y}}_1 \cdots \bar{\mathbf{y}}_L]$ to take advantage of the spatio-temporal smoothness of the background.

## 4. Numerical Experiments

### 4.1. Simulation Setting

We created an RI volume immersed in water ($\eta_b = 1.339$), fully included in the region $\Omega$ of size $(7.2 \times 7.2 \times 3.2)\mu\text{m}^3$ (Figure 2). This sample presents small features with RI values that are lower or higher than their surroundings. Then, we populated this sample with fluorophores randomly placed on a structure that is composed of an outer membrane as well as inner compartments. The smallest distance between two fluorophores is 20nm. We simulated $L = 1000$ SMLM acquisitions with a biplane modality, each corresponding to the activation of a single fluorophore. The two focal planes were set at $\pm 300\text{nm}$. The amplitude $a_l$ of each fluorophore emission was drawn from a Poisson distribution with mean $A = 1000$ and the wavelength of the emitted light is set at $\lambda = 647\text{nm}$. In addition, we simulated a pupil function for each focal plane with NA = 1.45 and 25 Zernike coefficients. Their values were drawn from the uniform distribution $U(-0.5, 0.5)$, except that the three first coefficients were set to 0 and that the fourth coefficient was drawn from $U(-0.1, 0.1)$ to better match the PSFs observed in real SMLM acquisitions. The background signals $\mathbf{b}_l$ for $l \in \{1, 100, 200, \dots, 1000\}$ were simulated by convolving a Gaussian kernel with a random image generated from a uniform distribution. We then scaled the obtained images so that their pixel values belong to the range [350, 450]. Backgrounds for intermediate frames were then obtained through interpolation. We set a large width for the Gaussian kernel so as to obtain a slowly varying background in both space and time. Finally, to control the noise, we scaled the noiseless measurements with a factor $r \in (0, 1]$ before applying the Poisson noise so that (4) writes as, $\forall l \in \{1, \dots, L\}$,

$$ \mathbf{y}_l = \text{Pois}(r(\mathbf{H}(\mathbf{f}, \mathbf{x}_l, a_l) + \mathbf{b}_l)). \quad (27) $$

By doing so, $r$ can be interpreted as the product between the excitation photon flux and the integration time. A

<footnote n="1">¹Matlab's command `regionfill` performs a smooth interpolation inward from the pixel values that are on the outer boundary of the mask.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: Simulation setup. The RI map, immersed in water ($n_b = 1.339$), is depicted on the left orthoviews. The sample is populated with fluorescent molecules that belong to the labeled region. They sequentially emit a spherical wave which is then propagated through the sample using the Lippmann-Schwinger model. Two focal planes (with pupil functions) are acquired. The widefield images are generated by summing all SMLM frames. The fluorescence images were saturated for visualization purpose. We display the labeled region with partial transparency so as to make the inner compartments visible. Scale bars: 500nm.</figcaption></figure>

small $r$ yields a higher level of noise, which increases the difficulty of the localization of molecules and the RI reconstruction.

We compare our joint-optimization framework with two baselines. They consist on the sole RI reconstruction with i) perfectly characterized molecules (i.e., true amplitudes and positions) or ii) the initial estimation of the amplitudes and positions obtained as described in Section 3.4.2. By doing so, we somehow obtain the worst-case and best-case scenarios. For each case, we obtained the best reconstruction by performing a grid search on the regularization parameter $\tau$.

For our joint-optimization framework, we set the parameters $T_{mol} = 4$, $T_{RI} = 1$, and $\alpha = 0.85$. In our implementation, $\mathbf{w}^0$ and $\mathbf{v}^0$ in Line 1 of Algorithm 3 are initialized from the previous call. We ran our optimization on a PowerEdge c4140 equipped with Intel Xeon Gold 6240 CPUs (2.60GHz) and a GPU NVIDIA Tesla V100 SXM3 (32 GB). An iteration of Algorithm 1 took 20 seconds on average. We used up to 2000 iterations, which corresponds to about 10 hours of computation.

## 4.2. Metrics and Visualization

To assess the quality of the reconstructed RI volume, we compute the relative error as well as the structural similarity index measure (SSIM) [43] with respect to the ground-truth. To assess the accuracy of the localization of the molecules, we compute the root-mean-square error (RMSE) with respect to the true positions. Note that we do not report detection metrics such as true/false detections as they are not really relevant in our setting where we consider only frames containing one molecule. Finally, given a list of molecule positions, we generate a 3D image through the Gaussian rendering technique [6]. To that

end, we represent the fluorophore positions as a sum of shifted Dirac

$$s(\mathbf{x}) = \sum_{l=1}^{L} \delta(\mathbf{x} - \mathbf{x}_l). \qquad (28)$$

Gaussian rendering then consists in convolving $s$ with an isotropic Gaussian kernel and sample the result on a grid. Here, we set the standard deviation of the Gaussian kernel to 10nm and the grid step to $h/10 = 10$nm.

## 4.3. Results

We first fix the noise level to $r = 1$ in (27).

### 4.3.1. Reconstructed Refractive Index

We display the RI volumes in Figure 3 and report there the relative errors and SSIM. When the positions and amplitudes of the molecules are perfectly known, we recover most of the details of the ground-truth. This is in line with our previous work [13]. On the contrary, the reconstruction obtained with the initial positions and amplitudes is unsuccessful. This highlights the importance of refining molecule positions and amplitudes jointly with RI reconstruction. We effectively see that our joint-optimization framework is able to recover an RI volume that is visually similar to the best-case scenario. The metrics confirm the visual assessment. Yet, one can observe some high frequency artifacts (ringing) on the reconstruction obtained with the joint optimization framework (plane $z = -0.2\mu m$). They are due to few badly refined molecule positions (outliers in Figure 8) that lead to a mismatch in the model. Finally, it should be noted that we could expect that the quality of the reconstruction varies with the axial position $z$. The reason is that an SMLM
<figure><img src="image_4.png" /><figcaption>Figure 3: **Reconstructions of the RI volume.** From left to right: Ground truth, initial guess, reconstruction with positions and amplitudes fixed to their initial values (see Section 3.4.2), reconstruction with the proposed joint-optimization framework, and reconstruction with positions and amplitudes fixed to their true values (gold-standard). The SSIM and relative errors are displayed in the first row at the left and right columns, respectively, of each corresponding reconstruction. Scale bar: 500nm.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: **Observable region of the sample.** Illustration of the region of the sample that is “illuminated” by a fluorophore far from (left) and close to (middle) the detection system. The right scheme illustrates the fact that SMLM data carry more information about RI regions that are close to the detection system (positive *z*).</figcaption></figure>

frame (from the activation of one molecule) carries infor-
mation about the part of sample that lies between the
activated molecule and the optical system (see Figure 4).
As such, there are more SMLM frames that carry informa-
tion about *z*-planes with positive *z* than frames that carry

information about *z*-planes with negative *z*. Moreover,
waves produced by fluorophores with negative *z*-positions
propagate through a larger layer of the sample, inducing
more scattering. These facts make that i) fluorophores
with negative *z*-positions are harder to localize, and ii)
*z*-planes of RI with negative *z* are harder to reconstruct.

4.3.2. Molecule Localization

It is noteworthy to recall that the primary objective of SMLM is to localize the fluorescent molecules with nano-metric precision. It follows that another benefit of our joint-optimization framework is an improvement of this localization. Indeed, our model accounts for sample-induced distortions that usually compromises the accurate localization of molecules [9].

We report in Table 1 the RMSE of the initial and re-
fined positions, as well as the RMSE of the initial and re-
fined amplitudes. In addition, we provide the RMSE of the
refined positions and amplitudes when the RI map is fixed
to the initial guess **f**^0 or the ground-truth (best-case sce-
<table><caption>Table 1: RMSE for the estimated positions and amplitudes. First row: Initial positions and amplitudes from our standard single-molecule localization. Second row: Positions and amplitudes from our joint-optimization framework. Third and fourth rows: Positions and amplitudes refined with the RI map fixed to the initial guess $\mathbf{f}^0$ and the ground-truth (GT), respectively. Amp.: Amplitude. Lat.: Lateral. Ax.: Axial.</caption><thead><tr><th></th><th>3D [nm]</th><th>Lat. [nm]</th><th>Ax. [nm]</th><th>Amp.</th></tr></thead><tbody><tr><td>Initial</td><td>163</td><td>69</td><td>148</td><td>109</td></tr><tr><td>Joint</td><td>74</td><td>15</td><td>72</td><td>76</td></tr><tr><td>with f<sup>0</sup></td><td>142</td><td>38</td><td>136</td><td>194</td></tr><tr><td>with f<sub>GT</sub></td><td>72</td><td>18</td><td>70</td><td>76</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 5: **Rendering of localized molecules (Y-projection).** Region-of-interest of the projection along Y of the rendered fluorescent volume. From top to bottom: Initial positions, positions refined with the joint-optimization framework, ground-truth. Field-of-view (XZ): (3600 × 600)nm². The images were saturated for visualization purpose. Scale bar: 200nm.</figcaption></figure>

nario). There is a gain of 89nm in the 3D RMSE for our joint-optimization framework. One sees that the lateral and axial RMSE are improved by 54nm and 76nm, respectively. Not only does our joint-optimization framework successfully recover the RI volume, but it also improves significantly the localization of the molecules. The proposed joint-optimization framework performs better than the refinement of the positions and amplitudes with the RI map fixed to $\mathbf{f}^0$. Moreover, it performs similarly to the refinement with the RI map fixed to the ground-truth. Those observations confirm that the joint-optimization framework is necessary to improve the localization and can even reach similar performance to the best-case scenario.

In Figure 5, we display a Y-projection of the fluorescent volume rendered from the molecule positions, where one can visually appreciate the gain in accuracy. The estimation of the amplitudes is improved as well. This can help to better estimate the uncertainty of localization [44, 45].

### 4.3.3. Influence of the Distribution of Fluorophores

From the phenomenon illustrated in Figure 4, one can expect that the quality of the reconstructed RI map is closely related to the spatial distribution of the fluorescent probes. In this section, we investigate this question

#### A: Distribution of fluorophores

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

#### B: RI reconstruction

<figure><img src="image_8.png" /><figcaption>Figure 6: **Reconstructions of the RI volume with different fluorophore distributions.** A: The four considered fluorophore distributions. B: Reconstructed RI maps. The SSIM and relative errors are displayed in the XZ view at the left and right corners of each reconstruction, respectively. The signal-to-noise ratio is displayed at the bottom-left corner for each noise level. Scale bar: 500nm.</figcaption></figure>

by comparing the reconstructions obtained with the four fluorophore distributions illustrated in Figure 6 (Panel A). These include the rather homogeneous distribution depicted in Figure 2 and a more concentrated distribution, both with two different numbers of molecules (i.e., numbers of frames).
<figure><img src="image_9.png" /><figcaption>Figure 7: **Reconstructions of the RI volume with different noise levels.** From left to right: Ground truth, reconstruction with the proposed joint optimization framework for $r = 0.1, 0.5$, and $1$. The SSIM and relative errors are displayed in the first row at the left and right corners of each reconstruction, respectively. The last row contains two examples of SMLM acquisition (ROI) for two different molecules at axial positions $z = 440$nm (top) and $z = 0$nm (bottom). The signal-to-noise ratio is displayed at the bottom-left corner for each noise level. Scale bar: 500nm.</figcaption></figure>

<table><caption>Table 2: RMSE of the estimated positions and amplitudes for the four distributions of fluorophores depicted in Figure 6 (Panel A). Amp.: Amplitude. Lat.: Lateral. Ax.: Axial. Dist.: Distribution.</caption><thead><tr><th></th><th>#Fluo</th><th>3D [nm]</th><th>Lat. [nm]</th><th>Ax. [nm]</th><th>Amp.</th></tr></thead><tbody><tr><td rowspan="2">Dist. 1</td><td>1000</td><td>74</td><td>15</td><td>72</td><td>76</td></tr><tr><td>100</td><td>80</td><td>13</td><td>79</td><td>76</td></tr><tr><td rowspan="2">Dist. 2</td><td>1000</td><td>85</td><td>3</td><td>85</td><td>77</td></tr><tr><td>100</td><td>82</td><td>3</td><td>82</td><td>55</td></tr></tbody></table>

As expected, the reconstructed RI map is significantly degraded when the distribution of fluorophores is more concentrated (Figure 6B, right column). Indeed, the emitted light has mainly propagated through a restricted area of the sample, limiting the information on the RI map carried by the measurements. On the contrary, the qual-

ity of the reconstructed RI map seems less sensitive to the number of fluorophores. Although some details are lost, the RI maps reconstructed with 100 frames remain qualitatively similar to their counterparts reconstructed from 1000 frames (Figure 6B).

Finally, we display in Table 2 the RMSE obtained after the joint optimization. Interestingly, the refinement of the positions and amplitudes of the molecules remains stable when reducing the number of frames. For the concentrated distribution, the axial RMSE is slightly degraded and, on the contrary, the lateral RMSE is drastically reduced, which might be due to the concentration of the distribution.

### 4.3.4. Robustness to Noise

Next, we are interested in the robustness of our framework to the measurement noise. To that end we vary the parameter $r$ in (27) from 0.1 to 1. Some examples of ob-
<figure><img src="image_10.png" /><figcaption>Figure 8: Box plots of the localization (left) and amplitude (right) error for different noise levels. Three noise levels are displayed with *r* = 0.1, 0.5, and 1. For each case, the left box plot (hatched) corresponds to the initial positions/amplitudes and the right box plot (solid) corresponds to the refined positions/amplitudes. Note that the estimated amplitudes were scaled by *r*⁻¹ to compare with the same ground-truth. For the box plots of the localization error after refinement (solid), we set the upper whiskers to 50nm so as to consider any larger error as outliers. This is in line with the expected 3D localization error in SMLM [6]. This bound is not relevant for the initial errors (hatched) as they are too large. For the hatched box plots, we thus set the default upper whiskers to *Q*3 + 1.5IQR, where IQR = ( *Q*3 − *Q*1) is the interquartile range and *Q*1, *Q*3 are the 25th and 75th percentile, respectively. Finally, the lower whiskers are always set to the smallest error among all molecules. Outliers are indicated by ×.</figcaption></figure>

tained measurements are shown in Figure 7 (last row),
where one can observe that the noise is stronger when *r* is
smaller. The RI reconstructions for each noise level are
displayed in Figure 7. Although the quality of reconstruc-
tion degrades when the noise increases, the shape and the
most prominent features are recognizable even for *r* = 0.1.
This suggests that our method is quite robust to noise.

The box plots of the localization errors are displayed in
Figure 8. For each noise level, we show the box plot for the
initial and refined positions to illustrate the improvement.
We again observe a certain robustness to noise, even for the
case r = 0.1 where the amplitudes were badly initialized.
It is noteworthy to mention the presence of outliers in the
displayed box plots even for r = 1. For some molecules, we
observed that the joint-optimization could not refine the
positions and amplitudes well. In few cases, the estimates
did even worsen. Fortunately, the number of such failures
is limited (e.g., 60 over 1000 molecules for r = 0.1).

5. Discussion

We have presented a joint-optimization framework to
estimate both the RI map and the position of fluores-
cent molecules from an SMLM acquisition stack. Our
method takes advantage of the sample-induced aberra-
tions to unveil the map of the refractive index of the sam-
ple. Such structural information complements fluorescence
imaging [10]. In addition to this unique feature, our frame-
work is able to improve the accuracy of molecule localiza-
tion. Our work shows that additional information about
the sample can be recovered from SMLM data. This is a

first step towards an exciting and new extension of SMLM.
As future refinement, we need to look into ways to identify
and correct the molecules that our framework can some-
times fail to refine. We shall then be in a good position to
apply our framework to real data.

Acknowledgments

T.-a.P. and M.U. were supported by the European Re-
search Council (ERC) under the European Unions Horizon
2020 research and innovation programme (Grant Agree-
ment No. 692726 GlobalBioIm). E.S. was supported by
the LabEx CIMI (ANR-11-LABX-0040).

References

[1] E. Betzig, G. H. Patterson, R. Sougrat, O. W. Lindwasser, S. Olenych, J. S. Bonifacino, M. W. Davidson, J. Lippincott-Schwartz, H. F. Hess, Imaging intracellular fluorescent proteins at nanometer resolution, Science 313 (2006) 1642–1645.

[2] S. T. Hess, T. P. Girirajan, M. D. Mason, Ultra-high reso-
lution imaging by fluorescence photoactivation localization mi-
croscopy, Biophysical Journal 91 (2006) 4258–4272.

[3] M. J. Rust, M. Bates, X. Zhuang, Sub-diffraction-limit imag-
ing by stochastic optical reconstruction microscopy (STORM),
Nature methods 3 (2006) 793–796.

[4] B. Huang, W. Wang, M. Bates, X. Zhuang, Three-dimensional super-resolution imaging by stochastic optical reconstruction microscopy, Science 319 (2008) 810–813.

[5] S. R. P. Pavani, M. A. Thompson, J. S. Biteen, S. J. Lord, N. Liu, R. J. Twieg, R. Piestun, W. Moerner, Three-dimensional, single-molecule fluorescence imaging beyond the diffraction limit by using a double-helix point spread function, Proceedings of the National Academy of Sciences 106 (2009) 2995–2999.
[6] D. Sage, T.-a. Pham, H. Babcock, T. Lukes, T. Pengo, J. Chao, R. Velmurugan, A. Herbert, A. Agrawal, S. Colabrese, A. Wheeler, A. Archetti, B. Rieger, R. Ober, G. M. Hagen, J.-B. Sibarita, J. Ries, R. Henriques, M. Unser, S. Holden, Super-resolution fight club: Assessment of 2D and 3D single-molecule localization microscopy software, Nature Methods 16 (2019) 387-395.

[7] H. P. Babcock, X. Zhuang, Analyzing single molecule localization microscopy data using cubic splines, Scientific Reports 7 (2017) 1-9.

[8] Y. Li, M. Mund, P. Hoess, J. Deschamps, U. Matti, B. Nijmeijer, V. J. Sabinina, J. Ellenberg, I. Schoen, J. Ries, Real-time 3D single-molecule localization using experimental point spread functions, Nature Methods 15 (2018) 367-369.

[9] F. Xu, D. Ma, K. P. MacPherson, S. Liu, Y. Bu, Y. Wang, Y. Tang, C. Bi, T. Kwok, A. A. Chubykin, P. Yin, S. Calve, G. E. Landreth, F. Huang, Three-dimensional nanoscopy of whole cells and tissues with in situ point spread function retrieval, Nature Methods 17 (2020) 531-540.

[10] Y. Park, C. Depeursinge, G. Popescu, Quantitative phase imaging in biomedicine, Nature Photonics 12 (2018) 578-589.

[11] Y. Xue, L. Waller, Computational hybrid imaging, arXiv preprint arXiv:2011.11574 (2020).

[12] A. C. Gilbert, H. W. Levinson, J. C. Schotland, Imaging from the inside out: Inverse scattering with photoactivated internal sources, Optics Letters 43 (2018) 3005-3008.

[13] T.-a. Pham, E. Soubies, F. Soulez, M. Unser, Diffraction tomography from single-molecule localization microscopy: Numerical feasibility, IEEE International Symposium on Biomedical Imaging (April 13-16, 2021).

[14] E. Wolf, Three-dimensional structure determination of semi-transparent objects from holographic data, Optics Communications 1 (1969) 153-156.

[15] J. W. Goodman, Introduction to Fourier Optics, Roberts and Company Publishers, 2005.

[16] A. Cornea, P. M. Conn, Fluorescence Microscopy: Super-Resolution and Other Novel Techniques, Elsevier, 2014.

[17] H.-Y. Liu, D. Liu, H. Mansour, P. T. Boufounos, L. Waller, U. S. Kamilov, Seagle: Sparsity-driven image reconstruction under multiple scattering, IEEE Transactions on Computational Imaging 4 (2017) 73-86.

[18] E. Soubies, T.-a. Pham, M. Unser, Efficient inversion of multiple-scattering model for optical diffraction tomography, Optics Express 25 (2017) 21786-21800.

[19] T.-a. Pham, E. Soubies, A. Ayoub, J. Lim, D. Psaltis, M. Unser, Three-dimensional optical diffraction tomography with Lippmann-Schwinger model, IEEE Transactions on Computational Imaging 6 (2020) 727-738.

[20] A. Aristov, B. Lelandais, E. Rensen, C. Zimmer, ZOLA-3D allows flexible 3D localization microscopy over an adjustable axial range, Nature Communications 9 (2018) 1-8.

[21] S. Liu, E. B. Kromann, W. D. Krueger, J. Bewersdorf, K. A. Lidke, Three dimensional single molecule localization using a phase retrieved pupil function, Optics Express 21 (2013) 29462-29487.

[22] M. F. Juette, T. J. Gould, M. D. Lessard, M. J. Mlodzianoski, B. S. Nagpure, B. T. Bennett, S. T. Hess, J. Bewersdorf, Three-dimensional sub-100 nm resolution fluorescence microscopy of thick samples, Nature Methods 5 (2008) 527-529.

[23] L. I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear Phenomena 60 (1992) 259-268.

[24] S. Lefkimiatis, J. P. Ward, M. Unser, Hessian Schatten-norm regularization for linear inverse problems, IEEE Transactions on Image Processing 22 (2013) 1873-1888.

[25] G. Ongie, A. Jalal, C. A. Metzler, R. G. Baraniuk, A. G. Dimakis, R. Willett, Deep learning techniques for inverse problems in imaging, IEEE Journal on Selected Areas in Information Theory 1 (2020) 39-56.

[26] F. Yang, T.-a. Pham, H. Gupta, M. Unser, J. Ma, Deep-learning projector for optical diffraction tomography, Optics Express 28 (2020) 3905-3921.

[27] H. Gupta, K. H. Jin, H. Q. Nguyen, M. T. McCann, M. Unser, CNN-based projected gradient descent for consistent ct image reconstruction, IEEE Transactions on Medical Imaging 37 (2018) 1440-1453.

[28] A. Goy, G. Rughoobur, S. Li, K. Arthur, A. I. Akinwande, G. Barbastathis, High-resolution limited-angle phase tomography of dense layered objects using deep neural networks, Proceedings of the National Academy of Sciences 116 (2019) 19848-19856.

[29] S. Kullback, R. A. Leibler, On information and sufficiency, The Annals of Mathematical Statistics 22 (1951) 79-86.

[30] M. Zehni, L. Donati, E. Soubies, Z. Zhao, M. Unser, Joint angular refinement and reconstruction for single-particle cryoem, IEEE Transactions on Image Processing 29 (2020) 6151-6163.

[31] M. Xie, Y. Sun, J. Liu, B. Wohlberg, U. S. Kamilov, Joint reconstruction and calibration using regularization by denoising, arXiv preprint arXiv:2011.13391 (2020).

[32] E. Soubies, F. Soulez, M. T. McCann, T.-a. Pham, L. Donati, T. Debarre, D. Sage, M. Unser, Pocket guide to solve inverse problems with GlobalBioIm, Inverse Problems 35 (2019) 1-20.

[33] P. Wolfe, Convergence conditions for ascent methods, SIAM Review 11 (1969) 226-235.

[34] L. Armijo, Minimization of functions having Lipschitz continuous first partial derivatives, Pacific Journal of Mathematics 16 (1966) 1-3.

[35] L. Tian, L. Waller, 3D intensity and phase imaging from light field measurements in an LED array microscope, Optica 2 (2015) 104-111.

[36] T.-a. Pham, E. Soubies, A. Goy, J. Lim, F. Soulez, D. Psaltis, M. Unser, Versatile reconstruction framework for diffraction tomography with intensity measurements and multiple scattering, Optics Express 26 (2018) 2749-2763.

[37] Y. Ma, H. Mansour, D. Liu, P. T. Boufounos, U. S. Kamilov, Accelerated image reconstruction for nonlinear diffractive imaging, in: 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Calgary AB, Canada, April 15-20, 2018, pp. 6473-6477.

[38] Y. Nesterov, Gradient methods for minimizing composite functions, Mathematical Programming 140 (2013) 125-161.

[39] A. Beck, M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM Journal on Imaging Sciences 2 (2009) 183-202.

[40] A. Beck, M. Teboulle, Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems, IEEE Transactions on Image Processing 18 (2009) 2419-2434.

[41] L. Tian, X. Li, K. Ramchandran, L. Waller, Multiplexed coded illumination for Fourier ptychography with an LED array microscope, Biomedical Optics Express 5 (2014) 2376-2389.

[42] J. Min, S. J. Holden, L. Carlini, M. Unser, S. Manley, J. C. Ye, 3D high-density localization microscopy using hybrid astigmatic/biplane imaging and sparse image reconstruction, Biomedical Optics Express 5 (2014) 3935-3948.

[43] Z. Wang, A. C. Bovik, H. R. Sheikh, E. P. Simoncelli, Image quality assessment: From error visibility to structural similarity, IEEE Transactions on Image Processing 13 (2004) 600-612.

[44] R. E. Thompson, D. R. Larson, W. W. Webb, Precise nanometer localization analysis for individual fluorescent probes, Biophysical Journal 82 (2002) 2775-2783.

[45] J. Chao, E. S. Ward, R. J. Ober, Fisher information theory for parameter estimation in single molecule microscopy: Tutorial, JOSA A 33 (2016) B36-B57.