Non-iterative determination of pattern phase in structured illumination microscopy using auto-correlations in Fourier space

Kai Wicker

doi:10.1364/OE.21.024692

1. Introduction

Structured illumination microscopy (SIM) is a technique for acquiring superresolution images in fluorescence microscopy [1–3]. This is achieved through the moiré effect: Projecting fine periodic patterns into the sample leads to a down-modulation of high spatial frequencies, which correspond to fine sample detail, i.e., these frequencies are shifted into the pass-band of the optical transfer function (OTF). The resolution is enhanced by separating this information from the unmodulated frequencies and computationally shifting it to its true location in Fourier space, thus increasing the support of the OTF.

This separation of superimposed information components requires the acquisition of several structured illumination raw images for different phases of the periodic illumination pattern. A good separation and thus an artefact-free reconstruction can only be achieved, if the pattern phases in the individual images are known with high precision. If this precision cannot be guaranteed by the experimental setup, the phases have to be determined a posteriori from the acquired data. Shroff et al. showed that the pattern phase of each individual raw image can be determined from the phase of the peak at the pattern frequency in Fourier space [4]. However, for very fine patterns, which are required for maximum resolution enhancement, this peak may be lost in noise, as the OTF transfer strength decreases with higher frequencies. Worse, due to the Stokes’ shift or for imaging modalities such as total internal reflection fluorescence (TIRF) the pattern frequency may lie outside the support of the detection OTF, making it inaccessible for this approach. In these cases an iterative optimisation using cross-correlations between separated information components can still yield accurate results [5]. Although the time-consuming re-calculation of correlations underlying this technique can be avoided through a mathematical trick, the iterative nature of the technique still leads to relatively long optimisation times compared to the more direct phase-of-peak technique by Shroff et al.

In this paper I present a method for calculating the pattern phase for an individual image through a single auto-correlation of its corresponding Fourier image, thereby combining the performance of the correlative approach with the speed of the non-iterative phase-of-peak method. Using simulations I analyse the performance of this approach and compare it to the previously published iterative technique [5]. For simplicity reasons I limit myself to two-dimensional samples and images, but the approach can be extended to three dimensions in a straight-forward way.

2. Image formation and reconstruction in SIM

For a harmonic pattern, the n^th SIM raw image can be generally written as

D_{n} (\vec{r}) = \sum_{m = - M}^{M} [S (\vec{r}) a_{m} exp {ι (2 π m \vec{p} \cdot \vec{r} + m ϕ_{n})}] \otimes h (\vec{r}),

where r⃗ denotes both sample as well as image coordinates (assuming a magnification of one). S(r⃗) denotes the sample fluorophore density, ⊗ the convolution operator, ⊗ the imaginary unit, p⃗ the base pattern frequency and M the number of harmonics in the periodic intensity pattern. a_m is the strength of the m^th harmonic and determines the pattern contrast. ϕ_n is the sought-after pattern phase of the n^th image. In Fourier space Eq. (1) becomes

{\tilde{D}}_{n} (\vec{k}) = \sum_{m = - M}^{M} exp {ι m ϕ_{n}} \underset{= : {\tilde{C}}_{m} (\vec{k})}{\underset{︸}{a_{m} \tilde{S} (\vec{k} - m \vec{p}) \tilde{h} (\vec{k}),}}

where ∼ denotes the Fourier transform of a function and C̃_m(k⃗) are the different information components.

Regarding the individual raw (Fourier) images D̃_n(k⃗) as elements of an image vector $\vec{\tilde{D}} (\vec{k})$ and likewise the individual components C̃_m(k⃗) as elements of a vector $\vec{\tilde{C}} (\vec{k})$ , we can rewrite Eq. (2) in matrix notation:

\vec{\tilde{D}} (\vec{k}) = M \vec{\tilde{C}} (\vec{k}),

where the elements of the component mixing matrix M are defined as

M_{n m} = exp {ι m ϕ_{n}} .

If the pattern phases are known, the components can be retrieved from the raw images by a simple inversion of Eq. (3):

