Efficient image reconstruction of high-density molecules with augmented Lagrangian method in super-resolution microscopy

Jia Li; Danni Chen; Junle Qu

doi:10.1364/OE.26.024329

1. Introduction

Super-resolution fluorescence microscopy has been considered a powerful technique to bring “optical microscopy into the nanodimension”. It effectively breaks the limitation in resolution capacity of conventional microscopy, and opens more promising opportunities accessible to life fundamental processes, such as DNA transcription, RNA translation and protein folding [1–3].

As an integral part of super-resolution microscopy methods, stochastic optical reconstruction microscopy (STORM) firstly takes advantage of photoswitchable fluorophores which can alternate between the activated state and the dark state to form raw images, and then constructs super-resolution image by localizing and merging the positions of sparsely distributed fluorescent molecules [4]. However, thousands of frames of raw images are generally required to achieve high spatial resolution, resulting in a long acquisition time and limiting the wide adoption of STORM for live-cell imaging.

To improve temporal resolution, one feasible way is to increase the density of activated fluorophores in each raw image such that fewer raw images are required [5]. This invalidates commonly-used single-molecule localization methods; accordingly, several algorithms have been developed to cater for localizing closely spaced molecules [6]. For fitting-based methods, represented by dominion astrophysical observatory stochastic optical reconstruction microscopy (DAOSTORM) [7] and multi-emitter fitting method [8], which simultaneously fit multiple overlapping point spread functions (PSFs) to the data, and determine the positions of molecules by minimizing mismatch between the data and the model [9]. Another category is deconvolution-based methods that typically pose localization as an image estimation problem. They give density profiles of fluorophores by employing tools from different fields, like Richardson-Lucy deconvolution [10], compressed sensing (CS) [11, 12], Bayesian analysis of the blinking and bleaching method [13]. Density is estimated on a sub-pixel discrete grid that is smaller than the camera pixel size, thereby reconstructing a higher resolution image. Various improvements of these algorithms have been achieved by more sophisticated modeling and optimization methods [14], including fast localization algorithm based on a continuous-space formulation (FALCON) [15] and 3D super-resolution imaging [16].

Unfortunately, the ever growing density of fluorophores caused by the demand for live-cell imaging may not fit the application scenario of compressed sensing analysis, necessitating algorithms that are not based on sparse recovery. To address the issue, we propose an efficient reconstruction algorithm based on augmented Lagrangian methods (ALMs) that does not have requirements for PSF or the degree of sparseness of fluorophore distribution in each frame. ALMs are a certain class of algorithms for solving constrained optimization problems, which replace the original optimization problem with a series of unconstrained problems. They can reduce the possibility of ill conditioning of the generated sub-problems by the inclusion of the Lagrange multiplier in the augmented Lagrangian function [17]. Recent developments that incorporate sparse matrix techniques and the use of partial updates have rekindled a lot of interest in this approach [18–20].

This work focuses on developing an image reconstruction algorithm suited to high-density fluorescent molecules, and the major contributions are twofold. First, the proposed scheme incorporates total variation (TV) norm through regularization term, resulting in the high-resolution image that is not only robust against noise, but also achieves a favorable tradeoff between localization accuracy and detection rate, thereby shortening the overall acquisition time for densely distributed fluorophores. We demonstrate the better algorithmic performance using fewer frames of simulated and real STORM data. Second, in terms of computation time, we achieve a higher speed compared to representative CS method. This is fulfilled by introducing quasi-Newton method and alternating direction method to solve sub-problems after constructing the augmented Lagrangian function.

2. Problem formulation

The imaging model in STORM can be described as a convolution of fluorophores distribution with the PSF of the optical system. We use the discrete form for image computation; furthermore, for ease of description, we turn an image of size $N \times N$ into vector of length $N^{2}$ through lexicographic ordering. In what follows, $s$ denotes the desired high-resolution grid of fluorophores locations. Mathematically, the vector notation of the raw camera image $b$ is given by

b = A s + n,

where

A

and

n

represent the convolution operator and the noise, respectively. In our work, the matrix

A

is a

p \times p

pixels symmetric Gaussian function that models the system PSF and is calculated as

A = \frac{N \times \exp [- 4 \ln 2 (\frac{x^{2} + y^{2}}{σ^{2}})]}{\int_{- \frac{p}{2}}^{\frac{p}{2}} \int_{- \frac{p}{2}}^{\frac{p}{2}} \exp [- 4 \ln 2 (\frac{x^{2} + y^{2}}{σ^{2}})] d x d y} .

