1. Introduction

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) [3], photo-activated localization microscopy (PALM) [4], 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) [16]. 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 $k_0$ and phase $\phi$ 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 $\phi$ of the fundamental, while the phase of the n-th harmonic in SSIM and non-linear SIM can be expressed as ${\phi }_{n-th}=n\phi$. 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 [21], 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 [20].

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 [24]).

2. Principles and algorithms

Assuming a three-step phase shifting strategy is applied for simplification, the image sequences $D_i$(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 $k_0$. 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 $\left|k_0\right|$ is smaller than the cutoff frequency $k_c$of the OTF. As the zero-frequency of the sample is shifted to $\pm k_0$, 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 $k_0$ as well as the estimated phase. The correlation based method can be applied when $\left|k_0\right|$ is approximately equal to or larger than k_c, 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 [22].

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 $M\left(1,1,1;{\gamma }_1,{\gamma }_2,{\gamma }_3\right)$. Here, ${\gamma }_i$can be selected casually in principle as long as ${\gamma }_i\ne {\gamma }_j\left(i\ne j\right)$. 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 $\tilde{S}\left(k\right)$, $\tilde{S}\left(k\pm k_0\right)$, and $\tilde{S}\left(k+2k_0\right)$ are uncorrelated.

A set of data was acquired using an homebuilt system [26], 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 $\pm k_0$, 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 $k_0$ 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 $k_0$, 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 $-k_0$ 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 [21].

 

Fig. 1 Comparison of different modulation frequency estimation methods. The peak inside the green boxes in (a)–(d) is located at $k_0$, while the peak inside the blue box is located at $-k_0$. The intensity of (a)–(d) was normalized by the corresponding maximum in each image. The correlation results were obtained when the phase shift was set as $2\pi /3$. Because the actual phase shift differed from the desired $2\pi /3$, the redundant component appeared in the cross-correlation and phase-only correlation results. A high-pass mask was applied in (a) to exclude components whose spatial frequency was smaller than $0.8k_c$. (e) The peak intensities are located at $k_0$(green square) and $-k_0$(blue circle). Abbreviations: AC, auto-correlation; CC, cross-correlation; CCwD, cross-correlation with deconvolution; PoC, phase-only correlation.

Download Full Size | PDF

2.2 Estimation of the phases

Once the modulation frequency has been determined, the phase of $C\left(k^{\prime }\right)$ can also be determined. As this phase is a function of the phase of the illumination pattern ${\psi }_{peak}=f\left({\phi }_1^{ill},{\phi }_2^{ill},{\phi }_3^{ill}\right)$, 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 ${\psi }_{peak}$ was calculated as follows:

To determine how ${\phi }_i^{ill}$ acts on ${\varphi }_1$, the product of the matrix (the inverse is calculated by the adjugate matrix method) is considered:

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 $b_i,{\beta }_i(i=1,2,3)$ originate from the inverse matrix that depends on the chosen phases ${\gamma }_i$, and ${\psi }_{peak}$ 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 ${\phi }_i^{ill}(i=1,2,3)$ 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 ${\phi }_i^{ill}\in \left[0,2\pi \right)$, 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: $2\sin A+3\sin B=\sqrt{3},2\cos A+\cos B=2,\sin B=0$. If the third equation is not known, solving the first two equations will result in the additional answer $\left(A,B\right)=\left(\frac{\pi }{3},0\right)$, 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 $M^{-1}\left(1,1,1;{\gamma }_1,{\gamma }_2,{\gamma }_3\right)$, 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 $M^{-1}\left(1,1,1;{\gamma }_1,{\gamma }_2,{\gamma }_3\right)$ has three “degrees of freedom.” In this way, three independent equations can be obtained that rely on only ${\gamma }_i$.

To obtain four or more independent equations, the inverse matrix can be changed to $M^{-1}\left({\sigma }_1,{\sigma }_2,{\sigma }_3;{\gamma }_1,{\gamma }_2,{\gamma }_3\right)$, where${\sigma }_i$ is positive and real. The product of this inverse matrix and the modulation matrix $M\left(1,1,1;{\phi }_1^{ill},{\phi }_2^{ill},{\phi }_3^{ill}\right)$ is:

To obtain a better approximation between ${\psi }_{peak}$ and ${\varphi }_1$, ${\tilde{R}}_2\left(k\right)$ can be updated before the phase-only correlation and estimation of ${\psi }_{peak}$ by eliminating the widefield component:

Once ${\tilde{R}}_2\left(k\right)$ 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 $M^{-1}\left({\sigma }_1,{\sigma }_2,{\sigma }_3;{\gamma }_1,{\gamma }_2,{\gamma }_3\right)$, 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 ${\psi }_{peak}$. 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 ${\phi }_s^{ill}$(s = 1/2/3), which was replaced by the above substitution, does not obey the Pythagorean theorem. For this reason, the phase ${\phi }_s^{ill}$ 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 $1.1k_c$, 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 $2\pi /3$ 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 $M_1^{-1}\left(1.1,1.1,0.9;0,\frac{2}{3}\pi ,\frac{4}{3}\pi \right)$, $M_2^{-1}\left(1.2,0.9,0.8;0,\frac{1}{3}\pi ,\frac{2}{3}\pi \right)$, $M_3^{-1}\left(1.15,0.85,1.3;0,\frac{4}{7}\pi ,\frac{10}{7}\pi \right)$, and $M_4^{-1}\left(0.8,1.2,1.3;0,\frac{4}{5}\pi ,\frac{8}{5}\pi \right)$ 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.

 

Fig. 2 Simulations of the proposed inverse matrix based phase estimation algorithm versus the auto-correlation and cross-correlation algorithms. The period of the illumination pattern was 195 nm and the minimum detectable period in this simulation was 214 nm. As the radial pattern consisted of certain frequencies, the auto-correlation algorithm failed to retrieve the correct phase under low modulation depth conditions even though the number of photons was large. The value of m in the lower left of each subfigure demotes the modulation depth. In this simulation, Poisson noise was introduced, the minimum detectable number of photons was 10, and 24 phases were simulated when calculating the mean phase error.

Download Full Size | PDF

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.

 

Fig. 3 (a) Simulation of the ground truth, (b) uncorrected reconstruction, and (c) phase corrected result using the inverse matrix based algorithm. The Fourier transform of the reconstructed images (b) and (c) are shown in (d). The green background pattern in (d)–(f) is the Fourier transform of the ground truth sample. (e) and (f) show the corresponding areas in (d), and the white dashed square in (e) highlights part of the redundant Fourier component that appears due to the error in the phase. The reconstruction results using (g) the auto-correlation, (h) cross-correlation, and (i) inverse matrix based method employed the same set of data that was also applied in Fig. 2 (m = 0.3, number of photons = 3162, three images were used).

Download Full Size | PDF

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 [26] 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.

 

Fig. 4 Experimental results: (a) the filtered widefield results, (b) the reconstructed results without correction, and (c) the results with phase correction via the inverse matrix based method (note the brightness was increased by 40% to improve the visibility). (d) The phase-only correlation results of the reconstructed image. Lower left: without correction, and upper right: with correction. The intensity profile along the cutline indicated by the arrows in (a)–(c) is shown in (e) using the results without brightness enhancement.

Download Full Size | PDF

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.

 

Fig. 5 Filtered widefield result (a), SIM results with the auto-correlation (b), cross-correlation (c), and inverse matrix (d) based phase estimation methods. The profile along the line segment in each subfigure is shown in (e). The profile was normalized by the corresponding maximum along each line segment. COS-7 cells were used in this experiment. Abbreviation: WF, widefield; AC, auto-correlation method; CC, cross-correlation method; IM, inverse matrix based method.

Download Full Size | PDF

4. Conclusions

In this paper, a phase-only correlation method was introduced to optimize the modulation spatial frequency vector $k_0$. As the estimation of the modulation frequency is the first step in the reconstruction process, the correct modulation frequency vector $k_0$ is the basis of all subsequent steps. This phase-only correlation method outperforms previous methods and significantly enhances the peak signal located at $k_0$.

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 $M^{-1}\left({\sigma }_1,{\sigma }_2,{\sigma }_3;{\gamma }_1,{\gamma }_2,{\gamma }_3\right)$ 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.