Analyzing engineered point spread functions using phasor-based single-molecule localization microscopy

The point spread function (PSF) of single molecule emitters can be engineered in the Fourier plane to encode three-dimensional localization information, creating double-helix, saddle-point or tetra-pod PSFs. Here, we describe and assess adaptations of the phasor-based single-molecule localization microscopy (pSMLM) algorithm to localize single molecules using these PSFs with sub-pixel accuracy. For double-helix, pSMLM identifies the two individual lobes and uses their relative rotation for obtaining z-resolved localizations, while for saddle-point or tetra-pod, a novel phasor-based deconvolution approach is used. The pSMLM software package delivers similar precision and recall rates to the best-in-class software package (SMAP) at signal-to-noise ratios typical for organic fluorophores. pSMLM substantially improves the localization rate by a factor of 2 - 4x on a standard CPU, with 1-1.5·104 (double-helix) or 2.5·105 (saddle-point/tetra-pod) localizations/second.


Introduction
Fluorescence microscopy is frequently employed in biological sciences due to its high selectivity and non-invasiveness. Conventionally, the obtainable optical resolution in fluorescence microscopy is given by Abbe's diffraction limit which is equal to the wavelength of the light divided by double the numerical aperture of the objective (~ 200 nm for visible light). A multitude of techniques summarized by the term superresolution (SR) microscopy or nanoscopy [1][2][3] have been developed, however, to obtain spatial information well below this limit. These techniques include (d)STORM (direct stochastic optical reconstruction microscopy) [4,5], PALM (photoactivatable localization microscopy) [6], SIM (structured illumination microscopy) [7], STED (stimulated emission depletion microscopy) [8], RESOLFT (reversible saturable optical fluorescence transitions) [9], SOFI (super-resolution optical fluctuation imaging) [10], SRRF (super-resolution radial fluctuations) [11] and MINFLUX (minimal photon fluxes localization microscopy) [12]. Single-molecule localization microscopy (SMLM) is the subcollection of super-resolution techniques in which the fluorescent emission profile, ordinarily referred to as a point spread function (PSF), of a single fluorophore is localized with a precision (~ 5 -40 nm) that can exceed the classical resolution limit by more than one order of magnitude [13][14][15][16]. SMLM is therefore an integral part of STORM and PALM, and has been extensively used in biological research [17][18][19], for example to study DNA transcription [20,21], CRISPR-Cas DNA screening [22][23][24], nuclear pore complexes [25,26], and microtubules [27]. In a conventional fluorescence microscope, a PSF from a single emitter in focus resembles an Airy pattern, which can be approximated by a 2-dimensional Gaussian function. This approach has been the basis of the earliest localization algorithms [13,16,28], which allow for determination of the emitter locations [16] as long as overlapping of PSFs is negligible. Besides Gaussian-based methods, these symmetrical PSFs have been analyzed via other mathematical frameworks, such as radial symmetry [29], cubic splines [30], or phasor (Fourier) analysis [31]. The shape of the PSF quickly deteriorates, however, if the emitter is out of focus (~100s of nm), leading to both a limited available axial range and inaccessibility of the absolute axial position [32]. Therefore, a variety of methods have been developed to modulate the shape of the PSF depending on the emitter's axial position [33]. Historically, the first method (astigmatism; AS) introduced a cylindrical lens in the emission pathway to create ellipsoid PSFs if the emitters are out of focus [34,35]. The extent of the deformation along with its orientation allows for determination of the axial position after a calibration procedure, and fitting of these PSFs could usually be performed by derivatized localization algorithms as the ones used for 2D PSFs [28,31,36]. However, the available axial range of astigmatism is limited to less than ~1 µm, which lead to the development of more advanced PSF shaping procedures that involve modulating the light in the pupil (Fourier) plane. Using a spatial light modulator (SLM), the principle was first employed to create a double helix (DH) pattern, in which the PSF is split in two separate lobes that non-degeneratively rotate around each other based on the emitters axial position, resulting in an usable axial range up to 2.5 µm [37]. Later, the same group theoretically maximized the information content of PSFs resulting in the Saddle-Point (SP) or Tetra-Pod (TP) designs, which are suitable for 3 µm (SP) or ≥ 6 µm (TP) axial ranges [38,39]. PSFs for both SP and TP are altered in the Fourier plane via a phase mask [38,39] or deformable mirror [40]. Determining the sub-pixel positions corresponding to the emitters via DH, SP, or TP PSFs, however, is more challenging than for symmetric or AS PSFs, as fitting with a single 2D Gaussian is insufficient. The current state-of-the-art fitting algorithms [41] rely on phase retrieving methods [40,42] or spline interpolation [26] to determine a PSF model based on calibration samples. A high-resolution PSF model can then be determined from these models which is fitted on experimental data. These methodologies can work with arbitrarily shaped PSFs, including DH, SP and TP. However, these methods are computationally expensive and thus time-consuming. Recently, real-time fitting localization of experimental PSFs have been achieved using graphical processing units (GPUs) [26], but this has not yet been achieved on computation processing units (CPUs), which would increase the accessibility and might allow implementations directly on the camera hardware. Here, we show fast retrieval of DH (1.5·10 4 loc/s) and SP/TP (2.5·10 5 loc/s) PSF localizations on a standard CPU via novel adaptations of the phasor-based single-molecule localization microscopy (pSMLM) algorithm [31]. We first explain the underlying methodology for DH and for SP/TP, termed circulartangent (ct-)pSMLM, and then explore the performance of the methods by analyzing simulated and experimental data. We have implemented all pSMLM versions (2D, AS, DH, SP/TP) in a recently published software package (SMALL-LABS [43]), resulting in user-friendly and open-source software to quickly perform sub-pixel localization including advanced background filtering options.