We take into account the pixelation noise in Eq. (2), and $N$ indicates the photon number. With the known wavelength $λ$ , the numerical aperture NA and the pixel width $WD$ , the full width at half maximum (FWHM) of PSF equals to $σ = λ / (2 \times NA \times WD)$ . Note $A$ also consists of an under-sampling operation. Compared to the measurement $A$ in [11], our convolution matrix $A$ offers mathematical simplicity and maintains good performance.

Given the $i th$ frame of the measured low-resolution images $b_{i}$ and the PSF of the imaging system, the corresponding unknown image $s_{i}$ containing information of the molecule positions can be reconstructed exactly or approximately by solving its regularization function $ϕ (\cdot)$

\begin{array}{l} minimize ϕ (s_{i}) \\ subject to b_{i} = A s_{i} + n . \end{array}

We solve the above problem for each frame independently, and then join all detected molecules $s = \sum_{i} s_{i}$ to obtain the reconstructed high-resolution image. In the following equations, we drop the indices $i$ for simplicity.

For the purpose of discussing the algorithms below, we define the following operators. The first one, $V (\cdot)$ , represents the vectorization of a matrix, which converts it into a column vector using a lexicographical order. Next are $D_{x}$ and $D_{y}$ , both denote the first-order forward finite difference operators, defined respectively as

D_{x} (s) = V (S_{x} - S) and D_{y} (s) = V (S_{y} - S),

in which

S_{x}

and

S_{y}

mean shifting

S

along the horizontal and vertical directions by one pixel, respectively so that

D = {[D_{x}, D_{y}]}^{T}

. The periodic boundary condition is used here.

3. Augmented Lagrangian method for super-resolution microscopy

In this section, we describe how to apply the augmented Lagrangian method to the image reconstruction of high-density fluorescent molecules. We first specify the original problem in Eq. (3) by designing a proper cost function and transform it into an unconstrained one. The resulting unconstrained optimization problem is then converted to a series of sub-problems by variable splitting. Thus, the solution of each sub-problem can be found by using corresponding algorithm.

3.1 The cost function

Image reconstruction in fluorescence microscopy is an inverse problem where the goal is to recover a high-resolution image from a series of blurry and noisy observation. Just like other inverse imaging methods, the design of a proper cost function is crucial. Generally speaking, the cost function consists of an image fidelity term and one or more regularization terms for specific purpose [21]. For the former, the deviations of the optimization result from the raw camera image is measured by the $l_{2}$ norm, i.e.,

F (s) = {‖ A s - b ‖}_{2}^{2} .

As for the regularization term, we constrain the total variation (TV) of the reconstructed image. Total variation is a popular and powerful technique in image reconstruction and image denoising [22, 23]. TV promotes piece-wise constant solutions, and it is suitable as long as the object to be reconstructed adheres to this prior condition. For molecule detection problem in fluorescent microscopy, the sharp features of image are composed of the activated fluorophores, while most dark regions are flat. We use TV norm due to its remarkable effectivity at simultaneously preserving the important image detail such as sharp features while smoothing away noise in flat regions, even at low signal-to-noise ratios. The molecule detections are hence promoted by reducing the total variation of the image subject closer to the original image. Additionally, current compressed sensing methods incorporate sparsity priors of the original image as TV regularization [24]. A similar approach has been previously proposed for laser scanning microscopy [25]. In this work, the anisotropic total variation norm is minimized to suppress the noise amplification, and we therefore have

R (s) = {‖ D (s) ‖}_{1} .

Combining Eqs. (5) and (6), the overall cost function in this paper is formulated as

C = F (s) + R (s) = \frac{μ}{2} {‖ A s - b ‖}_{2}^{2} + {‖ D (s) ‖}_{1} .

where

μ

is parameter for the tradeoff between the image fidelity term and the regularization term.

3.2 Optimization flow

We first transform Eq. (7) to an equivalent constrained problem by creating an intermediate variable $u$ . This procedure is often called variable splitting and has been recently used in several image processing applications [26, 27]. The rationale behind the variable splitting method is that it can be easier to solve the constrained problem than to solve its original unconstrained counterpart [18]. Moreover, since $s$ is non-negative, it is added as a constraint to the optimization. This leads to the following constrained problem

\begin{array}{l} \underset{s}{minimize} \frac{μ}{2} {‖ A s - b ‖}_{2}^{2} + {‖ u ‖}_{1} \\ subject to u = D (s), \\ s \geq 0 . \end{array}

In order to solve this optimization problem with the augmented Lagrangian scheme, we then derive the augmented Lagrangian function

