1. Introduction

Stochastic Optical Reconstruction Microscopy (STORM) and related techniques use stochastic switching and high-precision localization of single molecules to achieve image resolution beyond the diffraction limit [1–3]. In this imaging modality, sparse subsets of fluorescent emitters are activated stochastically, localized to sub-diffraction precision, and then inactivated. This process is repeated, and a final image is constructed from the individual localizations of emitters. Thus, the resolution of this image is no longer limited by diffraction but determined by the accuracy with which each emitter can be localized and the density of emitters [4–7].

However, because individual emitters are activated stochastically and observed as diffraction limited spots, the maximum density of resolvable emitters per imaging frame is limited. This emitter density limit in turn sets the number of camera frames that are required to construct a super-resolution image and, hence, the temporal resolution of STORM imaging. Thus, a lower permissible emitter density per frame may restrict the ability to follow processes in dynamic environments such as live cells [4,6,8]. In parallel, the inability to resolve overlapping fluorophores may also lead to underestimates of the number of emitters from regions of higher density, effectively limiting the dynamic range of the STORM image. For these reasons, better methods for handling overlapping emitters could substantially improve the performance of STORM imaging.

Recently, several methods have been introduced to resolve partially overlapping emitter images [9–15]. Among these methods, one powerful approach [15] allows the highest density of emitters to be localized to date by adopting techniques from the rapidly advancing field of compressed sensing. In the basic compressed sensing problem, an under-determined, sparse signal vector, i.e. a signal in which only a few basis elements are non-zero, is reconstructed from a noisy measurement in a basis in which the signal is not sparse. In the context of STORM, the sparse basis is a high resolution grid, in which fluorophore locations are presented—effectively an up-sampled, reconstructed image—while the noisy measurement basis is the lower resolution camera pixels, on which fluorescence signal are detected experimentally [15]. In this framework, the optimal reconstructed image is the one that contains the fewest number of fluorophores but reproduces the measured image on the camera to a given accuracy when convolved with the optical point-spread-function (PSF) (Fig. 1). Remarkably, this technique can accurately process emitter densities 4-6 fold higher than what could be handled by a sparse emitter fitting algorithm, which can only analyze non-overlapping images of emitters, and approximately 2-3 fold higher than what could be analyzed by other overlapping-emitter-fitting algorithms [15]. However, the interior point method used in the reported compressed sensing approach [15] is computationally intensive and, by our estimate, reconstruction of a typical STORM image of a 40 µm × 40 µm field of view using this approach would take a few months on a typical research computer.

 

Fig. 1 Schematic depiction of STORM imaging using compressed sensing. A subset of fluorescent emitters (red dots; top left) from a labeled sample (bottom left) are activated stochastically. The activated emitters are close enough that the individual emitters are not distinguishable in a 4 pixel × 4 pixel conventional image (top left). Using compressed sensing, a high resolution grid of fluorophore locations is reconstructed (top right) from the low resolution image (top left) under the constraint that this reconstruction contains the smallest number of emitters that can reproduce the measured conventional image up to a given accuracy. This process is repeated and the individual reconstructed frames are summed to produce the final high resolution reconstruction (lower right) of the original sample.

Download Full Size | PDF

Here, we introduce an alternative algorithm for reconstructing STORM images using the compressed sensing approach. This algorithm, referred to as L1-Homotopy (L1H) [16–18], is substantially faster than the previous approach which employed interior point methods, producing equivalent reconstructed STORM images in one to two orders of magnitude less time depending on the emitter density. Hence, the L1H algorithm should now allow the compressed sensing analysis of STORM data to become a routine practice.

2. Theory

In the compressed sensing approach to STORM image reconstruction, the reconstructed image is a high resolution grid of fluorophore locations that, when convolved with the optical PSF, reproduces the measured low resolution image on the camera to an accuracy set by the expected noise in the image. Mathematically, this reconstruction can be achieved by solving the following constrained-minimization problem:

