# Analysis of Super-resolution Single Molecule Localization Microscopy Data: a
tutorial

Mohamadreza Fazel1,2 and Michael J. Wester1,3
1Department of Physics and Astronomy, University of New Mexico,
Albuquerque, New Mexico 87106, USA
2Center for Biological Physics, Department of Physics, Arizona State
University,
Tempe, Arizona 85287, USA
3Department of Mathematics and Statistics, University of New Mexico,
Albuquerque, New Mexico 87106, USA

###### Abstract

The diffraction of light imposes a fundamental limit on the resolution of
light microscopes. This limit can be circumvented by creating and exploiting
independent behaviors of the sample at length scales below the diffraction
limit. In super-resolution single molecule localization microscopy (SMLM), the
independence arises from individual fluorescent labels stochastically
switching between dark and fluorescent states, which in turn allows the
pinpointing of fluorophores post experimentally using a sequence of acquired
sparse image frames. Finally, the resulting list of fluorophore coordinates is
utilized to produce high resolution images or to gain quantitative insight
into the underlying biological structures. Therefore, image processing and
post-processing are essential stages of SMLM. Here, we review the latest
progress on SMLM data processing and post-processing.

Author to whom correspondence should be addressed<EMAIL_ADDRESS>
##  1 Introduction

In the beginning of the 20th century, a major difficulty in the examination of
living organisms was their transparency, which allows transmission of a large
portion of light without effective scattering or attenuation [1]. Fluorescence
microscopy techniques were therefore developed which overcome this problem.
These techniques improve microscopy images by relying on light from
fluorescent probes that stain target structures within a sample instead of
weak scattered light from the entire sample. In spite of this accomplishment,
fluorescence microscopy still suffers from major issues: 1) light diffraction;
2) inefficient collection of fluorescent light by microscope lenses; 3) the
contribution of light from parts of the specimen that are out of the focal
point/plane; 4) photodamage to the sample due to high excitation laser
intensities or long data acquisition times.

Here, we discuss the fundamental resolution limit in microscopy images which
arises due to issues 1-2. In fluorescence microscopy under ideal conditions,
optimal resolution could be obtained if every fluorophore (point source of
light) was represented by a point in microscopy images. However, the
diffraction of light emitted by a fluorophore and inefficiently collected by a
microscope lens (objective) results in an extended intensity pattern called
the Airy disk. Furthermore, the actual pattern often differs from the Airy
disk due to optical aberrations, and is referred to as a point spread function
(PSF) or an impulse response.

The PSF plays a central role in determining the resolving power of
microscopes. For instance, two-point resolution is defined as the minimum
distance, $d$, between two point-like emitters in which their overlapping PSFs
can still be distinguished as two individual PSFs by the microscope user [2].
Rayleigh’s criterion for the two-point resolution expresses that two
noncoherent sources of light can be resolved when they are at least separated
by a distance equal to the separation between the maximum and the first
minimum of the diffraction pattern of one of the sources, Fig. 1a-b. Assuming
an Airy pattern, the Rayleigh’s resolution distance for point sources in a
lateral direction is given by [2, 3]

$d=0.61\frac{\lambda}{N_{a}}$ (1)

where $\lambda$ is the wavelength of the light. $N_{a}$ is the numerical
aperture, which quantifies the efficiency of an objective lens in collecting
light, Fig. 2.

Figure 1: Fundamental limit of resolution. (a) 1D diffraction pattern (Airy
pattern) from a circular aperture (lens). (b) The red curves are two Airy
patterns where their maxima are separated by $d$. The blue curve is the
intensity pattern from two noncoherent point sources of light with separation
$d$. (c) The absolute value of the Fourier transform of the Airy pattern in
panel a. $k_{\mathrm{max}}$ denotes the spatial frequency cutoff, which is the
largest spatial frequency passed by the lens.

On the other hand, the Abbe criterion provides a definition for resolution
based on the notion of the optical transfer function (OTF), which is simply
the Fourier transform of the PSF. The absolute value of the OTF is called the
modulation transfer function (MTF), which gradually decays with increasing
spatial frequency, going to zero at a maximum frequency, $k_{\mathrm{max}}$,
Fig. 1c. This frequency cutoff puts a fundamental limit on the ability of
microscopes (or any imaging system) to collect information, i.e., spatial
frequencies, from the sample [4]. This effect is called the Abbe diffraction
limit and results in the following lateral resolution [2]

$d=\frac{\lambda}{2N_{a}}$ (2)

which is related to the inverse of $k_{\mathrm{max}}$. For visible light,
$\lambda\sim 550$ nm, and for a typical objective, $N_{a}\sim 1.3$, so the
resolution limit is approximately $200$ nm.

We now proceed with a short description of resolution along the axial ($z$)
direction. Similar to the lateral resolution, the Abbe criterion can be used
to find the axial resolution. However, the larger size of the PSF in the axial
direction yields a smaller frequency cutoff compared to the lateral case. As a
result, the inverse of this cutoff leads to a larger (worse) resolution
distance, which can be approximated by [2, 5]

$d\sim\frac{2n\lambda}{(N_{a})^{2}},$ (3)

where $n$ is the index of refraction.

Figure 2: Numerical aperture. The numerical aperture ($N_{a}$) is related to
the fraction of light that can be collected by a lens. Here, the lens captures
light rays traveling within a cone with half angle $\theta$ given by
$\sin\theta\sim D/2f$ for an emitter sitting at the focal point. The numerical
aperture is defined by $N_{a}=n\sin\theta$, where the refractive index $n$
guarantees independence of $N_{a}$ from the medium by Snell’s law.

Now that we have introduced the fundamental limit of resolution as a
consequence of issues 1-2, we turn to issues 3-4. Although issues 3-4 do not
yield any fundamental limits, they often result in degradation in the quality
of microscope images. For example, undesired light from out-of-focus
structures stained by fluorophores, or autofluorescence from these structures,
contribute to the fluorescence light captured by the objective lens. Moreover,
some fluorescence microscopy techniques require high laser intensities or long
exposure times that might damage the sample.

To alleviate the issues associated with fluorescence microscopy, different
illumination and light collection techniques [6, 7, 8, 9, 10, 11] were
developed to reduce out-of-focus light and collect fluorescence emission more
efficiently. These techniques can be divided into two main categories: point
scanning and wide-field microscopes. The point scanning microscopes, such as
confocal [6] and 4-pi [9] microscopes, scan the specimen by illuminating a
single spot at a time. The main benefit of these techniques is the reduction
of out-of-focus light by embedding pinholes in the optical setup that only
allow transmission of the in-focus light. In addition to pinholes, the 4-pi
technique also uses two objective lenses to collect fluorescent light,
yielding a larger numerical aperture which, in turn, results in a higher
resolution. In wide-field microscopes, such as the total internal reflection
fluorescence (TIRF) microscope [7], the entire specimen is exposed to an
exciting laser, as opposed to scanning the sample point by point, allowing for
faster data acquisition.

Different fluorescence microscopy techniques pushed the Abbe diffraction
barrier to its very limits, however, it was not until the end of 20th century
that scientists were able to overcome this barrier and achieve resolutions
better than the diffraction limit [12]. Such techniques that surpass the
diffraction limit and achieve sub-diffraction resolutions (resolutions not
restricted by the cutoff in the frequency transmission, $k_{\mathrm{max}}$, by
an objective lens) are called super-resolution microscopy or nanoscopy [13,
14, 15, 16, 17, 18, 19]. These techniques accomplish sub-diffraction
resolution by creating independent behavior for individual fluorophores at
scales below the diffraction limit. Many super-resolution procedures use the
reversible switching property of some fluorescent probes transitioning between
a fluorescent state (emitting state) and a dark state (non-emitting state) to
obtain sub-diffraction resolution. These approaches can be classified into two
different groups based on how the probes are switched between the dark and
fluorescent states: targeted and stochastic switching procedures. The
microscopy techniques under the first category are STimulated Emission
Depletion (STED) microscopy [12, 20], Ground State Depletion (GSD) microscopy
[21], REversible Saturable Optically Linear Fluorescence Transition (RESOLFT)
microscopy [22] and Saturated Structured Illumination Microscopy (SSIM) [23].
These techniques utilize a laser to excite fluorophores within a region of
interest (ROI). Moreover, a second laser is used to deterministically switch
off the fluorophores in the diffraction limited vicinity of the target
fluorophore to accomplish sub-diffraction resolution. On the other hand,
Stochastic Optical Reconstruction Microscopy (STORM) [24], direct Stochastic
Optical Reconstruction Microscopy (dSTORM) [25], Photoactivated Localization
Microscopy (PALM) [26], Fluorescence Photoactivation Localization Microscopy
(FPALM) [27], Super-resolution Optical Fluctuation Imaging (SOFI) [28], DNA-
Point Accumulation for Imaging in Nanoscale Topography (DNA-PAINT) [29, 30],
and other variants introduced later [25, 29, 30, 31, 32, 33, 34, 35, 36, 37,
38, 39] can be categorized as stochastic switching techniques. These
techniques often make use of a wide-field microscopy method, such as TIRF, to
activate a small random subset of fluorophores to avoid simultaneous
fluorescent light from more than one fluorophore in a diffraction limited
region, which in turn allows localization of these fluorophores with a
precision better than the diffraction limit.

Most of the stochastic switching super-resolution techniques include an image
processing stage to localize fluorophores and are called single molecule
localization microscopy (SMLM), or just simply localization microscopy. In
what follows, we first discuss the SMLM techniques in more depth and next
focus on the processing and post-processing of the data produced from these
techniques.

##  2 Single Molecule Localization Microscopy

In SMLM, fluorophores can be stochastically switched between fluorescent and
dark states, then many images can be taken as a series in each of which only a
few fluorophores are emitting. The acquired images are next used to localize
isolated fluorophores with much higher precisions than the resolution limit,
which is given by

$\sigma\sim\frac{\sigma_{\mathrm{PSF}}}{\sqrt{I}}$ (4)

in the absence of background noise, where $\sigma$, $\sigma_{\mathrm{PSF}}$
and $I$ are the localization precision, size of the PSF and the number of
emitted photons by/intensity of the fluorophore, respectively. The resulting
precise localizations of these fluorophores can be built up into a super-
resolution reconstruction of the underlying structure, Fig. 3. However, major
detriments to accurate and precise localizations are: overlapping PSFs when
multiple fluorophores in close proximity are activated [40, 41]; small photon
counts per blinking event resulting in a low signal to noise ratio (SNR); poor
labeling of target structures; and low precision in axial location of
fluorophores.

Avoiding dense data regions with overlapping PSFs requires sparse activation
of fluorophores to prevent activation of more than one fluorophore in a
diffraction limited area [24, 25, 26, 27, 31, 32, 33]. Some other SMLM
approaches overcome the problem of overlapping PSFs by separating the
fluorescent signals from molecules with different emission spectra [34, 35,
36] or lifetimes [42]. Photo-bleaching of fluorophores has also been employed
to achieve low density images of PSFs with minimum overlaps [37, 38]. Another
promising SMLM technique is based on stochastic binding and unbinding of the
diffusing fluorescent emitters to the target, such as DNA-PAINT [29, 30, 39],
and lifeact [43, 44].

As mentioned, the sparsity constraint imposes activation of a small subset of
fluorophores and therefore a few localizations per frame. On the other hand,
enough localizations are required to obtain high-resolution images [45], which
in turn demands undesired long data acquisition times [46, 47]. This problem
can be alleviated by multiple-emitter fitting methods that are able to
localize emitters in denser regions of the data [48, 49, 50]; see section 4.3.

Here, we proceed with discussing the low photon counts and labeling
challenges. For various fluorescent probes, the number of emitted photons per
blinking event ranges from a few hundred to a few thousand with a small number
of blinking events per probe, where more photons and a larger number of
blinking events are desired for better localization precision and image
contrast [51]. Furthermore, optimal labeling of the sample is required to
achieve high resolution images. For optimal labeling, fluorophores have to
only stain target structures with a uniform and high labeling density.
Moreover, standard labeling procedures use antibodies to attach fluorophores
to a target structure. Some of these antibodies directly bind to the structure
and are called primary antibodies. On the other hand, some antibodies bind to
target structures via another antibody and are called secondary antibodies.
Therefore, the structures that are labeled using secondary antibodies often
look artificially larger [52].

Figure 3: SMLM concept. (a, c) Sample frames of raw SMLM data. (b) Image
reconstructed using localizations from a single frame. (d) Image reconstructed
using localizations from 50 frames. (e) Diffraction limited image, which is
the sum image of 5000 frames of raw SMLM data. (f) Super-resolution image
reconstructed using precise localizations from 5000 frames.

Cellular organelles are inherently 3D structures and therefore 3D microscopy
images are desirable. Super-resolution SMLM approaches are able to provide 3D
images of cellular structures by exploiting the variance in the PSF shape as a
function of axial distance from the focal plane. However, the PSF shape
changes slowly along the axial axis and is almost symmetric below and above
the focal plane, which hinder precise localization of fluorophres in the axial
direction. This problem can be overcome by PSF engineering, which
intentionally adds aberrations to the microscope to achieve faster changes in
PSF and thus higher axial precisions. 3D STORM microscopy was first
demonstrated by an engineered astigmatic PSF [53] and was later perfected
using more complex engineered PSFs [54, 55, 56, 57] such as the tetrapod PSF
[56] (see section 4.4).

##  3 Image Formation

Raw SMLM data is comprised of a sequence of image frames, where each frame is
a two dimensional array of pixels whose values are a function of the number of
photons captured by the camera over a fixed exposure time. Photons reaching
the camera often originate from multiple sources of light: 1) the fluorescent
light from in-focus fluorophores that label the target structure; 2) the
fluorescent light from out-of-focus fluorophores which might be bound to
undesired cellular structures, or autofluorescence from these structures; 3)
other sources of light that might exist in the area such as cosmic rays. The
contribution of light from sources other than the in-focus fluorophores is
unwanted and degrades the quality of images [45, 58]. The undesired light
reaching the camera gives rise to two types of background noise: 1) an
approximately uniform, homogeneous background; and 2) a heterogeneous
background called structured background [58, 59, 60, 61].

For each data frame, a subset of emitters is activated and the image frame is
described by

$\mathrm{Model}(x,y)=\sum_{i}I_{i}\delta(\xi-x_{i})\delta(\eta-
y_{i})\ast\mathrm{PSF}(x,y;\xi,\eta,z_{i})+\mathrm{b},$ (5)