L_{ρ} (s, u, d) = \frac{μ}{2} {‖ A s - b ‖}_{2}^{2} + {‖ u ‖}_{1} - d^{T} [u - D (s)] + \frac{ρ}{2} {‖ u - D (s) ‖}_{2}^{2} .

In Eq. (9), $d$ is the Lagrange multiplier associated with the constraint $u = D (s)$ , and $ρ > 0$ is the penalty parameter which determines the optimal step size to update $d$ . The Lagrange multiplier is introduced in mathematical optimization to find the local maxima or minima of a multivariable function $f (x_{i})$ subject to constraints $g_{i} (x) = c_{i}$ . ALM performs alternately between the primal minimization step and the dual ( $d$ ) updating step until the termination criterion is satisfied. Each step is accomplished by keeping the other one fixed. Here, because the primal is separated into variable $s$ and the auxiliary variable $u$ , we apply an algorithm known as the alternating direction method (ADM) to solve these sub-problems iteratively:

\begin{array}{l} s_{k + 1} = \arg \min_{s} \frac{μ}{2} {‖ A s - b ‖}_{2}^{2} - d_{k}^{T} [u_{k} - D (s)] + \frac{ρ}{2} {‖ u_{k} - D (s) ‖}_{2}^{2} \\ subject to s \geq 0, \end{array}

u_{k + 1} = \arg \min_{u} {‖ u ‖}_{1} - d_{k}^{T} [u - D (s_{k + 1})] + \frac{ρ}{2} {‖ u - D (s_{k + 1}) ‖}_{2}^{2},

d_{k + 1} = d_{k} - ρ [u_{k + 1} - D (s_{k + 1})],

where the subscript

k

denotes the

k th

iteration.

We now investigate the following three sub-problems.

