Inverse matrix based phase estimation algorithm for structured illumination microscopy

: 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.

technique is popular in biomedical research due to its lower phototoxicity and photobleaching effects compared to those of other methods.
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 0 k 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 n th n ϕ ϕ − = . 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][20][21][22][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 [24] (Code 1).

Principles and algorithms
Assuming a three-step phase shifting strategy is applied for simplification, the image sequences i D (i = 1, 2, 3) captured for a certain illumination pattern orientation can be expressed as: where * refers to convolution, S denotes the sample, m denotes the modulation depth, 0 k is the modulation light wave vector, i ϕ is the phase of the i-th illumination pattern, and de where the tilde (~) over S and D denotes the corresponding Fourier transform, i μ is real, ill i ϕ denotes the actual phase of the illumination pattern, de H  denotes the detection optical transfer function (OTF), and M is the modulation matrix.
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 0 k . 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 0 k is smaller than the cutoff frequency c k of the OTF. As the zero-frequency of the sample is shifted to 0 ±k , 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-ofpeak 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 0 k as well as the estimated phase. The correlation based method can be applied when 0 k 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].

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.
where conj is the complex conjugation operation and 1 ε is a positive regularization constant related to the noise that also serves to eliminate the zero denominator. Despite the fact that the correct phase of the illumination pattern is not known, the highfrequency and baseband components can be separated with the inverse matrix of ( ) where c 0 is a complex number and 1 2 1 2 , , , r r φ φ are non-negative real numbers. In the simulations and experiments in this study, 1,1,1; 0, , M π π was selected as it represents an ideal set of phase shifts. Also, for purposes of simplification, r 1 >r 2 was assumed. The separated results are: where i R  represents the separated result. Thus, 2 R  can be expressed as: where the first term on the right side of the equation is dominant as r 1 >r 2 and the second term is a redundant component caused by the incorrect phase. Hence, 2 R  can be considered an estimate of a particular high-frequency component.
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: where ⊗ denotes correlation, W  is the Fourier transform of the widefield, and 2 ε is a positive constant to avoid a zero denominator. The Fourier theorem can be applied to calculate the cross-correlation and estimate 0 k to sub-pixel accuracy. In the modulation frequency estimation, W  can be replaced by R  . This is because this estimation only requires the peak located at 0 k (or 0 −k ). The corresponding correlation result yields: Vol. 9, No. 10 | 1 Oct 2018 | BIOMEDICAL OPTICS EXPRESS 5040 The above approximation holds as ( ) 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 0 ±k , 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 0 k 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 0 k , 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 0 −k 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].

Estimation of the phases
Once the modulation frequency has been determined, the phase of ( ) C ′ k 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 peak ψ was calculated as follows: where ( ) arg ⋅ is the argument of a complex number. The first approximation holds because are considered to be uncorrelated; thus, the signal generated by ( ) is much smaller (typically 2-3 orders of magnitude) than the other terms in Eq. (9), and the second approximation is an assumption that appears in all correlation based algorithms [20,21]. This equation also holds when W  is replaced by are also uncorrelated.
To determine how ill i ϕ acts on 1 φ , the product of the matrix (the inverse is calculated by the adjugate matrix method) is considered: where 1 2 , x x are real numbers and , ( 1,2,3) 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: where . originate from the inverse matrix that depends on the chosen phases i γ , and 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 ( 1,2,3) ill i i ϕ = 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)  To obtain four or more independent equations, the inverse matrix can be changed to ( ) where C w is complex and is determined by the inverse matrix, and 3 x and 0 w are real positive numbers. Now,  To obtain a better approximation between peak ψ and 1 φ , ( ) 2 R k  can be updated before the phase-only correlation and estimation of peak ψ by eliminating the widefield component: 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 ( ) 1  1  2  3  1  2  3 , , ; , , 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 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): where the subscript t denotes the t-th inverse matrix where the subscript "est" denotes the estimated values of sin ill s ϕ and cos ill s ϕ , which were calculated from the other two phases.

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.1 c k , 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 / 3 The average phase errors of the inverse matrix based, cross-correlation, and autocorrelation 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 autocorrelation 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.
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][17][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 autocorrelation 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. 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.

Conclusions
In this paper, a phase-only correlation method was introduced to optimize the modulation spatial frequency vector 0 k . As the estimation of the modulation frequency is the first step in the reconstruction process, the correct modulation frequency vector 0 k is the basis of all subsequent steps. This phase-only correlation method outperforms previous methods and significantly enhances the peak signal located at 0 k .
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 ( )