which is the convolution of the emitter ${x}$ and ${y}$ locations with the PSF
plus a uniform background, Fig. 4. The sum is over the activated emitters.
$I_{i},x_{i},y_{i},z_{i}$ and $\delta$ represent the intensity and location of
the $i$th emitter and Dirac delta, respectively. $x$ and $y$ denote image
coordinates, and $\xi,\,\eta$ are auxiliary parameters.

Here, we explain different terms present in model (5) in more detail. First,
the additive term in eq. (5) models the homogeneous uniform background, while
the structured background, usually coming from the out-of-focus emitters, is
mixed with the in-focus emitters and is given by the convolution term.
Furthermore, note that while $x$ and $y$ coordinates only determine the
lateral location of the PSF, the PSF shape is a function of the ${z}$-location
(offset from the focal plane) of the emitter [59, 62]. Therefore, the out-of-
focus emitters have a different PSF compared to the in-focus emitters.
Moreover, some effects like dipole orientation [63, 64], sample movements [65,
66] and optical aberrations [64, 67] result in distortions of the PSF as well.
In the absence of these effects, the PSF has an approximate 2D Gaussian form
[40, 67, 68, 69, 70, 71, 72, 64]

$\mathrm{PSF}(x,y;x_{i},y_{i},z_{i})=\frac{1}{2\pi\sigma_{\mathrm{PSF}}(z_{i})^{2}}\exp\left[-\frac{(x-x_{i})^{2}+(y-y_{i})^{2}}{2\sigma_{\mathrm{PSF}}(z_{i})^{2}}\right].$
(6)

However, in cases where the Gaussian PSF is not an appropriate approximation,
theoretical PSF models [68, 73, 74] can be used or numerical PSFs can be
obtained from calibration experiments [75, 76, 77, 78]. Such calibrations
yield a numeric array of PSF samples which are then utilized to generate a
model [76, 78].

The pixel values recorded by the camera are not the same as the photon counts.
Rather, they are stochastic quantities due to noises introduced in the
detection process and the stochastic nature of photons. Here, we first explain
the main noises added by detection instruments (camera). The camera detectors
amplify the signal from the detected photons and a stochastic number of
photoelectrons are produced per photon which is used to generate the pixel
values [79, 80]. Moreover, the process of converting photoelectrons to pixel
values is also a stochastic process that introduces a random offset called the
read-out noise, which contributes to the detection noise. To obtain the
correct photon counts, therefore, these two effects have to be taken into
account.

We have so far discussed noises introduced by detection instruments, which
might be reduced or totally removed by future advances in the instrument
technology. However, another major source of stochasticity is shot noise (also
called Poisson noise), which comes from the particle nature of photons and
cannot be removed regardless of the imaging setup or the used equipment. The
expected photon count over a fixed exposure time for every individual pixel
can be calculated from model (5). However, the number of photons captured by a
detector over this period slightly deviates from the expected values and has a
Poisson distribution, which is called shot noise [81, 82, 83, 68], Fig. 4.

A complete noise model therefore has to consider both detection and shot
noises, which might be either pixel dependent (for sCMOS camers [82, 84]) or
pixel independent (for CCD and EMCCD cameras [79, 83, 80]). However, such
noise models are computationally expensive and thus appropriate approximate
forms are employed to formulate image formation for different cameras [83, 85,
82, 80, 84, 81, 68].

Figure 4: Image Formation. (a) Emitter locations. (b) Model is a convolution
of emitter locations with the PSF plus a uniform background. (c) The SMLM
image recorded by a detector is the model in panel b corrupted with the shot
(Poisson) noise. The emitters are assumed to be in-focus.

##  4 Image Processing

Image processing is a key step in super-resolution SMLM approaches. This step
is comprised of multiple stages: 1) pre-processing; 2) identification of
candidate emitters, and localization; 3) filtering and image rendering. Here,
we briefly describe each stage and more detailed discussions are postponed to
the subsequent sections.

The pre-processing step often alleviates the noise introduced during data
acquisition, such as camera and background noise. The estimation of emitter
parameters comes next. This step include finding locations and/or intensities
of emitters using the point spread function, which predates the advent of SMLM
approaches and had been employed in other scientific disciplines such as
electrical engineering [86, 87, 88], astronomy [60, 89, 83, 90], and particle
tracking [91, 92]. There are two major classes of localization approaches:
single-emitter localization algorithms, and multiple-emitter localization
algorithms. The single-emitter algorithms are only able to localize isolated
emitters where there are no overlapping PSFs. Single-emitter candidates are
usually identified by applying a threshold to detect local maxima and then
ROIs of a certain size including those local maxima are selected for further
analysis [93, 68, 69, 94]. Other common detection algorithms employ a wavelet
transform [95, 96, 97], different types of filters [98, 40] and other
techniques [99, 100] to identify candidate emitters. The performance of these
detection algorithms is highly correlated with SNR and signal to background
ratio. For a comparison of different detection algorithms see [101, 102].
After detecting the candidate emitters, single-emitter localization procedures
estimate their locations employing either a non-iterative or an iterative
algorithm.

However, a main challenge of the single-emitter algorithms is localizing
emitters in dense regions of data with multiple overlapping PSFs from closely
spaced activated emitters, Figs. 4-5. Such dense regions can be generated due
to either dense labeling of the sample or fast data collection. Multiple-
emitter fitting algorithms are capable of localizing emitters in such dense
regions of the data with overlapping PSFs. There are many such multiple-
emitter approaches that may be categorized based on their outputs [48] or
based on the algorithm itself [49, 50].

After the localization stage, there is often a rejection step that filters out
bad localizations to reduce artifacts in the final reconstructions [103, 104].
A popular filtering criteria is based on the found parameter values and their
uncertainties that removes the localizations with uncertainties larger than
given thresholds [24, 103, 61]. An additional filtering approach is based on
the nearest neighbor distances (NND) where localizations with less than $N$
neighbors within a certain distance are eliminated from the final list of
localizations [105, 106].

An example of bad localizations is fitting two or more overlapping PSFs as a
single emitter which gives rise to artifacts and degrades image quality. To
reduce this effect, the localization algorithms make use of different criteria
for recognizing these types of bad fits, Fig. 5. This can be done before
localization where ROIs with overlapping PSFs are recognized by a deviation in
shape of the identified bright blob from the given PSF [24]. Alternatively,
these types of bad fits can be recognized after the localization stage by
using a $\mathrm{p}$-value test. In this test, the p-values are calculated
assuming every individual localization results from a single PSF (the null
hypothesis). If the found p-values are smaller than a given threshold, the
fits are rejected [68]. Finally, the remaining localizations are used to
reconstruct a super-resolved image, Fig. 5.

###  4.1 Background Detection

The ultimate objective of SMLM techniques is reconstructing high resolution
images from precise and accurate [59] estimates of the emitter locations from
raw SMLM data. In order to accomplish this goal, a correct measure of
background noise is required as incorrect background leads to biased position
estimates [48]. The correction of uniform background noise is simple and
various approaches have been conceived to address this issue. These approaches
usually select a ROI and use that to compute a local uniform background noise
as the average of pixel values, median of pixel values, average of pixel
values after bleaching of the fluorophores, $X$th percentile of pixel values
or estimating the additive term in eq. (5) using an iterative approach [62,
68, 107, 108, 109, 110].

Structured background is significantly more complicated to remove and its
presence results in poor position estimates. A few methods have been put
forward to cope with this problem including an approach that uses a temporal
median filter to subtract structured background from the signal [111]. This
technique inspects the fluctuations of pixel values over time to find the
background value. The found background can be overestimated in dense regions
of the data where there is at least one active emitter at each time. An
alternative procedure detects all the emitters regardless of being signal or
background and then sets a threshold, removing the emitters with intensities
below that as structured background [112]. Wavelet decomposition has also been
employed to subtract structured background as well as uniform background prior
to emitter detection [113]. Recently, a deep learning method has been proposed
to detect structured background using the PSF shape [58]. Moreover, Faze1, et
al. used a Bayesian approach to model structured background with a collection
of PSF-sized dim emitters [61]. In the field of astronomy, methods such as
sigma clipping has been developed to deal with structured background in dense
data sets [60]. In the sigma clipping procedure, the brightness mean,
$m_{\mathrm{I}}$, and standard deviation, $\sigma_{\mathrm{I}}$, are
calculated and those intensities outside the range of
$[m_{\mathrm{I}}-\alpha\sigma_{\mathrm{I}},\,m_{\mathrm{I}}+\alpha\sigma_{\mathrm{I}}]$
are considered noise [60].

###  4.2 Single Emitter Fitting

The single-emitter localization algorithms can be classified into two major
categories: the algorithms that use non-iterative approaches to localize
emitters and the algorithms that use an iterative procedure. Studies show that
the iterative algorithms are more accurate than the non-iterative algorithms
[114]. However, iterative algorithms are computationally more demanding and
require a precise PSF model.

####  4.2.1 Non-iterative Algorithms

Non-iterative algorithms do not need any information about the PSF and are
usually fast and easy to implement. However, they are not as accurate as
iterative algorithms that utilize the PSF to generate a model of the data. The
lack of enough accuracy is often a consequence of different types of noise.

A few non-iterative approaches such as QuickPALM [115] calculate the emitter
locations as center of mass of ROIs containing single emitters [115, 54]. This
gives a good estimate of location, however, failure in background correction
results in biased localizations towards the center of the ROIs. Virtual Window
Center of Mass (VWCM) [116] ameliorates this issue by iteratively adjusting
the selected ROI to minimize the separation of the emitter location and the
center of the ROI.

FluoroBancroft borrows the Bancroft procedure from the satellite Global
Positioning System (GPS) to localize emitters [95, 117]. This approach uses
three pixel values within a PSF to draw three circles. The emitter is located
within the intersection of these three circles. The size of the intersection
region is a measure of the localization precision. A correct measure of
background is also of great importance in this approach to calculate accurate
radii.

Single emitters can also be localized by finding the gradient of the phase of
the Fourier transform of the ROIs. For a single emitter, eq. (5) reduces to

$I(m,n)=I_{0}\delta(\xi-x_{0})\delta(\eta-
y_{0})\ast\mathrm{PSF}(x_{m},y_{n};\xi,\eta,z_{0})+b$ (7)

where $m$ and $n$ count rows and columns of pixels in the ROI, and $x_{m}$ and
$y_{m}$ give the centers of those pixels. $(x_{0},y_{0},z_{0})$ is the
location of the single emitter. The Fourier transform of eq. (7) produces

$\tilde{I}(k,l)=H(k,l)\,\mathrm{exp}\left[-i2\pi\left(\frac{x_{0}}{M}k+\frac{y_{0}}{N}l\right)\right]+\tilde{b}$
(8)

where $k,l$ and $M,N$, respectively, denote pixel row and column indices and
the maximum number of pixels in rows and columns. Furthermore, $H$ is a real
quantity. For data sets with large SNR, the background term is negligible and
the Fourier Domain Localization Algorithm (FDLA) gives the emitter position as
the average of the gradient of the phase [118]:

$x_{0}=\mathrm{mean}\left(\left(\frac{\partial\phi}{\partial
k}\right)\frac{M}{2\pi}\right),\,\,\,\,y_{0}=\mathrm{mean}\left(\left(\frac{\partial\phi}{\partial
l}\right)\frac{N}{2\pi}\right)$ (9)

where
$\phi=\mathrm{arctan}\frac{\mathrm{Im}(\tilde{I})}{\mathrm{Re}(\tilde{I})}$.
The performance of this approach suffers from the presence of background noise
as well. Another approach localizes single emitters by calculating the first
Fourier coefficients in both $x$ and $y$ directions, and the phase of these
coefficients are then employed to find the emitter location [119].

Radial symmetry of the PSF has also been employed to calculate emitter
locations [98, 120]. Due to the radial symmetry of PSFs, intensity gradient
vectors for different pixels converge to the region with maximum intensity
where the emitter is located. This approach is robust in the presence of
uniform background noise and achieves precisions close to the theoretical
lower limit of precision given by the Cramer-Rao Lower Bound (CRLB).

####  4.2.2 Iterative Algorithms

Iterative algorithms are the most rigorous approaches for emitter
localization. In these approaches, the parameters in model (7) are adjusted in
an iterative manner to fulfill a certain criterion. In the localization
problem, the parameters are ($x_{0},y_{0},z_{0},I_{0},b$), the emitter
location, the number of photons from (intensity of) the emitter and a uniform
background noise. The criteria that are extensively utilized in the emitter
fitting literature are the Least Square (LS) difference between data and
model, and maximizing the likelihood function via a Maximum Likelihood
Estimate (MLE).

CRLB states that the fundamental limit of variance for estimating a parameter
from given data is obtained by the inverse of the Fisher Information [93, 81,
121]. Therefore, the fundamental limit of precision is given by the inverse of
the square root of the Fisher Information. Theoretically, MLE achieves the
best localization precision, equivalent to the CRLB [81, 59, 68, 67, 121, 122,
73, 123, 124]. LS performance is comparable to the MLE under certain
conditions described below [67, 73, 123, 124, 125].

The performance of weighted LS approaches that of MLE at high SNR, when the
Poisson noise (shot noise) can be well approximated by a Gaussian model, or
when read-out noise is dominant. Note that neither of these scenarios are
correct for SMLM data where the read-out noise is usually negligible in the
presence of Poisson noise and the SNR is not too high. In general, MLE yields
better localization accuracy and is more robust in the presence of PSF
mismatch, but is computationally more complex and requires an accurate model
of noise [67, 73].

Least Squares Fitting. The LS approaches iteratively vary the parameters of
the model to minimize the sum of differences between the pixel values from the
data and the model. The developed algorithms use a Gaussian PSF [93, 67],
theoretical PSFs [126, 127] or experimentally acquired PSFs [128, 129] to
calculate a model using eq. (7). The difference of the generated model and
data is then given by

$D=\sum_{\mathrm{pixel}}\,\frac{\left(\mathrm{data}-\mathrm{model}\right)^{2}}{\mathrm{expected\,variance}}$
(10)

where in weighted LS (WLS) the differences are scaled by the expected variance
of the noise, which scales the errors for individual pixels [54, 67, 128, 70,
126, 127]. A pixel with a high signal is expected to have a large noise
variance, and therefore it is allowed to have a larger error in the WLS
procedure. However, the scaling factor is replaced by one in the unweighted LS
algorithm, which we call LS hereafter, and does not accommodate noises [93,
94, 130, 129, 131]. The Levenberg-Marquardt iterative procedure [126, 129,
131], or other procedures [128], are then employed to iteratively adjust the
parameters of the model to reduce the difference (10).