1. $s$ -subproblem: We need to compute the derivative of $L_{ρ} (s, u, d)$ with respect to $s$ , finding the minimum when the derivative is equal to zero. This give rise to the formulation $μ A^{T} A s + ρ D^{T} D (s) = μ A^{T} b - D^{T} (d_{k}) + ρ D^{T} u_{k} .$
Since the convolution represented by $A$ in STORM may be a non-circulant matrix, and even for sometimes $A^{T} A$ might be circulant [28], it cannot be always diagonalized using the discrete Fourier transform. It should also be noted that Eq. (13) involves non-smooth terms and bound constraints, and as a result, it is not trivial to achieve an analytical formula for $s$ . As used in our previous work [29], we choose to solve Eq. (10) iteratively with an optimization technique called the L-BFGS-B algorithm, which stands for the Limited-memory Broyden-Fletcher-Goldfarb-Shanno method with simple Bounds on the variables [30]. This method is particularly suitable for optimization problems with a large amount of variables because of its limited memory requirement. More importantly, as a quasi-Newton optimization method, its superior convergence property and higher accuracy makes our proposed method promising in practical applications of super-resolution fluorescence microscopy for live-cell imaging. The reason for the choice of variable splitting and L-BFGS-B algorithm is that the second-order information is employed, and the convergence is generally faster than the first-order approaches such as fast iterative shrinkage/thresholding (FISTA) [31].
2. $u$ -subproblem: The solution of Eq. (11) can be obtained via the shrinkage formula [22] $u_{k + 1} = \max {| \frac{d_{k}}{ρ} + D (s_{k + 1}) | - \frac{1}{ρ}, 0} ⊙ sgn [\frac{d_{k}}{ρ} + D (s_{k + 1})],$
in which the multiplication $⊙$ and division are element-wise operations, and $sgn (\cdot)$ is the mathematical expression of the sign function, defined as follows:
$sgn (x) = {\begin{array}{l} 1 & x > 0 \\ 0 & x = 0 \\ - 1 & x < 0 . \end{array}$
3. $d$ -subproblem: The Lagrange multiplier $d$ is updated as described in Eq. (12). This is an obvious advantage over the penalty approach, namely that it is not necessary to enforce $ρ$ approach infinity to guarantee convergence for the original optimization problem. Instead, the existence of the Lagrange multiplier enables the penalty parameter to take a relatively smaller value, thereby avoiding ill-conditioning [17].

The pseudo-code in Table 1 describes the general procedure of this proposed algorithm.

Table 1. Pseudo-code of optimization procedure

View Table | View all tables in this article

4. Results

4.1 Simulated data analysis

We validate the performance of the proposed algorithm using both simulated and real STORM data. The images used for numerical simulations are generated by randomly distributing molecules in an area of $5.12 μm \times 5.12 μm$ , with the molecule density ranging from $1 {μm}^{- 2}$ to $20 {μm}^{- 2}$ . The pixel size is $80 nm \times 80 nm$ , and the high-resolution grid is set to $10 n m$ . A Gaussian model PSF with $λ = 670 nm$ and $NA = 1.4$ is adopted, assuming a photon number of 3000 per molecule. To make the simulation more realistic, the simulated images are then corrupted by Poisson shot noise and Gaussian readout noise with a mean value of 0.13 and a variance of 0.02. For the choice of optimal parameters $μ$ and $ρ$ , as illustrated in [22] and [29], the former value weighs the smoothness of the image against the noise level, and the latter plays an important role in controlling the convergence rate, so $μ = 1000$ and $ρ = 0.5$ are selected empirically in our simulations.

Figures 1(a) and 1(b) display true molecule positions at a density of $10 {μm}^{- 2}$ and the low-resolution raw image, respectively. The red square region has a size of $1 μm \times 1 μm$ and contains 11 molecules. The comparison of magnified reconstructed locations with ground-truth of red square region is in lower left corner of Fig. 1(a). Red dots and yellow dots are respectively true molecule positions and reconstructed positions, while two white dots indicate that the detected molecules are overlapped exactly with the ground-truth. It is evident that the result from the proposed algorithm matches well with true positions.

Fig. 1 Simulated STORM image analysis using the proposed algorithm. (a) True molecule positions with magnified reconstructed results for the red square region. (b) Low-resolution raw image.

Download Full Size | PDF

We then compare and analyze the simulation results in terms of the detection rate and localization accuracy. The detection rate use recall and Jaccard index (JAC) as assessment metrics, and the localization accuracy is measured by root-mean-square error (RMSE) which is the distance between reconstructed molecule locations $(x_{i}, y_{i})$ and their matched ground-truth $(x_{0}, y_{0})$ . Since detection rate and accuracy tend to be in opposition, we reach compromise between them by searching within $40 nm$ of the true positions. If the nearest simulated molecule is found, it is considered as true positive (TP); otherwise, it is counted as a false positive (FP). The remaining ground-truth molecules that are not paired with any detected molecule are categorized as false negative (FN) [32]. Recall, JAC and RMSE are hence respectively defined as

recall = \frac{TP}{FN+TP}, JAC = \frac{TP}{FN+FP+TP}, and RMSE = \sqrt{\frac{1}{TP} \sum_{i = 1}^{TP} {(x_{i} - x_{0})}^{2} + {(y_{i} - y_{0})}^{2}} .

For each molecule density, such as a density of $10 {μm}^{- 2}$ in Fig. 1, 100 frames of different molecule positions are generated by randomly distributing molecules. For each frame of molecule positions, we use only one frame of raw low-resolution image to reconstruct the molecule locations, and then we repeat the simulation for all $100$ frames to obtain average measurements. The reconstructed results using our proposed ALM framework are compared to that solved by compressed sensing (CS) method [33]. As illustrated in Fig. 2, the proposed algorithm and CS achieve similar recall and JAC values over the whole range of molecule densities, for instance, two methods produce a recall difference of 0.02 at the density of $19 {μm}^{- 2}$ . In other words, ALM identifies one molecule more than CS method per ${μm}^{2}$ in this density, which signifies that our method has a slightly higher capacity of identifying molecules than CS method. It is worth noting that ALM can provide noticeably improved location accuracy, particularly for densely distributed molecules. When molecule density changes from $5 {μm}^{- 2}$ to $20 {μm}^{- 2}$ , RMSE of the reconstructed image using CS method increases from $22 nm$ to $40 nm$ . In contrast, ALM exhibits smaller localization errors of $19 nm$ to $29 nm$ . This result is related to the fact that ALM can explore a larger solution space of the optimization problem, and therefore, the error distance between the reconstructed locations and their paired true positions is reduced. The RMSE here also accounts for the effect of TV regularization, which minimizes the error by reducing the noise. In addition, the computation time of CS method is more than 4 times longer than our algorithm when analyzing an image like that in Fig. 1. The effective decrease in computational complexity enables ALM to be applied practically in super-resolution fluorescence microscopy.

Fig. 2 Performance comparison of the proposed ALM and CS method on simulated STORM data. (a), (b) and (c) are comparisons of recall, JAC values and localization accuracy, respectively.

Download Full Size | PDF

Another simulation is conducted on a simulated STORM data with variable molecule densities, using our proposed ALM, single molecule localization performed by ThunderSTORM [34] and FALCON. FALCON is a fast sparsity based super-resolution fluorescence microscopy method. We run ThunderSTORM under the default options without additional post-processing, and set FALCON parameter of speed as normal. This can be regarded as further evidence that the reconstructed image quality is indeed improved by ALM. Figure 3(a) shows the ground-truth, and low-resolution raw images are generated in a way similar to the above, except that each molecule has a lower photon number of 2000 and current high-resolution grid is $8 n m$ . Figures 3(c)-3(e) are reconstructed images using 120 frames of raw images with average density of $3 {μm}^{- 2}$ . The correspondingly reconstructed results of yellow squares are given in lower left corner. It can be observed that three methods produce similar quality results in low-density regions, nevertheless, ALM resolves two closely spaced lines better in high-density regions, and reduced the spurious and noisy peaks at the positions between two lines. This is also in accordance with our expectations on TV regularization. Under the circumstance of 30 frames with average density of $12 {μm}^{- 2}$ , we notice that ThunderSTORM and FALCON detects fewer molecules in high-density regions in Figs. 3(f) and 3(g), which fails to reconstruct the line structure, the intersection point of two lines in particular. In comparison with the other two methods, our algorithm can localize these molecules and determine the structure with high accuracy, as shown in Fig. 3(h). This is also consistent with our observation in Fig. 3(b), which depicts the intensity profiles along the cyan line with three methods using different frames of low-resolution images. Although there is discontinuity in the reconstructed image, ALM achieves using fewer frames of raw image to resolve two lines with little noise, while ThunderSTORM and FALCON are unable to reduce the noise effectively during the process of reconstruction when using 120 frames of images and to localize the right molecule positions when using 30 frames.

Fig. 3 Comparison of ThunderSTORM, FACON and ALM on simulated STORM data. (a) True molecule positions. (b) Plots of intensity profiles along the cyan lines. (c) - (e) are reconstructed images using 120 frames of raw images with ThunderSTORM, FALCON and ALM respectively. (f) - (h) are reconstructed images using 30 frames of raw images with ThunderSTORM, FALCON and ALM, respectively. Yellow squares in lower left corners are correspondingly magnified reconstructed results.

Download Full Size | PDF

After evaluating the reconstructed results of different algorithms, we can now assess the impact of the proposed ALM by two image-based criteria: one is image quality, estimated by signal-to-noise ratio (SNR), and the other one is image resolution, measured by means of Fourier ring correlation (FRC). FRC computes the cross correlation between two halves of positions, and describes the length scale below which smaller details are not resolved in the image [35]. SNR is calculated according to the following formula in [32]

SNR = 10 \log_{10} \frac{{‖ s_{0} ‖}^{2}}{{‖ s_{0} - s ‖}^{2}},

where

s_{0}

and

s

are ground-truth image and reconstructed image, respectively.

Table 2 summarizes the measurements of JAC, RMSE, SNR, FRC and execution time on simulated STORM data in Figs. 3(c)-3(h). For different frames, we reconstruct high-resolution image from the three methods. Although they have “good enough” performance for using 120 frames of raw images, they lead to distinguishing outputs with 30 frames, where larger SNR and smaller FRC indicate that image quality and resolution are improved by our method. As listed in Table 2, the proposed ALM scheme is capable of localizing more molecules and improving JAC by 50% than ThunderSTORM, thus requiring fewer frames of raw images to achieve high resolution. These results are essentially in agreement with the conclusion in [32], i.e., a high SNR is often indicative of a successful tradeoff between detection rate and accuracy. The noise reduction capability of TV regularization is confirmed again by the reduced RMSE and increased SNR. Although FALCON takes shorter computation time than ALM, the latter exhibits a better performance in other aspects. It is noted that FALCON can spend more execution time in improving the performance by setting the parameter of speed as slow.

Table 2. Comparison of performance with different frames of raw images and methods.

View Table | View all tables in this article

4.2 Real data analysis

To verify the practicability of the algorithm, we use a publicly accessible data set from the single molecule localization challenge (EPFL) website [36]. The Tubulins high-density data set is composed of 500 raw images of $128 \times 128$ pixels (pixel resolution $100 nm$ ). Here an up-sampling factor of 10 is used. We combine every 10 images together to make one new raw image, and each of them contains more overlapped molecules, as illustrated in Fig. 4(c). Consequently, final reconstruction of STORM images is performed using 50 frames of low-resolution images.

Fig. 4 Comparison of FALCON, MUSICAL and ALM on experimental STORM data from EPFL website. (a) Reconstructed image using 500 frames of raw images with ThunderSTORM. (b) A frame of raw STORM image. (c), (d) and (e) are reconstructed images respectively using FALCON, MUSICAL and ALM with 500 frames. (f), (g) and (h) are reconstructed images respectively using FALCON, MUSICAL and ALM with 50 frames. White squares in upper right corners are correspondingly magnified reconstructed results.

Download Full Size | PDF

We perform single molecule localization reconstruction by ThunderSTORM with 500 frames as a comparison standard, as shown in Fig. 4(a). The result is compared to another high-density method in [37], which relies on correlations estimated form a movie stack by implementing the multiple signal classification algorithm (MUSICAL). Figures 4(c)-4(e) present the reconstructed images respectively using FALCON, MUSICAL and ALM with 500 frames. The corresponding results with 50 frames are given in the following row. We observe that the reconstructed image quality using the proposed ALM is improved, especially that the regions of complex microtubule structures are better reconstructed and look sharper, compared with that using FALCON, as indicated by the magnified results of white squares in the upper right corners of Figs. 4(c) and 4(e), also in Figs. 4(f) and 4(h). There are spurious peaks between adjacent microtubules in the reconstructed images using FALCON, and the width of microtubules appears wider. For the comparison to MUSICAL, it can be seen that only a fraction of overlapped molecules can be localized by this method, which results in obvious broken filaments, and the discontinuity is severer using 50 frames of raw images. This is reasonable because that the fewer frames are used the more molecules are contained in each frame. FALCON and MUSICAL suffers from noise and low detection rate, while ALM can achieve better performance in both detection rate and localization accuracy for high molecule density.

Although our algorithm resolves adjacent microtubules better, even for 50 frames with increased molecule density, it also exhibits discontinuity in the reconstructed image. This can be considered as a compromise between the number of frames and the image fidelity. Using fewer frames will shorten the overall acquisition time at the expense of possible broken features, but still can provide fine structure. This observation is similar to that in [28]. Also as suggested in [28], the proposed method has the potential to be firstly used as an alternative to reveal the fine structure of high-density data in short time. Then, other reconstruction methods supplement the broken regions.

The intensity profiles along the blue line and cyan line in Fig. 4(a) are measured and plotted in Fig. 5. It can be seen from Figs. 5(a)-5(f) that the dense microtubules cannot be resolved by FALCON and MUSICAL, while our algorithm generates clearly six microtubules and four microtubules and closer to the result from ThunderSTORM with 500 frames. We visualize that there are missing microtubules and spurious and noisy peaks at the positions between adjacent microtubules in reconstructed result from MUSICAL and FALCON, due to its failure in localizing dense molecules and denoising. However, ALM is capable of detecting Tubulins structures using 50 frames of raw images. The above results also affirm the potential of our method for live-cell imaging, using fewer frames to achieve an acceptable super-resolution image.

Fig. 5 (a) and (b) are plots of intensity profiles of FALCON reconstructed images along the blue and cyan lines, respectively. (c) and (d) are plots of intensity profiles of MUSICAL reconstructed images along the blue and cyan lines, respectively. (e) and (f) are plots of intensity profiles of ALM reconstructed images along the blue and cyan lines, respectively.

Download Full Size | PDF

5. Conclusion

In this paper, an ALM-based algorithm is developed for high-density data reconstruction in super-resolution fluorescence microscopy. We take advantage of the computational advances of ALM and the noise reduction capability of TV norm. The approach is applied to simulated data with noise levels as well as real experimental data. Evaluations against sparse recovery methods represented by CS method and FALCON method and single molecule localization method performed by ThunderSTORM illustrate its advantages such as shorter execution time, higher detection rate and image quality. These results allow the proposed algorithm to be a prime candidate for super-resolution microscopy. By modifying the matrix $A$ to fulfill variable PSFs, it is feasible to extend the algorithm to analyze three dimensional STORM data in the future.

Funding

National Key Research and Development Program of China (2017YFB0403804); National Basic Research Program of China (2015CB352005); National Natural Science Foundation of China (61605120, 11774242, 61525503, 61620106016); Natural Science Foundation of Guangdong Province (2017B020210006, 2014A030312008); Shenzhen Science and Technology R&D and Innovation Foundation (JCYJ20160422151611496, GRCK2017042110420047); Natural Science Foundation of Shenzhen University (827000161).

References

1. A. Lee, K. Tsekouras, C. Calderon, C. Bustamante, and S. Pressé, “Unraveling the thousand word picture: an introduction to super-resolution data analysis,” Chem. Rev. 117(11), 7276–7330 (2017). [CrossRef] [PubMed]

2. S. W. Hell and J. Wichmann, “Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy,” Opt. Lett. 19(11), 780–782 (1994). [CrossRef] [PubMed]

3. 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 313(5793), 1642–1645 (2006). [CrossRef] [PubMed]

4. M. J. Rust, M. Bates, and X. Zhuang, “Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM),” Nat. Methods 3(10), 793–795 (2006). [CrossRef] [PubMed]