Equation (1) was previously solved using interior point methods [15]—iterative algorithms broadly applicable to a large class of minimization problems. Unfortunately, at each iteration such methods must invert an N × N matrix—an operation with complexity that scales as O(N³). As N is typically of order 10³ – 10⁴, 10⁹ – 10¹² floating point operations are needed for each iteration. In practice, reconstructing a ~40 µm × 40 µm STORM image with ~20 nm resolution using Eq. (1) and interior point methods takes about 30 minutes per camera frame on a high-end desktop PC. Since the typical movie used to construct a single STORM image contains thousands to tens of thousands of frames, it would take thousands of hours to reconstruct a typical STORM image. The long computational time has, thus, limited the use of compressed sensing in the reconstruction of STORM images to small fields of view, for example ~3 µm × 3 µm [15].

However, interior point methods are not necessarily the most efficient way to solve this minimization problem, and there is a rapidly growing collection of alternative methods [19, 20]. Many of these alternative algorithms solve the following minimization problem:

However, Eq. (2) presents a challenge in the analysis of STORM data. Calculation of the appropriate noise parameter ε for the analysis of STORM data is simple. It can be done directly from the image data itself [15] (see Methods). However, identification of the value of λ for which Eq. (2) produces the equivalent solution to Eq. (1) is often not trivial since there is no known analytical connection between the noise parameter ε and λ for the general problem [21]. Thus, many of the alternative compressed sensing algorithms must first be modified to search for the appropriate value of λ before they can be applied to the analysis of STORM data.

Among the alternative algorithms that solve Eq. (2), we identified one algorithm for which the search for the appropriate λ may be particularly efficient for STORM imaging data: the L1-Homotopy (L1H) algorithm. This algorithm functions by exploiting a simple geometric property of the solution space of Eq. (2) to trace solutions of Eq. (2) from one value of λ to another value. In particular, as the value of λ is decreased, the balance shifts from favoring the sparsity of emitters to favoring the accuracy in the optimal solution, and this solution changes in a continuous, piece-wise-linear fashion as a function of λ. Importantly, the slope of this solution changes only at well-defined ‘break-points’ where individual emitters are either added to, removed from, or moved to a different location within the solution, and the value of λ at which each break-point occurs can be calculated from the solution at adjacent break-points. These properties are illustrated in Fig. 2 for a simple one-dimensional problem. The L1H algorithm begins at a large value of λ where the solution is obvious, $x=0$, and steps downward in λ from break-point to break-point until it reaches the desired value of λ.

 

Fig. 2 One-dimensional (1D) illustration of the properties of the solution space exploited by L1-Homotopy. (a) A simulated 1D image with two emitters before (blue) and after (red) convolution with a Gaussian point-spread-function. (b) Results of the L1H analysis of the image shown in (a) displayed as a kymograph of the amplitude of each up-sampled pixel, i.e. potential fluorophore location (rows), as a function of the homotopy parameter, λ. Note that as λ decreases, increasingly favoring accuracy over sparsity, the single initial peak splits into two peaks, representing two fluorophore localizations. (c) The amplitude of two adjacent up-sampled pixels, e.g. emitter locations, as a function of λ for the left emitter (top) and the right emitter (bottom). The amplitude of these pixels, depicted in red, blue, green and black symbols, corresponds to the pixels marked by red, blue, green and black arrows in (b). The amplitude of these pixels are piece-wise linear functions of λ. The break-points (dashed grey lines), where the slopes change, correspond to the addition, removal, or movement of a possible emitter, i.e. where the amplitude of a pixel changes from zero to non-zero or vice versa (marked by gray circles). Note that “movement” of an emitter from one up-sampled pixel to another is actually accomplished by removing it from the solution and then adding another emitter at a new location.

Download Full Size | PDF

While matrix inversions are needed at each break-point to calculate the position of the next break-point, these inversions typically involve a smaller matrix than those used by interior point methods. In particular, the required matrix contains a number of elements equal to the number non-zero up-sampled pixels. If this number—which scales as the number of reconstructed emitters, K—is small enough, we find that, in practice, the computational complexity of L1H is no longer dominated by the inversion of this matrix but rather by other matrix multiplications whose complexity scale roughly as O(N × K). By contrast, the matrix inversion typically dominates the cost of interior point methods, and its complexity scales as O(N³). Thus, when the solution vector is sparse, as it is in STORM images, K will be much less than N and the computational cost of L1H should be dramatically lower than the cost of interior point methods.

