Image Filtering in Structured Illumination Microscopy Using the Lukosz Bound

Various aspects of image filtering affect the final image quality in Structured Illumination Microscopy, in particular the regularization parameter and type of regularization function, the relative height of the side bands, and the shape of the apodization function. We propose an apodiza-tion filter without adjustable parameters based on the application of the Lukosz bound in order to guarantee a non-negative point spread function. Simulations of digital resolution charts and experimental data of chromatin structures and of actin filaments show artefact free reconstructions for a wide range of filter parameters. In general, a trade-off is observed between sharpness and noise suppression. A guide to super-resolution fluorescence microscopy, " J. Method of obtaining optical sectioning by using structured light in a conventional microscope, " Opt. Wide-field optically sectioning fluores-cence microscopy with laser illumination, " J. Microsc. 197, 1–4 (2000). 6. R. Heintzmann and C. Cremer, " Laterally modulated excitation microscopy: improvement of resolution by using a diffraction grating, " Proc. Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy, " J. True optical resolution beyond the Rayleigh limit achieved by standing wave illumination, " Proc. Saturated patterned excitation microscopy-a concept for optical resolution improvement, " J. Wide-field high-resolution structured illumination solid immersion fluorescence microscopy, " Opt. Three-dimensional resolution doubling in wide-field fluorescence microscopy by structured illumination, " Biophys. Super-resolution video microscopy of live cells by structured illumination, " Nat. Super-resolution 3d microscopy of live whole cells using structured illumination, " Nat. Structured illumination in total internal reflection fluorescence microscopy using a spatial light modulator, " Opt. Line scan-structured illumination microscopy super-resolution imaging in thick fluorescent samples, " Opt. A quantitative comparison of image restoration methods for confocal microscopy, " J. algorithm with total variation regularization for 3D confocal microscope deconvolution, " Microsc. Exact deconvolution for multiple convolution operators–an overview, plus performance characterizations for imaging sensors, " Proc. Phase-shift estimation in sinusoidally illuminated images for lateral superresolution, " J. Resolution enhancement in standing-wave total internal reflection microscopy: a point-spread-function engineering approach, " J. Resolution in structured illumination microscopy: a probabilistic approach, " J.image: a scientific image processing toolbox for MATLAB, " (1999-).'s disease derived cell lines HDLM-2 and L-428: comparison of morphology, immunological and isoenzyme profiles, " Leuk. The influence of the regularization parameter and the first estimate on the performance of Tikhonov regularized non-linear image restoration algorithms, " J.