5. S. A. Jones, S. H. Shim, J. He, and X. Zhuang, “Fast, three-dimensional super-resolution imaging of live cells,” Nat. Methods 8(6), 499–508 (2011). [CrossRef] [PubMed]

6. J. Huang, M. Sun, J. Ma, and Y. Chi, “Super-resolution image reconstruction for high-density 3D single-molecule microscopy,” IEEE Trans. Comput. Imag. 3(4), 763–773 (2017). [CrossRef]

7. S. J. Holden, S. Uphoff, and A. N. Kapanidis, “DAOSTORM: an algorithm for high- density super-resolution microscopy,” Nat. Methods 8(4), 279–280 (2011). [CrossRef] [PubMed]

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(5), 1377–1393 (2011). [CrossRef] [PubMed]

9. A. Small and S. Stahlheber, “Fluorophore localization algorithms for super-resolution microscopy,” Nat. Methods 11(3), 267–279 (2014). [CrossRef] [PubMed]

10. E. A. Mukamel, H. Babcock, and X. Zhuang, “Statistical deconvolution for superresolution fluorescence microscopy,” Biophys. J. 102(10), 2391–2400 (2012). [CrossRef] [PubMed]

11. L. Zhu, W. Zhang, D. Elnatan, and B. Huang, “Faster STORM using compressed sensing,” Nat. Methods 9(7), 721–723 (2012). [CrossRef] [PubMed]