One advantage of L1H over other alternative compressed sensing algorithms is that it naturally searches the space of possible λ values as part of its normal iteration, and L1H can be modified to identify the λ that corresponds to the appropriate value of ε. To accomplish this task, we introduced a simple stopping criterion into the L1H algorithm. After the identification of each break point, our algorithm computes the residual error between the reconstructed image at that break point and measured image, i.e. the L2 norm in Eq. (1). If that residual error is larger than the target error set by ε, then the algorithm continues onto the next break point. If the residual error is less than ε, then iteration stops, and the algorithm uses mid-point bisection to find the value of λ and the intensity of the localized emitters for which the residual error is equal to ε. This solution represents the final reconstructed image, and the final value of λ corresponds to the value for which the solutions of Eqs. (1) and (2) are equal. Because calculation of the residual error involves K × M calculations, it does not significantly increase the computational cost of the L1H algorithm.

Of the existing algorithms that solve Eq. (2), L1H is the only algorithm which has the property that identification of the correct λ is a byproduct of its natural iteration. Thus, we reasoned that L1H should be the optimal approach for STORM image reconstruction. However, it is possible that when combined with an efficient search algorithm for λ, other alternative compressed sensing algorithms might be comparable to L1H in terms of computational complexity. To test this possibility, we considered one of the more commonly used efficient compressed sensing algorithms, the Fast Iterative Shrinkage Thresholding Algorithm (FISTA) [22]. This algorithm also requires on the order of N × M calculations per iteration, and, thus has a computational cost per iteration that is similar to L1H. However, each iteration in FISTA does not produce an exact solution to Eq. (2), as is the case for L1H, but rather serves only to improve the accuracy of an approximation to Eq. (2). Thus, in order for FISTA to be competitive with L1H, it must produce an approximation of the solution with sufficient accuracy in a relatively small number of iterations.

We tested the number of iterations of FISTA required to produce a reasonable approximation to the solution of Eq. (2). Even when provided with the correct value of λ as identified by L1H, FISTA required 100 – 1000 times more iterations to converge to an approximate solution of a 7 pixel × 7 pixel data set than is required for L1H to solve the problem exactly. We expect similar limitations to arise for other iterative approaches such as gradient projection, proximal gradient, and augmented Lagrange multiplier methods [19, 20] which would also require an iterative search over λ in addition to their natural iterative improvement of the solution to Eq. (2). Thus, L1H appears to be particularly well suited for the analysis of STORM data using compressed sensing.

3. Methods

Algorithm and Image Analysis. The L1H algorithm used here was written in C with a Python interface and is based on the approach described by Donoho and Tsiag [18] with two modifications. First, because negative photon numbers are unphysical, we restrict the algorithm so that it only considers non-negative values for $x$. Second, we modified this algorithm so that it stops when the appropriate value of λ has been identified, as discussed above. The FISTA algorithm used here was also written in C with a Python interface following [22]. Source code and executables of the L1H and FISTA algorithms are available at http://zhuang.harvard.edu/software.html. We used the Matlab interior-point-method implementation of the compressed sensing algorithm provided in [15] for comparisons. The core of this code is the optimized interior point methods provided in the freeware package CVX [23]. All reported data were calculated with the default precision settings of this package; however, we observe little difference in the required computation time when these settings are optimized for speed. All analysis was conducted on a single core of an i7-3770 Sandy Bridge 3.4 GHz CPU.

To limit the size of the $A$ matrix, we divided the fluorescence images from the camera into sub-regions and analyzed these regions separately as was done previously with interior point methods [15]. Specifically, we divided the measured images into multiple, partially overlapping 7 pixel × 7 pixel regions. To analyze each of the 7 pixel × 7 pixel regions, we allowed emitters to be placed within an 8 pixel × 8 pixel space, which corresponds to a 1/2 pixel extension in each direction from the measured region, in order to take into account the fluorescence signal contributed by emitters outside of the measured region due to the finite PSF size, as shown in Fig. 3.To further limit errors due to edge effects, we only retained emitters located within a central 5 pixel × 5 pixel region. These 5 × 5-camera-pixel up-sampled images, with neighboring images touching but not overlapping with each other, were then stitched together to form the final STORM image. In addition, following [15], we included a background term in the $A$ matrix that allows for non-zero uniform background correction. As the regions that are analyzed are small (7 × 7 pixels) this background correction makes the solution fairly robust to the slowly varying background typically encountered in fluorescence images.

 