\vec{\tilde{C}} (\vec{k}) = M^{- 1} \vec{\tilde{D}} (\vec{k}) .

For isotropic resolution enhancement the process of image acquisition and component separation has to be repeated for a total of Q orientations of the pattern (indicated by pattern vectors p⃗_q) leading to separated components C̃_q_,_m(k⃗), where the additional index q indicates the orientation used.

All separated components are then shifted back to their correct frequencies in Fourier space and then recombined using a combination of weighted-averaging, generalised Wiener filtering and frequency apodisation [5, 6], leading to the sample estimation (i.e., reconstructed Fourier image)

{\tilde{S}}^{'} (\vec{k}) = \frac{\sum_{q, m} {a_{m} {\tilde{h}}^{*} (\vec{k} + m {\vec{p}}_{q}) {\tilde{C}}_{q, m} (\vec{k} + m \vec{p} q)}}{\sum_{q^{'}, m^{'}} {{| a_{m^{'}} \tilde{h} (\vec{k} + m^{'} {\vec{p}}_{q^{'}}) |}^{2}} + w} \tilde{A} (\vec{k}),

where w is the wiener filter parameter and Ã(k⃗) is the apodisation filter.

3. Current phase optimisation strategies

For the correct separation of components it is crucial that the pattern phase in the individual raw images be known with high precision. This phase can be retrieved from the phase of the peak at the pattern frequency p⃗ in Fourier space [4]. However, for fine patterns this method fails, as the peak is attenuated by the decreasing OTF strength and its phase is overshadowed by random phases stemming from Poisson noise.

For such fine patterns the pattern phases can be robustly determined by an iterative optimisation of the inverse mixing matrix M⁻¹ used for component separation [5]. Here the quality of separation is measured after each iteration by correlating the separated components. This cross-correlation based approach can still determine the pattern phase even if the pattern peaks themselves lie beyond the band limit of the OTF, as its effect is still detectable in the acquired data.

Although this approach is relatively fast its iterative nature inevitably leads to longer computation times than the more direct approach of Shroff et al. Furthermore, as with most minimisation problems, it can potentially yield wrong phases as a result of getting stuck in a local minimum.

4. Single-step correlative phase determination

Alternatively, we can analyse weighted auto-correlations of filtered Fourier images in order to retrieve the individual pattern phases. To do so, we first take the Fourier transformed image and filter (i.e., multiply) it with the complex conjugated OTF: D̃′_n(k⃗) = D̃_n(k⃗)h̃^*(k⃗). We then analyse the auto-correlation (denoted by ⍟) of this filtered Fourier image at frequency p⃗. Using Eq. (2) we get