12. 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(23), 28583–28596 (2013). [CrossRef] [PubMed]

13. S. Cox, E. Rosten, J. Monypenny, T. Jovanovic-Talisman, D. T. Burnette, J. Lippincott-Schwartz, G. E. Jones, and R. Heintzmann, “Bayesian localization microscopy reveals nanoscale podosome dynamics,” Nat. Methods 9(2), 195–200 (2011). [CrossRef] [PubMed]

14. U. Köthe, F. Herrmannsdörfer, I. Kats, and F. A. Hamprecht, “SimpleSTORM: a fast, self-calibrating reconstruction algorithm for localization microscopy,” Histochem. Cell Biol. 141(6), 613–627 (2014). [CrossRef] [PubMed]

15. J. Min, C. Vonesch, H. Kirshner, L. Carlini, N. Olivier, S. Holden, S. Manley, J. C. Ye, and M. Unser, “FALCON: fast and unbiased reconstruction of high-density super-resolution microscopy data,” Sci. Rep. 4(1), 4577 (2014). [CrossRef] [PubMed]

16. M. Ovesný, P. Křížek, Z. Švindrych, and G. M. Hagen, “High density 3D localization microscopy using sparse support recovery,” Opt. Express 22(25), 31263–31276 (2014). [CrossRef] [PubMed]