The WLS algorithm accomplishes accuracies close to the CRLB when the photon
count is high, but the noise variance needs to be known as well as an accurate
PSF model. The PSF mismatch, particularly in the tail of the PSF, results in
large errors when scaled by a small expected noise variance in the pixels far
from the emitter [48]. Therefore, the LS algorithm is more suitable when a
reasonable PSF model and/or noise model are not accessible.

Maximum Likelihood Estimator. To calculate the likelihood, we first need to
find expected photon counts per pixel. The contribution of an emitter located
at $(x_{0},y_{0},z_{0})$ to the $k$th pixel’s photon count is given by

$\Delta_{k}=\frac{I_{0}}{2\pi\sigma_{\mathrm{PSF}}(z_{0})^{2}}\int_{x_{k}-0.5}^{x_{k}+0.5}\int_{y_{k}-0.5}^{y_{k}+0.5}\mathrm{exp}\left[\frac{(x-x_{0})^{2}+(y-y_{0})^{2}}{2\sigma_{\mathrm{PSF}}(z_{0})^{2}}\right]dx\,dy$
(11)

where we used the Gaussian PSF (6).
$\Delta_{k},\sigma_{\mathrm{PSF}},I_{0},x_{k}$ and $y_{k}$, respectively,
denote the number of photons in the $k$th pixel from the emitter, the half
width of the Gaussian PSF, total number of photons from the emitter, the
emitter location and the center of the $k$th pixel.

Figure 5: Reconstructions from dense data with overlapping PSFs. (a) A frame
of dense raw SMLM data of a cross with parallel lines. The green circle shows
an example of two overlapping PSFs. (b) Reconstruction from a multiple-emitter
algorithm with no filtering. (c) Reconstruction from single-emitter algorithm
with no filtering. The localizations in the area between the lines are the
results of fitting overlapping PSFs with a single PSF. (d) Reconstruction from
single-emitter algorithm after filtering. The dense regions of data appear
sparse when processed with a single-emitter algorithm due to the inability to
localize overlapping PSFs, which is called the contrast inversion artifact.

The total photon count in the $k$th pixel is the sum of the photons from the
emitter and the uniform background noise

$\lambda_{k}=\Delta_{k}+b.$ (12)

Equation (12) yields the expected number of photons for pixel $k$ for a fixed
exposure time. Consequently, the number of detected photons in pixel $k$ has a
Poisson distribution, which gives the likelihood of the $k$th pixel

$P_{k}(D|\theta)=\frac{\lambda_{k}^{D_{k}}e^{-\lambda_{k}}}{D_{k}!}$ (13)

where $\theta$ stands for the set of parameters
($\theta=(x_{0},y_{0},z_{0},I_{0},b)$, where $z_{0}=0$ for 2D cases). $D$
represents data, which is a two dimensional array of pixels whose values are
related to the number of photons detected by the camera. $D_{k}$ selects the
$k$th pixel in $D$. Since the pixels are independent, the likelihood of the
ROI can be written as the product of the likelihoods of all the pixels in the
ROI

$P(D|\theta)=\prod_{k}P_{k}(D|\theta).$ (14)

The two main iterative algorithms employed in the literature to find the
parameters that optimize the above likelihood are variations of the Newton
method [67, 68, 71, 72, 74, 76, 77] and a modified version of Levenberg-
Marquardt [40, 70, 78, 132] adopted from LS procedures.

The Newton approach is employed to find the root of the derivative of the
likelihood function (14), and therefore one needs to calculate the second
derivative of the likelihood as well, which is computationally demanding. On
the other hand, the Levenberg-Marquardt algorithm only calculates the first
derivative of the likelihood, which makes it computationally less demanding in
comparison [70, 78]. Different strategies have been exploited to speed up the
Newton optimization algorithm including implementation on Graphical Processing
Units (GPUs), which allows parallel analysis of ROIs [68, 69, 71, 72, 77];
starting from better initial values [73]; and estimating $x$ and $y$ positions
individually utilizing the separability property of the Gaussian function
[72].

###  4.3 Multiple Emitter Fitting

Raw SMLM data is often acquired via a wide-field illumination procedure where
the whole sample is exposed to the excitation laser. Emitter activation is a
stochastic process, and hence there will always be activated emitters at close
proximity. Therefore, overlapping PSFs are unavoidable, even under sparse
activation conditions. The overlapping PSFs are eliminated in a filtering step
in single-emitter approaches, which results in losing information [122] as
well as the appearance of artifacts, for instance, contrast inversion, Fig. 5.

The inability of single-emitter algorithms to localize multiple activated
emitters in a diffraction limited vicinity enforces the sparse activation of
emitters. This is followed by a long acquisition time to build a full list of
localizations to reconstruct an image with high resolution. In some
experiments, such as studies of live or dynamic samples, fast data acquisition
is preferred and hence dense activation of emitters is inevitable. Therefore,
proper analysis of dense SMLM data with overlapping PSFs is necessary to
reduce data acquisition time, avoid artifacts and facilitate live sample
imaging.

Numerous multiple-emitter fitting algorithms have been devised to fit emitters
with overlapping PSFs, some borrowed from other areas, such as astronomy and
statistical analysis. The reported algorithms employ a very wide range of
approaches and have a broad spectrum of performance [49, 50]. These procedures
are often iterative algorithms or include an iterative training step.

####  4.3.1 Least Squares

The LS algorithm has been employed to fit multiple-emitters in the field of
astronomy [60]. It was modified for SMLM super-resolution microscopy, called
DAOSTORM in this context [133]. DAOSTORM uses isolated emitters in the raw
data to estimate the PSF, and then uses the found PSF to fit emitters in dense
data sets. The algorithm starts with an initial number of emitters located at
the brightest pixels of the ROI, and then uses least squares fitting to
localize them with sub-diffraction precision. Next, the residuum image is
calculated by subtracting the model from the data, and is used to detect new
emitters in the pixels brighter than a given threshold. The detected emitters
are then localized to obtain sub-diffraction precision. This step is repeated
until there is no pixel with intensity above the threshold in the residuum
image.

####  4.3.2 Maximum Likelihood

The MLE approach that was described before can be modified for multiple-
emitter fitting within ROIs with overlapping PSFs. The total photon counts in
the $k$th pixel is given by

$\lambda_{k}(N)=b+\sum_{i=1}^{N}\Delta_{k,i}$ (15)

where $\Delta_{k,i}$ is the number of photons received in the $k$th pixel from
the $i$th emitter and can be calculated using (11). $N$ and $b$ are the number
of emitters and the uniform background. The likelihood of the $k$th pixel is
then

$P_{k}(D|\theta)=\frac{\lambda_{k}(N)^{D_{k}}e^{-\lambda_{k}(N)}}{D_{k}!}.$
(16)

The likelihood of the ROI is obtained from the product of the likelihoods of
individual pixels (14).

The likelihood (16) has more than one emitter and therefore more parameters to
estimate, demanding more iterations and computational time. The MLE approach
is implemented in the same manner as single-emitter fitting to estimate the
parameters. Nevertheless, there is a new parameter, $N$, the number of
emitters, which cannot be directly estimated from the likelihood itself. This
is because adding more emitters tends to increase the likelihood, known as
over-fitting. The approaches that find the number of emitters are called model
selection algorithms. Several model selection algorithms have been reported
along with the MLE localization procedure, including thresholding of the
residuum image [134], p-value of the likelihood ratios [122], Bayesian
Information Criteria (BIC) [135, 136], PSF size and NND [137], and others
[138, 139, 140].

The 3D-DAOSTORM [134] is a 3D multiple emitter fitting approach. 3D fitting
procedures will be discussed in the next section and here we explain how this
procedure deals with overlapping PSFs. 3D-DAOSTORM fits overlapping PSFs by
fitting only choices from the brightest emitters in the ROI at the beginning.
It then subtracts the obtained model from the ROI intensities and uses the
residuum image to find pixels brighter than a given threshold to detect new
emitters. It employs MLE to localize the new emitters. The new emitters are
added to the list of detected emitters and this step is repeated until there
is no pixel brighter than the given threshold.

Simultaneous multiple-emitter fitting [122] starts from a model with one
emitter, $N=1$, and goes up to $N=N_{max}$. For each model, this method
localizes the emitters employing the MLE approach. The log-likelihood ratio
(LLR)

$\mathrm{LLR}=-2\log\left[\frac{P(D|\hat{\theta})}{P(D|D)}\right]$ (17)

has an approximate chi-square distribution [122], where $\hat{\theta}$ is the
parameters that maximize the likelihood and $P(D|D)$ gives the upper limit of
the likelihood. This LLR is then used to calculate p-values where the model
with lowest $N$ that meets a threshold p-value is accepted as the best fit.

As mentioned above, the MLE approach suffers from over-fitting, and adding
more parameters (emitters) tends to give larger values for likelihoods.
Bayesian Information Criteria is a model selection algorithm that penalizes
adding new parameters to an MLE problem. SSM-BIC [135] selects the model that
maximizes the following function

$\mathrm{BIC}(N)=\frac{(\mathrm{Data}-\mathrm{Model})^{2}}{\mathrm{Data}}+(3N+1)\log{(m)}$
(18)

where $m$ is the number of pixels in the given ROI. Note that there are $3N+1$
parameters in a 2D model with $N$ emitters and a known PSF. This approach has
also been implemented on GPUs with $\sim$100 times decrease in computational
time [136].

QC-STORM [137] uses a weighted likelihood

$P_{W}(D|\theta)=\prod_{k}W_{k}P_{k}(D|\theta)$ (19)

to localize emitters, where $W_{k}$ is the weight of the $k$th pixel and is
smaller for pixels closer to the edges of the ROI. The weighted likelihood
therefore suppresses the signal close to the edges of the ROI which makes it
an effective method to localize emitters within ROIs with signal
contaminations close to their edges. QC-STORM identifies ROIs with more than
one emitter based on the ROIs’ NNDs and the size of the PSF estimated using
weighted MLE. This algorithm has been implemented on GPUs and is capable of
processing very large fields of view.

Some other approaches use an iterative deconvolution algorithm to accomplish
maximum likelihood employing the Richardson and Lucy procedure [141, 142].
These approaches return a grid image with finer pixel sizes than the camera
pixel size, with non-zero pixel values at the emitter locations rather than
returning a list of localizations.

####  4.3.3 Bayesian Inference

The MLE algorithm suffers from overfitting as discussed above. Bayes’ formula
(20) provides an elegant way to include prior knowledge into the problem,
allowing the model to be restricted to reasonable number of parameters. It
however adds complications to the problem by including more distributions.
Furthermore, it has been shown that the Bayesian approach can achieve
localization uncertainties better than those from MLE by inclusion of
reasonable prior knowledge [143]. The Bayesian paradigm is equivalent to MLE
when there is no prior knowledge available. The posterior is given by

$P(\theta|D)=\frac{P(D|\theta)P(\theta)}{P(D)}$ (20)

where $P(\theta)$ and $P(D)$ are, respectively, the prior on the parameters
$\theta$ and the normalization constant,

$P(D)=\int P(D|\theta)P(\theta)d\theta$ (21)

called the evidence. Another difference of MLE and Bayesian approaches is that
MLE returns fixed emitter parameters (location and intensity), while the
Bayesian procedure returns a probability distribution, the posterior, for
every emitter parameter. A few fully Bayesian algorithms have been developed
for multiple-emitter fitting so far [61, 144, 145], where we discuss two of
them here.

The Bayesian multiple-emitter fitting (BAMF) algorithm [61] employs the
Reversible Jump Markov Chain Monte Carlo (RJMCMC) [146, 147] technique to
explore model spaces with different dimensions, or equivalently different
number of emitters, and returns a posterior distribution which is a weighted
average of different possible models. BAMF uses either a Gaussian PSF model or
an arbitrary input numerical PSF along with an empirical prior on emitter
intensities to attain precise and accurate localizations in dense regions of
data. This technique performs emitter fitting, model selection and structured
background modeling simultaneously, and therefore takes into account various
sources of uncertainties which are often ignored.

The 3B algorithm [144] analyzes the entire data set at the same time by
integrating over all possible positions and blinking events of emitters. This
algorithm makes inferences about the emitter locations, intensities, width of
the Gaussian PSF and blinking of emitters. The posterior of this problem is
given by

$P(a,b,N|D)=\frac{P(D|a,b,N)P(a)}{P(D)}$ (22)

where $a,b$ and $N$ represent the emitter parameters, blinking events and the
number of emitters, in turn. 3B uses a uniform prior for the locations and a
log-normal prior for other parameters. The discrete parameter, $b$, is then
integrated out using an MCMC [148] approach to obtain the posterior
distribution of $a$, $P(a,N|D)$. Next, the Maximum A Posteriori (MAP) is
computed using the conjugate gradient approach to obtain the emitter
locations, intensities and the PSF size. After that, the parameter $a$ is also
marginalized to get the model probability, $P(N|D)$, for model selection.

The 3B algorithm models the entire data set and is able to use all the
collected photons from multiple blinking events to achieve better localization
precision. The returned image is a probability map of the weighted average of
all possible models rather than a selected single model [149]. Moreover, it
needs a simple experimental setup for data collection [144]. However, this
technique is very slow because calculating the integrals to marginalize
parameters $a$ and $b$ is extremely computationally demanding. There have been
several attempts to speed up the algorithm including 3B implementation in
cloud computing [150], use of more informative priors [151], and initializing
the algorithm with better starting parameter values [152].

####  4.3.4 Compressed Sensing

A frame of SMLM data can be considered as a matrix $y$

$y=Ax+b$ (23)

where $x$ is the signal, which is an up-sampled discrete grid (image) with
non-zero elements at the emitter locations, $A$ is the PSF matrix, and $b$ is
the uniform background. The objective is to recover the non-zero elements of
the signal $x$ where most of the elements are zero due to sparse activation of
fluorescent emitters in SMLM microscopy. Compressed Sensing (CS) theory states
that a signal $x$ can be recovered from a noisy measurement $y$ if the signal
is sufficiently sparse [153]. This mathematically can be expressed as

$\displaystyle\mathrm{Minimize:}$ $\displaystyle\,||x||_{1}\hskip 79.6678pt$
(24) $\displaystyle\mathrm{Subject\,to:}$
$\displaystyle\,||y-(Ax+b)||_{2}\leq\epsilon$ (25)

where $||x||_{1}=\sum_{i}|x_{i}|$ is the L1-norm of the up-sampled image, and
the L2-norm of the residuum image is given by

