Three-dimensional structured illumination microscopy using Lukosz bound apodization reduces pixel negativity at no resolution cost

Christiaan H. Righolt; Sabine Mai; Lucas J. van Vliet; Sjoerd Stallinga

doi:10.1364/OE.22.011215

1. Introduction

Structured illumination microscopy (SIM) extends the resolution of fluorescence microscopy with a factor of up to two beyond the diffraction limit. It was first described theoretically in the 1960s [1] and demonstrated as both an optical sectioning and a superresolution method in the 1990s and 2000s [2–7]. Recent years have seen various improvements in the technology [8–15]. The standard way to reconstruct a SIM image has been through the inversion of the linear combination of frequency bands, followed by a Wiener style multi-band filter. In this Wiener style filter, the inverse signal-to-noise ratio has been replaced by a single regularization parameter. An apodization filter is then used as a final post-processing step, in order to suppress negative pixel value and ringing artifacts [11, 16–18]. It turns out, however, that apodization filters used in practice often suffer from severe negative pixel value artifacts, which are typically solved by clipping the negative pixel values to zero. This practice leads to unreliable final images for laymen users, and ultimately to biased biological results. Recently we described a method for SIM reconstruction in 2D [19] with the aim to reduce the number of free parameters and to ensure a non-negative Point Spread Function (PSF) by using the upper bound to the Modulation Transfer Function (MTF) as derived by Lukosz [20, 21] as apodization function.

In this paper we focus on the extension of the Lukosz bound apodization method to 3D-SIM. This method cannot be trivially extended from 2D to 3D due to the severe anistropy in the MTF of the microscope. In particular we will introduce, for the first time, the (numerical) parameterization of the anisotropic 3D cut-off of the SIM OTF, the generalization of the Lukosz bound to 3 (or any higher number of) dimensions, a method to compute the Lukosz bound for any arbitrary frequency cut-off surface with two possible ways to deal with anisotropic transfer functions, and the unification of all reconstruction steps in a generalized Tikhonov-Miller filter. This transforms SIM reconstruction from a patchwork of individual, unconnected steps to a unified, tunable optimization problem for quantitatively reliable super-resolution microscopy.

The content of the paper is as follows. In section 2 we will shortly summarize the main reconstruction steps of [19], derive the parametrization of the 3D SIM cut-off, give an extension of the Lukosz bound to 3D, and compare 3D Lukosz bound and triangular apodization. Experimental validation will be presented in section 3, were the methods developed throughout the paper will be applied to images of actin filaments and DAPI-stained DNA structure in a cell. The paper will conclude with a discussion in section 4.

2. Theory

2.1. Overview of 3D SIM reconstruction

The spatial coordinates (x, y, z)^T, with coordinate z along the optical axis, are normalized to u = (xn_med sinα/λ, yn_med sinα/λ, z4n_med sin² (α/2)/λ))^T, with n_med the refractive index of the medium in object space, λ the wavelength of the emitted light and α the marginal ray angle in object space. The spatial frequencies are denoted by v, for which we use the Fourier transform convention:

\hat{H} (v) = \int H (u) exp (- 2 π i u \cdot v) d u,

in which H (u) is the PSF and Ĥ (v) the optical transfer function (OTF). The illumination function W (u) is periodic and its spectrum consists of delta-peaks:

\hat{W} (v) = \sum_{m} w_{m} δ (v - q_{m}),

where w_m denotes the relative weight of the m^th frequency component. In 3D SIM, the zeroth and first orders of a diffraction grating are projected on the back focal plane of the objective. The interference pattern of the three resulting plane waves in object space leads to a periodic illumination pattern with a spectrum that has seven peaks, q_m ∈ (0, 0, 0)^T, (±2q_r, 0, 0)^T, (±q_r, 0, ±q_z)^T, where q_r = 1/p and

q_{z} = \frac{1}{4 {sin}^{2} (α / 2)} (1 - \sqrt{1 - \frac{{sin}^{2} α}{p^{2}}})

with p the normalized grating period. Images are recorded for a set of rotations R_l, (l = 1, 2,...N_r) and translations u_n = n(p/N_t, 0) (n = 1, 2,... N_t). For 3D SIM typically N_r ∈ {3, 5} rotations and N_t = 5 translations are used. In the remainder of the paper we will use N_r = 3.