\begin{array}{l} 𝒟_{n} & = [{\tilde{D}}_{n}^{'} ⍟ {\tilde{D}}_{n}^{'}] (\vec{p}) \\ = \int {\tilde{D}}_{n}^{'} (\vec{k}) {\tilde{D}}_{n}^{' *} (\vec{k} + \vec{p}) d^{3} k \\ = \sum_{m = - M}^{M} \sum_{m^{'} = - M}^{M} {exp {ι (m - m^{'}) ϕ_{n}} a_{m} a_{m^{'}} \\ \times \underset{= : 𝒞_{m, m^{'}}}{\underset{︸}{\int \tilde{S} (\vec{k} - m \vec{p}) {\tilde{S}}^{*} (\vec{k} - (m^{'} - 1) \vec{p}) {| \tilde{h} (\vec{k}) |}^{2} {| \tilde{h} (\vec{k} + \vec{p}) |}^{2} d^{3} k}}} . \end{array}

The multiplication of the images with the complex conjugated OTF has two effects. Firstly, this multiplication attenuates the white noise of the original Fourier images according to the signal strength of the OTF. This way noise will contribute less to the final auto-correlation, making it in effect a weighted auto-correlation as defined in [5] (ignoring the division by the constant sum of weights). Secondly, it neutralises any disturbing phase effects of the OTF, which only is real-valued if the corresponding PSF is symmetric. This is necessary for the retrieval of the pattern phase, as is shown below.

If we take a closer look at the last line of Eq. (7), we find that we can differentiate between two cases. The first one is m′ = m + 1. In this case we are correlating two shifted sample components which have a relative frequency distance of p⃗ and we are analysing this correlation also at frequency p⃗, where there is thus a strong correlation. The integral term becomes

𝒞_{m, m + 1} = \int {| \tilde{S} (\vec{k} - m \vec{p}) |}^{2} {| \tilde{h} (\vec{k}) |}^{2} {| \tilde{h} (\vec{k} + \vec{p}) |}^{2} d^{3} k,

which is purely real and positive. For the second case, m′ ≠ m + 1, this is different. We get

𝒞_{m, m^{'} \neq m + 1} = \int \tilde{S} (\vec{k} - m \vec{p}) {\tilde{S}}^{*} (\vec{k} - (m^{'} - 1) \vec{p}) {| \tilde{h} (\vec{k}) |}^{2} {| \tilde{h} (\vec{k} + \vec{p}) |}^{2} d^{3} k .

Here we are correlating shifted sample components which have a relative distance not equal to p⃗, but we are still looking at the correlation at p⃗, which should thus be low. Although its expectation value is not zero, but rather depends on the auto-correlation function of the Fourier transformed sample, it is low enough compared to the strong correlations for m′ = m + 1 so that we can ignore its contribution to Eq. (7). Considering hence only the m′ = m + 1 terms, this equation can be simplified to

𝒟_{n} \approx exp {- ι ϕ_{n}} \sum_{m = - M}^{M - 1} {a_{m} a_{m + 1} 𝒞_{m, m + 1}} .

As the sum is real and positive we can retrieve the pattern phase in the n^th image from the argument (or angle) of the complex valued 𝒟_n:

ϕ_{n} = - arg {𝒟_{n}} .

4.1. Periodic samples

The assumption that 𝒞_m_,_m′_≠_m₊₁ is negligible compared to 𝒞_m_,_m₊₁ may not always be fulfilled. This can be the case for samples which have strong periodic structures with frequencies corresponding to those present in the illumination pattern, i.e., mp⃗. Here most algorithms will have difficulties differentiating between illumination and sample structures and yield results biased towards the phase of the sample structure. This situation should ideally be avoided, e.g., by slightly rotating the illumination pattern or the sample. If this is not feasible, the impact of this periodic sample structure may be partially alleviated by estimating its contribution through an auto-correlation of the OTF-filtered wide-field image at frequency p⃗. This wide-field image can in turn be approximated as an average of all acquired raw images. This way the 𝒞_0,0 contribution can be avoided. However, for most samples this will not be necessary.

4.2. Determining changes in pattern contrast

Besides pattern phase this approach can in principle also determine fluctuations in the pattern contrast between raw images. For harmonic patterns (M = 1) Eq. (10) becomes

𝒟_{n} \approx exp {- ι ϕ_{n}} 2 a_{0} a_{1} 𝒞_{0, 1},

because a₋₁a₀ = a₀a₁ and 𝒞_−1,0 = 𝒞_0,1. Assuming a constant total image intensity a₀, any change in the pattern contrast c = 2a₁ can be detected by comparing the change in magnitude of 𝒟_n between different raw images:

a_{1}^{(n^{'})} / a_{1}^{(n)} = | 𝒟_{n^{'}} | / | 𝒟_{n} | .

Pattern contrast does not usually change for different pattern phases. There are however a few cases where the contrast can change. Examples are non-linear SIM using fluorescence saturation [7,8], where the effective contrast can change with fluctuations in the illumination intensity, or polarisation-coded SIM (picoSIM) [9], where the grating contrast may vary due to polarisation characteristics of the imaging optics. In these cases the determination of the relative pattern contrast can facilitate the correct separation of components and thus an artefact-free reconstruction.

5. Performance of the algorithm

I evaluated the performance of the single step phase optimisation using simulated data with known pattern phases. For comparability with our iterative approach [5] I used the same simulation parameters as well as pattern phases.

5.1. Simulation parameters

The experimental parameters used for the simulation were: excitation wavelength λ_ex = 488 nm; emission wavelength λ_em = 515 nm; numerical aperture NA = 1.4; refractive index of the embedding medium n_r = 1.52. This corresponds to a maximum detectable spatial frequency of (184 nm)⁻¹, defined by the detection OTF. The pixel size in the simulated raw image corresponds to 65 nm in sample space. The synthetic sample used for the simulations is shown in Fig. 1. No camera offset or other background was assumed in the simulation.

Fig. 1 The synthetic sample used for simulations.

Download Full Size | PDF

The illumination patterns were simulated for two-beam (i.e., sinusoidal, M = 1) illumination. As in [5], I simulated illumination with two different pattern frequencies |p⃗|: a fine pattern with a spatial frequency of |p⃗| = 2π (210 nm)⁻¹, or 87.6% of the maximum frequency supported by the OTF; and a very fine pattern with a spatial frequency of |p⃗| = 2π (185 nm)⁻¹, or 99.4% of the OTF support.

For each pattern frequency I simulated illumination patterns for three different pattern orientations (0°, 60°, 120°). The ideal pattern phases (0°, 120°, 240°) were varied randomly using a Gaussian distribution with a standard deviation of 10°. I generated raw images using twenty such sets of randomised pattern phases, {0° + ϕ_1,_p, 120° + ϕ_2,_p, 240° + ϕ_3,_p}, p = 1..20. For each of these, I simulated noisy images for 51 different signal-to-noise levels, using Poisson noise and an expectation value of 10^l^/10, l = 0..50, i.e., between 1 and 10⁵ photons, in the brightest pixel. The raw images generated this way were the same as the ones used in [5].

5.2. Analysis

Although this correlative method retrieves the phase of every individual raw image, I calculated the average phase error per dataset D̃_q_,_p_,_l = {D̃_1,_q_,_p_,_l, D̃_2,_q_,_p_,_l, D̃_3,_q_,_p_,_l} (i.e., three raw SIM images for one pattern orientation (index q), with one set of randomised phases (index p), for one particular photon level (index l)) rather than for individual images. This was done according to [5] for comparability with the iterative approach. For each of these datasets D̃_q_,_p_,_l I subtracted the dataset’s known real phases from the phases determined by the weighted autocorrelations, yielding the individual remaining phase errors Δϕ_1,_q_,_p_,_l, Δϕ_2,_q_,_p_,_l and Δϕ_3,_q_,_p_,_l. I then calculated the standard deviation of these three individual remaining phase errors, $ε_{q, p, l} = {〈 {(Δ ϕ_{n, q, p, l} - {〈 Δ ϕ_{m, q, p, l} 〉}_{m})}^{2} 〉}_{n}^{1 / 2}$ , where 〈..〉_n denotes the mean of an expression over the index n. This ε_q_,_p_,_l is called the phase error of the dataset D̃_q_,_p_,_l. This approach disregards any global pattern phase offsets, which is only relevant for the iterative approach, as it determines relative phase steps between image rather than the absolute phase of every individual image. This phase error is the measure of how well the phases were found for the particular dataset D̃_q_,_p_,_l. For each of the 51 photon levels (l) I then calculated the mean of the datasets’ phase errors over all 20 different phase variations (p) and three orientations (q), E_l = 〈ε_q_,_p_,_l〉_q_,_p, calling this the average phase error for a particular photon level (l). As a measure of robustness, I furthermore calculated the standard deviation from this average phase error, $Δ E_{l} = {〈 {(ε_{q, p, l} - E_{l})}^{2} 〉}_{q, p}^{1 / 2}$ .

5.3. Results

Figure 2 shows the result of this analysis. The red line shows the average phase error E_l of the new non-iterative (single step) approach, with the shaded red area marking its standard deviation ΔE_l. For comparison the dark blue line and shaded blue area show the results for the iterative approach, corresponding to the data published in [5].

Fig. 2 Performance of phase optimisation. The red line shows the average phase error over the photon level (expectation value of the brightest pixel in a set of three raw images) of the non-iterative (single step) auto-correlation technique and as such is a measure for the accuracy for the technique. The red shaded area shows the standard deviation around this average phase error and is a measure for the robustness of the technique. The blue line and shaded area show the corresponding average phase errors for the iterative approach [5]. (a) shows the results for the coarser pattern with a period of 210 nm, (b) shows the finer pattern with a period of 185 nm. For realistic photon levels (> 100) the single step approach yields more accurate and more robust results than the iterative approach.

Download Full Size | PDF

Figure 2(a) shows the results for the coarser of the two gratings. For very low signal-to-noise levels (less than 3 photons expected in the brightest pixel) the iterative approach manages to retrieve the phase better than the single step approach. However, for these extremely low light situations, both techniques yield phase errors which are too high to guarantee artefact-free reconstructions. For more realistic values of 10 photons or more both methods yield similar average phase errors, reaching values below 1° for higher photon numbers. This corresponds to an error of less than 0.6 nm or λ/800. Furthermore, its significantly lower standard deviation indicates a much stronger robustness of the single step approach, which unlike the iterative approach never produced any outliers.

A similar tendency can be observed for the finer grating in Fig. 2(b). For very low signal-to-noise levels the single step approach fails to retrieve the phase correctly. The iterative approach performs slightly better, but even here average phase errors of above 5° would likely lead to some artefacts in the reconstructed images. For photon numbers of 100 and higher the single step approach yields better results than the iterative one, reaching average phase errors of less than 1°, which corresponds to an error of about 0.5nmor λ/1000. Here, too, the much smaller standard deviation indicates a greater robustness of the single step approach, which again did not produce any outliers.

6. Reconstruction of experimental data

I applied single step phase optimisation to the reconstruction of experimental SIM images of autofluorescent lipofuscin accumulating in the retinal pigment epithelial (RPE). Experimental parameters were [10]: excitation wavelength λ_ex = 488 nm; emission wavelength λ_em > 500 nm; numerical aperture NA = 1.4, oil immersion. SIM patterns were generated using two-beam illumination with frequencies |p⃗| = {(310.6 nm)⁻¹, (311.4 nm)⁻¹, (301.5 nm)⁻¹} for the three orientations α = {0.9°, 120.8°, 59.7°}. For comparison I reconstructed the acquired data using a wrong assumption of equidistant phase steps and no optimisation, optimisation using the iterative approach and optimisation using the single step approach.

Figure 3 shows the results of this reconstruction. Figure 3(a) shows the reconstruction under the wrong assumption of equidistant steps of the pattern phase between the raw images. These wrong phases were not optimised; consequently the reconstructed images exhibits strong artefacts, such as the honeycomb-like structures visible in the close-up in Fig. 3(a). In the corresponding Fourier image, Fig. 3(b), residual pattern peaks stemming from an imperfect separation of components are visible. Images reconstructed using iterative phase optimisation, Fig. 3(c), as well as the new single step phase optimisation, Fig. 3(e), both yield artefact-free reconstructions, with no visible differences discernible between them. Their respective Fourier images, Figs. 3(d) and 3(f), no longer exhibit any unwanted residual peaks.

Fig. 3 SIM reconstructions of retinal pigment epithelial (RPE). (a) shows the reconstruction with the (obviously) wrong assumption of equidistant phase stepping, leading to strong artefacts. In the corresponding Fourier image (b) residual pattern peaks can be observed, which stem from a bad separation of components. Both iterative phase optimisation (c) and single step phase optimisation (e) lead to an artefact-free reconstruction with no visibly discernible differences. Their respective Fourier images (d) and (f) are void of any unwanted residual pattern peaks. Images reconstructed from data provided by Gerrit Best from the group of Christoph Cremer.

Download Full Size | PDF

7. Discussion

Imprecise knowledge of pattern phases is one of the main causes for artefacts in SIM image reconstructions. This is a significant problem in systems which cannot guarantee the necessary control of pattern phase. But also in more robust systems, bleaching of the illumination structure into the sample may lead to a shift in the effective position of illumination patterns in subsequent images. It can therefore be necessary to determine the pattern phases a posteriori from the data. Through correlative approaches the pattern phase can be determined even in situations where the pattern is no longer visibly detectable in the acquired raw images (e.g., for extremely fine gratings or low signal-to-noise levels). I have presented a non-iterative approach which determines the pattern phase through weighted auto-correlations of individual raw Fourier images. Analysis of simulated SIM datasets shows that this approach yields accurate results with high reliability under realistic signal-to-noise conditions. Images reconstructed using this technique could not be visibly discerned from images reconstructed using the iterative approach [5]. In addition to improved robustness, this technique is also significantly faster than the iterative approach. The main reason for this is the non-iterative nature of the optimisation. Furthermore, phase determination using weighted auto-correlations requires the computation of only N autocorrelations, whereas the optimised iterative approach requires the calculation of N²L auto- and cross-correlations before the start of the iterative optimisation, N being the number of raw images and L the largest whole number for which the shifted OTF h̃(k⃗ − Lp⃗) still has overlapping support with the unshifted OTF h̃(k⃗).

The improved accuracy, robustness and speed over previous approaches opens the possibility of on-the-fly image reconstructions with phase corrections. A further practical advantage towards this goal is that the phase in one raw image can already be calculated while the following raw images are still being acquired, as the single step is based on auto-correlations and therefore does not simultaneously require all raw images for the optimisation – unlike the iterative approach, which is based on cross-correlations between different images. The technique has thus the potential to significantly improve live-cell SIM applied under imperfect imaging conditions.

Acknowledgments

I would like to thank Rainer Heintzmann for interesting discussions about this technique, as well as Gerrit Best and Christoph Cremer for providing the RPE data used for the experimental reconstructions.

References and links

1. R. Heintzmann and C. Cremer, “Laterally modulated excitation microscopy: improvement of resolution by using a diffraction grating,” Proc. SPIE 3568, 185–196 (1999). [CrossRef]

2. M. G. L. Gustafsson, “Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy,” J. Microsc. 198, 82–87 (2000). [CrossRef] [PubMed]

3. J. T. Frohn, “Super-resolution fluorescence microscopy by structured light illumination,” Ph.D. thesis, Eidgenössische Technische Hochschule Zürich, Switzerland (2000).

4. S. A. Shroff, J. R. Fienup, and D. R. Williams, “Phase-shift estimation in sinusoidally illuminated images for lateral superresolution,” JOSA A 26, 413–424 (2009). [CrossRef] [PubMed]

5. K. Wicker, O. Mandula, G. Best, R. Fiolka, and R. Heintzmann, “Phase optimisation for structured illumination microscopy,” Opt. Express 21, 2032–2049 (2013). [CrossRef] [PubMed]

6. 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, 4957–4970 (2008). [CrossRef] [PubMed]

7. R. Heintzmann, T. M. Jovin, and C. Cremer, “Saturated patterned excitation microscopy - a concept for optical resolution improvement,” JOSA A 19, 1599–1609 (2002). [CrossRef]

8. M. G. L. Gustafsson, “Nonlinear structured-illumination microscopy: wide-field fluorescence imaging with theoretically unlimited resolution,” PNAS 102, 13081–13086 (2005). [CrossRef] [PubMed]

9. K. Wicker and R. Heintzmann, “Single-shot optical sectioning using polarization-coded structured illumination,” J. Opt. 12, 084010 (2010). [CrossRef]

10. G. Best, R. Amberger, D. Baddeley, T. Ach, S. Dithmar, R. Heintzmann, and C. Cremer, “Structured illumination microscopy of autofluorescent aggregations in human tissue,” Micron 42, 330–335 (2011). [CrossRef]

Non-iterative determination of pattern phase in structured illumination microscopy using auto-correlations in Fourier space

Abstract

1. Introduction

2. Image formation and reconstruction in SIM

3. Current phase optimisation strategies

4. Single-step correlative phase determination

4.1. Periodic samples

4.2. Determining changes in pattern contrast

5. Performance of the algorithm

5.1. Simulation parameters

5.2. Analysis

5.3. Results

6. Reconstruction of experimental data

7. Discussion

Acknowledgments

References and links

Cited By

Figures (3)

Equations (13)

Optics Express