$||y-(Ax+b)||_{2}=\sqrt{\sum_{i}\left(y_{i}-(Ax+b)_{i}\right)^{2}}.$ (26)

The inequality allows fluctuations from the expected model due to different
types of noise.

Various optimization approaches have been utilized to minimize the L1-norm in
the presence of the given restriction in (25) including a convex optimization
algorithm [154], L1-Hotomopy [155], gradient descent [113], and others [156,
157]. These algorithms are able to detect and localize the emitters in very
dense regions of the data. However, due to the large size of the up-sampled
image, the CS algorithms are slow and the resolution cannot be better than the
grid size of the up-sampled image.

The issues mentioned above have been addressed in later literature using
different approaches. FALCON [113] accelerates the algorithm by implementing
CS on GPUs. It also ameliorates the grid-size problem by refining the found
locations in the subsequent steps after the deconvolution stage. CS has
recently been implemented over continuous parameter spaces to remove the
limits imposed by the up-sampled grid [158, 159]. To lower the computational
cost of the CS algorithm, a recent paper models the entire sequence at the
same time rather than employing a frame by frame analysis of the data [160].
Another approach implements the CS algorithm in the correlation domain by
calculating the frame cross-correlations [161].

#### Singular Value Decomposition

Assuming $A$ as an $n\times n$ square matrix with non-zero determinant, it can
be factorized into

$A=U\Sigma U^{-1}$ (27)

where U is an $n\times n$ matrix whose columns are the eigenvectors of the
decomposition, and $\Sigma$ is a diagonal $n\times n$ matrix where the
diagonal elements are the eigenvalues. A non-square matrix $B$ can also be
decomposed in a similar fashion, called the singular value decomposition (SVD)

$B_{n\times m}=V_{n\times n}\Lambda_{n\times m}W_{m\times m}$ (28)

where $V$ and $W$ are, respectively, $n\times n$ and $m\times m$ matrices.
$\Lambda$ is a diagonal $n\times m$ matrix with diagonal elements the square
roots of the eigenvalues of $BB^{T}$, where $T$ stands for transpose [162].

The MUltiple SIgnal Classification ALgorithm (MUSICAL) for super-resolution
fluorescence microscopy [163] takes $B$ as a collection of frames of super-
resolution images where each column of $B$ is a frame of the raw data. $B$ is
then a non-square matrix that can be factorized into a diagonal matrix and
unitary matrices where the eigenvectors are eigenimages. The eigenimages are
next classified into signal and noise based on the corresponding eigenvalues
using a given threshold. Finally, MUSICAL calculates the projection of the PSF
at different locations in the eigenimages to identify and localize the
emitters. An alternative method makes use of the SVD in the Fourier domain,
analyzing the sequence of raw data frame by frame to localize the emitters
[164].

#### Deep Learning

Deep learning approaches are non-iterative optimization algorithms that
perform calculations in a parallel manner and hence are very fast. Theses
approaches are based on Artificial Neural Networks (ANN) that are inspired by
animal brains and neural systems. The building blocks of the brain are neurons
equivalent to perceptrons or sigmoid neurons in ANNs [165], which are
explained in the following.

Figure 6: Perceptron takes several binary inputs and gives a binary outcome.
$w_{i}$s are the weights of the inputs.

Perceptrons take a few binary inputs and generate a binary output, Fig. 6. The
output of the perceptron is one if the sum of inputs times their weights is
larger than a threshold and is zero otherwise, eq. (29).

$\mathrm{Output}=\begin{cases}0\hskip 14.22636pt\mathrm{if}\hskip
14.22636pt\sum_{i}w_{i}\,\,\mathrm{input}_{i}<\mathrm{threshold}\\\ 1\hskip
14.22636pt\mathrm{if}\hskip
14.22636pt\sum_{i}w_{i}\,\,\mathrm{input}_{i}>\mathrm{threshold}\end{cases}$
(29)

It can be shown that certain combinations of perceptrons produce logical
operations such as AND, OR and NAND, which are the underlying bases of
computation, and any mathematical function could be generated using them
[165]. Therefore, ANNs are capable of producing any mathematical function
using perceptrons. A network of perceptrons can be trained to perform
different tasks by adjusting the weights, $w_{i}$. On the other hand, sigmoid
neurons are more sophisticated versions of perceptrons where the output is a
value within the interval of $[0,1]$ rather than a binary output. Sigmoid
neurons are more flexible in the training process and are used in neural
networks.

ANNs have been employed to attack various problems in the field of biomedical
imaging [166, 167], specifically for SMLM image processing [168]. ANNs
designed to localize emitters in raw SMLM data can be categorized into three
different types based on their training approach:
1) ANNs are trained using simulated data where the ground truth is available
or using localizations found by a standard iterative algorithm [169, 170,
171]. In the training stage, the ANN learns by minimizing the sum of distances
between the found localizations and the given locations. Using synthesized
data, there will always be adequate training data.
2) ANNs are also trained using two well-registered sets of data acquired from
the same sample where one of them is used as ground truth. The ground truth
image has high SNR that can be acquired employing different procedures such as
confocal microscopy [172, 173], using an objective with high numerical
aperture [174], or using a sparse data set to reconstruct a super-resolved
image with high SNR [175] to train the network. In the training stage, the
network learns by minimizing the difference between the output image and the
acquired image with high SNR. The minimization of differences can be
implemented via a standard iterative optimization algorithm [172], or by using
a sub-network in the training stage, which takes the output of the ANN along
with the ground truth and labels the output as real or fake [174, 175].
3) In an alternative training approach, there is no ground truth used in the
training step. Instead, the network is trained by reproducing the input data
from the list of found emitters and minimizing the difference between the
original input and the synthesized image [176]. This training procedure is
called unsupervised learning.

The deep learning procedures have multiple advantages including: fast data
analysis, no required input parameters or thresholds, and performance
comparable to the MLE algorithms [171, 177]. However, the training process is
very sensitive and has to be done very carefully. Two major pitfalls of
training are the hallucination and generalization problems [166]. Deep
learning algorithms might make mistakes in identifying patterns from random
inputs when there is not adequate training, which is called the hallucination
problem. If there are new patterns that are not seen by the algorithm before,
it fits these new data by the old patterns, this is called the generalization
problem.

#### Others

WindStorm [178] uses a temporal filter to estimate the background and remove
it from the raw data. It then implements deconvolution by dividing the Fourier
transform of the clean image with the Fourier transform of the PSF. It next
performs frequency truncation. The recovered locations are given by the peaks
of the deconvolved image in the spatial domain. This is a fast and non-
iterative approach and the found locations can be used as initial values for
iterative algorithms.

Wedge Template Matching (WTM) [179] identifies and localizes emitters by
matching an entire or partial template of the PSF to the regions of the data
with overlapping emitters. The WTM algorithm picks either an entire or partial
PSF template based on the degree of overlapping between PSFs. It finds the
candidate pixels containing emitters using cross correlation of the template
with the image, and then finds the locations of the detected emitters with
sub-pixel accuracies.

Other approaches employ machine learning algorithms [180], independent
component analysis and a shape matching approach in the frequency domain
[181], and other algorithms [182] to identify and localize emitters in dense
regions of the data.

###  4.4 3D Emitter Fitting

Biological samples are 3D in nature and 3D microscopy approaches are required
to gain better insight into biological structures. A standard super-resolution
microscope cannot provide precise axial location of an emitter due to slow
changes in the 3D PSF as a function of axial position ($z$), and the symmetry
of the PSF below and above the focal plane [125], see the Airy PSF in Fig. 7.
Super-resolution microscopes can however provide precise axial location of an
emitter by some alterations in the optical setup. Several different
modifications have been reported in the literature encompassing multi-focal
methods [128, 183, 184, 185] which image multiple focal planes at the same
time; engineered PSFs [53, 54, 55, 56, 57, 186, 187] where the axial
information is encoded into the PSF shape; and other procedures [188, 189,
190, 191]. Here, we briefly discuss the multi-focal and PSF engineering
techniques.

The multi-focal approaches allow precise axial emitter localization by
providing multiple $z$-slices of the PSF. This can be achieved by splitting
the emission beam into multiple light-paths with different optical path
lengths [183, 184, 76]. The lateral locations of the PSF are often registered
to be the same on all paths and the differences along the $z$-axis are known.
Moreover, the intensity is often uniformly split between different planes.
Therefore, the localization task is finding emitter parameters (coordinates
and intensity) in only one plane using data from all the planes, which in turn
breaks the PSF symmetry and provides high axial precisions.

Figure 7: Different engineered PSFs at multiple $z$-locations (corrupted with
Poisson noise and a uniform background). The Airy pattern has a symmetric
shape above and below the focal plane. The astigmatic, double-helix and
tetrapod PSFs are engineered to break the symmetry below and above the focal
plane and yield better localization precision along the axial direction.

The PSF engineering approaches attain precise axial localization of emitters
by inscribing axial information in the shape of the PSF by simple wavefront
modifications. The modification in wavefront can be achieved by inserting
extra optical components [53, 55, 192] or computer controlled phase modulators
[54, 56, 187, 193] in the optical path. Several engineered PSFs have been
devised for 3D super-resolution microscopy including an astigmatic PSF [53],
double-helix PSF [54], phase ramp PSF [55], corkscrew PSF [186], self bending
PSF [187], and tetrapod PSF [56, 193], Fig. 7. These PSFs have special
characteristics, such as relative motions of different parts, width variations
in different directions that quickly changes as a function of axial location,
Fig. 7. The engineered PSFs cover different ranges of axial locations. The
tetrapod PSF was designed to accomplish precise lateral emitter locations over
the largest axial range compared to the other PSFs [56, 125, 50].

Most of the 2D emitter fitting algorithms discussed in the previous sections
can be adapted for 3D emitter fitting [50, 194]. These approaches employ
either a theoretical 3D PSF [53, 68, 92, 159, 195, 196, 197] or an
experimentally acquired numeric PSF to calculate the model [115, 130, 75, 76,
77, 78, 198, 199, 200, 201, 202, 203, 204]. More complicated PSFs such as
tetrapod or double-helix have very intricate features that cannot be exactly
expressed by an analytical function, so the PSF is generated using empirical
data from calibration experiments. In the following, we discuss the PSF
calibration procedure.

Data for PSF calibration is often a stack of images acquired at different
$z$-locations by placing a fluorescent emitter at the origin of the lateral
plane, $(x_{0}=0,y_{0}=0)$, and incrementing its axial location below and
above the focal plane. At far distances from the emitter, the PSF at the
$z_{i}$-plane in the stack is given by the scalar diffraction theory [3, 76]

$\mathrm{PSF}(x,y;z_{i})=\Bigg{|}\iint\left(\mathcal{P}(k_{x},k_{y})e^{2\pi
ik_{z}(k_{x},k_{y})z_{i}}\right)e^{2\pi
i\left(k_{x}x+k_{y}y\right)}dk_{x}\,dk_{y}\Bigg{|}^{2}$ (30)

where the frequency cutoff, $k_{\mathrm{max}}$, due to the microscope
objective (see Fig. 1), has to be considered in performing the integration.
The above integral gives the inverse Fourier transform of the term inside the
parentheses, which quantifies two effects that contribute to the wavefront
distortion and thus deformation in the PSF shape including: 1) the pupil
function, $\mathcal{P}=\mathcal{A}e^{i\Phi}$, is due to the scattering and
attenuation of light traveling via different media such as through the sample
itself ($\Phi$ can also be due to the phase intentionally added for PSF
engineering); 2) the phase due to the axial ($z$) location of the emitter
given by $k_{z}z_{i}$, where
$k_{z}^{2}={n^{2}}/{\lambda^{2}}-(k_{x}^{2}+k_{y}^{2})$. $n$ and $\lambda$
are, respectively, the refractive index of the medium and the wavelength of
light in vacuum.

The camera only records the square magnitude of the above integral (which is
the light intensity), while the phase information is mostly lost. However,
given the stack of intensity PSF images acquired at multiple $z$-locations,
the phase $\Phi$ can still be found employing phase retrieval techniques such
as the Gerchberg-Saxton algorithm [76, 205, 206] and the MLE approach [140].
Substituting the obtained phase in eq. (30), one is able to obtain the PSF at
any given $z$-location. However, the integral (30) is computationally
expensive and therefore the PSF model is generated at just a few $z$-slices.
The PSF at a desired location can then be numerically produced via either
linear or spline interpolation of appropriate generated $z$-slices [76, 78].

###  4.5 Drift Correction

Drift is a common problem in super-resolution microscopy procedures, in which
the sample alters its location overtime, resulting in distortion and
degradation in the quality of final images. SMLM microscopy reconstructs high-
quality images from a list of localizations collected over the course of an
experiment, where even slight disturbances in the experiment lead to serious
defects in the results. For instance, mechanical vibrations of the microscope
stage or fluctuations in temperature result in rotational or translational
movements of the sample during data acquisition. Since such disturbances are
unavoidable, algorithms are required to measure and correct for drift in the
image processing step.

Multiple drift correction algorithms have been employed in super-resolution
microscopy, such as use of fixed fiducial markers during the experiment [26,
207, 208, 209, 210, 211, 212], image cross-correlation between frames at
different times [213, 214, 215, 216, 217, 218, 219, 220], computing drift
directly from the list of found emitter coordinates [221, 222, 223, 224, 225]
and other procedures [226, 227, 228]. In what follows, we briefly introduce a
few of these drift correction methods.

Fiducial markers are fixed point sources of light or structures during data
acquisition that are used as reference points for drift correction. Since
fiducial markers are fixed, their movements can be measured in the
localization step and used to eliminate drift errors. Fiducial markers can
introduce light corruption to the sample and some fiducial markers emit
photons for a limited time, which restricts the data acquisition time. Cross-
correlation approaches are the most common procedure for removal of drift
error and different variants have been reported in the literature. Image noise
deteriorates performance of the cross-correlation approaches to calculate
drift errors. One algorithm maximizes the cross-correlation of the first frame
with the rest of the sequence to calculate drift [215, 218]. A fast
implementation of cross-correlation is achieved in the Fourier domain [216,
217, 220], and some other approaches employ the cross-correlation of the sum
images to reduce the effect of noise [214, 219]. Some approaches have been
developed using the list of found localizations to estimate sample drifts,
such as the NND distribution of locations [223, 224, 225].

Wester, et al. [224] made use of a combination of the image registration
approach and NND distribution of localizations to measure drift in the sample.
The algorithm employs periodic 3D registration of the sample using brightfield
images to remove drift errors. However, brightfield registration is only
accurate to around $10\,\mathrm{nm}$ and there might be still residual drift
remaining. This approach thus uses post-processing of the localizations to
extract residual drift from the NND distribution of emitter coordinates. This
procedure is robust and also capable of calculating axial drift.