Fig. 3 Analysis of a sub-region of the image by L1-Homotopy. The wide-field fluorescence image detected by the camera is first divided into partially overlapping, sub-regions of 7 × 7 camera pixels in size (a camera pixel is represented by a black box). One such region is shown here. In the reconstruction of this region, we allow emitters to be placed within an 8 × 8 camera pixel space, corresponding to a 1/2 pixel extension in each direction from the 7 × 7 camera pixel region, and the emitter positions are allowed on a finer, 8-fold up-sampled grid (red dots). To further limit errors that might arise due to edge effects, only emitter localizations located within a central 5 pixel × 5 pixel block defined by the thicker black line are kept. Multiple of these non-overlapping, 5-pixel wide, up-sampled images are then stitched together to form the final STORM image.

Download Full Size | PDF

For a noise-free detector, the noise parameter $\epsilon$ would be theoretically set by the Poisson noise of the photons, i.e. $\epsilon =\sqrt{{\displaystyle {\sum }_i{b}_i}}$. However, it was previously shown that a slightly larger noise parameter produced higher quality reconstructions [15]. For simulated data, the optimized noise parameter $\epsilon$ was shown to be 1.5 times the theoretical noise parameter. For real data collected with an electron-multiplying CCD camera, the optimized noise parameter $\epsilon$ was shown to be 2.1 times the theoretical noise parameter. We also used these optimized $\epsilon$values in this work.

Simulated Data. For all simulated data, we randomly distributed emitters on an 8-fold up-sampled grid and generated the number of photons for each emitter from a lognormal distribution with a mean and standard deviation of 3000 and 1700 photons, respectively. A Gaussian PSF of width equal to one camera pixel was used to convolve the signal from individual emitters, and the photons from all up-sampled pixels within a low-resolution, camera pixel were summed. Finally, a uniform background of 70 photons was added to each camera pixel and the image was corrupted with Poisson distributed noise with a mean set by the number of photons in each pixel. In all simulations, the camera pixel is assumed to be 167 nm in size. These values were chosen to match our typical experimental data. Using this approach, 7 × 7 camera pixel images were generated, and these simulated images were analyzed by both the L1H and CVX methods to compare the accuracy of the reconstructions. We also generated 256 × 256 camera pixel images to compare the computation time of these two methods. Reported computation times are the average of three independent trials.

Algorithm Comparison. To compare the quality of the reconstructions between different algorithms, we used a single-emitter fitting algorithm [24, 25], the multi-emitter fitting algorithm 3d-DAOSTORM [13], L1H, and CVX to analyze 256 × 256 cameral pixel, five frame data simulated with the same parameters as above. We used the 2d option of 3d-DAOSTORM which produces similar results to the original DAOSTORM algorithm [9]. We converted the up-sampled reconstructed images produced by L1H and CVX to molecule lists by combining all non-zero pixels with at least one of their 8 adjacent neighbors also with a non-zero value into a single molecule whose location is determined by the weighted centroid of these connected pixels as described previously [15]. Only emitters reconstructed within a central 150 × 150 camera pixel region of interest were included in the analysis to minimize edge effects. The reconstructed density is the number of localized emitters in the region of interest. The XY localization error is defined as the median Cartesian distance between each reconstructed molecule and the nearest real molecule. This distance is related to the standard deviation of a Gaussian distributed position error,$\sigma$, via $\sqrt{2\ln (2)}\sigma$.

Real Data. To test the L1H algorithm on real data, we collected images of immunostained microtubules. BS-C-1 cells were fixed with 3% paraformaldehyde and 0.1% glutaraldehyde in phosphate buffered saline (PBS) for 10 minutes at room temperature and then were treated with 0.1% sodium borohydride for 7 minutes to reduce auto-fluorescence. After blocking with 3% w/v bovine serum albumin (BSA) and 0.1% v/v Triton X-100 in PBS, the microtubules were labeled with a mouse monoclonal anti-beta-tubulin primary antibody (Sigma, T5201) followed by a donkey anti-mouse secondary antibody (Jackson Labs, 715-005-150) conjugated with Alexa Fluor 405 and Alexa Fluor 647 [13]. Alexa Fluor 647 is a photoswitchable dye and Alexa Fluor 405 facilitates the photoactivation of Alexa Fluor 647. STORM images were collected on a Olympus IX-71 inverted optical microscope as described previously [24, 25]. Briefly, a 647 nm laser was used to both image Alexa Fluor 647 and switch it off, and a 405 nm laser was used to reactivate the Alexa Fluor 647 to the fluorescent state. The 405 nm laser was used at low intensities such that only a subset of the Alexa Fluor 647 molecules were activated, allowing the locations of these molecules to be determined using either L1H or CVX.