Introduction
The diffraction limit to resolution is broken in optical microscopy by numerous techniques developed in recent years [1,2].Localization microscopy (PALM, STORM) circumvents the diffraction limit by localizing switchable fluorophores with nanometer precision, STED improves resolution by reducing the fluorescence emission to a sub-diffraction sized region by de-exciting fluorophores near the boundary of the excitation spot with a ring shaped depletion spot.The technique of Structured Illumination Microscopy (SIM) extends the diffraction limit by a factor of two by illuminating the sample with a series of periodic patterns of closely spaced lines of different orientations and phases.The final super-resolution image is reconstructed from this set of recorded images.
The idea to use heterodyne detection to enhance the effective spatial resolution in microscopic imaging was introduced in the 1960s by Lukosz [3].The application of these techniques to achieve optical sectioning rather than resolution improvement were further developed in the mid-1990s [4,5].Implementations to increase the effective lateral resolution were developed around the turn of the century [6][7][8][9].Two-dimensional SIM was further extended by exploiting the non-linear fluorescence response [10,11] and by combination with solid immersion imaging [12].The principles of SIM were implemented in 3D several years ago [13].More recently the temporal resolution of the technology has improved in both two- [14] and three- [15] dimensional SIM, based on the use of Spatial Light Modulators (SLMs) [16].A reduction of out-of-focus light can be achieved using a combination with line scanning [17].
Image reconstruction in SIM boils down to first disentangling the different spatial frequency bands, and subsequently stitching these bands together to extend the bandwidth of the reconstructed image.This step is often integrated with a filtering approach for suppressing noise and for tailoring the effective response as a function of spatial frequency.Commonly used filters in SIM image reconstruction are borrowed from the field of image restoration, i.e. methods to extend resolution beyond the diffraction limit by incorporating prior knowledge such as nonnegativity [18], noise models [19], background level [20], and image smoothness [21].The default filtering technique in SIM is (generalized) Wiener filtering [22][23][24], with a signal-tonoise ratio (SNR) that is assumed constant across the entire spatial frequency range.The free parameter in this type of filter is the so-called regularization parameter, which measures the 'strength' of the prior knowledge.It is not a priori clear what value this parameter should have.Moreover, if the parameter takes a small value, several artefacts appear, in particular edge ringing, halos surrounding bright objects, and negative pixel values.The origin of the artefacts is an effective transfer function which is flat up to a point close to the extended diffraction limit.A suitable counter measure is therefore the application of a so-called apodization filter, which suppresses the response at high spatial frequencies compared to the response at low spatial frequencies [8,13].Apodization filters are applied ad-hoc and often introduce further free parameters which have optimal settings that are a priori unclear.An additional image artefact common to SIM is the enhancement of low-frequency noise structures in the reconstructed image.The effect of the different free parameters on this artefact have not been explored so far.
The goal of the current paper is to investigate the role of filtering and filter parameters on the image quality of SIM.The issues of edge ringing and negative Point Spread Function (PSF) values were treated systematically by Lukosz half a century ago [25,26].He derived a set of conditions on the Optical Transfer Function (OTF) of the system that must be satisfied in order to guarantee a non-negative PSF.This Lukosz bound to the OTF thus provides a necessary (but not a sufficient) condition for having a non-negative PSF.Here, we propose to use the Lukosz bound for apodization in SIM image reconstruction and filtering.We will present and evaluate different ways to incorporate the Lukosz bound.The advantage of our proposal is that it provides a systematic approach as opposed to ad-hoc procedures and that it explicitly rules out free apodization parameters.
The content of the paper is as follows.The key results of Lukosz are summarized in section 2, image formation and filtering for SIM is presented in section 3. The application of the Lukosz bound in image reconstruction is described in section 4. Simulations are presented in section 5, experiments in section 6.The paper concludes with a discussion of the results and an outlook to further work in section 7.7

The Lukosz bound
Lukosz considered linear shift-invariant band-limited imaging systems, characterized by a PSF H (u), or equivalently its Fourier transform, the OTF Ĥ (v).Here, and thoughout the paper we use the designation u for the spatial coordinates, and v for the spatial frequencies.The PSF and OTF are a Fourier transform pair, for which we use the convention: (1) We will first consider 1D band-limited imaging systems, with a cut-off spatial frequency q c .Lukosz showed in his papers [25,26] that the Modulation Transfer Function (MTF), the absolute value of the OTF, must be smaller than the Lukosz bound in order to have a non-negative PSF, i.e. a PSF that is positive-definite everywhere: This function describes a staircase with infinite, progressively smaller, steps toward the origin, as shown in Fig. 1.It is a necessary condition that the MTF should be below the 'Lukosz stairs', but not a sufficient condition.An additional condition is that the average of the OTF over its support must be less than 1/2.A suitable curve in practice is the continuous curve connecting the lower corners of the stairs: This practical Lukosz bound OTF is shown in Fig. 1, as well as the incoherent OTF and the corresponding Lukosz bound PSF and incoherent PSF.It turns out that indeed a positive definite PSF is obtained, at least for spatial frequencies up to the second diffraction minimum.The achieved sharpness is only marginally less than the sharpness of the incoherent PSF.
Next, consider a 2D imaging system.Now, the cut-off spatial frequency q c (φ ) can depend on the azimuthal angle φ .According to Lukosz the MTF for a spatial frequency vector (v cos ψ, v sin ψ) should be below the product of the two 1D Lukosz bounds along the orthogonal axes (cos φ , sin φ ) and (− sin φ , cos φ ) for any value of φ .The 2D bound for azimuthal angle φ is thus: where the spatial frequencies in the rotated coordinate frame are v = (v cos (ψ − φ ) , v sin (ψ − φ )).The MTF should be below this bound for any value of φ giving the overall 2D Lukosz bound as: The cut-off spatial frequency is constant as a function of azimuthal angle in case the OTF is rotationally symmetric.In that case the recipe for the 2D Lukosz bound reduces to: and 0 elsewhere, i.e. the minimum is taken between the product of the two 1D-Lukosz bounds for coordinate frames rotated over angles 0 and π/4.Fig. 1 shows the 2D OTF and PSF according to the Lukosz bound of the rotationally symmetric imaging system, as well as the conventional incoherent OTF and PSF.The Lukosz bound PSF is positive-definite up to the second dark diffraction ring (where it drops slightly below zero to a value around -0.006) and is a bit sharper than the incoherent PSF (full-width at half-maximum (FWHM) is around 5% smaller).The possibly negative pixel values arising from the small negative PSF values around the second dark diffraction ring are not expected to give rise to significant edge ringing or halo artefacts.Also, the effects are mitigated by a background signal which is always non-zero in practice and by blurring through the finite pixel size.
In case the OTF of an imaging system is equal to the Lukosz bound, it is intuitively clear that this provides the most compact non-negative PSF that can be achieved.In a sense the Lukosz bound, therefore, describes the 'ideal' imaging system.It is further noted that the Lukosz bound does not improve much over the standard incoherent response.