The images are linear combinations of shifted band-passed images. These band-passed images can be recovered by an inversion algorithm [22–24]:

{\hat{I}}_{bands, l m} (v) = w_{m} \hat{H} (v) \hat{T} (v - R_{l} \cdot q_{m}),

where the object spectrum is denoted by T̂ (v). A final reconstructed image can now be obtained by moving each band in frequency space and filtering it with kernel F̂_lm (v):

{\hat{I}}_{gen} (v) = \sum_{l m} s_{m} {\hat{F}}_{l m} (v) {\hat{I}}_{band, l m} (v + R_{l} \cdot q_{m}),

where s_m is an additional weight factor. The generalized Tikhonov-Miller functional to be minimized in SIM reconstruction [19] is given by:

ℒ_{gen} = \int [\sum_{l m} \hat{g} (v_{l m}) {| {\hat{B}}_{l m} (v) {\hat{I}}_{gen} (v) - s_{m} \hat{C} (v) {\hat{I}}_{band, l m} (v_{l m}) |}^{2} + κ {| \hat{A} (v) {\hat{I}}_{gen} (v) |}^{2}] d v .

In Eq. (6), v_lm = v + R_l · q_m is the band-shifted spatial frequency vector, B̂_lm (v) = s_mw_mĤ (v_lm) denotes the ideal or desired response of the imaging system, κ the regularization parameter, Â (v) the regularization weight function, Ĉ (v) the apodization function, and ĝ (v) a data misfit weight function, that can be used to suppress the hexagonal imprint artifact that is typical for SIM [23]. A Gaussian data misfit function can be chosen:

\hat{g} (v) = 1 - α exp [- \frac{v_{x}^{2} + v_{y}^{2}}{2 σ_{l}^{2}} - \frac{v_{z}^{2}}{2 σ_{a}^{2}}],

where α is the amplitude, σ_l the lateral frequency width and σ_a the axial frequency width. The functional is minimized with respect to the reconstructed image Î_gen (v), leading to the filter kernels:

{\hat{F}}_{l m} (v) = \frac{\hat{C} (v) \hat{g} (v + R_{l} \cdot q_{m}) {\hat{B}}_{l m} {(v)}^{*}}{κ \hat{A} {(v)}^{2} + \sum_{l m} \hat{g} (v + R_{l} \cdot q_{m}) {| {\hat{B}}_{l m} (v) |}^{2}}

and the effective OTF for the reconstructed image:

{\hat{H}}_{gen} (v) = \frac{\hat{C} (v) \sum_{l m} s_{m} w_{m} \hat{g} (v + R_{l} \cdot q_{m}) {\hat{B}}_{l m} {(v)}^{*} \hat{H} (v + R_{l} \cdot q_{m})}{κ \hat{A} {(v)}^{2} + \sum_{l m} \hat{g} (v + R_{l} \cdot q_{m}) {| {\hat{B}}_{l m} (v) |}^{2}} .

2.2. The 3D SIM cut-off

In order to compute the Lukosz-bound for 3D SIM, the anisotropic SIM cut-off needs to be derived first. This cut-off will be derived in spherical coordinates (ρ, θ, ϕ), such that v_x = ρ cosϕ sinθ, v_y = ρ sinϕ sinθ and v_z = ρ cosθ. The SIM OTF of Eq. (9) is non-zero when ∑_lm Ĥ (v + R_l · q_m) is non-zero. The cut-off of the SIM OTF is, therefore, given by the outer boundary segments of the incoherent OTFs moved to positions {R_l · q_m}. The cut-off of the incoherent OTF has been derived by Sheppard [25]. The cross-section of the incoherent OTF with a meridional plane has a two-lobed shape, each lobe consisting of two elliptical arcs. The SIM OTF cut-off ρ(θ, ϕ), is thus given by the intersection of the radial line in direction (θ, ϕ) and one of the moved elliptical arcs of the incoherent OTF. Using the translation q_m = (q_m,x, q_m,y, q_m,z), see Eq. (3), and using our coordinate normalization, we then find:

{sin}^{2} α {(\sqrt{{(v_{x} - q_{m, x})}^{2} + {(v_{y} - q_{m, y})}^{2}} - 1)}^{2} + 16 {sin}^{4} (\frac{α}{2}) {(v_{z} - q_{m, z} - \frac{cos α}{4 {sin}^{2} (α / 2)})}^{2} = 1 .