In a recent work [225], drift correction is performed using an entropy
minimization scheme. The SMLM reconstruction is considered to be a
probablility distribution about the true emitter locations, and an entropy
metric expressing the uncertainty in the position estimation is used as an
optimization cost. To simplify the calculation, the cost is replaced by a
closed form upper bound and points beyond a fixed distance threshold in an
estimated covariance matrix are excluded. The drift model is taken to be a
piecewise cubic spline over a series of dataset segments. This algorithm
outperforms the cross-correlation procedures.

###  4.6 Fitting Quality and Image Quality

The quality of a fit is defined as how close the recovered location of an
emitter is to its true location. The fitting quality depends on the emitter
detection, emitter fitting and thresholding stages. Each of these steps were
discussed in previous sections. There are a few metrics to assess the fitting
quality including the Jaccard index (JAC), root mean square errors (RMSE) or
accuracy, precision [49, 50], a vectorial model for localization uncertainties
[229], and a fit confidence test that examines the stability of a localization
by exploring the model in the vicinity of the maximum likelihood [230].

Image resolution is often used as a reporter of image quality. However, there
is no widely accepted definition for resolution in super-resolution microscopy
literature [59, 124, 231]. There are however a few common metrics that are
used to assess microscopy image resolution such as localization precision
[59], SNR, Fourier ring correlation (FRC) [45], and parameter-free resolution
estimate [232]. In what follows, we discuss some of these metrics.

#### Fitting Quality

The uncertainty in the location of an isolated emitter is characterized by
precision and accuracy. Assuming that an emitter with location $x_{0}$ blinks
multiple times, there are then multiple corresponding found locations,
$x_{f}$, for the emitter. For an optimal fitting algorithm and under ideal
conditions, the found locations are still slightly different and spread over a
small region around the emitter location, $x_{0}$. The localization precision,
$\sigma_{x}$, is the standard deviation of the distribution of $x_{f\,}$. This
is a fundamental limit for localization precision imposed by the random nature
of photons, and is not a result of instrument imperfections or flaws in the
experiment design [93, 81, 59]. Due to the stochastic nature of photons, each
time that the emitter turns on, the number of photons reaching a certain pixel
deviates by a small amount, which is called shot noise, resulting in small
variations in the found emitter locations, Fig. 8. The fundamental limit of
precision for an optimal fitting algorithm is given by the square root of the
inverse of the diagonal elements of the Fisher information matrix [59]

$\sigma_{\theta_{i}}=\frac{1}{\sqrt{I\int_{-\infty}^{+\infty}\frac{1}{L(D|\theta)}\left(\frac{\partial
L(D|\theta)}{\partial\theta_{i}}\right)^{2}d\theta}}$ (31)

Figure 8: Precision and accuracy concepts. Blue dots are the true emitter
position and the orange dots stand for found locations. (a) Distribution of
the found locations from an optimal fitting algorithm. Precision is
proportional to the size of the distribution of the found locations. (b)
Distribution of the found locations from a biased fitting algorithm. Accuracy
is described as the deviation of the center of the distribution of found
locations from the true position.

where $L(D|\theta)$ is the likelihood given in (14), $\theta,\theta_{i}$ and
$I$ stand for the collection of all the parameters, the $i$th parameter, and
the number of photons, respectively. The fundamental limit of localization
precision for lateral coordinates due to shot noise is given by eq. (4).

Accuracy describes the deviation of the average of the found locations from
the true emitter location, Fig. 8. Unlike localization precision, localization
accuracy does not have a fundamental limit and can be zero [59]. However, due
to instrument imperfections, optical aberrations and noise, the resultant
localizations are often a little biased [233, 234, 235, 236]. A set of
acceptable localizations of an emitter have similar spreads of precisions and
accuracies.

In practice, it is often challenging to achieve unbiased localizations and
localization precisions close to the square root of the CRLB variance derived
from a too simplistic model that only takes into account Poisson noise. This
is due to other existing factors in super-resolution experiments such as PSF
mismatch, equipment (for instance cameras), and background noise. In what
follows, we explain each of these effects:

PSF mismatch. Error-free PSF input to the iterative approaches is of key
importance for accurate and precise emitter localization [64]. There are
several factors that lead to major deviations between the expected PSF and the
resultant PSF, including dipole orientation, emitter motions, and changes in
the optical properties of the environment like the refractive index [48, 59,
67, 73, 236, 237]. Fluorophore emitters are electric dipoles in nature and do
not emit photons in an isotropic manner, despite the assumption of spatial
uniform distribution of photons. A rigid emitter therefore has a nonuniform
photon distribution which degrades the fitting quality when using an isotropic
approximation of the PSF [63, 64, 238]. However, these emitters are usually
free to randomly rotate, which provides an isotropic distribution of photons.
Another cause of PSF mismatch is emitter movements that exist in some
experiments due to the nature of the sample, like live cell imaging or
particle tracking, or due to bad fixation of the sample. This can result in a
smearing pattern in the PSF known as the motion blur effect, specifically for
large exposure times and fast diffusions, when the emitter emits photons while
moving [65, 66, 239]. Furthermore, imaging cellular structures within their
native environment often requires light traveling through thick tissues of the
specimen. The optical properties of the sample, such as the refractive index,
changes as a function of location which leads to a disturbance in the light
wavefront and therefore distortion of the PSF [236, 237, 240, 241].

Camera. Several properties of cameras, such as noise introduced in the
detection process and pixelation, affect the localization precision and
accuracy. Two major detection noises introduced in cameras are read-out noise
and the noise added in the amplification process that can affect localization
uncertainty, see section 3. Model (13) does not take into account the
detection noises, but can be modified to account for these noises by
subtracting the read-out offset from the pixel values (input data) and
rescaling them by the amplification factor in EMCCD cameras, where the read-
out offset and amplification are uniform across the field of view. However,
the detection noises require a more careful consideration in sCMOS cameras
where these noises are unique for each pixel [242, 84]. The pixel dependent
noise in sCMOC cameras gives rise to a complicated and computationally
expensive noise model. As such, in an effort to reduce the computational
complexity, Huang et al. showed that this model can be approximately by a
Poisson distribution with an additional source of photons [84, 243].

Scientific cameras are not able to record the exact location of photons
reaching them. However, they report the number of photons reached on a region
within a certain area called a pixel. This causes loss of accurate locations
of photons from an isolated emitter, and therefore results in localization
uncertainty. This effect is called pixelation and is often negligible for
relatively high SNR [93, 85]. Moreover, sensitivity of cameras can be
nonuniform over the field of view, which is another potential source of
deterioration of localization precision [59].

Background. Another limiting factor for precision and accuracy in localization
microscopy is the background noise. Failure in correct background estimation
can lead to biased localization estimates as was discussed in section 4.1. The
homogeneous background is usually modeled by an additive offset to the
intensity described in eq. (12). However, the situation is worse in the
presence of structured background, which is more difficult to model and can
give rise to more uncertainty in localizations.

We have so far discussed fitting quality for isolated emitters. However,
overlapping PSFs from multiple active emitters in a diffraction limited area
are unavoidable in SMLM. JAC is therefore introduced as a standard metric to
measure the quality of fitting for a group of emitters in a dense region with
overlapping PSFs [49, 50, 61]. JAC evaluates the rate of emitter detection
versus the density of emitters. This is particularly important for multiple-
emitter fitting algorithms which are supposed to detect and localize emitters
in dense regions of data. To perform this test, sets of SMLM data with
different emitter densities are simulated and processed with the given fitting
algorithm. Next, the pairs of matching emitters within the sets of found and
true locations are discovered by minimizing the cost metric between the two
sets. JAC is then calculated as the ratio of the number of matched emitters
(ME) to the total number of localizations within both sets

$\mathrm{JAC}=\frac{\mathrm{ME}}{\mathrm{FE+TE}}$ (32)

where $\mathrm{FE}$ and $\mathrm{TE}$ refer to the number of found emitters
and true emitters, respectively. Furthermore, Ram, et al. provided an approach
to assess the quality of fitting for two nearby emitters employing Fisher
information theory [244].

Finally, theoretical inspection of fitting quality considering all the
important factors is a difficult task, however statistical and experimental
techniques have been proposed to calculate localization uncertainties taking
into account all those factors [105, 245, 246]. An approach has been proposed
that prescribes a simple way to compute localization precision from NND
distribution of found localizations [105]. Furthermore, Lee, et al. used well
isolated emitters in sparsely labeled biological samples to experimentally
measure localization precision under realistic experimental conditions [245].

#### Image Quality

After conducting an experiment and performing the image analysis, the final
result is a super-resolved image reconstructed from the list of found
localizations. What a scientist eventually cares about is the
quality/resolution of this image and the amount of details that it reveals.
The quality of the list of fits used to reconstruct an image directly affects
the image resolution. For instance, the localization precisions must be
smaller than a desired resolution. However, although the fitting quality is an
essential factor in assessing the resolution of SMLM reconstructed images, it
does not provide a comprehensive assessment of resolution for images of
complex biological structures. For instance, an important factor in evaluating
the resolution of reconstructed images of such complex structures is labeling
density [103, 247]. Furthermore, resolution highly depends on the type of
structure being examined and differs from experiment to experiment [59, 124,
248]. For example, visualization of fine actin filaments requires uniform and
high labeling densities [43, 249], while in imaging nuclear pore complexes
(NPCs), the localization precision is of more importance since the sample
itself has a discrete structure [250]. To obtain a desired resolution, the
labeling density should result in a NND of at least half of the desired
resolution between the localized emitters according to the Nyquist-Shannon
sampling theorem [59, 251, 252, 253, 254]. Therefore, more advanced
theoretical and experimental methods for evaluation of resolution in
localization microscopy have been devised. Here, we first provide a concise
discussion on the labeling issue and then introduce some of the theoretical
and experimental metrics for assessing the resolution.

Low labeling density can be due to several problems in super-resolution
microscopy, such as insufficient labeling, or broken labels or filtering
localizations in the thresholding stage because of inadequate fitting quality.
For some cases, such as experiments where there are multiple identical
structures, this problem can be overcome by aligning and fusing the images of
similar structures [255, 256, 257]. Several alignment and fusion algorithms
have been reported in the literature, including a few that require a template
or structural assumptions in the alignment step [258, 259, 260, 261, 262] as
well as template-free particle fusion algorithms [106, 263, 264, 265].

A range of theoretical methods for evaluation of resolution in SMLM
reconstructed images have been developed such as the information transfer
function (ITF) [266], Fourier ring correlation (FRC) [45, 267], super-
resolution quantitative image rating and reporting of error locations
(SQUIRREL) [268], parameter-free resolution estimation [232], and others [269,
270]. In the following, we describe a few of these methods.

ITF. Mukamel, et al. introduced the notion of ITF that accounts for photon
statistics in localization microscopy as well as the resolution limit due to
objective lenses. ITF is a generalization of the concept of MTF (see the
Introduction) and reduces to the square of the MTF for conventional microscopy
[266]. The ITF formalism gives image resolution bounds by providing a minimum
on precision for image frequency estimates

$\Delta_{I(k)}\leq F^{-1}(k)$ (33)

where $\Delta_{I},k$ and $F$ stand for uncertainties, spatial frequency, and
the ITF, respectively. This idea of resolution equally applies to both
conventional and stochastic microscopy and unifies the concept of resolution
[266]. Although ITF provides a generalized concept of resolution, it requires
models of the target structure which makes it limited in practice [59].

SQUIRREL. SQUIRREL builds an intensity image by convolving an input super-
resolution image with a resolution scaling function. It then compares the
obtained intensity image with the wide-field counterpart of the input super-
resolution image to identify artifacts and disappearance of details [268].
SQUIRREL does not require any prior knowledge of the sample and is able to
identify common super-resolution artifacts. However, out-of-focus light
affects the performance of SQUIRREL and moreover, it cannot recognize small-
scale artifacts. This approach returns a quantitative map of image anomalies
and artifacts rather than reporting a resolution measure. The returned map can
be used to optimize the experiment and image analysis, and improve the
reconstructed image qualities. SQUIRREL was used to show that high image
resolution does not necessarily correlate with artifact-free reconstructions
[268].

FRC. The FRC algorithm was adopted from electron microscopy [271] and adapted
for resolution measurement in stochastic fluorescent microscopy [45, 267].
This approach uses two images from two identical structures or reconstruction
of two images of the same structure by splitting the list of localizations
into two halves. It then calculates the cross-correlations of the two images
at different image frequencies, where components with frequency $k$ are given
by a ring with radius $|k|$ in the Fourier domain. The higher frequencies
represent finer details in the image. It is there the two images start to
deviate and hence their cross-correlation rapidly declines. The frequency at
which the cross-correlation falls below a given threshold can be utilized to
calculate the image resolution. The method takes into account both
localization precision and labeling density. It however has to be used with
care for structures with more complicated geometries and discontinuous
boundaries [59, 124].

Parameter-free resolution estimate. Parameter-free resolution estimation [232]
is based on the idea of correlation in the frequency domain similar to FRC. It
however only uses a single image and furthermore does not require any
threshold. The core idea is that structure information is mostly preserved in
lower frequencies while high frequencies mostly represent noise. Therefore
applying a frequency mask with a decreasing size will remove the contribution
of noise and the correlation will gradually increase to a maximum in a certain
frequency cutoff. After this frequency, the mask eliminates the structure
information and the correlations falls. The frequency at which the maximum
occurs is then employed to calculate the resolution [232].

Now we turn to experimental approaches for characterization of resolution
reported in the literature. For instance, DNA-rulers, DNA-origami structures,
or standard biological structures such as NPCs with known spacing between
emitters, have been employed as references to benchmark the ability of
microscopy techniques in resolving closely spaced emitters [250, 272, 273].
Some other well characterized biological structures such as microtubules and
synthesized SMLM data for different structures, such as the Siemens star and
crossing lines, have also been employed to inspect the resolution of SMLM
reconstructed images [49, 50, 61, 122, 274].

Finally, resolution depends on multiple experimental and image analysis
factors. A proper choice of labels, buffers and fixation protocols can lead to
better localization uncertainties, less artifacts and therefore better image
resolution. The choice of a localization algorithm also helps in optimizing
both localization precision and localization density, reducing the image
artifacts. The image rendering and visualization techniques also influence the
image resolution. It has been demonstrated that different visualization
procedures yield different resolutions for images reconstructed from the same
list of emitter coordinates [275].