Software and hardware
All software was written and ran in MATLAB (MathWorks, UK) version 2018b on a 64-bit Windows 10 computer equipped with an Intel i5-8600 CPU @ 3.10 GHz, 16 GB RAM.

SMALL-LABS software
Our software package expands the original SMALL-LABS software [43] in several ways. Firstly, we added the original pSMLM-3D algorithm for 2D or astigmatism PSF sub-pixel localization, as well as the novel variations discussed in this manuscript. Next, a custom GUI was written to increase user accessibility. Lastly, the pre-and post-processing options are expanded with wavelet filtering [44], cross-correlation drift correction in three dimensions [45], and average shifted histogram result image generation [46,47].

Saddle-point PSF simulations
PSF simulations have been performed as described earlier [16,31] with NA = 1.25, emission light at 500 nm, 100 nm/pixel camera acquisition and 1000 PSFs for every intensity/noise combination. We used a full vectorial model of the PSF needed to describe the high NA case typically used in fluorescent superresolution imaging. The center of the PSF is located within ± 1 pixel of the center of the image. Zernike polynomials (primary astigmatism) and (secondary astigmatism) are introduced in a 0.5:-0.65 ratio [40], and z-positions were chosen randomly between -1.5 and +1.5 µm away from the focal plane.

Sub-pixel localization of single-molecule data
For double-helix (DH) sub-pixel localization, the four datasets from the 2016 SMLM challenge [41] were analyzed, which use experimental PSF models. These datasets differ in signal-tonoise (SNR) values ('high SNR', which mimics Alexa647 fluorophores, and 'low SNR', which mimics fluorescent proteins) and in emitter density ('low density' (LD) at 0.2 loc/µm 2 and 'high density' (HD) at 2 loc/µm 2 ). For double-helix localization, the following settings were used. For SMALL-LABS-pSMLM-DH, the temporal window length and the minimum duration of fluorophore on-time before it is discarded were both set to 150 frames. Filtering for region-ofinterest (ROI) finding was performed with a β-Spline wavelet filter with the threshold set to 1.9 times the standard deviation of a filtered frame. Single-lobe DH location was performed with a phasor radius of 4 pixels (low density) or 2 pixels (high density). The z-position was calculated via a calibration with identical phasor radius. For the SMAP software with fit3dSpline sub-pixel localization [26], a calibration was performed with a 33 x 33 pixel ROI. Then, localizations were identified via a mean calibrated PSF, with a 2.9 pixel Gaussian blur, using all calibrated z-positions. A threshold set to an absolute cutoff value of 86 (high SNR), 76 (low SNR, LD), or 29 (low SNR, HD) photons was used. The calibrated spline PSF was fitted with a 15 x 15 pixel ROI. Then, localizations with relative log-likelihood lower than -2 (low density) or -5 (high density) were discarded. For the localization of simulated saddle-point PSFs, we note that in order to prevent localization artefacts at specific z positions for ct-pSMLM (SI fig 1), a large ROI (> ~ 1.2x the maximum distance between the lobes) had to be used to determine lateral localization accurately, while a smaller ROI (~ 2.3 x 2.3 µm) was required for accurate axial localization, as ct-pSMLM with a large ROI cannot accurately describe the axial position around the focus. Then, all PSFs were localized directly with ct-pSMLM as described in the result section or via the SMAP software with fit3dSpline sub-pixel localization [26]. For ct-pSMLM, we used a 23 x 23 pixel ROI to calculate the z-position and a 43 x 43 pixel ROI to calculate the x and y position, and a 4 th -order polynomial was used to fit the calibration curve. For SMAP, calibration was performed with a 15-px Gaussian blur to find PSFs, a 51 x 51 pixel ROI, and 10 nm axial distance between every localization. For localization, a 10-px Gaussian blur was used to find PSFs, with a threshold of 20. A 41 x 41 pixel ROI spline fitting with 100 iterations based on the calibrated data was used to localize the PSFs. Localization of experimental saddle-point data was performed with the SMALL-LABS-pSMLM software package. A median background subtraction with temporal window length of 150 frames and minimum duration of fluorophores to be discarded of 100 frames was used. Localizations were identified via a bandpass filter with a 95 threshold percentile. Potential lobes of saddle-point point spread functions were identified with a 3 pixel radius ROI 2D phasor fitting routine. Ct-pSMLM fitting was then performed with an 11 pixel radius ROI around the center of localizations. Calibration was performed using simulated point spread functions at varying z positions, consisting of 5000 photons on a noiseless background, with deformations similar to experimental data. Three-dimensional cross correlative drift correction was performed via the SMALL-LABS-JH software, with 10 lateral subpixels and 10 temporal bins. The average shifted histogram image was created using ThunderSTORM [46], using 50 nm axial bins and 10 lateral subpixels.