17. J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed. (Springer, 2006), Chap. 17.

18. M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, “Fast image recovery using variable splitting and constrained optimization,” IEEE Trans. Image Process. 19(9), 2345–2356 (2010). [CrossRef] [PubMed]

19. M. J. Allison, S. Ramani, and J. A. Fessler, “Accelerated regularized estimation of MR coil sensitivities using augmented Lagrangian methods,” IEEE Trans. Med. Imaging 32(3), 556–564 (2013). [CrossRef] [PubMed]

20. W. Hao, R. Feng, Z.-F. Han, and C.-S. Leung, “ADMM-based algorithm for training fault tolerant RBF networks and selecting centers,” IEEE Trans. Neural Netw. Learn. Syst. 29(8), 3870–3878 (2017). [PubMed]

21. J. Li, Y. Shen, and E. Y. Lam, “Hotspot-aware fast source and mask optimization,” Opt. Express 20(19), 21792–21804 (2012). [CrossRef] [PubMed]

22. S. H. Chan, R. Khoshabeh, K. B. Gibson, P. E. Gill, and T. Q. Nguyen, “An augmented Lagrangian method for total variation video restoration,” IEEE Trans. Image Process. 20(11), 3097–3111 (2011). [CrossRef] [PubMed]

23. C. Li, W. Yin, H. Jiang, and Y. Zhang, “An efficient augmented Lagrangian method with applications to total variation minimization,” Comput. Optim. Appl. 56(3), 507–530 (2013). [CrossRef]