##  5 Applications of Single Molecule Localization Microscopy

Super-resolution microscopy has had significant contributions in shedding
light on numerous biological problems since its advent in the 1990s. The major
improvement in resolution of microscopy images helped with unraveling details
of many biological structures and processes, such as nuclear pore complexes
[250, 264, 276], actin structures [43], cell membrane [277, 278], cell
division [279], neurons [280, 281], plant cell biology [282, 283], live cell
imaging [253, 284, 285], and many more [14, 17, 18, 19, 286, 287].

In addition to visualization of intricate structures, a tremendous
quantitative insight into cell biology can be gained through post-processing
of resulting localizations from SMLM techniques. These post-processing
algorithms include single particle tracking (SPT) [125, 288, 289], cluster
analysis [290, 291], molecule counting [292, 293, 294, 295, 296, 297, 298,
299, 300] and particle fusion [106, 255, 256, 264, 262, 257, 258, 260, 261,
259, 263]. In the following, SPT and clustering approaches are described in
more detail.

###  5.1 Single Particle Tracking

While super-resolution microscopy only deals with spatial resolution, Single
Particle Tracking (SPT) is concerned with spatial as well as temporal
resolution of dynamic samples. The microscopic world of biology is a live and
vibrant place where many phenomena and interactive processes take place over
time. An effective examination of these types of processes and interactions
requires a comprehensive knowledge of dynamics and the motion of particles in
those environments.

SPT was first demonstrated by optical visualization of gold nano-particles on
cell membranes [301, 302], which paved the path for use of SPT in the study of
cell membranes [303, 304]. Since the first demonstration of SPT in cell
biology, there have been a burst of research using SPT that have
revolutionized our understanding of the cell membrane [305, 306, 307, 308]. An
example of such research is the use of SPT in shedding light on the confined
motion of receptors on the cell membrane imposed by the actin cytoskeleton
[306, 309]. SPT has also been employed to explore intracellular environments,
including influenza and HIV virus internalization [310, 311, 312], and cargo
transport along microtubules [313].

SPT is comprised of three main steps: 1) detection and localization; 2)
linking; and 3) post processing. The detection and localization of probes in
SPT is similar to emitter localization in SMLM techniques and most of the
fitting approaches discussed before can be employed for this stage of SPT. As
described previously, more photons from the probes lead to more precise
localizations. On the other hand, to increase the temporal resolution, a
higher frame rate is desired during data acquisition. Higher frame rate means
shorter exposure time and less photons and hence larger localization
uncertainties. Therefore, bright probes are necessary for both high temporal
resolution and precise localizations. Other important elements are photo-
stability and size of the probes [314]. Probes with long-lasting photo-
stability and small sizes are desired for obtaining longer trajectories and
less perturbation in the mobility and dynamics of the particles under study.
There are three major categories of probes used in SPT: nano-particles,
organic dyes and quantum dots (QDs) [289]. Nano-particles reflect light rather
than emitting fluorescent light and have high photo-stability, but their size
is rather large compared to typical proteins. Organic dyes are very small
compared to typical proteins, but they are really dim and bleach rapidly
[289]. The size of QDs are in-between organic dyes and nano-particles. They
are relatively bright and photo-stable, which make them a perfect probe for
most SPT experiments. QDs are manufactured with different emission wavelengths
and are also suitable for multi-color SPT [315, 316, 317].

After detection and localization, the detected probes are linked together
across consecutive frames to construct particle trajectories. Linking is a
non-trivial problem, even in the presence of only a few probes, Fig. 9. Yet
this is the case where there is an equal number of probes in every given
frame, in the absence of blinking, bleaching, probes entering or leaving the
field of view, and missed probes in the detection stage. The reported linking
algorithms in SPT can be categorized in two classes: 1) deterministic; 2)
probabilistic; and 3) deep learning approaches [318, 319, 320, 321]. In what
follows, we only describe deterministic and probabilistic techniques.
Deterministic approaches give the same results every time, but probabilistic
procedures use a random number generator at some point in the analysis and
depending on the sequence of random numbers produced, they give slightly
different answers each time.

Given a set of localizations, deterministic approaches construct trajectories
by linking the localizations so that the resulting trajectories minimize a
certain cost function [322]. The reported cost functions include distance
[323, 324, 325, 326, 327], as well as concepts such as minimal energy path
[328, 329, 330] or inertia [331] borrowed from classical mechanics, and ad hoc
cost functions that can be a function of parameters like probe location,
intensity, and motion direction [332, 333, 334, 335]. A range of different
optimization algorithms have been employed to minimize the aforementioned cost
functions, including the Hungarian [326, 333, 336] and the greedy [331, 332,
335] methods.

Probabilistic SPT approaches make use of a Bayesian algorithm to obtain the
most probable tracks in the linking problem. SPT algorithms that implement a
Bayesian procedure have demonstrated better performances compared to
deterministic approaches due to inclusion of prior information [337].

Figure 9: Three possible links in the presence of two probes. The Red circles
and green curves, respectively, represent detected probes and possible
trajectories.

Here, we briefly describe the probabilistic approach. Assume a system with
state $x_{t}$ at time $t$, which leads to an observation $z_{t}$, where
$x_{t}\neq z_{t}$ due to noise. The state at time $t$ is a stochastic quantity
related to the state at time $t-1$ as follows

$x_{t}\sim f_{t}(x_{t-1},u_{t})$ (34)

where $f_{t}$ and $u_{t}$ are, respectively, a function that characterizes the
motion type and a variance. The observation also deviates from the state of
the system due to noise and is given by

$z_{t}\sim g(x_{t},v_{t})$ (35)

where $v_{t}$ is the noise variance. The objective is finding the state of the
system at time $t$, $x_{t}$, given the list of observations up to this time,
$z_{1:t}$. Using Bayes theorem, it can be shown that [338]

$P(x_{t}|z_{1:t})=\frac{P(z_{t}|x_{t})P(x_{t}|z_{1:t-1})}{P(z_{t}|z_{1:t-1})}$
(36)

where $P(z_{t}|x_{t})$ is the likelihood given by the noise model (35), and
$P(x_{t}|z_{1:t})$ is the posterior. The prior is given by

$P(x_{t}|z_{1:t-1})=\int P(x_{t}|x_{t-1})P(x_{t-1}|z_{1:t-1})dx_{t-1}$ (37)

where $P(x_{t}|x_{t-1})$ is defined by eq. (34). Equations (36) and (37)
provide a recursive formula to obtain the current state of the system from the
list of observations up to this time along with the previous state of the
system. However, these equations do not have an analytical solution for a
general case. However, an analytical solution can be derived if the posterior
has a Gaussian form at any time point, $f_{t}$ is a linear function and $g$
describes Gaussian noise, known as the Kalman filter [339, 340, 341, 342, 343,
344, 345]. For instance, diffusion is a linear process and can be
characterized by a Kalman filter. Another common filter is the particle filter
[99, 338, 346, 347, 348, 349, 350], which is suitable for non-linear processes
with non-Gaussian noise.

The recursive method described above provides an approach that links new
probes to the trajectories on a frame by frame basis. However, the Multiple
Hypothesis Tracking (MHT) algorithm evaluates the possibility of different
tracks, taking into account localized probes over either the entire data set
or over a time window [351, 352, 353, 354, 355, 356]. While this procedure
accomplishes the most optimal solution to the linking problem, it is
computationally extremely expensive and impractical. Greedy algorithms
accomplish approximations to the MHT optimal solution with tremendously less
computational cost [353, 357].

Jaqaman et al. [353] proposed a relatively computationally inexpensive
approximation to MHT using a greedy algorithm. This procedure uses a cost
matrix approach to incorporate appearance (birth) and disappearance (death) of
the probes in the model. They introduced ghost probes to consider three
different linking scenarios when there is not an equal number of detected
probes in consecutive frames, Fig. 10. The real probes are connected together
using a cost calculated from the localization precisions and kinetics of the
probes

$P(x_{2},t_{2}|x_{1},t_{1})=\mathcal{N}(x_{1},x_{2},\sigma^{2})$ (38)

where $x$ stands for the probe locations and $\sigma^{2}$ is the variance due
to diffusion and localization precisions, $\sigma_{x}$, defined as

$\sigma^{2}=2D\Delta t+\sum_{i=1}^{2}\sigma_{xi}^{2}$ (39)

where $D$ and $\Delta t$ are the diffusion constant and exposure time, in
turn. A real probe can be connected to a ghost probe when a probe goes either
off or on in a subsequent frame and the penalty is computed from the blinking
statistics of the probes. The blinking process is assumed to be memoryless
with $k_{\mathrm{on}}$ and $k_{\mathrm{off}}$ being the rate of the probes,
respectively, going from off to on, and on to off

$\begin{split}P(\mathrm{\mathrm{off}}\rightarrow\mathrm{on}|k_{\mathrm{on}},\Delta
t)=k_{\mathrm{on}}\Delta t\\\
P(\mathrm{on}\rightarrow\mathrm{off}|k_{\mathrm{off}},\Delta
t)=k_{\mathrm{off}}\Delta t.\end{split}$ (40)

Additionally, the probabilities for a probe staying on or off for $N$ frames
is given as

$\begin{split}P(\mathrm{Stay\,off}|N\Delta
t)=P(\mathrm{Not\,going\,on}|N\Delta t)=(1-k_{\mathrm{on}}\Delta t)^{N}\\\
P(\mathrm{Stay\,on}|N\Delta t)=P(\mathrm{Not\,going\,off}|N\Delta
t)=(1-k_{\mathrm{off}}\Delta t)^{N}.\end{split}$ (41)

Eventually, ghost probes can also be linked together where the costs are given
by the transpose cost matrix of the real probes. The Hungarian procedure is
utilized to find the connections with the least cost.

The first stage only connects probes across consecutive frames and produces
short trajectories, Fig. 10. In a second stage, the resulting short
trajectories are connected together to obtain probe trajectories across the
entire given sequence of the frames, which is called gap-closing. The gap-
closing cost function is also a function of blinking statistics and
localization precisions, except this time the probes are allowed to be linked
together in non-consecutive frames.

Figure 10: An instance of probe linking in two consecutive frames in the
presence of ghost probes. The blank and filled circles represent ghost and
real probes.

Once the trajectories are constructed, they can be employed to extract
information about the underlying dynamics of the biological sample under study
[289, 358, 359, 239, 288]. Mean Square Displacement (MSD) is the most common
post-analysis quantity of trajectories employed for this purpose. It is
defined as the average of spatial distances between all pairs of localizations
that are temporally separated by certain times, called the time lag. Given a
trajectory comprised of $N$ localizations, the MSD can be calculated as
follows:

$\mathrm{MSD}(t_{\mathrm{lag}}=m\Delta
t)=\frac{1}{N-m}\sum_{i=1}^{N-m}\left[x(t_{i}+m\Delta t)-x(t_{i})\right]^{2}$
(42)

where $x(t_{i})$ is the location of the probe at time $t_{i}$ and
$x(t_{i}+m\Delta t)$ is the probe location at $t_{\mathrm{lag}}=m\Delta t$
later. The MSD is a function of the time lag, $t_{\mathrm{lag}}$, and this
function has been extensively investigated and well characterized for
different motion types. By comparing the resulting MSD function from the
obtained trajectories and the theoretical MSD forms, one is able to determine
the underlying dynamic of the system under study. For diffusion or Brownian
motion, the MSD is proportional to the time lag

$\mathrm{MSD}(t_{\mathrm{lag}})=2dDt_{\mathrm{lag}}$ (43)

where $d$ is the spatial dimensionality of the motion.

###  5.2 Clustering

Clustering is the task of classifying a given set of data points into subsets
in such a way that data points in a subset are more similar to each other than
those in other subsets, based on a given property [360]. Clustering has been
employed in many scientific spheres, including localization microscopy,
recognizing patterns and detecting communities.

In SMLM, there are different levels of clustering for a list of localizations
resulting from an experiment [290, 291]. In SMLM, there are multiple observed
localizations from an emitter over the course of data acquisition that form a
cluster, which we will call here an intra-cluster. A higher level of
clustering is clustering between emitters/molecules that signifies
interactions between the molecules in a cluster at the biological level. We
call this type of cluster inter-cluster. There is also co-clustering that
takes place among different types of molecules imaged via multi-color super-
resolution microscopy [361, 362, 363, 364].

Multiple algorithms have been developed or adapted to identify all sorts of
clusters within SMLM localizations, including algorithms that inspect the
extent of clustering within a data set [365, 366, 367, 368, 369, 370, 371],
density based clustering algorithms [371, 372, 373, 374, 375, 376] and
Bayesian clustering approaches [147, 290, 371, 377, 378, 379]. These
approaches have been popular in post-processing of SMLM results of cell
membrane samples, where it is conjectured that molecular clusters indicate
underlying membrane structure or cell signaling [380, 381, 382, 383, 384,
385].

Some methods examine the degree of grouping or group properties of clusters,
such as the shape or average density within a data set, Fig. 11. A few of
these approaches are NNDs [365], Hopkins statistics [366], and Ripley’s
function and pair correlation [367, 368, 369, 370]. The NND probability
distribution for uniform randomly distributed data is given by

$f(r)=2\pi\rho re^{-\pi\rho r^{2}}$ (44)

where $\rho$ is the density of the data points. The deviation of the NND
distribution for a given data set from (44) indicates either clustering or a
more regularly spaced data set, Fig. 11. Hopkins statistics ($H$) tests for
spatial randomness of a point pattern by comparing NNDs between a given data
set and a uniform randomly distributed data set. Values of $H$ near $0.5$
imply randomly distributed data, while values near 1 indicate highly clustered
data, and values near zero signify more regularly spaced data, Fig 11.

The Ripley’s K-function, $K(r)$, is another statistical test that explores the
extent of clustering within a given data set. Defining $k_{i}(r)$ as the
number of points within distance $r$ from the $i$th point, $K(r)$ is given by
the average of the $k_{i}$ for all data points divided by the mean density
over all data points. The clustering behavior of the given data set is then
given by

$\begin{split}K(r)<\pi r^{2}\hskip 28.45274pt\mathrm{regularly\,spaced}\hskip
8.53581pt\\\ K(r)=\pi r^{2}\hskip 28.45274pt\mathrm{random}\hskip 43.2483pt\\\
K(r)>\pi r^{2}\hskip 28.45274pt\mathrm{clustered}\hskip 38.41139pt\end{split}$
(45)

Pair correlation analysis is the same idea as Ripley’s K-function extended to
images via pixels rather than a set of points.