Substituting the lateral coordinates in Eq. (10) with v_x = v_z cosϕ tanθ and v_y = v_z sinϕ tanθ and squaring this equation leads to a quartic equation in v_z after several rearrangements:

a_{4} {v_{z}}^{4} + a_{3} {v_{z}}^{3} + a_{2} {v_{z}}^{2} + a_{1} v_{z} + a_{0} = 0 .

The polynomial coefficients are:

a_{0} = Q_{C}^{2} - {4 q_{m, x}}^{2} - {4 q_{m, y}}^{2},

a_{1} = 4 Q_{B} Q_{C} + 8 q_{m, x} cos ϕ tan θ + 8 q_{m, y} sin ϕ tan θ,

a_{2} = - 2 Q_{A} Q_{C} + 4 Q_{B}^{2} - 4 {tan}^{2} θ,

a_{3} = - 4 Q_{A} Q_{B},

a_{4} = Q_{A}^{2},

with the constants:

Q_{A} = 4 {tan}^{2} (\frac{α}{2}) + {tan}^{2} θ,

Q_{B} = q_{m, x} cos ϕ tan θ + q_{m, y} sin ϕ tan θ + (1 + 4 q_{m, z}) {tan}^{2} (\frac{α}{2}) - 1,

Q_{C} = \frac{1}{{sin}^{2} α} - {(\frac{1}{tan α} - 2 q_{m, z} tan (\frac{α}{2}))}^{2} - {q_{m, x}}^{2} - {q_{m, y}}^{2} - 1 .

Equation (11) has four solutions v_z,m,j (θ, ϕ) for each peak in the illumination pattern. The z-solutions can be transformed to spherical coordinates leading to:

ρ_{m, j} (θ, ϕ) = v_{z, m, j} (θ, ϕ) \sqrt{1 + {tan}^{2} θ} .

The jth solution represents an actual intersection when it is positive. Negative solutions match positive solutions for a different m because of symmetry in the system, and imaginary solutions occur when there is no intersection. The maximum of all solutions, excluding imaginary and negative numbers yields the SIM cut-off:

ρ (θ, ϕ) = max_{m, j} {ρ_{m, j} (θ, ϕ) > 0} .

The above derivation is valid, except when θ ∈ {0, π/2, π}. For θ = π/2, v_z = 0; for standard 3D SIM this leads to:

ρ (\frac{π}{2}, ϕ) = 2 q_{r} cos ϕ_{r} + 2 \sqrt{1 - q_{r}^{2} {sin}^{2} ϕ_{r}} .

For θ ∈ {0, π}, v_x = v_y = 0 and Eq. (10) can be solved directly as

ρ (0, ϕ) = ρ (π, ϕ) = q_{z} + \frac{\sqrt{1 - {(q_{r} - 1)}^{2} {sin}^{2} α} - cos α}{4 {sin}^{2} (α / 2)} .

The resulting parametrization of the SIM cut-off is shown in Fig. 1.

Fig. 1 Plot of the frequency cut-off in normalized frequency coordinates for standard 3D SIM in which the pattern is created by the zeroth and first order diffraction orders with three equally spaced rotation angles for q_r = 0.7 and q_z = 0.2.

Download Full Size | PDF

2.3. The 3D SIM Lukosz-bound

Lukosz presented several conditions a transfer function should satisfy for its corresponding PSF to be non-negative in 1D and 2D [20, 21]. Note that these conditions are necessary but not sufficient. Previously we discussed how this set of rules can be used to derive the Lukosz-bound for 2D SIM in order to get a non-negative SIM PSF [19].

The 1D bound is given by connecting the tips of Lukosz’ “stair case”:

{\hat{Λ}}_{1} (v, q_{c}) = cos (\frac{π | v |}{| v | + q_{c}}) for | v | \leq q_{c} and {\hat{Λ}}_{1} (v) = 0 elsewhere,

where q_c is the cut-off frequency. For the 2D case, the product of two 1D bounds along two orthogonal directions has to be taken. The actual bound is the minimum over all possible orientations of the 2D coordinate system. This formalism can be extended into 3D:

{\hat{Λ}}_{3} (v) = min {\prod_{j = 1}^{3} {\hat{Λ}}_{1} (v_{j}^{'}, q_{c, j} (ℛ)) | ℛ \in 𝒢},

where the minimum is taken over a set 𝒢 of rotations ℛ of the 3D coordinate frame, where v′ are the spatial frequencies in the rotated frame and q_c,j (ℛ) are the cut-offs along the rotated cardinal directions. The axial and lateral directions are fundamentally different, which requires care in implementing this recipe. We have devised two implementations that work well in practice. In the first, only rotations around the optical axis are considered; in the second, the full group of 3D rotations, SO(3), is considered in combination with an anisotropic pre-strechting of the spatial frequency space.

The first option of dealing with the anisotropy between the lateral and axial directions is to ensure that one of the cardinal directions remains aligned with the optical axis. This leaves one degree of freedom, the rotation ϕ around the v_z-axis. The coordinate system is now described by orthogonal axes (cosϕ, sinϕ, 0), (−sinϕ, cosϕ, 0) and (0, 0, 1). A spatial frequency vector v = (vsinξ cosψ, vsinξ sinψ, vcosξ) will have coordinates v′ = (vsinξ cos(ψ − ϕ), vsinξ sin(ψ − ϕ), cosξ) in the rotated frame. The minimum for all rotations for this bound is then given by:

{\hat{Λ}}_{3, ϕ} (v) = min {{\hat{Λ}}_{1} (v_{x}^{'}, ρ (\frac{π}{2}, ϕ)) {\hat{Λ}}_{1} (v_{y}^{'}, ρ (\frac{π}{2}, ϕ + \frac{π}{2})) {\hat{Λ}}_{1} (v_{z}^{'}, ρ (0, ϕ)) | ϕ \in [0, 2 π)},

where the respective cut-offs are described by Eqs. (22) and (23). The 3D bound should still adhere to the 1D bound in the direction (ξ, ψ). The 3D Lukosz-bound for rotations around the optical axis, Λ̂_OAR (v), therefore becomes:

{\hat{Λ}}_{OAR} (v) = min {{\hat{Λ}}_{3, ϕ} (v), {\hat{Λ}}_{1} (v, ρ (ξ, ψ))} .

Note that Λ̂_OAR (v_x, v_y, 0) is the 2D bound we described previously in [19].

The second option is to use all possible 3D rotations:

{\hat{Λ}}_{F R} (v) = min {\prod_{j = 1}^{3} {\hat{Λ}}_{1} (v_{j}^{'}, q_{c, j} (ℛ)) | ℛ \in SO (3)} .

Straightforward application of this recipe leads to a bound that is too pessimistic (too low) because of the severe anisotropy of the underlying spatial frequency cut-off. It appears that this can largely be solved by first applying an anisotropic stretch such that the bounding box of the cut-off surface has a cubic shape. To achieve this we rescale the v_z-axis by a factor:

ζ = \frac{q_{max, r}}{q_{max, z}} = \frac{q_{r} + 2}{q_{z} + 1 / 2} .

Note that this scale factor can be changed accordingly when a different spatial normalization is chosen; 2 and 1/2 are the lateral and axial frequency cut-offs of the lens OTF for our normalization. This scaling leads to several changed coordinates, marked with a tilde:

{\tilde{v}}_{z} = ζ v_{z}, \tilde{v} = \sqrt{v_{x}^{2} + v_{y}^{2} + ζ^{2} v_{z}^{2}} and \tilde{θ} = arccos (\frac{ζ v_{z}}{\sqrt{v_{x}^{2} + v_{y}^{2} + ζ^{2} v_{z}^{2}}}) .

The cut-off can also be obtained in the rescaled coordinate system:

\tilde{ρ} (\tilde{θ}, ϕ) = ρ (θ, ϕ) \sqrt{{sin}^{2} θ + ζ^{2} {cos}^{2} θ},

All of the cardinal direction cut-offs q_c,j (ℛ) and spatial frequency components v′_j in the rotated frame can be calculated using a parametrization of the rotation matrices ℛ with the three Euler angles, and can then be plugged into Eq. (28). This then yields Λ̂_FR (v), the 3D Lukosz-bound based on the full set of SO(3) rotations.

2.4. Comparison of apodization functions

The two options presented for the 3D SIM Lukosz-bound can be used as apodization function. A triangular apodization is typically used in practice, which is commonly created using a distance transform [11, 23]. In this case the apodization function becomes:

{\hat{C}}_{D T} (v) = \frac{𝒟 {\sum_{l m} \hat{H} (v + R_{l} \cdot q_{m}) > 0}}{max 𝒟 {\sum_{l m} \hat{H} (v + R_{l} \cdot q_{m}) > 0}},

in which 𝒟 denotes the distance transform and ∑_lm Ĥ (v + R_l · q_m) > 0 is the footprint of the SIM OTF. The triangle can alternatively be formed over every line originating in the origin to the cut-off in that direction. This leads to an apodization function:

{\hat{C}}_{line} (v) = 1 - \frac{v}{ρ (θ, ϕ)} for v \leq ρ (θ, ϕ) and {\hat{C}}_{line} (v) = 0 elsewhere .

We calculated these four apodization functions, the two triangle-based functions and the two Lukosz-bound functions, in the frequency domain for a 256 × 256 × 128 grid using Matlab (MathWorks, Natick, USA). We chose the pixel size such that the footprint of the SIM OTF fitted exactly. We used q_r = 0.65 and q_z = 0.16 for three equidistant rotation angles. The corresponding PSF’s, the inverse Fourier transforms of the respective functions were computed using the chirp z-transform algorithm [26]. The relative minimum values of these PSF’s are summarized in Table 1 together with the full width at half maximum (FWHM) in the cardinal directions. Lines through the PSF’s are plotted for the distance transform and SO(3)-rotated Lukosz-bound in Fig. 2.

Table 1. Comparison of the PSF’s corresponding to the four different apodization options for 3D-SIM with three rotation angles. In addition to the absolute values of the measures, a normalized measure with respect to the value for C_DT (u) is included as well. The relative minimum value of the PSF (compared to the maximum) is listed as min(PSF)/max(PSF). The lateral and axial full width at half maximum (FWHM lat. and FWHM ax.) are tabulated in the last four columns. The lateral anisotropy leads to a lateral FWHM difference between the R_l and (R_l + R_l₊₁)/2 directions of a half percent or less in all four cases. Note that the lateral and axial coordinates are normalized differently.

View Table

Fig. 2 Cross-sectional plots of the PSF’s corresponding to C_DT (u) in blue and Λ_FR (u) in red for 3D-SIM with three rotation angles. The PSF’s are the inverse Fourier transforms of Ĉ_DT (v) and Λ̂_FR (v). The left plot depicts the line through the origin in the lateral v_z = 0 plane in one of the pattern directions of R_l. The middle plot shows a similar line in the v_z = 0 plane in the direction midway between the directions of R_l and R_l₊₁. The right plot shows the values of the PSF’s along the optical axis.

Download Full Size | PDF

The traditional apodization filter, C_DT (u), yields negative values of −9% of the value at the origin. Negative PSF values lead to ringing artifacts in the reconstructed SIM image. The minimum occurs around u_x ≈ u_y ≈ 0.25, which causes these artifacts to be at the length scale of interest for this type of microscopy. Using a different method to create the triangle, namely by parametrizing the SIM cut-off, and drawing a triangle in each direction leads to a twofold reduction of the negativity. An order of magnitude improvement can, however, be gained by switching to the 3D SIM Lukosz-bound with the full SO(3) rotation group. This comes at the price of a small (around ≈ 10%) increase in the lateral FWHM. The axial FWHM is hardly affected. Note that, as for the 2D case, these 3D Lukosz-bounds have a small residual amount of negativity, because the conditions on which they are based are necessary conditions, not sufficient conditions.

All performance measures are a bit worse for the Lukosz-bound implementation with only rotations around the optical axis. The computational load of finding the Lukosz-bound for only one rotational degree of freedom is, however, significantly lower. The calculation of Λ̂_FR (v) is computationally intensive because it needs to be performed in a 6D-space (three spatial frequency components, three Euler angles parametrizing all rotations). Fortunately, this function does not change from image to image when the imaging system parameters, in particular the pitches and rotation angles of the illumination pattern, remain constant.

3. Experiments

In order to investigate the performance of the 3D-Lukosz bound as apodization function we recorded raw SIM image series of two test samples. The first sample is of tubulin from bovine pulmonary artery endothelial (BPEA) cells on a Molecular Probes test slide (F14781), the second sample contains the DAPI (4,6-diamidino-2-phenylindole)-stained DNA structure of a fixed, complex, multinucleated Reed-Sternberg cell from the Hodgkin’s lymphoma HDLM-2 cell line [27]. Images were acquired using a Zeiss Elyra SIM microscope with a 63x/1.40 NA Plan-Apochromat objective lens and n_med = 1.518 immersion oil equipped with an Andor EM-CCD iXon 885 camera with 8 μm pixel size. The back-projected pixel size was 79 nm, because a 1.6x tube lens was used as well. The tubulin image was acquired using 488 nm laser excitation, a projected grating period of 444 nm and an emission passband of 495–575 nm. The DAPI image was acquired using 405 nm laser excitation, a projected grating period of 365 nm and a 420–480 nm passband for the emitted light.

Image reconstruction of the 3D image was performed in two steps. First, we extracted the spatial frequency bands using the methods developed by Wicker and Heintzmann [22–24]. Second, we obtained the final image using the filters described in section 2.1 with the 3D Lukosz-bound based on SO(3) rotations as apodization function Ĉ (v). We used conventional total energy regularization, Â (v) = 1 and the same amplitude α = 0.95 for the Gaussian data misfit function as for the 2D reconstruction. We fit a Gaussian to the lens PSF, Fourier transformed this Gaussian and took half the widths in Fourier space as the widths for the data misfit function: σ_l = 0.9 and σ_a = 0.3. The images were reconstructed for a range of regularization parameters κ without boosting the side-bands, i.e. s_m = 1 for all m. A slice of the 3D volume depicting tubulin filaments is shown for κ = 1 × 10⁻⁵ in Fig. 3. A slice depicting the DNA structure in the multinucleated cell for κ = 1×10⁻³ is shown in Fig. 4. We choose these values because they visually led to the best trade-off between noise suppression and image sharpness and contrast. The influence of the regularization parameter on these slices is shown in Media 1 and Media 3 for both images. A z-stack movie, going through the slices along the optical axis is visible in Media 2 and Media 4 for both samples. Note that no clipping of negative pixel values has been applied to any of the reconstructions.

Fig. 3 3D SIM image of tubulin filaments from a BPEA cell sample reconstructed using the full 3D SIM Lukosz-bound as apodization function. The intensity of the reconstructed image is linearly stretched, without clipping, between gray values 0 and 255 for display purposes. Shown is the result for a lateral z-slice (A) and an axial y-slice (B) for regularization parameter κ = 10⁻⁵. The lateral slice is shown for a range of values for the regularization parameters in Media 1 and for all z-positions in Media 2. The scale bar is 5 μm.

Download Full Size | PDF

Fig. 4 3D SIM image of the DNA structure (DAPI) in a multinucleated Reed-Sternberg cell reconstructed using the 3D SIM Lukosz-bound as apodization function. The intensity of the reconstructed image is linearly stretched, without clipping, between gray values 0 and 255 for display purposes. Shown is the result for a lateral z-slice (A) and an axial y-slice (B) for regularization parameter κ = ×10⁻³. The lateral slice is shown for a range of values for the regularization parameters in Media 3 and for all z-positions in Media 4. The scale bar is 5 μm.

Download Full Size | PDF

We compared these images with the reconstruction obtained by Zeiss’ ZEN software (release version 8.0) using triangular apodization, also without clipping negative pixel values to zero. The regularization parameter was set to the same values, 10⁻⁵ and 10⁻³ respectively, as for Lukosz-bound apodization. The ratio of the minimum to the maximum pixel value and the percentage of negative pixels were determined for both reconstructions as follows. First, we estimated the zero-level in the ZEN reconstructions from the background, because the ZEN software maps the images to 16-bit files. The foreground object was subsequently found by a coarse segmentation – an isodata threshold followed by a binary dilation – for the same slices as shown in Figs. 3 and 4. We took the minimum as the 1st and the maximum as the 99th percentile over the region of interest (ROI) thus found in order to be independent of outliers. We also estimated the number of pixels in the ROI with value below the background level. The distribution of normalized pixel values over the ROI’s are plotted in Fig. 5.

Fig. 5 Probability density function (pdf) of the normalized pixel values in reconstructed 3D SIM images. For each image, the object pixels are normalized by linearly scaling the background level to 0 and maximum pixel value to 1. The blue lines denote the pdf’s, the bin-widths are 0.01, the red lines indicate the zero-intensity levels. This is done for: the ZEN reconstruction with triangular apodization of the tubulin sample (a); our reconstruction with Lukosz-bound apodization of the tubulin sample (b); the triangular apodized ZEN reconstruction of the DAPI cell sample (c); and our Lukosz-bound apodization in the reconstruction of the DAPI cell sample (d).

Download Full Size | PDF

For the tubulin sample, the ZEN reconstruction led to a minimum pixel value of −7.8% of the maximum and 20% of the pixels in the ROI had a negative value; our reconstruction (Fig. 3) led to a minimum pixel value of −1.8% of the maximum and 3.9% negative pixels in the ROI. For the DAPI image, the ZEN reconstruction led to a minimum pixel value of −46% of the maximum with 42% negative pixels in the ROI. Our reconstruction (Fig. 4) had a minimum value of −18% the maximum; 9.3% of the pixels were negative in the ROI. Note that the depths of the minima do not solely depend on the choice of apodization function. The data misfit function, regularization term, noise properties and frequency content of the sample play a role as well. The experiments confirm our simulation results and show that using the Lukosz-bound as apodization function substantially reduces the negativity in the reconstructed image. For our images this led to a two to fivefold reduction in negativity. The sharpness of the Lukosz-bound reconstruction is visually the same as the sharpness of the ZEN reconstruction (data not shown). The reduction of negativity does, however, improve image contrast when seen on a stretched linear scale.

4. Discussion

In summary, we presented an extension of the Lukosz-bound to 3D-SIM and incorporated this bound as apodization function in 3D-SIM image reconstruction. The reconstruction is formulated in the framework of Tikohonov-Miller filtering, an approach that unifies the linear combination of bands shifted in spatial frequency space with filtering steps for noise suppression, regularization and apodization. We presented a formalism to calculate the 3D Lukosz-bound for any given surface that describes the cut-off of a 3D OTF in spatial frequency space. An analytical description of this cut-off surface for 3D SIM, valid for unequally spaced rotation angles and rotation-dependent grating periods, was derived as well. We compared the 3D SIM Lukosz-bound to other choices for the apodization function that are currently used in SIM reconstructions. The Lukosz-bound results in a PSF with an order of magnitude reduction of negativity at the expense of a moderate increase in width when compared to the distance transform based triangle function. Switching to Lukosz-bound based apodization reduces ringing artifacts and avoids the need for clipping the final image as shown for the tubulin and DAPI results in Figs. 3 and 4.

Apodization with the 3D Lukosz-bound does not require fine-tuning of any free parameters in the apodization function, such as taking it to a power other than one, as is sometimes done with the triangular apodization function [23]. Note that raising the 3D Lukosz-bound to a power smaller than one would violate the non-negativity criteria. Raising it to a power larger than one causes an unnecessary reduction of the “effective” width in frequency space, leading to a wider PSF.

References and links

1. W. Lukosz, “Optical systems with resolving powers exceeding the classical limit,” J. Opt. Soc. Am. 56, 1463–1471 (1966). [CrossRef]

2. M. A. A. Neil, R. Juskaitis, and T. Wilson, “Method of obtaining optical sectioning by using structured light in a conventional microscope,” Opt. Lett. 22, 1905–1907 (1997). [CrossRef]

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

4. G. E. Cragg and P. T. So, “Lateral resolution enhancement with standing evanescent waves,” Opt. Lett. 25, 46–48 (2000). [CrossRef]

5. M. A. A. Neil, A. Squire, R. Juskaitis, P. I. H. Bastiaens, and T. Wilson, “Wide-field optically sectioning fluorescence microscopy with laser illumination,” J. Microsc. 197, 1–4 (2000). [CrossRef] [PubMed]

6. 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]

7. J. T. Frohn, H. F. Knapp, and A. Stemmer, “True optical resolution beyond the Rayleigh limit achieved by standing wave illumination,” Proc. Natl. Acad. Sci. U. S. A. 97, 7232–7236 (2000). [CrossRef] [PubMed]

8. R. Heintzmann, T. Jovin, and C. Cremer, “Saturated patterned excitation microscopy - a concept for optical resolution improvement,” J. Opt. Soc. Am. B 19, 1599–1609 (2002). [CrossRef]

9. M. G. L. Gustafsson, “Nonlinear structured-illumination microscopy: Wide-field fluorescence imaging with theoretically unlimited resolution,” Proc. Natl. Acad. Sci. U. S. A. 102, 13081–13086 (2005). [CrossRef] [PubMed]

10. R. Fiolka, M. Beck, and A. Stemmer, “Structured illumination in total internal reflection fluorescence microscopy using a spatial light modulator,” Opt. Lett. 33, 1629–1631 (2008). [CrossRef] [PubMed]

11. 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]

12. P. Kner, B. B. Chhun, E. R. Griffis, L. Winoto, and M. G. L. Gustafsson, “Super-resolution video microscopy of live cells by structured illumination,” Nat. Methods 6, 339–342 (2009). [CrossRef] [PubMed]

13. L. Shao, P. Kner, E. H. Rego, and M. G. L. Gustafsson, “Super-resolution 3d microscopy of live whole cells using structured illumination,” Nat. Methods 12, 1044–1046 (2011). [CrossRef]

14. L. Wang, M. C. Pitter, and M. G. Somekh, “Wide-field high-resolution structured illumination solid immersion fluorescence microscopy,” Opt. Lett. 36, 2794–2796 (2011). [CrossRef] [PubMed]

15. O. Mandula, M. Kielhorn, K. Wicker, G. Krampert, I. Kleppe, and R. Heintzmann, “Line scan - structured illumination microscopy super-resolution imaging in thick fluorescent samples,” Opt. Express 20, 24167–24174 (2012). [CrossRef] [PubMed]

16. C. Berenstein and E. Patrick, “Exact deconvolution for multiple convolution operators–An overview, plus performance characterizations for imaging sensors,” Proc. IEEE 78, 723–734 (1990). [CrossRef]

17. L. P. Yaroslavsky and H. J. Caulfield, “Deconvolution of multiple images of the same object,” Appl. Opt. 33, 2157–2162 (1994). [CrossRef] [PubMed]

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

19. C. H. Righolt, J. A. Slotman, I. T. Young, S. Mai, L. J. van Vliet, and S. Stallinga, “Image filtering in structured illumination microscopy using the lukosz bound,” Opt. Express 21, 24431–24451 (2013). [CrossRef] [PubMed]

20. W. Lukosz, “Übertragung Nicht-negativer Signale Durch Lineare Filter,” J. Mod. Opt. 9, 335–364 (1962).

21. W. Lukosz, “Properties of linear low-pass filters for nonnegative signals,” J. Opt. Soc. Am. 52, 827–829 (1962). [CrossRef]

22. K. Wicker, “Increasing resolution and light efficiency in fluorescence microscopy,” Ph.D. thesis, King’s College, London (2010).

23. 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]

24. K. Wicker, “Non-iterative determination of pattern phase in structured illumination microscopy using autocorrelations in Fourier space,” Opt. Express 21, 24692–24701 (2013). [CrossRef] [PubMed]

25. C. J. R. Sheppard, M. Gu, Y. Kawata, and S. Kawata, “Three-dimensional transfer functions for high-aperture systems,” J. Opt. Soc. Am. A 11, 593–598 (1994). [CrossRef]

26. J. L. Bakx, “Efficient computation of optical disk readout by use of the chirp z transform,” Appl. Opt. 41, 4879–4903 (2002). [CrossRef]

27. H. G. Drexler, G. Gaedicke, M. S. Lok, V. Diehl, and J. Minowada, “Hodgkin’s disease derived cell lines HDLM-2 and L-428: comparison of morphology, immunological and isoenzyme profiles,” Leukemia Res. 10, 487–500 (1986). [CrossRef]

	$\frac{min (PSF)}{max (PSF)}$	÷C_DT	FWHM lat.	÷C_DT	FWHM ax.	÷C_DT

C_DT (u)	−0.0905	1	0.286	1	1.44	1
C_line (u)	−0.0483	0.534	0.308	1.077	1.35	0.935
Λ_OAR (u)	−0.0158	0.175	0.321	1.120	1.53	1.062
Λ_FR (u)	−0.0109	0.120	0.320	1.118	1.44	0.997

Three-dimensional structured illumination microscopy using Lukosz bound apodization reduces pixel negativity at no resolution cost

Abstract

1. Introduction

2. Theory

2.1. Overview of 3D SIM reconstruction

2.2. The 3D SIM cut-off

2.3. The 3D SIM Lukosz-bound

2.4. Comparison of apodization functions

3. Experiments

4. Discussion

References and links

Supplementary Material (4)

Cited By

Figures (5)

Tables (1)

Equations (33)

Optics Express