The fast imaging speed and low-intensity requirement of structured illumination microscopy (SIM) have made it one of the most widely used imaging tools in live cell imaging. In order to obtain a high fidelity reconstructed image, a precise estimation of the phase of the illumination pattern is required, especially in those structured illumination based techniques that rely on high-order harmonics to improve the resolution. This can be achieved in one of two fundamental ways. The first is to build a high-end control system capable of shifting a sinusoidal pattern with high precision, while the second is to apply estimation algorithms to determine how patterns shift during post-processing. The latter method is preferred in low-cost super-resolution imaging systems; however, existing algorithms are either time-consuming or fail due to noise and a low modulation depth. In this paper, we introduce additional matrixes into the phase estimation algorithm and propose an inverse matrix based phase estimation method with which analytical solutions of the phases can be determined without iteration. The proposed algorithm was validated via simulation and experiments using a home-made total internal reflection fluorescent SIM system (TIRF-SIM). When tested, the method obtained the true phase even when the modulation depth was low. The source code is now available for download by researchers and others.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
When compared to electron microscopy, optical microscopy is less damaging to live cells; however, conventional optical microscopy is unable to resolve details smaller than 200 nm. Recent developments in fluorescence microscopy have enabled nanometer resolution, and super-resolution microscopy has great potential in imaging live cells. Due to the chemical specificity of proteins, it can also be applied to fluorescence labeling. A variety of fluorescence super-resolution techniques can be employed to uncover interactions between different organelles, and other techniques, such as stimulated emission depletion microscopy (STED) [1,2], stochastic optical reconstruction microscopy (STORM) , photo-activated localization microscopy (PALM) , and structured illumination microscopy (SIM) [5–7] can support higher resolutions than conventional optical microscopy. Among these, SIM is fast and can be adapted to different kinds of fluorophores. Despite the fact that the resolution enhancement in SIM is limited to a factor of two, in addition to its speed and adaptability, this technique is popular in biomedical research due to its lower phototoxicity and photobleaching effects compared to those of other methods.
Various methods have been proposed to improve the resolution in SIM, such as total internal reflection fluorescent microscopy (TIRF) [8–10] based structured illumination, plasmonic SIM (PSIM) [11,12], 3D-SIM [13–15], and saturated SIM (SSIM) . While SSIM requires a greater intensity than SIM and may cause more photodamage to the cell, photo-switchable proteins can be used to reduce the illumination intensity requirement [17,18].
A key step in the post-processing of images is in determining the modulation frequency and estimating the phase. A high fidelity reconstruction typically requires the correct modulation frequency and phase of each illumination pattern, and errors in the estimated parameters degrade the image quality and introduce artifacts. The phase estimation algorithms are employed to estimate the phase of the fundamental, while the phase of the n-th harmonic in SSIM and non-linear SIM can be expressed as . Consequently, the error in the estimation step is also magnified by a factor of n if the n-th harmonic is applied. For this reason, SSIM and non-linear SIM is extremely sensitive to estimation errors.
Hence, many methods have been proposed to obtain the correct phase of each sinusoidal pattern [19–23]. For example, iterative cross-correlation based algorithms provide a precise phase estimation in most cases; however, they are time-consuming and may get trapped in a local minimum , which is the case for most optimization algorithms and limits their application. The auto-correlation based phase estimation algorithm can be used to directly calculate the phase of the illumination pattern without optimization, but its performance is generally poor compared to that of iterative methods when the SNR is poor, and assumes the modulation vector differs from the spatial frequency vector for some periodic samples .
In this paper, we propose a new algorithm based on an inverse matrix that improves on existing algorithms. The proposed algorithm solves equation sets describing the phases of the illumination pattern analytically without requiring optimization. In this way, the algorithm avoids being trapped in a local minimum, which is a risk in iteration based algorithms. Compared with the auto-correlation method, this approach performs better in low SNR and low modulation depth conditions. The phase-only correlation method is introduced in this paper to enhance the amplitude of the peak resulting from the modulation of the illumination field. The source code for the inverse matrix based method is available for download (Code 1 ).
2. Principles and algorithms
Assuming a three-step phase shifting strategy is applied for simplification, the image sequences (i = 1, 2, 3) captured for a certain illumination pattern orientation can be expressed as:
In the Fourier domain, this can be represented as:
In SIM, the phase optimization step in the reconstruction process comes after estimating the modulation frequency. Therefore, the prerequisite for phase estimation is to obtain the correct modulation frequency vector . The phase-of-peak method proposed in the past estimates the modulation frequency vector and phases of the illumination pattern using the Fourier transform of the acquired images. This method assumes that the modulation frequency is smaller than the cutoff frequency of the OTF. As the zero-frequency of the sample is shifted to , a peak appears at the corresponding position in the Fourier transform and the phase-of-peak method regards the phase of this peak as the phase of the illumination pattern. However, this assumption is not always satisfied. To address this, the cross-correlation method was proposed to correct the modulation matrix. Unlike the phase-of-peak method, which estimates the phase based on a single point, the cross-correlation method calculates the correlation between different components. All of the detectable frequency components in the Fourier domain contribute to the “peak” located at as well as the estimated phase. The correlation based method can be applied when is approximately equal to or larger than kc, as in TIRF-SIM [20,21,25]. In addition to the above, other approaches that enhance the contrast of the peak based on the cross-correlation include the deconvolution method .
2.1 Estimation of the modulation frequency vector
In the estimation process, a low-pass filter similar to the Wiener filter is first applied to avoid the high-frequency noise beyond the OTF.
Despite the fact that the correct phase of the illumination pattern is not known, the high-frequency and baseband components can be separated with the inverse matrix of . Here, can be selected casually in principle as long as . The reason the inverse matrix can separate different frequency components is illustrated by the following equation, which was derived by calculating the product on the left side:
At the conclusion of the above steps, the noise has been reduced and the different frequency components have been roughly separated. Next, the modulation frequency vector is estimated via the phase-only correlation, which yields:
The above approximation holds as , , and are uncorrelated.
A set of data was acquired using an homebuilt system , and the results after applying the existing and phase-only algorithms are shown in Fig. 1. The corresponding reconstruction results are shown later in this paper. The peak intensity, which was located at , is shown in Fig. 1(e). Even though a high-pass mask was applied to the auto-correlation results, it was still difficult to determine the modulation frequency as the amplitude at served as a local maximum and there was a large amount of background noise. When the deconvolution method was used to enhance the signal located at , a high-pass mask was still required in order to detect the correct modulation frequency. In contrast, when the phase-only algorithm was used, the intensity became a global maximum, which obviated the need for a mask, as implied in Fig. 1(e). The signal amplitude located at was significantly enhanced, which demonstrates the potential of the phase-only estimation algorithm compared to existing methods. The phase-only correlation algorithm may increase the sensitivity and improve the performance of the cross-correlation method, which can be used to optimize the phase based on the conjugated peak .
2.2 Estimation of the phases
Once the modulation frequency has been determined, the phase of can also be determined. As this phase is a function of the phase of the illumination pattern , it is possible to compute the phases analytically. With this in mind, an inverse matrix based algorithm is proposed, the principle of which is as follows.
To mitigate the influence of the noise, the phase was calculated as follows:Eq. (9), and the second approximation is an assumption that appears in all correlation based algorithms [20,21]. This equation also holds when is replaced by , as will be dominant in. And , are also uncorrelated.
To determine how acts on , the product of the matrix (the inverse is calculated by the adjugate matrix method) is considered:Eq. (10):
As the ratios of the imaginary and real parts in the right and left sides of Eq. (11) are equal, we have:
The simplified expression is:
Note that originate from the inverse matrix that depends on the chosen phases , and can be determined in the cross-correlation. Thus, the remaining unknown variables represent the phases of interest. Also, note that the sign of the imaginary and real parts of Eq. (11) should be the same, as this is useful in the inverse matrix based phase estimation method explained below. As are the only unknowns, theoretically, the above equations can be solved via three independent equations. While there will be multiple solutions in Eq. (13) for , unique solutions that match the correct phases can be found using Eq. (11). As the constraints associated with Eq. (11) are more robust than those for Eq. (13), such as the sign issue mentioned earlier, solutions in this paper that satisfy Eq. (13) but do not hold for Eq. (11) are considered redundant.
However, as it difficult to solve trigonometric equations involving three variables, Eq. (13) is solved by introducing additional independent equations to avoid applying the Pythagorean trigonometric identity of sine and cosine. For example, consider the following set of equations: . If the third equation is not known, solving the first two equations will result in the additional answer , despite the fact that the other answers were also obtained. In contrast, if all three equations are known, it is straightforward to solve them as the Pythagorean identity need only be applied once.
The classical inverse (demodulation) matrix, which is , has only three independent arguments because m can be considered as a weighting factor of the modulation matrix, and therefore does not influence the phase of Eq. (11). Hence, the inverse matrix has three “degrees of freedom.” In this way, three independent equations can be obtained that rely on only .
To obtain four or more independent equations, the inverse matrix can be changed to , where is positive and real. The product of this inverse matrix and the modulation matrix is:
To obtain a better approximation between and , can be updated before the phase-only correlation and estimation of by eliminating the widefield component:
Once has been updated, Eqs. (8) and (13) can be applied to solve Eq. (13). Since there are up to six independent arguments in the inverse matrix , it is feasible to obtain six independent equations in the form given by Eq. (13). In addition, the equation set can be considered to be a linear equation set that is consistent with the Pythagorean theorem. In this case, no redundant answer is obtained as the constraints are sufficiently robust. However, it may be computationally expensive and time-consuming to manipulate six equations because this requires calculating additional correlations in order to obtain the phase . Here, four independent equations are used to convert the original problem into two trigonometric equations with only two unknown angles (phases). To eliminate the redundant answers from the answer set, a cost function is introduced based on Eq. (11):
Due to the noise in the captured image, it is likely that the sine and cosine values of the phase (s = 1/2/3), which was replaced by the above substitution, does not obey the Pythagorean theorem. For this reason, the phase was estimated via the following equation:
3. Simulations and experiments
In order to verify the above algorithm, a simulation was conducted in which different numbers of photons were emitted from a particular object, and the resulting phase was estimated using the proposed inverse matrix based algorithm, the cross-correlation and the auto-correlation algorithm. The object in question was a Siemens star (see the ground truth discussion below) and the spatial frequency of the illumination pattern was , which is larger than the detection cutoff frequency and feasible due to the Stokes shift of the fluorescence light or the surface plasmonic wave. The designated phase shift was and the uniformly distributed phase error was in the range [-0.3, 0.3].
The inverse matrixes employed in the simulations and experiments in this paper were , , , and as these parameters were the most promising among the six different sets evaluated.
The average phase errors of the inverse matrix based, cross-correlation, and auto-correlation algorithms for different modulation depths are shown in Fig. 2. The program used to implement these algorithms was written in MATLAB and executed on a personal computer with an Intel Core i5-4200H 2.8 GHz processor. The run-times of the inverse matrix based, cross-correlation, and auto-correlation methods were about 0.27, 0.33, and 0.18 s, respectively, when estimating three phases for an image size of 513 × 513 pixels. In an ideal system, the performance the three algorithms is quite similar, as shown in Fig. 2(a), while the auto-correlation method consumes less computation time. Thus, in this case, the auto-correlation method is a good choice for reconstruction.
However, when the modulation depth is lower than expected, the inverse matrix and iterative cross-correlation based algorithms are better than the auto-correlation method. In this case, the execution time of the inverse matrix based phase estimation method was 20% shorter than that of the iterative method.
If the non-linear effect is utilized, 2N + 1 images are required to cover the N-th harmonic [16–18]. In this case, the images can be divided into small groups of three images each (there may be duplicate images in the different subgroups). Once the true phases of the illumination patterns have been estimated, the modulation depth estimation step can be improved because it is also based on the correlation result in the correlation-based estimation algorithms [21,25].
The parameters that are required for SIM reconstruction are obtained in the above step. The image process was simulated using a standard object that was the same as that in Fig. 2, in which the maximum phase error was 0.6 rad. As shown in Fig. 3(f), there was no redundant component as the phases of the illumination pattern were the correct ones as per the proposed method. The modulation depth in Figs. 3(a)–(f) was 0.6. Apart from the estimated phases, the other aspects of the reconstruction method remained the same and the reconstruction result shown in Fig. 3(c) is similar to the ground truth image illustrated in Fig. 3(a). However, the stripes in Fig. 3(b) appear broken as a result of the existing redundant component.
To demonstrate the reliability of the inverse matrix based algorithm, the image process was simulated with a low modulation depth, after which the image was reconstructed using the above three methods. The results are shown in Figs. 3(g)–(i). When the modulation depth was m = 0.3, one set of phases (the second pattern orientation) estimated using the auto-correlation method was completely incorrect. Thus, the corresponding reconstruction result, which is shown in Fig. 3(g), is worse than the results of the other two phase estimation methods.
A home-made TIRF-SIM system  was constructed to verify the inverse matrix based phase estimation algorithm. In the experiment, 100 nm diameter yellow-green fluorescent beads (excitation wavelength: 505 nm/emission wavelength: 515 nm, Model F8803, Life Technologies) were imaged and the results are shown in Fig. 4. The images were captured using a high numerical aperture TIRF objective (100 × / 1.49 TIRF, Nikon) and the results were reconstructed using the inverse matrix based phase estimation algorithm. All aspects of the code used for the simulation remained unchanged, except for the phase estimation method. The data used to obtain Fig. 1 is from the same data set as that in Fig. 4. As shown in the figure, the redundant component that appeared in Fig. 1 disappeared in this figure. This indicates that the phases estimated by the inverse matrix algorithm were the correct phases. It can be seen in the image and subfigures in Fig. 4(c) that the beads are less distorted than those in Fig. 4(b). The contrast in Fig. 4(c) is also better. The brightness of each bead in the phase corrected results are similar, which is highly consistent with the prior knowledge, while the brightness fluctuates in the results without phase correction.
Microtubule was also imaged and reconstructed via the three phase estimation algorithm. The sample material was Alexa Fluor 647 (excitation wavelength: 650 nm / emission wavelength: 665 nm). The sample preparation and staining procedures were as follows.
COS-7 (African green monkey kidney fibroblast-like cell line) cells were purchased from Boster Biological Technology Co., Ltd., Wuhan, China, and cultured in Dulbecco's Modified Eagle's Medium (DMEM; Thermo Fisher Scientific, Inc.). All media were supplemented with 10% (v/v) fetal bovine serum (FBS; Thermo Fisher Scientific, Inc.), and were kept at 37°C.
COS-7 cells were seeded in Nunc Glass Bottom Dishes (Φ 12 mm, Thermo Fisher Scientific, Inc.) at a density of 1.5-2.0 × 104 per well in growth medium (150 µL).
Antibodies were labeled using the manufacturers’ protocols. In brief, 10 µl of 1 M aqueous NaHCO3 (pH 8.0) was added to 50–100 µg of antibody to a final antibody concentration of ~2 µM. Alexa Fluor 647 N-hydroxysuccinimidyl ester (Invitrogen) were dissolved in dimethyl sulfoxide and added to antibody solutions at a final concentration of ~4 µM to achieve the desired dye to protein ratios. Labeled antibodies were purified by gel filtration using NAP-5 columns (GE Healthcare).
The immunostaining procedure consisted of fixation for 10 min with 3% paraformaldehyde and 0.1% glutaraldehyde in PBS, washing with PBS, reduction for 5 min with 0.1% sodium borohydride in PBS to reduce background fluorescence, washing with PBS, blocking for 30 min with 3% bovine serum albumin and 0.25% (v/v) Triton X-100 in PBS (blocking buffer), staining for 45 min with rat primary antibody to tubulin (Abcam) diluted in blocking buffer to a concentration of 10 µg ml−1, washing with PBS, incubation for 45 min with secondary antibodies (~1–2 dyes per antibody; Jackson ImmunoResearch Laboratories) at a concentration of ~2.5 µg ml−1 in blocking buffer, washing with PBS, fixation after labeling for 10 min with 3% paraformaldehyde and 0.1% glutaraldehyde in PBS, and finally washing with PBS.
As shown in Fig. 5(e), reconstructed results using the cross-correlation and inverse matrix based phase estimation methods were better than those using the auto-correlation method due to the accuracy of the estimated phases.
In this paper, a phase-only correlation method was introduced to optimize the modulation spatial frequency vector . As the estimation of the modulation frequency is the first step in the reconstruction process, the correct modulation frequency vector is the basis of all subsequent steps. This phase-only correlation method outperforms previous methods and significantly enhances the peak signal located at .
By analyzing the inverse matrix and the outcomes of the separation steps, an inverse matrix based phase estimation algorithm was developed. Unlike the cross-correlation phase optimization algorithm, which uses iteration when determining the correct phase, the inverse matrix based algorithm does not require iteration and cannot be trapped in a local minimum. Thus, by way of comparison, the inverse matrix algorithm was demonstrated to be a strong competitor to the auto-correlation algorithm. While the performance of the inverse matrix method is similar to the cross-correlation method, the proposed method is faster, and the two methods are both better than the auto-correlation method when the modulation depth is low.
Although the inverse matrix based phase estimation method has the potential to be widely used, more work is required. While a correct set of phases is essential for high quality reconstruction, it should be noted that the quality of the reconstruction is not solely dependent on the phases as the SNR of the image is also important. Thus, the use of an improved noise reduction algorithm would improve the reliability of the reconstruction results and reduce the phase estimation errors. In addition, as the parameters of the inverse matrix affect the errors of the phase estimation algorithm, the selection of an appropriate set of parameters has the potential to reduce the average phase error in the inverse matrix based method. These topics can be explored in future research.
National Basic Research Program of China (973Program) (2015CB352003); National Natural Science Foundation of China (NSFC) (61735017, 61335003, 61427818, 61827825); Natural Science Foundation of Zhejiang province (LR16F050001); and the Fundamental Research Funds for the Central Universities (2017FZA5004).
The authors declare that there are no conflicts of interest related to this article.
1. 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]
2. T. A. Klar, S. Jakobs, M. Dyba, A. Egner, and S. W. Hell, “Fluorescence microscopy with diffraction resolution barrier broken by stimulated emission,” Proc. Natl. Acad. Sci. U.S.A. 97(15), 8206–8210 (2000). [CrossRef] [PubMed]
4. 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]
6. P. Kner, B. B. Chhun, E. R. Griffis, L. Winoto, and M. G. L. Gustafsson, “Super-resolution video microscopy of live cells by structured illumination,” Nat. Methods 6(5), 339–342 (2009). [CrossRef] [PubMed]
7. R. Heintzmann and C. G. Cremer, “Laterally modulated excitation microscopy: improvement of resolution by using a diffraction grating,” Proc. SPIE 3568, 185–196 (1999). [CrossRef]
9. E. Chung, D. Kim, Y. Cui, Y.-H. Kim, and P. T. C. So, “Two-Dimensional Standing Wave Total Internal Reflection Fluorescence Microscopy: Superresolution Imaging of Single Molecular and Biological Specimens,” Biophys. J. 93(5), 1747–1757 (2007). [CrossRef] [PubMed]
10. R. Fiolka, M. Beck, and A. Stemmer, “Structured illumination in total internal reflection fluorescence microscopy using a spatial light modulator,” Opt. Lett. 33(14), 1629–1631 (2008). [CrossRef] [PubMed]
12. F. Wei, D. Lu, H. Shen, W. Wan, J. L. Ponsetto, E. Huang, and Z. Liu, “Wide Field Super-Resolution Surface Imaging through Plasmonic Structured Illumination Microscopy,” Nano Lett. 14(8), 4634–4639 (2014). [CrossRef] [PubMed]
13. M. G. L. Gustafsson, L. Shao, P. M. Carlton, C. J. R. Wang, I. N. Golubovskaya, W. Z. Cande, D. A. Agard, and J. W. Sedat, “Three-Dimensional Resolution Doubling in Wide-Field Fluorescence Microscopy by Structured Illumination,” Biophys. J. 94(12), 4957–4970 (2008). [CrossRef] [PubMed]
14. L. Shao, P. Kner, E. H. Rego, and M. G. L. Gustafsson, “Super-resolution 3D microscopy of live whole cells using structured illumination,” Nat. Methods 8(12), 1044–1046 (2011). [CrossRef] [PubMed]
15. L. Schermelleh, P. M. Carlton, S. Haase, L. Shao, L. Winoto, P. Kner, B. Burke, M. C. Cardoso, D. A. Agard, M. G. L. Gustafsson, H. Leonhardt, and J. W. Sedat, “Subdiffraction Multicolor Imaging of the Nuclear Periphery with 3D Structured Illumination Microscopy,” Science 320(5881), 1332–1336 (2008). [CrossRef] [PubMed]
16. M. G. Gustafsson, “Nonlinear structured-illumination microscopy: wide-field fluorescence imaging with theoretically unlimited resolution,” Proc. Natl. Acad. Sci. U.S.A. 102(37), 13081–13086 (2005). [CrossRef] [PubMed]
17. D. Li, L. Shao, B.-C. Chen, X. Zhang, M. Zhang, B. Moses, D. E. Milkie, J. R. Beach, J. A. Hammer 3rd, M. Pasham, T. Kirchhausen, M. A. Baird, M. W. Davidson, P. Xu, and E. Betzig, “ADVANCED IMAGING. Extended-resolution structured illumination imaging of endocytic and cytoskeletal dynamics,” Science 349(6251), aab3500 (2015). [CrossRef] [PubMed]
18. E. H. Rego, L. Shao, J. J. Macklin, L. Winoto, G. A. Johansson, N. Kamps-Hughes, M. W. Davidson, and M. G. Gustafsson, “Nonlinear structured-illumination microscopy with a photoswitchable protein reveals cellular structures at 50-nm resolution,” Proc. Natl. Acad. Sci. U.S.A. 109(3), E135–E143 (2012). [CrossRef] [PubMed]
19. S. A. Shroff, J. R. Fienup, and D. R. Williams, “Phase-shift estimation in sinusoidally illuminated images for lateral superresolution,” J. Opt. Soc. Am. A 26(2), 413–424 (2009). [CrossRef] [PubMed]
20. K. Wicker, “Non-iterative determination of pattern phase in structured illumination microscopy using auto-correlations in Fourier space,” Opt. Express 21(21), 24692–24701 (2013). [CrossRef] [PubMed]
23. A. Lal, C. Shan, and P. Xi, “Structured Illumination Microscopy Image Reconstruction Algorithm,” IEEE J. Sel. Top. Quantum Electron. 22(4), 50–63 (2016). [CrossRef]
24. R. Cao, “Inverse matrix phase algorithm: SIM reconstruction software using inverse matrix based phase estimation method”, Zenodo (2018) [retrieved 17 September 2018], https://doi.org/10.5281/zenodo.1314351.
25. M. Müller, V. Mönkemöller, S. Hennig, W. Hübner, and T. Huser, “Open-source image reconstruction of super-resolution structured illumination microscopy data in ImageJ,” Nat. Commun. 7, 10980 (2016). [CrossRef] [PubMed]
26. Y. Chen, R. Cao, W. Liu, D. Zhu, Z. Zhang, C. Kuang, and X. Liu, “Widefield and total internal reflection fluorescent structured illumination microscopy with scanning galvo mirrors,” J. Biomed. Opt. 23(4), 1–9 (2018). [CrossRef] [PubMed]