A range of clustering algorithms have been developed that take advantage of
density fluctuations within a data set to identify clusters. DBSCAN [372] and
Voronoi tessellation [375, 376] algorithms are among the most popular density
based approaches in SMLM data post-analysis. DBSCAN has two inputs. One of
them is the maximum distance between the points within a cluster and the other
is the minimum number of data points within a cluster. The algorithm starts
from a random point and adds points to the set until there is no more points
within the given maximum distance. If the number of points inside the set is
more than the given threshold, it will be returned as a cluster, otherwise it
is considered as an outlier. Voronoi tessellation draws boundary line segments
between every pair of data points in the given data set so that the distance
from every member of a pair to the line between them is equal. The final
geometry is comprised of polygons where each polygon contains only one data
point, called a Voronoi diagram. This algorithm identifies clusters by
inspecting the number of vertices and the area of the polygons at different
neighboring levels.

Bayesian clustering algorithms for SMLM data are superior to density based
approaches because they take advantage of prior knowledge. Available prior
information can be included in Bayesian algorithms via prior distributions.
Rubin-Delanchy, et al. developed an algorithm to identify inter-clusters
within SMLM localizations, called Bayesian cluster identification [290, 377].
This algorithm classifies the localizations into either background or signal
by means of a binomial distribution. A large number of clustering proposals
are then put forward using a procedure similar to Ripley’s K-function. Scores
for different clustering proposals are calculated using the posterior, where
the Dirichlet distribution is the prior on different arrangements. At the end,
the configuration with the highest score is returned as the solution.

Figure 11: NND distribution and Hopkins statistics. First and second rows
represent the NND distribution and Hopkins statistics, respectively. (a)
randomly distributed, (b) clustered, (c) more regularly spaced data.

Post, et al. proposed a Bayesian clustering algorithm to identify intra-
clusters (emitters) on the surface of spherical nano-particles [378]. Their
algorithm takes full advantage of localization precisions, but it does not
return the number of emitters on the surface of every individual nano-
particle. Otherwise, they find the mean number of emitters on the surface of
nano-particles by fitting a log-normal distribution to the number of
localizations per nano-particle. They next employ this piece of information to
localize those emitters and then inspect their inter-clustering behavior using
the NND distribution.

Fazel et al. developed Bayesian Grouping of Localizations (BaGoL) that takes
advantage of localization precisions as well as mean number of
blinking/binding events per emitter to identify and localize emitters with
sub-nanometer precision in dense regions [379]. This algorithm uses NND
distributions and an intensity filter to recognize and remove outliers.

Lin, et al. developed a non-Bayesian intra-clustering algorithm based on a
hierarchical clustering paradigm, which allows exploitation of localization
precisions [383]. This algorithm proposes different clustering configurations
for a given set of localizations similar to the hierarchical procedure, and
next performs a hypothesis test to pick one of the proposed clustering models.

## Conclusion and Future Outlook

Two decades after the advent of super-resolution microscopy, these techniques
discussed above have become standard tools in direct and non-invasive
examination of sub-cellular environments with unprecedented resolution. Among
different super-resolution techniques, the SMLM approach involves data
processing and post-processing stages that are reviewed in this survey. Yet,
despite all the progress in this area, the existing SMLM data analysis
techniques face multiple challenges. These include the existence of sample
induced optical aberrations, and the lack of common metrics in assessing the
robustness of image analysis and image resolution. Furthermore, extracting
meaningful biological information from the resulting list of SMLM emitter
coordinates are complex and difficult tasks. Therefore, the area of SMLM data
post-processing is a fast growing field opening many promising venues of
research in this area which, in turn, demands new computational tools.

The 3D nature of biological structures has been a long time motivation for
developing 3D SMLM microscopy techniques. Particularly, it has been of
interest to observe subcellular structures within their native environments.
However, fluorescent light traveling within such complex biological
environments is significantly perturbed, degrading the quality of the final
SMLM images. Adaptive optics [236] is an established experimental technique to
compensate for such sample induced aberrations, but the potential of
computational methods to deal with this major problem is still under-explored
[241]. Another main challenge in 3D SMLM is visualization of 3D structures,
which demands the development of more sophisticated visualization tool-boxes
able to robustly and rapidly handle millions of localizations.

Another major challenge in SMLM data analysis is assessment of fitting and
image qualities. In spite of numerous reported metrics and methods in the
literature (covered in section 4.6), there is still no standardized universal
set of metrics for comparing analysis results. For example, there can be
disagreement between different existing metrics in evaluating the robustness
of a fitting algorithm or the resolution of an image [268].

An obstacle in obtaining meaning from SMLM and SPT data is the development of
multi-color and multi-target analysis techniques that can identify spatial
relationships between multiple species of proteins, and techniques for
inferring protein-protein interactions [36, 386, 387, 287]. The main issues in
this area are image registration and colocalization between different
channels. These are typically performed using different correlation techniques
[388, 389], but such methods are often hindered by noise [390, 391]. The use
of machine learning techniques seems to be a promising solution [392].

Solutions to the current SMLM data analysis challenges discussed here as well
as many others might be found in leveraging tools from cutting-edge machine
learning methods such as Bayesian techniques and neural networks. However,
each of these paradigms comes with its own benefits and shortcomings. For
example, Bayesian techniques provide robust analysis tools for data with low
SNR, and takes into account different existing sources of uncertainty.
Nonetheless, Bayesian techniques are often slow and require multiple input
parameters. On the other hand, deep learning techniques are fast and
parameter-free, but they do not provide any uncertainty estimate over the
learned quantities and require careful training. A promising emerging data
analysis avenue is Bayesian neural networks which addresses some of the
mentioned issues by combining Bayesian techniques and neural networks.
Finally, an important guideline in developing SMLM data analysis methods is
the accessibility of these advanced machine learning tools to non-expert users
by developing user friendly image processing packages and providing thorough
algorithmic and implementation oriented (user) documentation.

## Acknowledgement

We thank Keith A. Lidke and Steve Pressé for useful comments.

## Data Availability

Data sharing is not applicable to this article as no new data were created or
analyzed in this work.

## Conflict of Interest

Authors declare no conflict of interests.