We find that compressed sensing is sensitive to differences between the assumed PSF and the actual PSF; thus, we used an experimentally determined PSF to analyze the real STORM data with the L1H and CVX algorithms. The PSF of the STORM microscope was determined by analyzing the images of well separated single Alexa Fluor 647 molecules. The pixel intensity profile of the images of these molecules was fit with a cubic spline. Five hundred splines were aligned to each other at the centroid positions, and then averaged to create a single cubic spline estimate of the PSF. This spline PSF estimate was then used to create the $A$ matrix used for this data.

4. Results

In theory, interior point methods and L1H should produce identical solutions. Nonetheless, as both algorithms are subject to stability issues and round-off errors, we first tested computationally whether the two algorithms produced identical reconstructions of STORM data. We constructed a set of simulated 7 × 7 camera pixel STORM data and compared the results derived from our L1H algorithm and the algorithm for the interior point method (CVX) [15, 23]. Figure 4(a)-4(c) shows that the two algorithms produce reconstructions in which the vast majority of fluorophores are in the same locations. There are occasionally differences between the two reconstructions, as depicted in Fig. 4(c). However, as shown in Fig. 4(d)-4(f), these differences are exceedingly rare and when they do occur, the reconstructed positions of the fluorophores differ only by one up-sampled pixel (1/8 of the camera pixel). Moreover, we find that the L1H and CVX produce solutions that differ in the L1 norm and in the residual image error by much less than 0.1%, as shown in Fig. 4(g)-4(h). Thus, we conclude that the two algorithms produce functionally equivalent reconstructions.

 

Fig. 4 L1-Homotopy produces nearly identical reconstructions to interior point methods. (a-c; left) Simulated, high-density, 7 × 7 camera pixel single-frame conventional image. The locations of individual emitters are marked with the red crosses. (a-c; middle) The CVX reconstruction of the emitter locations. (a-c; right) The L1H reconstruction of the same images. The arrows in (c) highlight a rare difference between the two solutions. (d) Histogram of the distance between an emitter in the CVX reconstruction and the nearest emitter in the L1H reconstruction for the same simulated data. The histogram is plotted for three different emitter densities. (e) The average distance between every emitter in the CVX solution to the nearest emitter in the L1H solution as a function of emitter density. (f) The average fractional difference in the number of emitters found in the CVX solution and the L1H solution as a function of emitter density. (g) The average percent difference between the L1 norm of the CVX solution and the L1H solution. (f) The average percent difference between the residual image error of the CVX and the L1H solution. (e) – (f) Both the mean (blue) and median (red) are provided. Error bars represent standard deviation measured directly (blue) or estimated from the inter-quartile range (red).

Download Full Size | PDF

To further confirm that L1H produces functionally equivalent reconstructions to CVX, we generated a series of simulated 256 × 256 camera pixel STORM movies across a range of emitter densities and analyzed them with CVX and with L1H. Figure 5 reveals that L1H and CVX produce reconstructions that are effectively identical in the density and localization precision of reconstructed emitters. However, Fig. 6 reveals that L1H produced these reconstructions in roughly 10 – 250 fold less time than CVX. The speed improvement with L1H decreases as the emitter density increases. This density dependence reflects the fact that the computation time of L1H scales with the number of emitters while that of CVX scales primarily with the number of up-sampled pixels. Based on the above results, we conclude that L1H produces the same reconstructions as CVX but in a fraction of the time.

 

Fig. 5 Comparison of the reconstructed emitter density and localization error derived from L1H, CVX, a single-emitter fitting algorithm, and a multi-emitter fitting algorithm. (a) The density of reconstructed emitters as a function of the density of simulated emitters for a single-emitter fitting algorithm (green squares), a multi-emitter fitting algorithm (blue pluses), L1H (red crosses), and CVX (black diamonds). The dashed black line has the slope of 1. (b) The XY localization error for each algorithm labeled as in panel (a). The two panels in (a) and (b) cover different density ranges.

