## Abstract

The artefact-free reconstruction of structured illumination microscopy images requires precise knowledge of the pattern phases in the raw images. If this parameter cannot be controlled precisely enough in an experimental setup, the phases have to be determined *a posteriori* from the acquired data. While an iterative optimisation based on cross-correlations between individual Fourier images yields accurate results, it is rather time-consuming. Here I present a fast non-iterative technique which determines each pattern phase from an auto-correlation of the respective Fourier image. In addition to improving the speed of the reconstruction, simulations show that this method is also more robust, yielding errors of typically less than *λ*/500 under realistic signal-to-noise levels.

© 2013 OSA

## 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

*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*is the strength of the

_{m}*m*harmonic and determines the pattern contrast.

^{th}*ϕ*is the sought-after pattern phase of the

_{n}*n*

^{th}image. In Fourier space Eq. (1) becomes

*C̃*(

_{m}*k⃗*) are the different information components.

Regarding the individual raw (Fourier) images *D̃ _{n}*(

*k⃗*) as elements of an image vector $\overrightarrow{\tilde{D}}\left(\overrightarrow{k}\right)$ and likewise the individual components

*C̃*(

_{m}*k⃗*) as elements of a vector $\overrightarrow{\tilde{C}}\left(\overrightarrow{k}\right)$, we can rewrite Eq. (2) in matrix notation:

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

*Q*orientations of the pattern (indicated by pattern vectors

*p⃗*) leading to separated components

_{q}*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)

*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**^{−1} 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

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, this is different. We get

*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*

^{th}image from the argument (or angle) of the complex valued

*𝒟*:

_{n}#### 4.1. Periodic samples

The assumption that *𝒞 _{m}*

_{,}

_{m′}_{≠}

_{m}_{+1}is negligible compared to

*𝒞*

_{m}_{,}

_{m}_{+1}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

*a*

_{−1}

*a*

_{0}=

*a*

_{0}

*a*

_{1}and

*𝒞*

_{−1,0}=

*𝒞*

_{0,1}. Assuming a constant total image intensity

*a*

_{0}, any change in the pattern contrast

*c*= 2

*a*

_{1}can be detected by comparing the change in magnitude of

*𝒟*between different raw images:

_{n}## 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

*λ*= 515 nm; numerical aperture

_{em}*NA*= 1.4; refractive index of the embedding medium

*n*= 1.52. This corresponds to a maximum detectable spatial frequency of (184 nm)

_{r}^{−1}, 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.

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)^{−1}, 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)^{−1}, 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,}

*, 240° +*

_{p}*ϕ*

_{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

^{5}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}_{,}

*} (i.e., three raw SIM images for one pattern orientation (index*

_{l}*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}_{,}

*I subtracted the dataset’s known real phases from the phases determined by the weighted autocorrelations, yielding the individual remaining phase errors Δ*

_{l}*ϕ*

_{1,}

_{q}_{,}

_{p}_{,}

*, Δ*

_{l}*ϕ*

_{2,}

_{q}_{,}

_{p}_{,}

*and Δ*

_{l}*ϕ*

_{3,}

_{q}_{,}

_{p}_{,}

*. I then calculated the standard deviation of these three individual remaining phase errors, ${\epsilon}_{q,p,l}={\u3008{\left(\mathrm{\Delta}{\varphi}_{n,q,p,l}-{\u3008\mathrm{\Delta}{\varphi}_{m,q,p,l}\u3009}_{m}\right)}^{2}\u3009}_{n}^{1/2}$, where 〈..〉*

_{l}*denotes the mean of an expression over the index*

_{n}*n*. This

*ε*

_{q}_{,}

_{p}_{,}

*is called the phase error of the dataset*

_{l}**D̃**

_{q}_{,}

_{p}_{,}

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

_{l}**D̃**

_{q}_{,}

_{p}_{,}

*. For each of the 51 photon levels (*

_{l}*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}_{,}

*, calling this the*

_{p}*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, $\mathrm{\Delta}{E}_{l}={\u3008{\left({\epsilon}_{q,p,l}-{E}_{l}\right)}^{2}\u3009}_{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*. For comparison the dark blue line and shaded blue area show the results for the

_{l}*iterative*approach, corresponding to the data published in [5].

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

*λ*> 500 nm; numerical aperture

_{em}*NA*= 1.4, oil immersion. SIM patterns were generated using two-beam illumination with frequencies |

*p⃗*| = {(310.6 nm)

^{−1}, (311.4 nm)

^{−1}, (301.5 nm)

^{−1}} 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.

## 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*^{2}*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]