## References

  * [1] Malte Renz.  Fluorescence microscopy – a historical and technical perspective.  Cytometry Part A, 83(9):767–779, 2013.
  * [2] Max Born and Emil Wolf.  Principles of optics: electromagnetic theory of propagation, interference and diffraction of light.  Elsevier, 2013.
  * [3] Joseph W Goodman.  Introduction to fourier optics. 3rd.  Roberts and Company Publishers, 2005.
  * [4] Sjoerd Stallinga, Antoine Delon, and Jörg Enderlein.  Optimal transfer functions for bandwidth-limited imaging.  arXiv preprint arXiv:2110.04839, 2021.
  * [5] James Pawley.  Handbook of biological confocal microscopy.  Springer Science & Business Media, 2010.
  * [6] Minsky Marvin.  Microscopy apparatus, December 19 1961.  US Patent 3,013,467.
  * [7] Daniel Axelrod.  Cell-substrate contacts illuminated by total internal reflection fluorescence.  The Journal of cell biology, 89(1):141–145, 1981.
  * [8] Winfried Denk, James H Strickler, and Watt W Webb.  Two-photon laser scanning fluorescence microscopy.  Science, 248(4951):73–76, 1990.
  * [9] Stefan Hell and Ernst HK Stelzer.  Fundamental improvement of resolution with a 4pi-confocal fluorescence microscope using two-photon excitation.  Optics Communications, 93(5-6):277–282, 1992.
  * [10] Brent Bailey, Daniel L Farkas, D Lansing Taylor, and Frederick Lanni.  Enhancement of axial resolution in fluorescence microscopy by standing-wave excitation.  Nature, 366(6450):44, 1993.
  * [11] Arne H Voie, DH Burns, and FA Spelman.  Orthogonal-plane fluorescence optical sectioning: three-dimensional imaging of macroscopic biological specimens.  Journal of microscopy, 170(3):229–236, 1993.
  * [12] Stefan W Hell and Jan Wichmann.  Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy.  Optics letters, 19(11):780–782, 1994.
  * [13] Bo Huang, Mark Bates, and Xiaowei Zhuang.  Super-resolution fluorescence microscopy.  Annual review of biochemistry, 78:993–1016, 2009.
  * [14] Eugenio F Fornasiero and Felipe Opazo.  Super-resolution imaging for cell biologists: Concepts, applications, current challenges and developments.  Bioessays, 37(4):436–451, 2015.
  * [15] Christian Eggeling, Katrin I Willig, Steffen J Sahl, and Stefan W Hell.  Lens-based fluorescence nanoscopy.  Quarterly reviews of biophysics, 48(2):178–243, 2015.
  * [16] Bartosz Turkowyd, David Virant, and Ulrike Endesfelder.  From single molecules to life: microscopy at the nanoscale.  Analytical and bioanalytical chemistry, 408(25):6885–6911, 2016\.
  * [17] Steffen J Sahl, Stefan W Hell, and Stefan Jakobs.  Fluorescence nanoscopy in cell biology.  Nature reviews Molecular cell biology, 18(11):685, 2017.
  * [18] Lothar Schermelleh, Alexia Ferrand, Thomas Huser, Christian Eggeling, Markus Sauer, Oliver Biehlmaier, and Gregor PC Drummen.  Super-resolution microscopy demystified.  Nature cell biology, 21(1):72, 2019.
  * [19] Jeroen Vangindertael, R Camacho, W Sempels, H Mizuno, P Dedecker, and KPF Janssen.  An introduction to optical super-resolution microscopy for the adventurous biologist.  Methods and applications in fluorescence, 6(2):022003, 2018.
  * [20] Thomas A Klar and Stefan W Hell.  Subdiffraction resolution in far-field fluorescence microscopy.  Optics letters, 24(14):954–956, 1999.
  * [21] Stefan W Hell and Matthias Kroug.  Ground-state-depletion fluorscence microscopy: A concept for breaking the diffraction resolution limit.  Applied Physics B, 60(5):495–497, 1995.
  * [22] Stefan W Hell.  Far-field optical nanoscopy.  science, 316(5828):1153–1158, 2007.
  * [23] Mats GL Gustafsson.  Nonlinear structured-illumination microscopy: wide-field fluorescence imaging with theoretically unlimited resolution.  Proceedings of the National Academy of Sciences, 102(37):13081–13086, 2005.
  * [24] Michael J Rust, Mark Bates, and Xiaowei Zhuang.  Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (storm).  Nature methods, 3(10):793, 2006.
  * [25] Sebastian Van de Linde, Anna Löschberger, Teresa Klein, Meike Heidbreder, Steve Wolter, Mike Heilemann, and Markus Sauer.  Direct stochastic optical reconstruction microscopy with standard fluorescent probes.  Nature protocols, 6(7):991, 2011.
  * [26] Eric Betzig, George H Patterson, Rachid Sougrat, O Wolf Lindwasser, Scott Olenych, Juan S Bonifacino, Michael W Davidson, Jennifer Lippincott-Schwartz, and Harald F Hess.  Imaging intracellular fluorescent proteins at nanometer resolution.  Science, 313(5793):1642–1645, 2006.
  * [27] Samuel T Hess, Thanu PK Girirajan, and Michael D Mason.  Ultra-high resolution imaging by fluorescence photoactivation localization microscopy.  Biophysical journal, 91(11):4258–4272, 2006.
  * [28] Thomas Dertinger, Ryan Colyer, Gopal Iyer, Shimon Weiss, and Jörg Enderlein.  Fast, background-free, 3d super-resolution optical fluctuation imaging (sofi).  Proceedings of the National Academy of Sciences, 106(52):22287–22292, 2009.
  * [29] Joerg Schnitzbauer, Maximilian T Strauss, Thomas Schlichthaerle, Florian Schueder, and Ralf Jungmann.  Super-resolution microscopy with dna-paint.  Nature protocols, 12(6):1198, 2017.
  * [30] Sandeep Pallikkuth, Cheyenne Martin, Farzin Farzam, Jeremy S Edwards, Matthew R Lakin, Diane S Lidke, and Keith A Lidke.  Sequential super-resolution imaging using dna strand displacement.  PloS one, 13(8):e0203291, 2018.
  * [31] Keith A Lidke, Bernd Rieger, Thomas M Jovin, and Rainer Heintzmann.  Superresolution by localization of quantum dots using blinking statistics.  Optics express, 13(18):7052–7062, 2005.
  * [32] B Christoffer Lagerholm, Laurel Averett, Gabriel E Weinreb, Ken Jacobson, and Nancy L Thompson.  Analysis method for measuring submicroscopic distances with blinking quantum dots.  Biophysical Journal, 91(8):3050–3060, 2006.
  * [33] Alexander Egner, Claudia Geisler, Claas Von Middendorff, Hannes Bock, Dirk Wenzel, Rebecca Medda, Martin Andresen, Andre C Stiel, Stefan Jakobs, Christian Eggeling, et al.  Fluorescence nanoscopy in whole cells by asynchronous localization of photoswitching emitters.  Biophysical journal, 93(9):3285–3290, 2007.
  * [34] L Stirling Churchman, Zeynep Ökten, Ronald S Rock, John F Dawson, and James A Spudich.  Single molecule high-resolution colocalization of cy3 and cy5 attached to macromolecules measures intramolecular distances through time.  Proceedings of the National Academy of Sciences, 102(5):1419–1423, 2005.
  * [35] Mark Bates, Bo Huang, Graham T Dempsey, and Xiaowei Zhuang.  Multicolor super-resolution imaging with photo-switchable fluorescent probes.  Science, 317(5845):1749–1753, 2007.
  * [36] Shalini T Low-Nam, Keith A Lidke, Patrick J Cutler, Rob C Roovers, Paul MP van Bergen en Henegouwen, Bridget S Wilson, and Diane S Lidke.  Erbb1 dimerization is promoted by domain co-confinement and stabilized by ligand binding.  Nature structural & molecular biology, 18(11):1244, 2011.
  * [37] Matthew P Gordon, Taekjip Ha, and Paul R Selvin.  Single-molecule high-resolution imaging with photobleaching.  Proceedings of the National Academy of Sciences, 101(17):6462–6465, 2004.
  * [38] Xiaohui Qu, David Wu, Laurens Mets, and Norbert F Scherer.  Nanometer-localized multiple single-molecule fluorescence microscopy.  Proceedings of the National Academy of Sciences, 101(31):11298–11303, 2004.
  * [39] Alexey Sharonov and Robin M Hochstrasser.  Wide-field subdiffraction imaging by accumulated binding of diffusing probes.  Proceedings of the National Academy of Sciences, 103(50):18911–18916, 2006.
  * [40] Steve Wolter, Ulrike Endesfelder, Sebastian van de Linde, Mike Heilemann, and Markus Sauer.  Measuring localization performance of super-resolution algorithms on very active samples.  Optics express, 19(8):7020–7033, 2011.
  * [41] Alex Small.  Multifluorophore localization as a percolation problem: limits to density and precision.  JOSA A, 33(7):B21–B30, 2016.
  * [42] Jan Christoph Thiele, Dominic A Helmerich, Nazar Oleksiievets, Roman Tsukanov, Eugenia Butkevich, Markus Sauer, Oleksii Nevskyi, and Jorg Enderlein.  Confocal fluorescence-lifetime single-molecule localization microscopy.  ACS nano, 14(10):14190–14200, 2020.
  * [43] Julia Riedl, Alvaro H Crevenna, Kai Kessenbrock, Jerry Haochen Yu, Dorothee Neukirchen, Michal Bista, Frank Bradke, Dieter Jenne, Tad A Holak, Zena Werb, et al.  Lifeact: a versatile marker to visualize f-actin.  Nature methods, 5(7):605, 2008.
  * [44] Hanieh Mazloom-Farsibaf, Farzin Farzam, Mohamadreza Fazel, Michael J Wester, Marjolein BM Meddens, and Keith A Lidke.  Comparing lifeact and phalloidin for super-resolution imaging of actin in fixed cells.  Plos one, 16(1):e0246138, 2021.
  * [45] Robert PJ Nieuwenhuizen, Keith A Lidke, Mark Bates, Daniela Leyton Puig, David Grünwald, Sjoerd Stallinga, and Bernd Rieger.  Measuring image resolution in optical nanoscopy.  Nature methods, 10(6):557, 2013.
  * [46] Alexander R Small.  Theoretical limits on errors and acquisition rates in localizing switchable fluorophores.  Biophysical journal, 96(2):L16–L18, 2009.
  * [47] Patrick Fox-Roberts, Richard Marsh, Karin Pfisterer, Asier Jayo, Maddy Parsons, and Susan Cox.  Local dimensionality determines imaging speed in localization microscopy.  Nature communications, 8:13558, 2017.
  * [48] Alex Small and Shane Stahlheber.  Fluorophore localization algorithms for super-resolution microscopy.  Nature methods, 11(3):267, 2014.
  * [49] Daniel Sage, Hagai Kirshner, Thomas Pengo, Nico Stuurman, Junhong Min, Suliana Manley, and Michael Unser.  Quantitative evaluation of software packages for single-molecule localization microscopy.  Nature methods, 12(8):717, 2015.
  * [50] Daniel Sage, Thanh-An Pham, Hazen Babcock, Tomas Lukes, Thomas Pengo, Jerry Chao, Ramraj Velmurugan, Alex Herbert, Anurag Agrawal, Silvia Colabrese, et al.  Super-resolution fight club: assessment of 2d and 3d single-molecule localization microscopy software.  Nature methods, 16(5):387, 2019.
  * [51] Graham T Dempsey, Joshua C Vaughan, Kok Hao Chen, Mark Bates, and Xiaowei Zhuang.  Evaluation of fluorophores for optimal performance in localization-based super-resolution imaging.  Nature methods, 8(12):1027, 2011.
  * [52] Fabian U Zwettler, Sebastian Reinhard, Davide Gambarotto, Toby DM Bell, Virginie Hamel, Paul Guichard, and Markus Sauer.  Molecular resolution imaging by post-labeling expansion single-molecule localization microscopy (ex-smlm).  Nature communications, 11(1):1–11, 2020.
  * [53] Bo Huang, Wenqin Wang, Mark Bates, and Xiaowei Zhuang.  Three-dimensional super-resolution imaging by stochastic optical reconstruction microscopy.  Science, 319(5864):810–813, 2008.
  * [54] Sri Rama Prasanna Pavani, Michael A Thompson, Julie S Biteen, Samuel J Lord, Na Liu, Robert J Twieg, Rafael Piestun, and WE 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(9):2995–2999, 2009.
  * [55] David Baddeley, Mark B Cannell, and Christian Soeller.  Three-dimensional sub-100 nm super-resolution imaging of biological samples using a phase ramp in the objective pupil.  Nano Research, 4(6):589–598, 2011.
  * [56] Yoav Shechtman, Steffen J Sahl, Adam S Backer, and WE Moerner.  Optimal point spread function design for 3d imaging.  Physical review letters, 113(13):133902, 2014.
  * [57] Sudhakar Prasad.  Rotating point spread function via pupil-phase engineering.  Optics letters, 38(4):585–587, 2013.
  * [58] Leonhard Möckl, Anish R Roy, Petar N Petrov, and WE Moerner.  Bgnet: Accurate and rapid background estimation in single-molecule localization microscopy with deep neural nets.  arXiv preprint arXiv:1909.08151, 2019.
  * [59] Hendrik Deschout, Francesca Cella Zanacchi, Michael Mlodzianoski, Alberto Diaspro, Joerg Bewersdorf, Samuel T Hess, and Kevin Braeckmans.  Precisely and accurately localizing single emitters in fluorescence microscopy.  Nature methods, 11(3):253, 2014.
  * [60] Peter B Stetson.  Daophot: A computer program for crowded-field stellar photometry.  Publications of the Astronomical Society of the Pacific, 99(613):191, 1987.
  * [61] Mohamadreza Fazel, Michael J. Wester, Hanieh Mazloom-Farsibaf, Marjolein B. M. Meddens, Alexandra S. Eklund, Thomas Schlichthaerle, Florian Schueder, Ralf Jungmann, and Keith A. Lidke.  Bayesian multiple emitter fitting using reversible jump markov chain monte carlo.  Scientific Reports, 9:13791, 2019.
  * [62] Christian Franke, Markus Sauer, and Sebastian van de Linde.  Photometry unlocks 3d information from 2d localization microscopy data.  Nature methods, 14(1):41, 2017.
  * [63] Johann Engelhardt, Jan Keller, Patrick Hoyer, Matthias Reuss, Thorsten Staudt, and Stefan W Hell.  Molecular orientation affects localization accuracy in superresolution far-field fluorescence microscopy.  Nano letters, 11(1):209–213, 2010.
  * [64] Sjoerd Stallinga and Bernd Rieger.  Accuracy of the gaussian point spread function model in 2d localization microscopy.  Optics express, 18(24):24461–24476, 2010.
  * [65] Yau Wong, Zhiping Lin, and Raimund J Ober.  Limit of the accuracy of parameter estimation for moving single molecules imaged by fluorescence microscopy.  IEEE Transactions on Signal Processing, 59(3):895–911, 2010.
  * [66] Hendrik Deschout, Kristiaan Neyts, and Kevin Braeckmans.  The influence of movement on the localization precision of sub-resolution particles in fluorescence microscopy.  Journal of biophotonics, 5(1):97–109, 2012.
  * [67] Anish V Abraham, Sripad Ram, Jerry Chao, ES Ward, and Raimund J Ober.  Quantitative study of single molecule location estimation techniques.  Optics express, 17(26):23352–23373, 2009.
  * [68] Carlas S Smith, Nikolai Joseph, Bernd Rieger, and Keith A Lidke.  Fast, single-molecule localization that achieves theoretically minimum uncertainty.  Nature methods, 7(5):373, 2010.
  * [69] Tingwei Quan, Pengcheng Li, Fan Long, Shaoqun Zeng, Qingming Luo, Per Niklas Hedde, Gerd Ulrich Nienhaus, and Zhen-Li Huang.  Ultra-fast, high-precision image analysis for localization-based super resolution microscopy.  Optics Express, 18(11):11867–11876, 2010.
  * [70] Steve Wolter, Anna Löschberger, Thorge Holm, Sarah Aufmkolk, Marie-Christine Dabauvalle, Sebastian Van De Linde, and Markus Sauer.  rapidstorm: accurate, fast open-source software for localization microscopy.  Nature methods, 9(11):1040, 2012.
  * [71] Norman Brede and Melike Lakadamyali.  Graspj: an open source, real-time analysis package for super-resolution imaging.  Optical Nanoscopy, 1(1):11, 2012.
  * [72] Rebecca Starr, Shane Stahlheber, and Alex Small.  Fast maximum likelihood algorithm for localization of fluorescent molecules.  Optics letters, 37(3):413–415, 2012.
  * [73] Kim I Mortensen, L Stirling Churchman, James A Spudich, and Henrik Flyvbjerg.  Optimized localization analysis for single-molecule tracking and super-resolution microscopy.  Nature methods, 7(5):377, 2010.
  * [74] François Aguet, Dimitri Van De Ville, and Michael Unser.  A maximum-likelihood formalism for sub-resolution axial localization of fluorescent nanoparticles.  Optics Express, 13(26):10503–10522, 2005.
  * [75] Ginni Grover, Keith DeLuca, Sean Quirin, Jennifer DeLuca, and Rafael Piestun.  Super-resolution photon-efficient imaging by nanometric double-helix point spread function localization of emitters (spindle).  Optics express, 20(24):26681–26695, 2012.
  * [76] Sheng Liu, Emil B Kromann, Wesley D Krueger, Joerg Bewersdorf, and Keith A Lidke.  Three dimensional single molecule localization using a phase retrieved pupil function.  Optics express, 21(24):29462–29487, 2013.
  * [77] Hazen P Babcock and Xiaowei Zhuang.  Analyzing single molecule localization microscopy data using cubic splines.  Scientific reports, 7(1):552, 2017.
  * [78] Yiming Li, Markus Mund, Philipp Hoess, Joran Deschamps, Ulf Matti, Bianca Nijmeijer, Vilma Jimenez Sabinina, Jan Ellenberg, Ingmar Schoen, and Jonas Ries.  Real-time 3d single-molecule localization using experimental point spread functions.  Nature methods, 15(5):367, 2018.
  * [79] Sudhir K Madan, BASABI Bhaumik, and JM Vasi.  Experimental observation of avalanche multiplication in charge-coupled devices.  IEEE Transactions on Electron Devices, 30(6):694–699, 1983.
  * [80] Rainer Heintzmann, Peter K Relich, Robert PJ Nieuwenhuizen, Keith A Lidke, and Bernd Rieger.  Calibrating photon counts from a single image.  arXiv preprint arXiv:1611.05654, 2016.
  * [81] Raimund J Ober, Sripad Ram, and E Sally Ward.  Localization accuracy in single-molecule microscopy.  Biophysical journal, 86(2):1185–1200, 2004.
  * [82] Tingwei Quan, Shaoqun Zeng, and Zhenli Huang.  Localization capability and limitation of electron-multiplying charge-coupled, scientific complementary metal-oxide semiconductor, and charge-coupled devices for superresolution imaging.  Journal of biomedical optics, 15(6):066005, 2010.
  * [83] Donald L Snyder, Abed M Hammoud, and Richard L White.  Image recovery from data acquired with a charge-coupled-device camera.  JOSA A, 10(5):1014–1023, 1993.
  * [84] Fang Huang, Tobias MP Hartwich, Felix E Rivera-Molina, Yu Lin, Whitney C Duim, Jane J Long, Pradeep D Uchil, Jordan R Myers, Michelle A Baird, Walther Mothes, et al.  Video-rate nanoscopy using scmos camera–specific single-molecule localization algorithms.  Nature methods, 10(7):653, 2013.
  * [85] Jerry Chao, Sripad Ram, E Sally Ward, and Raimund J Ober.  Ultrahigh accuracy imaging modality for super-localization microscopy.  nature methods, 10(4):335, 2013.
  * [86] William Hadley Richardson.  Bayesian-based iterative method of image restoration.  JoSA, 62(1):55–59, 1972.
  * [87] SS Reddi.  Multiple source location-a digital approach.  IEEE Transactions on Aerospace and Electronic Systems, (1):95–105, 1979.
  * [88] Norman Bobroff.  Position measurement with a resolution and noise-limited instrument.  Review of Scientific Instruments, 57(6):1152–1157, 1986.
  * [89] Leon B Lucy.  An iterative technique for the rectification of observed distributions.  The astronomical journal, 79:745, 1974.
  * [90] Mohammad Faisal, Aaron D Lanterman, Donald L Snyder, and Richard L White.  Implementation of a modified richardson-lucy method for image restoration on a massively parallel computer to compensate for space-variant point spread of a charge-coupled-device camera.  JOSA A, 12(12):2593–2603, 1995.
  * [91] Mark Goulian and Sanford M Simon.  Tracking single proteins within cells.  Biophysical journal, 79(4):2188–2198, 2000.
  * [92] D Thomann, DR Rines, PK Sorger, and G Danuser.  Automatic fluorescent tag detection in 3d with super-resolution: application to the analysis of chromosome movement.  Journal of microscopy, 208(1):49–64, 2002.
  * [93] Russell E Thompson, Daniel R Larson, and Watt W Webb.  Precise nanometer localization analysis for individual fluorescent probes.  Biophysical journal, 82(5):2775–2783, 2002.
  * [94] Steve Wolter, Mark Schüttpelz, Marko Tscherepanow, Sebastian Van de Linde, Mike Heilemann, and Markus Sauer.  Real-time computation of subdiffraction-resolution fluorescence images.  Journal of microscopy, 237(1):12–22, 2010.
  * [95] Sean B Andersson.  Localization of a fluorescent source without numerical fitting.  Optics express, 16(23):18714–18724, 2008.
  * [96] Adel Kechkar, Deepak Nair, Mike Heilemann, Daniel Choquet, and Jean-Baptiste Sibarita.  Real-time analysis and visualization for single-molecule based super-resolution microscopy.  PLoS One, 8(4):e62918, 2013.
  * [97] I Izeddin, J Boulanger, V Racine, CG Specht, A Kechkar, D Nair, A Triller, D Choquet, M Dahan, and JB Sibarita.  Wavelet analysis for single molecule localization microscopy.  Optics express, 20(3):2081–2095, 2012.
  * [98] Hongqiang Ma, Fan Long, Shaoqun Zeng, and Zhen-Li Huang.  Fast and precise algorithm based on maximum radial symmetry for single molecule localization.  Optics letters, 37(13):2481–2483, 2012.
  * [99] Ihor Smal, Erik Meijering, Katharina Draegestein, Niels Galjart, Ilya Grigoriev, Anna Akhmanova, ME Van Royen, Adriaan B Houtsmuller, and Wiro Niessen.  Multiple object tracking in molecular bioimaging by rao-blackwellized marginal particle filtering.  Medical Image Analysis, 12(6):764–777, 2008.