Image formation in SIM
Image formation in SIM consists of two distinct steps, both of which will be discussed in this section; 1) the actual physical image acquisition through the microscope, and 2) the postprocessing that yields the final reconstructed image.Here, the spatial coordinates in object space are normalized with λ /NA ob with NA ob the objective numerical aperture, while the co-ordinates in image space are normalized with λ /NA tube with NA tube the tube lens numerical aperture.The imaging system has unit magnification in these normalized coordinates.The spatial frequencies are normalized accordingly with factors NA/λ .The stop of the imaging system thus corresponds to the unit circle, and the normalized cut-off spatial frequency of a conventional widefield fluorescence microscope is then equal to two.The following functions will be used throughout: T (u) , T (v) for the object function/spectrum; W (u) , Ŵ (v) for the illumination function/spectrum; and H (u) , Ĥ (v) for the PSF/OTF of the widefield fluorescence microscope.The symbol ⊗ will denote convolution.

Image acquisition and phase shifting
The premise of SIM is that a periodic illumination function is used that is rotated and shifted with respect to the sample.In particular, a periodic pattern of stripes is used: with W av the average illumination intensity, w a parameter characterizing the modulation depth (optimum value w = 1/2), and where p = 1/(2q) is the spacing of the lines.Such an illumination pattern can be made from the interference pattern of two plane waves traveling at opposite angles with respect to the optical axis.The illumination pattern in the back focal plane of the objective lens (assuming epi-illumination is used) then consists of two dots at positions R(q, 0) and R(−q, 0), with R the pupil radius, from which it follows that q ≤ 1.The Fourier transform of the illumination function consists of a set of discrete delta-peaks: with q m = m(q, 0) and with non-zero Fourier coefficients Typically N r = 3 rotations and N t = 3 translations are used (for 2D-SIM).The recorded image for rotation l and translation n is given by: This can be rewritten in the spatial frequency domain as: with phase shifting matrix: The Fourier transform of each recorded image is a linear sum of bands in spatial frequency space.These bands can be distilled by an inversion procedure: This inversion works provided the number of translations N t ≥ K, with K the number of independent Fourier components of the illumination function.In case the translations are not equidistant, the Moore-Penrose pseudo-inverse can be used to obtain the band images.

Weighted sum of shifted bands reconstruction
The most simple and straightforward method to obtain a final super-resolution image is to take a weighted linear combination of the bands shifted in spatial frequency space [6,7,27,28]: Here (s −2 , s 0 , s 2 ) = S av (s, 1, s) are a set of weight coefficients, which can in principle be chosen at will.The parameter S av can be used to scale the overall effective OTF: such that Ĥlin (0) = 1.The parameter s can be used to tune the relative height of the side bands compared to the central band.It is common to pick this side band height parameter a priori -often as s = 1 or s = 1/w -and not estimate it as part of a filtering process [8,13,24].Interestingly, the band stitching process can be expressed as a convolution operation in spatial frequency space.It follows that this is equivalent to point multiplication of the recorded images with a suitable mask function in real space, and then adding all contributions: with mask function: So, the mask function for each recorded image has the same period, orientation and phase as the illumination pattern for the recorded image, only the modulation depth is characterized by the parameter s instead of w.