Assessment of localization performance
For double-helix (DH), localizations between ground-truth (GT) and software 1 (S1; SMALL-LABS-pSMLM-DH) and between GT and software 2 (S2; SMAP with fit3dSpline) are linked on a frame-by-frame basis, with a maximum allowed lateral distance of 250 nm, and a maximum allowed axial distance of 500 nm. The median offset between GT and S1 and between GT and S2 is calculated and subtracted from the S1 and S2 datasets, to avoid introducing consistent offset errors in the RMSE calculations. The linking of localizations between GT and S1/S2 is repeated, as localizations can be shifted in/out of the maximum linking distance due to the median offset. Of this linked dataset, the Jaccard index is calculated as follows: = + + where TPo, FP, and FN are the true positive, false positive, and false negative localizations, respectively. Then, only localizations that are present in all three datasets (GT, S1 and S2) are selected, and of these localizations, the root mean square error (RMSE) in a single dimension is calculated as follows: where pi indicates the position of localization i in any dimension, and S indicates S1 or S2. For saddle-point (SP), 1000 PSFs for every signal-to-noise combination were simulated (section 2.3), after which a calibration curve was created in SMALL-LABS-ct-pSMLM or SMAP with PSFs containing 3·10 4 photons on a 1 photon/pixel background. For SMAP localizations, obtained localizations well outside the expected regime (10 pixels or further removed from center) were discarded, and frames containing multiple or no localizations were fully discarded. Note, localizations obtained with SMAP that were clearly misfitted (an offset in z by at least 3 times the average z offset calculated by ct-pSMLM) were discarded; no such discarding was performed for ct-pSMLM. For both ct-pSMLM and SMAP, the x, y and z positions were compared with the ground-truth, and the standard deviation of this offset was calculated for every intensity and background combination and is shown in the results. We note that the mean of the offset was centered around 0 for every tested intensity/noise/software combination.