Download Full Size | PDF

 

Fig. 6 L1-Homotopy reconstructs images with a substantially higher speed than interior point methods, but is slower than the emitter fitting algorithms. (a) The average analysis time for a 256 × 256 camera pixel image as a function of emitter density for the two compressed sensing algorithms CVX (black diamonds) and L1H (red crosses) as well as a multi-emitter fitting algorithm (blue pluses) and a single-emitter fitting algorithm (green squares). (b) The ratio of the analysis time per frame for CVX to L1H as a function of emitter density.

Download Full Size | PDF

To determine the computational speed improvement for real data, where emitter density is non-uniform but can vary dramatically within a single image, we analyzed STORM images of immunostained microtubules in mammalian cells. Figure 7(a) shows that for the real STORM data we again find that CVX and L1H produce essentially identical solutions. However using the time it took CVX to analyze the first 10 frames of this 5000 frame movie (which is 256 × 256 camera pixels per frame), we estimate that a complete analysis of all 5000 frames would take roughly two months using a single core of a high end PC. In contrast, L1H was able to analyze the entire 5000 frame movie in just under 3 hours using the same computer. The analysis time for the first 10 frames using L1H was 340-fold faster that using CVX, as shown in Fig. 7(b). The average emitter density in this field of view is approximately 1 μm⁻²; however, the speed enhancement is substantially larger than that observed for the simulated data at a uniform density of 1 μm⁻² because of the large variation in the local emitter density observed in the real STORM data. Since the typical STORM image does not have a uniform emitter density, we anticipate this experimentally-derived speed enhancement to be representative of what will be achieved in the analysis of a wide range of real STORM data.

 

Fig. 7 L1-Homotopy analysis of experimental STORM data. (a) A sub-area of a single frame of a high-emitter-density STORM data set acquired from Alexa-647-labeled microtubules in BS-C-1 cells. Individual molecules found with L1H (green circles), CVX (blue crosses), and a single-emitter localization algorithm (red crosses) are plotted. (b) Average analysis time per frame (256 × 256 camera pixel) for L1H and CVX estimated from the analysis time for the first 10 frames of a 5000 frame STORM movie of microtubules in a BSC-1 cell. CVX takes 340-fold longer than L1H to analyze these frames. (c) The full reconstructed STORM image of this data set using a single-emitter localization algorithm. (d) The same data set reconstructed with L1H. The red arrows and arrow heads indicated two regions with high and low microtubule density, respectively. (e) A zoom-in of the area outlined by the red box in (c). (f) A zoom-in of the same area outlined by the red box in (d). The image reconstructed by the L1H compressed sensing algorithm is much smoother than that by the single-emitter localization algorithm because roughly 4-fold more fluorophores are localized by the L1H algorithm. Scale bars in (c,d) are 10 µm and 1 µm in (a, e, f).

Download Full Size | PDF

The results in Fig. 7 also confirmed the previous conclusion that analysis of relatively high emitter density data by compressed sensing yields superior results compared to analysis by a single-emitter fitting algorithm [24], where overlapping emitters are discarded. Because of the relatively high emitter density of this data set, the single-emitter localization algorithm only identified 23% of the localizations that L1H identified. Figure 7(c)-7(f) show that the single-emitter algorithm produced an image with reduced sharpness, continuity, and contrast as compared to the reconstruction produced with the compressed sensing approach using the L1H algorithm. Moreover, the single-emitter fitting algorithm produces an image with a lower dynamic range, which largely evened out the intensity differences between regions of high and low microtubule densities whereas such difference are obvious in the images generated by the L1H algorithm. (For example, compare the regions highlighted by the red arrows and arrow heads). Thus, by resolving partially overlapping emitters, the compressed sensing approach to STORM analysis can produce reconstructed images with improved sharpness, contrast, and dynamic range. It is worth noting that the results from compressed sensing based analysis are sensitive to differences between the assumed PSF and the actual PSF. In real samples, the PSF from out of focus emitters will naturally be different than that of in focus emitters and hence can cause localization errors when a fixed PSF is used in the analysis. Despite this concern, Fig. 7 reveals that the compressed sensing reconstruction of real data with some natural variation in the focus of the emitters is still superior to that produced by the single-emitter fitting algorithm.