Filtering approaches
A next step in sophistication is to apply a noise suppression filter with kernel F (v) to the image obtained by the weighted sum of shifted bands reconstruction method, producing a superresolution image Îapo (v) = F (v) Îlin (v).It urns out that the level of noise suppression for this post filtering approach can be improved by integrating the filtering operation into the process of stitching together the bands in spatial frequency space.This is done by making the filter kernel dependent on the spatial frequency band.This so-called generalized filtering approach gives rise to: The filter kernels can be derived from minimizing a Tikohonov-Miller (TM) type of functional [29]: where the filtered super-resolution image Îgen (v) is considered as the variable in the minimization procedure.The first term of the TM-functional is the 'data misfit term', i.e. the summed squared difference between the data and the forward model, and the second term is the 'regularization term', representing the prior knowledge of the reconstructed object.Here the functions Blm (v) describe the ideal or desired response of the imaging system.They are usually identified with the weighted and shifted microscope OTF s m w m Ĥ (v + R l • q m ), but can in principle be chosen differently.The parameter κ is the regularization parameter and Â (v) is the regularization weight function.We consider: with p an integer, so for p = 0 the total signal energy is regularized, for p = 1 the total signal gradient energy, for p = 2 the total signal Laplacian energy, etc. Minimization of the TMfunctional gives: and an OTF: The conventional choice for the functions Blm (v) = s m w m Ĥ (v + R l • q m ) gives rise to an OTF that is flat up to the cut-off in the limit of a small regularization parameter κ.Such OTF's do not satisfy the Lukosz bound and thus give rise to edge ringing, halo artefacts and negative pixel values.
We point out that these TM-filters correspond to (generalized) Wiener-filters if the signalto-noise power spectrum is used as the regularization function [30].In the literature on SIM these filters with total signal energy regularization (p = 0) are, however, usually referred to as (generalized) Wiener-filters [13].We will follow this nomenclature, although we believe it is basically incorrect.

Data misfit weight function
Light originating from out-of-focus layers gives rise to a blurred background and hence a reduced contrast in the captured images.Moreover it may lead to honeycomb artefacts in the reconstruction because the low frequency blur is transfered to the regions close to the side band peaks in spatial frequency space.Suppression of out-of-focus light can be done by subtracting a low-pass filtered version of the acquired images before feeding the acquired images into the image reconstruction procedure.An alternative is to modify the filters with a weight function for suppressing the effect of the data misfit function on the final filter close to the side band peaks in spatial frequency space [31].It follows that the regularization dominates over the data misfit term for the spatial frequencies for which the out-of-focus light has a significant influence.Here we show that this modification of the filtering step can be incorporated into the framework of TM-reconstruction.The data misfit term in the TM-functional must be changed in such a way that the relative weight of that term compared to the regularization term decreases with decreasing spatial frequency.The TM-functional for generalized filtering then becomes: with ĝ (v) a function that increases with increasing spatial frequency.For example, a Gaussian data misfit weight gives rise to: where the magnitude α and width σ are in principle adjustable parameters.Minimization of the TM-functional leads to: and an OTF: A consequence of relatively large values of the Gaussian magnitude α (close to one) is that dips in the effective OTF arise at the spatial frequencies around the center of the main band and of the side bands.The dips can thus lead to violations of the Lukosz bound, as that is a monotonically decreasing function.Although the resulting OTF does not lead to an imaging system that gives rise to a non-negative output for any conceivable object structure, it may be expected to do so for the particular case where a large out-of-focus contribution is present, i.e. for thick samples.The Gaussian function should, therefore, be tailored to the thickness and density of the sample in question by varying α and/or σ for each case.

The Lukosz bound in SIM filtering
The filtering approaches presented so far are clearly inappropriate and some form of apodization is in order.According to Lukosz the undesired artefacts are eliminated if the final OTF is below the Lukosz bound.In this section, we will describe how to incorporate the Lukosz bound in the filtering step.

2D SIM Lukosz bound
The bound Λ (v) for (2D) SIM can be derived along the lines presented in sect. 2 provided the spatial frequency cut-off is known.For the case of SIM this corresponds to the combined support of the central band and the side bands.For the typical case of N r = 3 rotations this results in a 'flower'-shaped support parametrized by: where q ≤ 1 determines the pitch of the illumination pattern and where: is the reduced azimuth angle corresponding to each circle segment of the flower-shaped support.Fig. 2 shows the OTF and PSF according to the 2D-SIM Lukosz bound for a value q = 0.9.The support of the OTF is reduced somewhat to a hexagonal shape, the spatial frequency cut-off in the x-direction is 2(1 + q), the spatial frequency cut-off in the y-direction is √ 3(1 + q).The PSF resulting from the Lukosz bound does not differ much from the incoherent PSF for a spatial frequency cut-off 2(1 + q) with a FWHM of 0.28 (in units λ /NA) in both the x and y-direction.The first dark diffraction ring has a hexagonal imprint with 6 'pits', where the PSF drops slightly below zero (minimum value −1 × 10 −2 ).This small violation of positive definiteness is not expected to give rise to highly visible artefacts.However, the hexagonal imprint may be reduced by using an anisotropically stretched version of the Lukosz bound: where μ is a stretching parameter.The anisotropically stretched Lukosz bound OTF gives rise to a non-negative and isotropic PSF (at least up to the second dark diffraction ring) for a value μ = 0.0656.In the following we will use the non-stretched version of the Lukosz bound for the sake of simplicity.