Single-molecule microscopy
For SMLM experiments, we used a home-built super-resolution microscope similar to one reported previously [24]. Briefly, light from a fiber-coupled 642 nm laser (Omicron, Germany) was collimated using an achromatic lens (f = 30mm, Thorlabs) and conducted to a parabolic mirror (RC12APC-P01, Thorlabs). The laser light was then focused using an achromat lens (f = 150mm, Thorlabs) in front of a polychroic mirror (ZT532/640rpc, Chroma) into the backfocal plane of an 100x oil-immersion objective (CFI Plan Apo, NA = 1.45, Nikon Japan) such that a highly inclined illumination (HiLo) profile with a total laser power of ~70 mW was achieved. Emitted fluorescence passing the objective, the polychroic mirror and a bandpass filter (ZET532/640m-TRF, Chroma) was then guided into a 4f geometry using the following lenses (

STORM experiment
A SAFe sample containing immobilized Cos7 fibroblasts from Green African Monkeys (ATCC) with Alexa Fluor 647 labeled tubulins was purchased from Abbelight (Paris, France). A nitrogen-flushed buffer containing 50 mM TRIS pH8, 10 mM NaCl, 10% glucose, 50 mM 2-mercaptoethanol, 68 µg/mL catalase, and 200 µg/mL glucose oxidase [27] was added to the sample chamber which was sealed off before the measurements. 60.000 frames of 20 ms length were recorded using the setup described in section 2.6. Analysis of the singlemolecule data was performed as specified in section 2.4.

Principles of engineered PSF localization with pSMLM: Double-helix: DH-pSMLM
To localize double-helix (DH) PSFs, we rely on the fact that pSMLM-2D provides accurate lateral localization even when using a relatively small ROI around the center of an emitter [31]. Therefore, the two lobes rotating around each other (Fig. 1a) can be localized separately. During calibration, the distance and rotation between the two lobes is plotted against the axial position (Fig. 1b). The rotation is fitted with a third-order polynomial. This polynomial is weighted on the inverse of the standard deviation of each axial position if more than one calibration bead has been used. The lateral position is calculated as being the average lateral position of the two lobes, corrected for a 'wobble' factor. This wobble factor is determined in x and y as function of the emitters axial position (Fig. 1c) by comparing the lateral localization at all axial positions with the lateral localization at the axial center of the calibration dataset. The average of this wobble effect over an axial sliding window (user-defined, default value is set to 5 axial positions) is determined during calibration and stored for future correction of lateral localization calculation (Fig. 1d). To extract positional information, first a standard pSMLM-2D fitting is performed [31]. The localizations in each frame are compared with each other to find pairs within the expected distance regime (determined during calibration; minimum and maximum of distance between lobe centers, with a ~10% error margin), and are discarded if no pair can be found. During the linking of the lobes, priority is given to lobes that only have a single possible counter-lobe over those that have multiple options to reduce mis-fitting of closely positioned DH PSFs. The axial position is then determined from the rotation of the two lobes via the calibration curve. The obtained distance between lobes is checked against the distance determined during calibration at the found axial position, and the localization is discarded if these values differ more than ~100 nm (user defined). Lastly, the lateral position is determined from the mean of the 2D-determined position of the two lobes, and corrected for the wobble determined during calibration (Fig 1d).
3.2 Principles of engineered PSF localization with pSMLM: Saddle-point and tetra-pod: ct-pSMLM We analyze saddle-point (SP) and tetra-pod (TP) PSFs with an adapted phasor-based localization methodology. SP and TP have similar characteristics and show separation of a single point when in focus into two lobes above and below the focus in perpendicular directions [38,39]. Moreover, they are based on similar PSF deformations introduced by primary and secondary astigmatism Zernike coefficients [40]. We modified a spectral phasor-based approach [48] in which the convolution of arbitrary profiles in real space is a linear combination of their respective phasor representations in phasor space. In this approach, the normalized intensity ratio between the original profiles in the convoluted profile (real space) is represented as the distance of the original phasor profiles to the convoluted phasor profile (phasor space). This entails that if two profiles are combined with a 1:1 ratio, the convoluted phasor representation is on the mid-point of the line between the phasor representations of the original profiles. In SP and TP PSFs, the final spatial representation of the PSF is a convolution of two separated lobes of identical intensity. Thus, SP and TP PSFs can be treated as a 1:1 ratio of arbitrary profiles that are separated at a varying distance, depending on the axial position of the emitters. Note that by orientating the respective optical components correctly, this separation can be achieved perfectly on the xor y-axis. Therefore, the value for the separation dlobes, along with the orientation of this separation provides suitable information for calibration of SP and TP PSFs.
To determine the separation dlobes with phasor-based singlemolecule localization microscopy (pSMLM), we assume that the width of the individual lobes in the direction of the convolution is identical to the width of the convoluted PSF in the other, unconvoluted spatial direction. For illustration, we show a combination of two 2-dimensional Gaussian distributions (Fig.  2a,b,c). The phasor representation of the individual Gaussian distributions is represented by a single phasor for both dimensions, each having a certain, but different, angle representing the emitter's position in real space [31]. Then, if reasoned from the convoluted PSF (Fig. 2c) to obtain the individual lobes, the tangent at the magnitude circle in the convoluted spatial dimension (broad spatial dimension; small phasor magnitude; black cross located on the red circle in Fig.  2c) will intersect the magnitude of the smaller spatial dimension (large phasor magnitude; represented as a blue circle) at two points ( Fig. 1c; magenta and orange dots). These points are a measure for original arbitrary profiles with identical spatial sizes in both dimensions that combine in a 1:1 intensity ratio to result in the convoluted profile. The angle between these two obtained intersectional points (θlobes) in phasor space is a direct normalized value for the distance dlobes in real space (Fig. 2c). We call this method circular-tangent pSMLM (ct-pSMLM). The obtained dlobes is used to create well-defined calibration curves for both SP and TP (SP shown in Fig. 2d,) that can be fitted with arbitrary functions (e.g. a fourth-order polynomial) to deduce axial positional information from experimental PSFs (Fig. 2e). The lateral localization information of the SP or TP PSFs is still inherently present in the original phasor-representation of the complete PSF.
Determining localization of the SP or TP PSFs in the SMALL-LABS-pSMLM software consists of two parts: finding the central positions and further analysis with ct-pSMLM. The mid-points of single PSFs are determined by first checking whether two detected emitters that could represent two lobes of a single SP/TP PSF belong to the same PSF. If these emitters have little deviation in one dimension (<0.5 px) and are slightly separated in the other dimension (less than the calibrated maximum distance), the mid-point of these emitters is calculated and stored. If no other lobe can be found, it is assumed the located emitter is the mid-point of the SP/TP PSF. Then, ct-pSMLM is performed around the central point with a reasonably large region of interest (> 2 µm) to obtain dlobes and to calculate the axial position.

Double-helix
To evaluate the performance of DH-pSMLM, we performed fitting of simulated datasets [44] via the full pSMLM-updated SMALL-LABS software package and compared with the currently best performing non-machine learned localization algorithm (experimental PSF spline fitting methodology incorporated in SMAP [26]). As the ground truth of these datasets is publicly available, we were able to extract (Table 1) quantitative performance parameters such as the expected deviation of localization accuracy in all three dimensions (root mean squared error, RMSE) and the Jaccard index JACC (a measure for correctly and incorrectly localized particles [41]). These performance parameters are calculated from localizations that were found in both software packages and in the Ground-Truth datasets. We observe that both SMALL-LABS-pSMLM and SMAP have comparable RMSE errors (Table 1) in the order of 10-25 nm for the low density (LD), high signal to noise (SNR) dataset which are similar to the ones reported previously for SMAP [26,41]. However, at low signal to noise levels, SMAP outperforms SMALL-LABS-pSMLM on all performance indicators. This is presumably due to SMAP using the full PSF at once, while SMALL-LABS-pSMLM splits localization in two steps. This results in SMALL-LABS-pSMLM working with a lower apparent signal to noise level, causing a lower localization accuracy. We note that the reported RMSE values for SMAP analysis of the LD, low SNR dataset are counter-intuitively better than those of SMAP analysis of the LD, high SNR dataset. This is a result of the RMSE calculation methodology used (Material and methods), as only localizations that are found in both software analyses as well as in the ground-truth are used for RMSE calculations. We observe that SMALL-LABS-pSMLM outperforms SMAP in terms of localization recall rates (Jaccard index, Table 1) at high SNR (23% increased), but not at low SNR (4% decrease). The Jaccard values for SMAP are slightly lower than reported earlier [41] (Material and methods), but can be compared directly with the Jaccard values for SMALL-LABS-pSMLM reported here. We note that no background subtraction is performed in SMAP, while SMALL-LABS-pSMLM subtracts the background based on foreground temporal variations (SMALL-LABS [43]). Both software packages are not capable of recognizing HD PSFs with a good recall rate, although SMALL-LABS-pSMLM outperforms SMAP in all conditions, as single DH lobes are localized with only 5 x 5 pixel ROIs, decreasing the influence of the other nearby emitters. SMALL-LABS-pSMLM is ~3 -4 x faster compared to SMAP for low density datasets, and ~2x faster for high density datasets. However, most analysis time for SMALL-LABS-pSMLM (~65%) is spend on the background correction and format conversion rather than approximate localization or sub-pixel DH-pSMLM localization (Supplementary Table 1). The localization procedure itself can achieve 1 -1.5·10 4 localizations per second on a standard CPU.

Saddle-point and tetra-pod
The performance of the localization of saddle-point (SP) PSFs was assessed and compared to experimental PSF spline fitting [26]. As ct-pSMLM is a non-iterative method, high localization rates of up to 2.5·10 5 localizations per second were achieved on standard CPUs (Fig. 3a). This is an order of magnitude lower than traditional pSMLM-3D [31], mostly due to the large required region of interest around the PSF (>2 µm; here 23 x 23 px), and partly due to the additional computations required for ct-pSMLM. Taken alone, the additional computations of ct-pSMLM compared to pSMLM-3D only result in a 10 -40% decrease in localization rates (~10% for at a large region of interests 23 x 23 px; ~40% decrease for 7 x 7 px). The lateral localization accuracy of ct-pSMLM is in line with that of experimental PSF spline fitting (Fig 3b), and decreases from ~100 nm (~1 pixel) at typical (~200 -1300) photon values for fluorescent proteins to ~10 nm (~0.1 pixels) at typical (~(2 -11) x 10 3 ) photon values for organic fluorophores. The localisation accuracy is roughly one order of magnitude lower than the lateral localization accuracy of non-engineered PSFs at high photon values (~0.08 pixels and ~0.01 pixels, respectively [31]), and roughly 1.5x worse than AS PSFs (~0.05 pixels [31]) caused by lower effective signal to noise ratio due to the expanded PSF. We observed a lower limit in lateral localization accuracy for SMAP fitting of ~10 nm (~0.1 pixel), which has an unknown origin.
The axial localization accuracy of ct-pSMLM increases with increasing photon values as well (Fig 3c). The average axial accuracy at typical photon values for organic fluorophores is ~40 nm, which is ~3x worse than AS PSFs with similar total photon counts and photons/pixel background [31]. The best obtainable axial accuracy is limited by the sub-optimal fitting of the calibration curve to around 11 nm, which is similar to AS PSFs (Fig 3c, SI fig 2, [31]). We attribute this lower axial accuracy of SP PSFs compared to AS PSFs again to lower effective signal to noise ratios due to expanded PSFs. Up to 20% of SMAPlocalized emitters had to be discarded from calculating the z offset, as these were substantially misfitted (> 3 times the corresponding ct-pSMLM z offset, see Methods). We furthermore demonstrate the implementation of ct-pSMLM in SMALL-LABS-pSMLM by analysing an experimental STORM experiment showing labelled microtubule of a Monkey cell line (Fig 3d, Material and methods). The total analysis time for the SMALL-LABS-pSMLM analysis for the ~8 GB, 60.000 frames dataset containing 1.5 million localizations was ~15 minutes on a standard CPU, including file conversion (Supplementary Table  1), of which the ct-pSMLM sub-pixel fitting routine comprised just 77 seconds.

Discussion
Here we present additions on the phasor-based single-molecule localization microscopy (pSMLM) framework to localize double helix (DH), saddle-point (SP), and tetra-pod (TP) PSFs with very good accuracy and speed on standard CPUs. In the current implementation, DH-pSMLM can achieve up to 1.5·10 4 localizations/second on a 3.10 GHz processing unit, while ct-pSMLM, the basis for SP and TP localization, can achieve up to 2.5·10 5 localizations/second. Specifically, ct-pSMLM is designated for real-time localization methods, combined with computationally inexpensive filtering and background subtraction methods, to better enable (automated) feedbackoriented SMLM instrumentation. Possibly a pSMLM-based methodology could be implemented on the integrated circuits of cameras to further increase end-user accessibility of advanced single-molecule techniques.
The DH-pSMLM implementation in the SMALL-LABS software has similar performance as the current state-of-the-art methods when using organic fluorophores, while decreasing the overall analysis time when ran on a CPU. We note that our implementation is not particularly sensitive to overfitting, as stringent constraints on the pair-finding are used. This allows for some false positives during the initial localization steps which are later discarded. Ct-pSMLM improves on our previous phasor implementation [31] by offering a direct way of determining the distance between two emission peaks of a single PSF, well-suited for the quantification of the axial position in SP and TP PSFs. Here, perfect horizontal and vertical elongation of the PSFs is a requirement for ct-pSMLM to perform. Our algorithm is capable of retrieving the emitters location with a precision similar to current best non-machine learning localization algorithm [41], and is mostly limited by the fitting of the calibration curve. Naturally, for ct-pSMLM to work correctly, no other emitters or highly inhomogeneous background should be present in the fitting region. As SP and TP PSFs require large ROIs (~23 x 23px), this results in substantially lower accessible emitter density compared to approaches using standard and astigmatic PSFs. For high-density engineered PSF localization approaches, we point to alternative approaches such as deep learning [49,50] or matching pursuit [51]. We incorporated the novel pSMLM-derivative localization methodologies in the SMALL-LABS software [43]. The updated SMALL-LABS-pSMLM software package expands the original work with a user-friendly GUI, wavelet filtering, drift-correction in 3D, and result image generation. We believe that the software package strikes an excellent balance between fast analysis, accurate results, experimental freedom, good expandability, and hassle-free installation and operation. The software is freely available at https://github.com/HohlbeinLab/SMALL-LABS-pSMLM.