24. G.-H. Chen, J. Tang, and S. Leng, “Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets,” Med. Phys. 35(2), 660–663 (2008). [CrossRef] [PubMed]

25. N. Pavillon and N. I. Smith, “Compressed sensing laser scanning microscopy,” Opt. Express 24(26), 30038–30052 (2016). [CrossRef] [PubMed]

26. T. Goldstein and S. Osher, “The split Bregman algorithm for $l_{1}$ ,” SIAM J. Imaging Sci. 2(2), 323–343 (2009). [CrossRef]

27. H. Bingsheng, H. Liu, Z. Wang, and X. Yuan, “A strictly contractive Peaceman-Rachford splitting method for convex programming,” SIAM J. Optim. 24(3), 1011–1040 (2014). [CrossRef] [PubMed]

28. O. Solomon, M. Mutzafi, M. Segev, and Y. C. Eldar, “Sparsity-based super-resolution microscopy from correlation information,” Opt. Express 26(14), 18238–18269 (2018). [CrossRef] [PubMed]

29. J. Li, S. Liu, and E. Y. Lam, “Efficient source and mask optimization with augmented Lagrangian methods in optical lithography,” Opt. Express 21(7), 8076–8090 (2013). [CrossRef] [PubMed]

30. J. L. Morales and J. Nocedal, “Remark on ‘algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound constrained optimization’,” ACM Trans. Math. Softw. 38(1), 1–4 (2011). [CrossRef]

31. A. Beck and M. Teboulle, “Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems,” IEEE Trans. Image Process. 18(11), 2419–2434 (2009). [CrossRef] [PubMed]

32. D. Sage, H. Kirshner, T. Pengo, N. Stuurman, J. Min, S. Manley, and M. Unser, “Quantitative evaluation of software packages for single-molecule localization microscopy,” Nat. Methods 12(8), 717–724 (2015). [CrossRef] [PubMed]

33. T. Cheng, D. Chen, B. Yu, and H. Niu, “Reconstruction of super-resolution STORM images using compressed sensing based on low-resolution raw images and interpolation,” Biomed. Opt. Express 8(5), 2445–2457 (2017). [CrossRef] [PubMed]

34. M. Ovesný, P. Křížek, J. Borkovec, Z. Švindrych, and G. M. Hagen, “ThunderSTORM: a comprehensive ImageJ plug-in for PALM and STORM data analysis and super-resolution imaging,” Bioinformatics 30(16), 2389–2390 (2014). [CrossRef] [PubMed]

35. R. P. J. Nieuwenhuizen, K. A. Lidke, M. Bates, D. L. Puig, D. Grünwald, S. Stallinga, and B. Rieger, “Measuring image resolution in optical nanoscopy,” Nat. Methods 10(6), 557–562 (2013). [CrossRef] [PubMed]

36. Biomedical Imaging Group, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, “Benchmarking of Single-Molecule Localization Microscopy Software,” http://bigwww.epfl.ch/smlm/.

37. K. Agarwal and R. Macháň, “Multiple signal classification algorithm for super-resolution fluorescence microscopy,” Nat. Commun. 7, 13752 (2016). [CrossRef] [PubMed]

	Method	JAC %	RMSE (nm)	SNR (dB)	FRC (nm)	Time (s)
120 frames	ALM	0.94	13.56	5.1921	32.27	660
	ThunderSTORM	0.65	28.12	3.1603	50.80	1208
	FALCON	0.92	23.35	3.3575	39.13	120
30 frames	ALM	0.91	22.47	3.0612	40.11	120
	ThunderSTORM	0.49	38.77	2.0168	55.43	234
	FALCON	0.86	27.44	2.2580	50.72	30

	Method	JAC %	RMSE (nm)	SNR (dB)	FRC (nm)	Time (s)
120 frames	ALM	0.94	13.56	5.1921	32.27	660
	ThunderSTORM	0.65	28.12	3.1603	50.80	1208
	FALCON	0.92	23.35	3.3575	39.13	120
30 frames	ALM	0.91	22.47	3.0612	40.11	120
	ThunderSTORM	0.49	38.77	2.0168	55.43	234
	FALCON	0.86	27.44	2.2580	50.72	30

Efficient image reconstruction of high-density molecules with augmented Lagrangian method in super-resolution microscopy

Abstract

1. Introduction

2. Problem formulation

3. Augmented Lagrangian method for super-resolution microscopy

3.1 The cost function

3.2 Optimization flow

4. Results

4.1 Simulated data analysis

4.2 Real data analysis

5. Conclusion

Funding

References

Cited By

Figures (5)

Tables (2)

Equations (17)

Optics Express