This paper proposes a cascading algorithm (CSR) based on compressed sensing, which aims to reduce intensive computations in super-resolution imaging of fluorescence microscopy. Performance of existing algorithms such as CVX and L1H drop sharply when applied to obtain finer images with high density molecules. CSR fully exploits the extreme sparsity property of molecules in the compressed sensing model and progressively restricts solution space stage by stage. We perform a comprehensive study of existing algorithms and the proposed algorithm under different resolutions and molecules’ densities. Simulation and experimental results confirm the performance advantage of CSR when applied to recover dense molecules.
© 2015 Optical Society of America
Light microscopy, as a tool of observing the micro-world, gradually cannot satisfy the need in certain biological fields because of light diffraction. The true scale of molecules often ranges around a few nanometers  while microscopy can only probe several hundred nanometers. The scale gap between captured images and actual images leads to inaccurate judgement on the observed samples with a lot of useful information missing. Mathematically, point sources of molecules are diffracted into airy disks by point-spread functions (PSFs) . These airy disks often overlap, especially for high density molecules. It’s a big challenge to distinguish molecules among overlapped disks. In addition, background and photon counting noise further complicate camera images. The combination of diffraction, overlapping and noise makes super-resolution imaging of molecules an even more challenging task.
STORM (Stochastic Optical Reconstruction Microscopy)  and (F) PALM ((Fluorescence) Photo-Activated Localization Microscopy) [4, 5] solve the overlapping problem to a large extent. They exploit the photo-switchable characteristics of the fluorescent material such as Alexa Fluor 647 and stochastically control that only a small fraction of molecules are activated into the light state at one point of time. In order to reconstruct super-resolution images, single molecule localization (SML) studies diffraction patterns and localizes molecules’ positions. Various SML algorithms related with STORM/ (F) PALM are studied for faster or more accurate reconstruction. Fitting-based algorithms [6–8] identify low density molecules with very fast speed but fail to recover high density molecules. CSSTORM  proposes a compressed sensing [10, 11] based model and solves the problem with the MATLAB package CVX , which recalls at least three times as many molecules as fitting-based algorithms but with long computation time. L1H  adopts an efficient algorithm L1 homotopy (L1H)  to reduce the computation time of the proposed model. However, it is still not satisfactory for live cell imaging, especially when applied to recover finer-resolution images. In order to improve localization accuracy, Barsic et al. present a two-stage algorithm MP + CO, which uses the matching pursuit algorithm (MP) to localize molecules’ coarse positions and convex programming(CO) for fine images . Although effective, little has been discussed regarding to algorithmic details of the two-stage technique.
With the purpose of faster super-resolution imaging of dense molecules, this paper proposes a cascading multi-stage algorithm, Cascading Super-Resolution (CSR). We observe that considerable superfluous cost has been wasted on computing no-molecule areas because of ultra imbalance between extremely sparse molecules and large solution space in the compressed sensing model . Here, the sparsity means that the number of molecules is much smaller than the dimension of solution space. Even for high density molecules, this sparsity still stands. Meanwhile, there is no prior information about molecules’ positions available. The main idea of CSR lies on removing the no-molecule areas as many as possible and providing coarse positions of molecules as prior information. The direct super-resolution reconstruction with heavy computational cost is substituted by multiple stages with significantly reduced cost.
The CVX  and L1 homotopy  algorithms are integrated into the proposed CSR framework. The CVX algorithm initially computes the coarsest-resolution images and L1H further solves fine molecule locations in cascading stages. The cascade between adjacent stages is realized with two efficient techniques. CSR and existing algorithms are compared with different molecule densities and resolutions. The simulation and real data experiments verify that CSR has significantly better performance compared with CVX and L1H and higher reconstruction quality than L1H in super-resolution imaging of dense molecules.
This section first presents the problem models used for different stages in the proposed cascading algorithm CSR and then explains the algorithm with an example. Finally, details of the proposed algorithm are presented.
2.1. Problem models
In the model of CSSTORM , the relationship between true images and sampled images is supposed to be linear. Formally, it solves the following optimization problem with inequality constrains:Equation 1 can be solved by the interior-point method in CVX  but with high computational cost.
L1H traces a parameter-controlled solution path, which can be seen as a piecewise function with several discrete breakpoints. The parameter λ helps determine where these breakpoints are. Then, we just need to check whether solutions in these breakpoints satisfy the given accuracy. It maintains a k-step solution property when the density k is not so high . Hence, it is efficient when k-step property holds but still computationally intensive when this property fails. High density molecules are likely to cause failure of the k-step property.
2.2. An example
CSR is briefly explained with a simple example in the following. RF is assumed to be 8, which means that each raw pixel should be expanded into 8×8 grids. CSR decomposes 8=2×2×2 by prime factorization and substitutes three stages with incremental subRFs of 2, 4, and 8 respectively for direct reconstruction.
Figure 1 illustrates the iterative process of one raw pixel in CSR. There are three molecules inside the raw pixel. Initially, the pixel is divided into 2×2=4 grids with subRF equal to 2. Two pixels in diagonal (light grey boxes in Fig. 1(b)) are identified by the pre-reduction algorithm, and considered as the prior information of molecules’ approximate locations with coarse resolution. The solution space is reduced in half. With subRF equal to 4, three pixels (dark grey boxes in Fig. 1(c)) are further located by the first cascading control and segmented into 12 grids for the next stage with RF equal to 8. Now, the whole space with 64 pixels are reduced to 12 pixels with more than 80% eliminated. Finally, locations of the three molecules (dark boxes in Fig. 1(d)) with the required resolution are found by the second cascading control.
2.3. The cascading algorithm CSR
Details of the cascading algorithm CSR are presented in Algorithm 1. In general, it contains three stages: prime factorization, pre-reduction and cascading control. Initially, prime factorization decomposes the refining factor to multiple factors and computes the sub-factors (sub-RFs) with incremental multiplications, as described in Line 1. For example, the refining factor 8 can be decomposed into multiplication of 2, 2 and 2 and then subRFs for multiple stages are 2, 4 and 8 respectively. The smallest subRF is chosen as the refining factor to perform pre-reduction, which primarily computes the coarsest solution, as shown in Line 4. Items of the solution with larger values than the threshold T are refined and mapped into sub-areas Θ, which indicates possible appearance of molecules, as shown in Lines 5–17. Cascading control proceeds to search molecules’ locations within known sub-areas stage by stage until the required resolution is obtained, as shown in Lines 18–21.
Two efficient methods are derived to map the large items of the coarse solution into fine-resolution sub-areas: 2D space enlargement and constrained fine-resolution reconstruction. 2D space enlargement first transforms the coarse vector back to a two-dimensional image. Then, each pixel represents a block of pixels with fine resolution. For example, one coarse pixel with RF equal to 2 represents a block of 2×2 pixels with RF equal to 4. Sub-areas enlarged from the large items are merged with overlapped areas removed. When image borders are encountered, 2D space enlargement makes sure that the enlarged areas do not cross the borders. Constrained fine-resolution reconstruction involves reconstructing images within prior constrained areas. The point-spread functions (PSFs) are computed separately for each stage. Based on the known constrained items obtained from the last stage, only the corresponding columns of the PSF matrix are computed in the cascading control stages, as shown in Line 19 of Algorithm 1. The size reduction of PSF matrix significantly decreases the computational cost on matrix multiplications and reversions. In addition, an extra raw pixel margin is eliminated to avoid false identifications on marginal areas, as shown in Line 24 of Algorithm 1.
3. Simulations and experiments
The cascading algorithm CSR makes use of CVX to obtain the initial approximation of molecules’ locations and progressively reduces time-consuming computations involved in L1 homotopy algorithm . Through theoretical analysis, CSR can perform faster reconstruction of raw images than L1H and CVX. In this section, the effectiveness and efficiency of the proposed CSR technique are verified with simulations and real-data experiments.
A series of STORM raw images are simulated with 22×22 camera pixels. One raw pixel is assumed to be 106 nanometers. Molecules are randomly distributed in fine images. Then, diffracted raw images are obtained by applying Gaussian convolution function with one-raw-pixel width. Uniform 100 photons background and poisson noise are added.
Three algorithms CVX , L1H  and the proposed CSR are implemented in MATLAB. Specifically, since the pixel intensity in STORM images is nonnegative, the L1 homotopy algorithm is modified to solve nonnegative values based on the Python code version in . All simulations and experiments are conducted on a PC platform with an i7-4770 3.4GHz CPU. For the recovery of simulated and real images, similar methods as  are adopted. The whole image frame is computed window by window. Each window contains 7×7 raw pixels with a one-raw-pixel margin. Reconstructed pixels on the margin are removed to avoid inaccurate localizations. Then, only reconstructions of 5×5 raw pixels are valid. There is a two-raw-pixel overlap between adjacent windows. For example, the first window has the overlap of the 6th and 7th pixels with the second window in x axis. There are 4×4 windows to compute for each simulated raw image. The noise weight ε in Eq. 1 is set as ε = 1.5 in simulations and ε = 2.1 in real-data experiments. The performance evaluation involves two metrics: 1) reconstruction accuracy, represented by recall rate and 2) computation time . It is a key measurement for super-resolution imaging whether molecules’ precise locations are found. Recall rate reflects this measurement. It means a successful recall to the molecule when the corresponding location has larger values than the given threshold. Recall rate is calculated as the ratio of recalled molecules within all molecules. The computation time of each algorithm is recorded for each 7×7 window averaged for each frame.
Two groups of raw images are recovered with different densities and refining factors. Given refining factors of 8, 12 and 16, the first group evaluates performance of the three algorithms under densities ranging from 1 molecule/μm2 to 12 molecules/μm2 with an interval of 1 molecule/μm2. The corresponding reconstructed scales are approximately 13 nanometers, 9 nanometers and 7 nanometers. Average results of 30 independent trials are shown in Fig. 2. It is first observed that the computation time of CVX (the blue curve with dots) is almost in parallel with the density axis while that of the other two algorithms, L1H (the black curve with pluses) and CSR (the red curve with stars) increases linearly with molecule density. Although L1H and CSR both take much less time than CVX, the CSR curve has much smaller slope than L1H curve. L1H has better performance when densities and refining factors are low, as shown in Fig. 2(a). CSR begins to catch up with L1H as the density is higher. The crossing point happens at the density of 4 molecules/μm2 for RF equal to 8. There are almost no crossing points for higher refining factors 12 and 16. CSR clearly surpasses L1H and CVX, as shown in Figs. 2(c) and 2(e). In addition, the performance advantage of L1H against CVX gradually drops with the increase of density while CSR maintains a stable superiority. CSR has more performance advantages than L1H and CVX when applied to reconstruct dense molecules with high refining factors.
The recall rate of all algorithms decreases with the increase of density, as shown in Figs. 2(b), 2(d) and 2(f). CVX has the best performance compared with L1H and CSR. L1H presents slightly higher recall rate than CSR for low density molecules but performs inferior for high density molecules. When molecule density grows to 9 molecules/μm2 with RF equal to 12, CSR identifies more molecules than L1H, as shown in Fig. 2(d). The crossing points of CSR and L1H differ for refining factors 8 and 16, as shown in Figs. 2(b) and 2(f). One instance of the reconstructed results is chosen from the simulation and shown in Fig. 3(a). The molecule density is approximately 5 molecules/μm2. Compared with the real image, CVX reconstructs the proximal image while the other two algorithms have differences with CVX in some points. Although CSR does not recall the point indicated by the blue arrow, it shows the ability to control the false positive identifications because of the CVX algorithm applied in the pre-reduction stage. As shown in Fig. 3(a), the point in the blue circle is identified as two points by L1H while CVX and CSR recover it as one.
The performance of CSR under certain refining factors has been demonstrated and analyzed. Next, in order to comprehensively study the effect of refining factors, the other group of simulated raw images are recovered with refining factors ranging from 2 to 18, given a density of 5 molecules/μm2, as shown in Figs. 3(b) and 3(c). Higher refining factors are not chosen in order to trade off the recall rate. When refining factors are too high, it may sacrifice the reconstruction quality [18, 19]. As shown in the figures, the performance of CSR under some refining factors such as 8, 12, 16 and 18 shows the same results to Fig. 2. Under the other refining factors such as 5, 7, and 11, behaviors of the CVX algorithm applied in the pre-reduction stage are more or less maintained by CSR. This reflects the nature of the cascading algorithm, which will be discussed later.
In order to visually display the benefit of the proposed cascading algorithm, a three-hollow-ring pattern image is simulated with high density in the crossing regions of rings and low density in the other regions. Two configurations are set: relatively low average density with RF equal to 8 and high average density with RF equal to 16. The first configuration includes 50 frames with 32×32 camera pixels. For each frame, 50 molecules are randomly placed onto the areas of three hollow rings. Assuming 106nm camera pixel size, the radius and width of each ring are 715nm and 80nm respectively. Molecule density for each frame is estimated as 5 molecules/μm2 on average. There are approximately 2500 molecules in total. The second configuration contains 100 frames with 17×17 camera pixels. Each frame has 20 random molecules on the ring pattern and density around 8.8 molecules/μm2 on average. There are approximately 2000 molecules in total. The radius and width of each ring are 357nm and 40nm respectively.
One raw frame and the real image are shown in Fig. 4(a). The reconstructed images by algorithms CVX, L1H, CSR and single-molecule fitting algorithm levenberg-marquardt (LM)  are displayed in Fig. 4(b). Here, the open source tool RapidSTORM  is used to calculate results of the LM algorithm. From the result shown in the first row of Fig. 4(b), three hollow rings by algorithms CVX, L1H and CSR can be clearly identified with RF equal to 8 while RapidSTORM only gives a pattern of three circles with vague hollow areas under the same RF. CSR and L1H performs slightly weaker than CVX in the recovery of three-ring crossing areas indicated by the green arrows and two-ring crossing areas indicated by the white arrows. Under the higher resolution with RF equal to 16, shape of the pattern recovered by RapidSTORM is completely blurry with imperfect circles, as shown in the second row of Fig. 4(b). With CVX, L1H and CSR, the whole ring pattern can be identified although with fuzzy hollow areas. The crossing areas of three rings with very high molecule density are vague. As to the areas of two ring crossing, CVX provides the clearest image while the other two algorithms partially identify the details. However, the image by CSR has slightly sharper edges than L1H in these areas.
The reconstruction time and recall rate of the simulated three-hollow-ring pattern are also recorded for each algorithm in Table 1. Speed gain ratio represents the ratio of the algorithm’s speed gain against CVX. The time of RapidSTORM is not compared because these two methods will not likely be applied to the same data. RapidSTORM is much faster  (several milliseconds to reconstruct the simulated dataset), but it does not yield accurate results in high-density experiments, as shown in Fig. 4(b). Similar conclusions can be derived from the table. CVX takes much longer time than CSR and L1H. The performance advantage of L1H against CVX decreases while advantages of CSR become stronger with the increase of the refining factor. CSR surpasses L1H for higher density and RF. As to the recall rate, similar conclusions can be drawn that RapidSTORM recalls just a small fraction of molecules while CVX, L1H and CSR identify the majority of them. CVX recovers the most molecules compared with L1H and CSR. Although CSR has lower recall rate than L1H for low density configuration, it performs better reconstruction for dense molecules with higher refining factor.
3.2. Experiments with real data
To test the performance of the proposed CSR algorithm on real data, we collected images of immunostained mitochondrial sub-organelle structures in cos-7 cells. Cos-7 cells were grown on 18mm ♯1 coverslips (Marienfeld). After fixation with 3% paraformaldehyde and 0.05% glutaraldehyde in phosphate buffered saline (PBS) for 10 minutes at room temperature, cells were treated with 100 mM glycine in PBS for 7 minutes to reduce the unreacted aldehyde groups formed during fixation. After 15-minute incubation in blocking buffer (3% w/v normal donkey serum, 0.2% v/v Triton X-100 in PBS), the mitochondrial outer membrane protein Tom20 was labeled with a rabbit polyclonal anti-Tom20 primary antibody (Santa Cruz). The samples were then rinsed with PBS three times. Corresponding secondary antibodies conjugated with Alexa Fluor 647 (Jackson Labs) were added to the samples in blocking buffer and incubated for 30 minutes, followed by rinsing thoroughly with PBS.
Cell imaging was performed in a standard imaging buffer that contains 50 mMTris, pH 8.0, 10 mMNaCl, 0.56 mg/mL glucose oxidase, 40μg/mL catalase, 10% (w/v) glucose, 2 mMcyclooctatetraene, 50 mM β-mercaptoethanol. The STORM microscope was built upon a Nikon Ti-E inverted microscope with perfect focus system (PFS). Laser excitation for Alexa 647 channel was provided by a 656.5nm DPSS laser (CNI MXL-N-656.5). The fluorescent signal was collected by a 100× TIRF objective (Nikon), and was recorded using EMCCD (AndoriXon-Ultra) with 100× EM gain with 30 milliseconds exposure time and at 30 frames per second.
The whole dataset contains 5000 frames with the size of 325×145 camera pixels. For each frame, the molecule density is high for each mitochondria. 300 frames of one image patch with 100×100 pixels are reconstructed using RapidSTORM , CVX and the proposed algorithm CSR with refining factors of 8 and 16. The reconstructed pixel scales are around 13nm and 7nm for the raw pixel size of 106nm. With the purpose of reducing effect of background, each pixel is subtracted with 200 photons. Image results and analysis are shown in Fig. 5 and Fig. 6 respectively. One raw frame is shown in Fig. 5(a), from which no details can be seen. Recovered images are displayed in Figs. 5(b) and 5(c). Although the speed to reconstruct dense molecules is improved by CSR, it still takes more than 10 days to recover the whole 5000 frames while RapidSTORM reconstruction just costs less than 10 seconds. The full dataset of 5000 frames is recovered with RapidSTORM and shown in Fig. 5(d).
The time for CVX, L1H and CSR to reconstruct 13nm resolution images are 165s, 110s and 31s per frame respectively. L1H takes around 67% time of CVX while CSR further improves 72% timing cost of L1H, as shown in Fig. 6(a). When they are performed to recover 7nm resolution images, the timing costs of CVX, L1H and CSR are 651s, 379s and 41s for each frame respectively. Obviously, finer-resolution reconstruction requires higher timing cost. The improvements of L1H and CSR against CVX are 62% and 94% respectively. The speed advantage of L1H decreases while advantage of CSR performs better for recovering higher resolution images. This result reconfirms the conclusion in simulations. CVX and CSR recover more continuous images than RapidSTORM. The profile of an image line, indicated by the white line in Fig. 5(c) for each algorithm has been plotted in Fig. 6(b). These results show that slightly degraded images can be reconstructed by CSR compared with CVX but with much faster speed.
We discuss two aspects of the cascading algorithm. 1) Decomposability of the refining factor. In the prime factorization of CSR, the required refining factor is decomposed into a number of subRFs. The number of subRFs is defined as decomposability. Under refining factors with lower decomposability, the advantage of the cascading algorithm cannot be fully exploited.
Hence, the results of CSR in these RFs show no significant difference with CVX because they still largely reflect performance of CVX used in pre-reduction. The choice of proper resolutions involves many factors, such as photon counts of molecules, imaging systems and so on [18, 19]. In practice, they can be adjusted accordingly to improve the decomposability. 2) Effect of algorithms in pre-reduction and cascading control stages. The performance of CSR depends on the algorithms chosen in each stage. The disadvantage of the algorithms also results in performance degradation of CSR. We make use of the speed advantage of L1H, but cannot avoid the quality degradation of reconstructed images. We have considered the other choices as well. Initially, we tried to apply CVX in both pre-reduction and cascading control stages. However, when applied to the cascading control stage, CVX sometimes leads to infeasible values within partial matrix columns. Hence, we adopt L1H for cascading control stages. The performance degradation induced by L1H is acceptable in practice.
This paper proposes an efficient algorithm CSR to speed up super-resolution imaging of dense molecules in fluorescence microscopy. CSR combines the advantages of CVX and L1H. It progressively increments the refining factor and restricts the large solution space in multiple stages. The experimental results show better performance and reconstruction quality of CSR against L1H for higher density molecules and finer resolutions.
We thank Nanobioimaging Limited for assistance with collecting real STORM mitochondrial data.
References and links
1. Nature Education, “The relative scale of biological molecules and structures,” http://www.nature.com/scitable/content/the-relative-scale-of-biological-molecules-and-14704956.
2. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (CUP Archive, (1999). [CrossRef]
4. S. T. Hess, T. P. K. Girirajan, and M. D. Mason, “Ultra-high resolution imaging by fluorescence photoactivation localization microscopy,” Biophys. J. 11, 4258–4272 (2006). [CrossRef]
5. E. Betzig, G. H. Patterson, R. Sougrat, O. W. Lindwasser, S. Olenych, J. S. Bonifacino, M.W. Davidson, J. Lippincott-Schwartz, and H. F. Hess, “Imaging intracellular fluorescent proteins at nanometer resolution,” Science 5793, 1642–1645 (2006). [CrossRef]
6. A. Yildiz, J. N. Forkey, S. A. McKinney, T. Ha, Y. E. Goldman, and P. R. Selvin, “Myosin V walks hand-over-hand: single fluorophore imaging with 1.5-nm localization,” Science 5628, 2061–2065 (2003). [CrossRef]
8. F. Huang, S. L. Schwartz, J. M. Byars, and K. A. Lidke, “Simultaneous multiple-emitter fitting for single molecule super-resolution imaging,” Biomed. Opt. Express 2, 1377–1393 (2011) [CrossRef] [PubMed]
10. E. J. Candés, Y. C. Eldar, D. Needell, and P. Randall, “Compressed sensing with coherent and redundant dictionaries,” Appl. Comput. Harmon. A. 31, 59–73 (2011). [CrossRef]
11. E. J. Candés, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE T. Inform. Theory 52, 489–509 (2006). [CrossRef]
12. M. Grant, S. Boyd, and Y. Ye, “CVX: Matlab software for disciplined convex programming,” http://cvxr.com/cvx/.
13. H. P. Babcock, J. R. Moffitt, Y. Cao, and X. Zhuang, “Fast compressed sensing analysis for super-resolution imaging using L1-homotopy,” Opt. Express 21, 28583–28596 (2013). [CrossRef]
14. D. L. Donoho and Y. Tsaig, “Fast solution of L1-norm minimization problems when the solution may be sparse,” IEEE T. Inform. Theory 54, 4789–4812 (2008). [CrossRef]
16. M. S. Asif and J. Romberg, “L1 Homotopy,” http://users.ece.gatech.edu/sasif/homotopy/
17. S. Wolter, U. Endesfelder, S. Van de Linde, and M. Heilemann, “Measuring localization performance of super-resolution algorithms on very active samples,” Opt. Express 19, 7020–7033 (2011). [CrossRef] [PubMed]
18. S. Ram, E. S. Ward, and R. J. Ober, “Beyond Rayleigh’s criterion: a resolution measure with application to single-molecule microscopy,” in Proceedings of the National Academy of Sciences of the United States of America (2006), pp. 4457–4462. [CrossRef]
20. D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” J. Soc. Ind. Appl. Math. 11, 431–441 (1963). [CrossRef]
21. S. Wolter, A. Löschberger, T. Holm, S. Aufmkolk, M. C. Dabauvalle, S. van de Linde, and M. Sauer, “Rapid-STORM: accurate, fast open-source software for localization microscopy,” Nat. Methods 9, 1040–1041 (2012). [CrossRef] [PubMed]