To further compare the performance between the compressed sensing approaches and other STORM image reconstruction algorithms, we also analyzed the 256 pixel × 256 camera pixel simulated data with a single-emitter localization algorithm that can only handle sparse emitter densities and with the multi-emitter fitting algorithm which can handle moderate-to-high emitter densities [9, 13]. The comparison of the performance of these algorithms as a function of emitter density in Figs. 5 and 6 illustrates the benefits and potential limitations of a compressed sensing approach with respect to other algorithms. Above a density of ~0.3 emitters per μm⁻², the single-emitter localization algorithm begins to miss a substantial fraction of emitters, and above an intermediate density of ~3 emitters per μm⁻² the multi-emitter fitting algorithm also begins to miss a substantial number of emitters. However, L1H and CVX both reconstruct the vast majority of the emitters up to a density of 5-6 emitters per μm⁻², as shown in Fig. 5(a). Thus, for the highest emitter densities, a compressed sensing approach is required to reconstruct most emitters.

We also quantified the accuracy with which these emitters are reconstructed, as seen in Fig. 5(b). At low densities, the multi-emitter fitting algorithm has the lowest localization error because it employs maximum likelihood estimation as opposed to the least squares minimization [27] employed by our single-emitter fitting algorithm. This algorithm continues to outperform compressed sensing in terms of localization error up to an intermediate density of 3-4 emitters per μm⁻². Beyond this density, a compressed sensing approach is superior. These results are consistent with those reported previously [15]. However it is important to note that this improvement in high-density image reconstruction comes at a cost. Compressed sensing requires slightly more analysis time than multi-emitter fitting, even when using the significantly faster L1H algorithm. See Fig. 6(a).

Finally, to qualitatively capture and summarize these benefits, we present the results on a simulated STORM data with variable emitter density derived from several different STORM image analysis methods, including a single-emitter fitting algorithm, DAOSTORM (a multi-emitter fitting algorithm), deconSTORM (a recently reported deconvolution analysis method [14],) and L1H. As seen in Fig. 8, in the regions of low density all algorithms produce similar quality images, but in regions of high emitter densities, the single- and multi-emitter fitting algorithms miss most of the emitters. deconSTORM and compressed sensing retain these emitters in the high-density regions and, thus, more faithfully reflect the dynamic range of the image. In addition, compressed sensing produces the sharpest reconstructed images within the high-density regions. This ability to accurately process data with a large range of emitter densities is a key advantage of compressed sensing approaches. Thus, while compressed sensing is still a computationally expensive analysis method, with the use of L1H it should now be practical to routinely apply this method. Moreover, compressed sensing is a versatile analysis method, and it should be possible to modify the A matrix to include PSFs of variable size and shape, which will in turn extend compressed sensing to the analysis of three dimensional data. The significant decrease in computational complexity afforded by L1H will partially offset the substantial increase in computational complexity of 3D implementations of compressed sensing. Finally, it may be possible to further increase computational speed by implementing L1H on GPUs or clusters of CPUs.

 

Fig. 8 A comparison of different analysis methods on simulated STORM data. (a) A single frame of the STORM movie. (b) The true locations of the emitters used for the simulation. (c) The STORM image reconstructed using a single-emitter fitting algorithm that can only handle sparse emitter densites (d) The STORM image reconstructed using a multi-emitter fitting algorithm, DAOSTORM, that can handle moderate-to-high emitter densities. (e) The STORM image reconstructed using the deconSTORM algorithm. (f) The STORM image reconstructed using the L1H algorithm. All scale bars are 500 nm. The simulated STORM movie had an average density of 5 emitters/μm² and was 200 frames long.

Download Full Size | PDF

5. Conclusions

In conclusion, the speed improvement enabled by the L1H algorithm makes it more practical to apply compressed-sensing-based image analysis to high emitter density data. The high density of fluorescent emitters that can be localized should substantially reduce the number of imaging frames required to reconstruct a super-resolution image, increasing the time resolution of STORM and related super-resolution imaging methods, and benefit the study of ultra-structural dynamics inside cells.