Lukosz bound filtering
There are two ways in which the Lukosz bound can be used in apodization.The first way is simply to apply it as a post-processing apodization filter Λ (v) after a super-resolution image has been generated with the generalized filtering approach.This results in an OTF: The second way to use the Lukosz bound is to use modified versions of the functions that characterize the ideal or desired response of the imaging system.If we use Blm (v) = s m w m Ĥ (v + R l • q m ) / Λ (v) then the OTF becomes: Clearly, this alternative route of incorporating the Lukosz bound gives rise to an OTF that is very similar to the one obtained by applying the Lukosz bound as apodization filter after the filtering operation.The difference is that now the regularization term is effectively damped with a factor Λ (v) 2 resulting in a an effectively smaller regularization effect at higher spatial frequencies.In the limit κ → 0 the effective OTF approaches the Lukosz bound, just as for the generalized Wiener filtering with Lukosz apodization.This integrated Lukosz filtering approach gives rise to an effective OTF that satisfies the Lukosz bound for all values of the regularization parameter, with the exception of the signal energy regularization case (p = 0).See the Appendix for more details about ensuring the Lukosz bound is not violated in this case.

Simulation results
The reconstruction methods presented in the previous sections of this paper are demonstrated here using simulations.We simulated a virtual resolution chart consisting of 512×512 pixels with pixel size 1/16 (in normalized units of λ /NA).The virtual resolution chart, shown in Fig. 3(a) contains various structures (points, crossing lines, bar patterns, uniform blocks) that can be used to judge image quality.Several groups of bar patterns can be identified, the number indicates the pitch in units of λ /16NA.Pattern "4" is thus at the cut-off for widefield imaging, pattern "2" near the cut-off for SIM.The final images were binned 2×2 giving a 256×256 image with pixel size corresponding to Nyquist sampling at the maximum extended cut-off spatial frequency of SIM (equal to 4 in normalized units of NA/λ ).The simulated recorded images were corrupted by shot noise, which is taken to be the dominant source of noise.The reconstructions were performed with an implementation of Eqs.(17,20) in which either the Lukosz bound was used as apodization filter Eq. (29) or in which the Lukosz bound was incorporated in the Tikhonov-Miller functional Eq.(30).All simulations were performed in MATLAB (Math-Works, Natick, USA) using the DIPimage toolbox [32].
Fig. 3(a) shows an image of the simulated resolution chart, Fig. 3(b) of the simulated widefield image where on average the intensity corresponded to 800 photons per (λ /8NA) 2 area.Key SIM image reconstruction artefacts are shown in Fig. 3(c) (Media 1), which shows the effect of the side band height parameter s on the results obtained by the weighted sum of shifted bands reconstruction, and Fig. 3(d) (Media 2), which shows the effect of the regularization parameter κ on the results obtained by the generalized Wiener filtering in the image reconstruction without any apodization.The simple and straightforward weighted sum of shifted bands reconstruction shows resolution improvement at the expense of noise enhancement for increasing values of the side band height parameter.The optimum value for s seems to be in the range 1.0-1.5.Generalized Wiener filtering without apodization results in grave image artefacts.Significant effects of negative pixel values, edge ringing and an enhanced noise structure overlaying the genuine image can be observed, especially for low values of the regularization parameter.
Figure 4  aforementioned image artefacts.Notice that residual noise structures are still clearly visible in the regime of small regularization parameter values, in particular in the non-sparse parts of the image.
Additional simulations were performed in order to quantify the sharpness and the noise enhancement effect.The signal-to-noise ratio (SNR) of uniform objects, as well as the full-width at half-maximum (FWHM) and the peak signal-to-noise ratio (PSNR) of point objects, have been computed as a function of the regularization parameter κ, for different values of the side band height parameter (s = 1, s = 2, and s = 3) for both signal energy regularization (p = 0) and signal gradient energy (p = 1) regularization for generalized Wiener filtering with Lukosz bound apodization and for integrated Lukosz filtering.The SNR was measured over the large uniform rectangle in the top left of resolution chart and the PSNR was computed over 64 reconstructed peak signals of point sources for each case.
Figure 5 shows the resulting graphs, where the SNR and PSNR are normalized to the values obtained for widefield imaging with the same number of collected photons.The FWHM for low values of the regularization parameter is about 0.28λ /NA, which is a factor of about 1.8 smaller than the widefield FWHM-value 0.51λ /NA.The SNR in that regime of regularization parameters is lower than the widefield value, implying that the resolution improvement is obtained at the expense of noise enhancement.The PSNR starts out at a higher value than the widefield case because the narrower spot already implies an increase in peak signal.The FWHM and both the SNR and PSNR increase with increasing regularization parameter, indicating a trade-off between resolution improvement and noise suppression.The exception is the case of integrated Lukosz-filtering with p = 0 regularization, for which violations of the Lukosz bound may occur for side band height parameters larger than approximately 1.6, as discussed in teh appendix.The same trade-off is apparent from the effect of the side band height parameter s.The FWHM and both the SNR and PSNR increase with decreasing s.Signal gradient energy regularization (p = 1) is less sensitive to variations in the side band height parameter and shows flatter curves  in the regime of low regularization parameters than signal energy regularization (p = 0), and is therefore the more favourable regularization from the point of view of robustness.Generalized Wiener filtering with Lukosz apodization is to be prefered over integrated Lukosz filtering for p = 0 regularization due to the Lukosz bound violation for large side band height parameters for the integrated Lukosz filtering method.The differences between the two filtering approaches for p = 1 regularization are much smaller.A small shift in the regime of regularization parameter values where the resolutionnoise-suppression trade-off occurs can be recognized.In addition, integrated Lukosz filtering gives rise to a slightly lower SNR and PSNR curves, due to an effective lower regularization at high spatial frequencies (the regularization term is effectively weighted by the square of the Lukosz bound OTF).This drawback is compensated by a more natural noise spectrum, i.e. the small additional amount of high frequency noise visually masks the low frequency noise structure (see also Fig. 4).
Different apodization functions than the Lukosz bound have been in use so far.A tri-angular apodization function [13], for example, results in a bit sharper image but also in more visible noise artefacts, as it does not comply with the Lukosz bound close to the cut-off.An alternative is the approach based on the distance transform aplied to the support of the OTF [31].Here the apodization has the form d (v) ζ , where d (v) is the (normalized) Euclidean distance from the point v to the cutoff of the SIM-OTF and where ζ is an exponent.This approach may be expected to give similar results as the Lukosz bound apodization provided the exponent ζ is a bit higher than one.
The data misfit term of the TM-functional is often interpreted as the logarithm of the likelihood for the reconstructed ground truth image given the measured image data corrupted by noise.A least squares type of data misfit term then corresponds to uniform Gaussian noise.Shot noise, however, is governed by the Poisson distribution and should therefore give rise to a data misfit term in the TM-functional that is different from least squares.It also raises the question if the root cause of the observed noise structure lies in the mismatch between the nature of the noise source and the character of the data misfit term in the TM-functional.To that end we have performed additional simulations where uniform Gaussian noise, e.g.camera readout noise, is added to the raw images.It turns out that the same noise structures are observed as for the shot noise simulations, indicating that the precise nature of the noise source is not particularly relevant for the structured noise artefact.
Finally, we note that the post-filtering technique, in which a filtering operation is applied to the weighted sum of shifted bands reconstruction gives rise to an amount of noise enhancement intermediate between the non-filtered weighted sum of shifted bands reconstruction and the full-blown generalized filtering approach, both for the generalized Wiener-filtering approach followed by Lukosz bound apodization and for the integrated Lukosz bound filtering approach.

Experiments
The effects of regularization and apodization have also been investigated using experimental data.We have imaged the chromatin distribution in HDLM-2, which is a Hodgkins lymphoma cell line [33], using a Zeiss Elyra SIM microscope.The cells were fixed on the microscope slide in formaldehyde and the chromatin was counterstained with DAPI (4',6-diamidino-2phenylindole).An excitation wavelength of 405 nm was used, the passband of the filtercube was 420-480 nm for the emitted light.A 63X objective lens with numerical aparture of 1.40 was used with immersion oil with refractive index of 1.518.The images were captured on an Andor EM-CCD iXon 885 camera with a pixel size of 8 μm.A tube lens with additional zoom factor 1.6 was used giving a projected pixel size of 79 nm in the object plane and a projected pixel size in the reconstructed SIM images of 43 nm.The projected pixel size of the captured images corresponds to Nyquist sampling for an emission wavelength of 444 nm.The reciprocal pitch of the illumination pattern was calculated from the data using the methods of refs.[31,34] to be 0.74 times the cut-off frequency of the lens.We have performed additional experiments on actin filaments in bovine pulmonary artery endothelial (BPAE) cells.The Factin in these cells was stained with Alexa488, the slide was a Molecular Probes prepared slide (F36924).Imaging and reconstruction settings were the same as for the HDLM-2 cells, except the excitation wavelength was 488 nm and the passband of the filtercube was 495-575 nm.Image reconstruction of a single 2D-slice was done by first applying the phase shift estimation methods developed by Wicker and Heintzmann [31,34] for extracting the spatial frequency bands and then combining and filtering the bands using the Tikhonov-Miller filtering approach with the currently proposed Lukosz bound as apodization function.In the phase shift estimation only the 0th and ±2nd order bands were extracted so as to emulate a true 2D-SIM data acquisition.Fig. 6 and Fig. 7 shows the resulting reconstructions for the two experiments as a function of the regularization parameter for the different filtering and regularization methods also shown in Fig. 4. The other reconstruction parameters were the data misfit weighting parameters α = 0.95 and σ = 1, which is at 50% of the lateral cut-off frequency of the lens, and side band height parameter s = 3.The filters with this Gaussian data misfit term are described in Eq. (24).The relatively high value for the side band height is used to compensate for the loss of high-frequency image content due to the data misfit weighting.The experimental data indicates that Lukosz bound apodization with p = 0 (signal energy) regularization seems stepping jitter.Finally, different ways of regularization, better suited to suppressing typical SIM artefacts, may be explored.

Appendix
This integrated Lukosz filtering approach gives rise to an effective OTF that satisfies the Lukosz bound for all values of the regularization parameter, with the exception of the signal energy regularization case (p = 0).The normalization of the OTF (must be equal to one for v = 0) implies that the r.h.s. of Eq. ( 30) must be multiplied with a normalization factor: Now the effective OTF is no longer necessarily a decreasing function of the regularization parameter κ.It follows that violations of the Lukosz bound can then occur, as the effective OTF is equal to the Lukosz bound for κ = 0.The requirement that the effective OTF must satisfy the Lukosz bound leads to the constraint that the function: must be smaller than its value for v = 0 for all v.It may be expected that, whenever this constraint is not satisfied, this occurs for spatial frequencies at the side band peaks, as the function z (v) has local maxima there.For these spatial frequencies it may be deduced that the constraint is satisfied provided the side band height parameter is sufficiently small: with H s and Λ s equal to the widefield OTF and to the Lukosz bound at the side-peak spatial frequency, respectively.For a side band peak at 90% of the wide field cut-off (H s = 0.037 and Λ s = 0.45) and a maximum modulation depth of the illumination (w = 1/2) we find s max = 1.56.

Fig. 1 .
Fig. 1.(a) Plot of the Lukosz bound OTF and incoherent OTF in 1D.(b) Plot of the Lukosz bound PSF and incoherent PSF in 1D.(c) Plot of the Lukosz bound OTF and incoherent OTF in 2D.(d) Plot of the Lukosz bound PSF and incoherent PSF in 2D.

#Fig. 2 .
Fig. 2. (a) Plot of the 2D SIM Lukosz bound OTF.(b) Plot of the 2D SIM Lukosz bound PSF.(c) Cross-sections of the 2D SIM Lukosz bound OTF along the x and y-direction, and the incoherent OTF for cutoff 2 + 2q as reference.(d) Cross-sections of the 2D SIM Lukosz bound PSF along the x and y-direction, and the incoherent PSF for cut-off 2 + 2q as reference.

Fig. 3 .
Fig. 3. (a) Digital resolution chart used in the simulations.The numbers with the bar patterns indicate the pitch in pixel units, i.e. in units of λ /16NA.(b) Simulated widefield image, (c) Screenshot of a series of SIM image reconstructions with different side band height parameters using a weighted sum of bands method (Media 1).The screenshot is the reconstruction with high side band height parameter s = 3 showing noise amplification.(d) Screenshot of a series of SIM image reconstructions with generalized Wiener filtering without apodization for different regularization parameter values (Media 2).Shown is the reconstruction at low regularization parameter κ = 10 −4 , which demonstrates halo artefacts and low-frequency noise structures.

Fig. 4 .
Fig. 4. (a) Screenshot of a series of SIM image reconstructions with generalized Wiener filtering with signal energy (p = 0) regularization and Lukosz bound apodization (Media 3).Shown is the reconstruction for κ = 5.62 × 10 −4 .(b) Screenshot of a series of SIM image reconstructions with generalized Wiener filtering with signal gradient energy (p = 1) regularization and Lukosz bound apodization (Media 4).Shown is the reconstruction for κ = 10 −3 .(c) Screenshot of a series of SIM image reconstructions with integrated Lukosz filtering with signal energy (p = 0) regularization (Media 5).Shown is the reconstruction for κ = 10 −1 .(d) Screenshot of a series of SIM image reconstructions with integrated Lukosz filtering with signal gradient energy (p = 1) regularization and Lukosz bound apodization (Media 6).Shown is the reconstruction for κ = 1.78×10 −1 .The regularization parameters for the different screenshots have been manually optimized.

1 Fig. 5 .
Fig. 5. Simulated image quality measures as a function of regularization parameter for different side band heights.The columns show the FWHM, the SNR of uniform image patches, and the peak SNR of single spots for generalized Wiener filtering with Lukosz bound apodization with p = 0 (signal energy) regularization ((a)-(c)), for Lukosz apodization with p = 1 (signal gradient energy) regularization ((d)-(f)), integrated Lukosz filtering with p = 0 regularization ((g)-(i)), and for integrated Lukosz filtering with p = 1 regularization ((j)-(l)).The SNR and PSNR are normalized by the values for the corresponding widefield images.

Fig. 6 .
Fig. 6.Reconstructions from experimental data of the chromatin distribution in a HDLM-2 cell.(a) Screenshot of a series of reconstructions with generalized Wiener filtering with p = 0 (signal energy) regularization and Lukosz bound apodization (Media 7).Shown is the reconstruction for κ = 3.16.(b) Screenshot of a series of SIM image reconstructions with generalized Wiener filtering with p = 1 (signal gradient energy) regularization and Lukosz bound apodization (Media 8).Shown is the reconstruction for κ = 1.78 × 10 −1 .(c) Screenshot of a series of SIM image reconstructions with integrated Lukosz filtering with p = 0 (signal energy) regularization (Media 9).Shown is the reconstruction for κ = 3.16 × 10 −2 .(d) Screenshot of a SIM image reconstruction with integrated Lukosz filtering with p = 1 (signal gradient energy) regularization and Lukosz bound apodization (Media 10).Shown is the reconstruction for κ = 5.62 × 10 −3 .The regularization parameters for the different screenshots have been manually optimized.The intensities are in arbitrary units, the scale bars are 5 μm.

#Fig. 7 .
Fig. 7. Reconstructions from experimental data of actin in BPAE cells.(a) Screenshot of a series of reconstructions with generalized Wiener filtering with p = 0 (signal energy) regularization and Lukosz bound apodization (Media 11).Shown is the reconstruction for κ = 1 × 10 −1 .(b) Screenshot of a series of SIM image reconstructions with generalized Wiener filtering with p = 1 (signal gradient energy) regularization and Lukosz bound apodization (Media 12).Shown is the reconstruction for κ = 1.78 × 10 −1 .(c) Screenshot of a series of SIM image reconstructions with integrated Lukosz filtering with p = 0 (signal energy) regularization (Media 13).Shown is the reconstruction for κ = 3.16 × 10 −2 .(d) Screenshot of a SIM image reconstruction with integrated Lukosz filtering with p = 1 (signal gradient energy) regularization and Lukosz bound apodization (Media 14).Shown is the reconstruction for κ = 5.62 × 10 −1 .The regularization parameters for the different screenshots have been manually optimized.The intensities are in arbitrary units, the scale bars are 5 μm.The cutouts have been linearly stretched for easier visual inspection.