Image reconstruction for structured-illumination microscopy with low signal level

Kaiqin Chu; Paul J. McMillan; Zachary J. Smith; Jie Yin; Jeniffer Atkins; Paul Goodwin; Sebastian Wachsmann-Hogiu; Stephen Lane

doi:10.1364/OE.22.008687

1. Introduction

Super-resolution fluorescence microscopy is an active research field with many techniques developed over the past decade such as STED [1], PALM [2], structured illumination microscopy (abbreviated as SIM) [3], STORM [4], and SOFI [5], among others. No matter what the technology is, the super-resolution is achieved through a two-step process. In the first step, the data is acquired by a diffraction limited imaging system and in the second step, the data is processed using some prior information about the way the image is recorded. Among current techniques, SIM is attractive due to its wide field nature and ability to use general-purpose florescence labels. However just as with any fluorescence microscopy, the applicability of SIM suffers due to the photobleaching of fluorophores and phototoxicity to living biological samples. So far, most efforts to address this problem are focused on the first step (data acquisition) where a light sheet (either through a cylindrical lens or total internal reflection) is used to limit the exposure of the sample, thus alleviating the photobleaching issue [6,7]. Here we are interested in addressing the problem by improving the data processing step. In particular, we propose to use a total variation (abbreviated as TV) constrained image reconstruction method for data processing in SIM. The current standard reconstruction method, which is a Wiener-type filter, requires that each image is taken with a high signal level [3]. This accelerates the photobleaching of the sample and limits the number of time points of live cell images. In this paper, we will demonstrate that the proposed TV method requires much less signal and slows down the photobleaching of the sample. The sample can therefore be observed for a longer time, and additionally, the shorter exposure times may lead to faster acquisition speeds.

The constraint of total variation (abberivated as TV) was first proposed to perform denoising [8], while currently it has been extended to image deconvolution problems. TV has been used in photography as well as 2D and 3D microscopy [9–11]. However, there is no reported work on using TV for SIM, partly because there are multiple steps involved in the processing of SIM data. As explained in [3], the raw data is a mixture of low and high frequency bands with high bands being scaled and down-shifted into the passing band of the imaging system. In particular, the scaling factor between bands depends on the contrast and the starting phase of the illumination pattern. As the pattern is typically formed through the interference of 2 or 3 beams, the presence of the biological sample in the beam path will affect the contrast and initial phase of the interference pattern seen by the microscope. Thus the contrast and initial phase vary from sample to sample and cannot be calibrated before the experiment. Here we absorb those two parameters into one and call it the scaling factor since it was used to scale the higher order optical transfer function with respect to the central one. A complete description of SIM processing is given in [3]. Briefly, we first need to separate the bands from the raw data and then find the shift frequency and the scaling factor for the high frequency bands. Lastly we can assemble the bands properly with a Wiener-type filtering method as follows (Eq. (11) of [3]):

S (k) = \frac{\sum_{d, m} a_{d, m}^{*} O_{m}^{*} (k + m p_{d}) I_{d, m} (k + m p_{d})}{\sum_{d, m} {| a_{d, m} O_{m} (k + m p_{d}) |}^{2} + w^{2}} A (k) .

Here S is the estimate of the sample information in the Fourier space, O_m is the optical transfer function for the m^th order band, I_d,m is the measured band for the m^th order and d^th direction, and a_d,m is the scaling factor corresponding to m^th order and d^th direction. In the original equation, a_d,m was absorbed in the optical transfer function but is explicitly expressed here. The superscript ^* implies the complex conjugate. p_d is the spatial frequency of the d^th direction illumination pattern and itself is a vector. w is the Wiener parameter and A is an apodization function to suppress the ringing artifacts often associated with the Wiener filter. The experimental parameters such as a_d,m and p_d are estimated from the raw data. When the signal is high, we can follow standard steps to find the scaling factors and shift frequencies. However, when the signal is low, the image shot-noise and the electronic detector noise lowers the contrast of the interference pattern seen by the camera, thus making the estimation of the scaling factor difficult and usually resulting in a low amplitude estimate. In this paper, we suggest a new way of estimation which reduces the effect of the noise. Combining this new method of preprocessing and TV constrained image reconstruction yields much improved results with lower exposure levels, thus prolonging the lifetime of the fluorophores within the sample for study.

2. Estimation of experimental parameters through pre-filtering

In 3D SIM, 15 images are acquired by shifting the illumination pattern in 5 phase steps per orientation of the illumination pattern and repeating for 3 orientations. If the phase steps are not identical, a phase estimation should be performed to achieve better image quality [12, 13]. In this paper we assume that the phase steps are identical. We also assume that the distortions to the excitation pattern, due to samples or other experimental factors, can be ignored. For each direction of the illumination pattern, we will have a central band and a number of side bands. Before the actual deconvolution, the experimental parameters such as the scaling factors and the shift frequencies of the side bands with respect to the central band have to be determined. Typically they are found through examination of the overlap areas between the side bands and central band [3]. The recovered sample information from the overlap area obtained with, for example, inverse-filtering, should be identical, i.e.,

\frac{I_{d, m} (k + m p_{d})}{a_{d, m} O_{m} (k + m p_{d})} = \frac{I_{d, 0} (k)}{O_{0} (k)},

where k and k + mp_d are confined to the overlap area. Thus, the scaling factor between bands (central band and the side band), a_d,m, can be estimated through

a_{d, m} = \frac{I_{d, m} (k + m p_{d}) O_{0} (k)}{I_{d, 0} (k) O_{m} (k + m p_{d})}

A more robust way to estimate the scaling factor is by using OTF compensation and averaging over the overlap areas, i.e.,

a_{d, m} = \frac{\int S_{d, 0}^{ol, *} (k) S_{d, m}^{ol} (k) d k}{\int {| S_{d, 0}^{ol} (k) |}^{2} d k},

where the superscript ol means overlap and

\begin{array}{l} S_{d, m}^{ol} (k) & = & \frac{I_{d, m} (k + m p_{d}) O_{0} (k)}{{| O_{m} (k + m p_{d}) |}^{2} + {| O_{0} (k) |}^{2}} \\ S_{d, 0}^{ol} (k) & = & \frac{I_{d, 0} (k) O_{m} (k + m p_{d})}{{| O_{m} (k + m p_{d}) |}^{2} + {| O_{0} (k) |}^{2}} \end{array}

When the signal is high or when the noise can be ignored, Eq. (4) is a good estimate of the scaling factor of the side band with respect to the central band. However, with low exposures, the noise comes from both detection electronics and inherent photon noise and a simple background subtraction would not remove the noise completely or produce a white noise spectrum. The noise is aggregated in the central band, using Eq. (4) directly will give a low valued result in general because of the extra noise spectrum in the central band.

In order to demonstrate the noise effect on the central band, we first illuminate 100 nm sized fluorescent beads (TetraSpeck green fluorescent microspheres, Life Technologies) with high and low exposure (the transmission filter for the illumination is set to be T = 10% and 0.1%, respectively, while the exposure time is t = 5 ms for both cases. The illumination laser power at the sample plane is about 2.5 mW when the transmission filter is 100%, its wavelength is 488 nm and the center emission wavelength of the beads is 512 nm). The data was obtained with a commercial structured illumination microscope(DeltaVision OMX, API-GE Healthcare, with an Olympus Plan Apo N 60x 1.42NA objective, the resolution is 120 nm in XY and 300 nm in Z). Figures 1(a) and 1(b) show the diffraction limited images which are the average of 5 phase-shifted images. The Fourier transforms of these images are the central bands, amplitudes of which at the f_z = 0 plane are shown in Figs. 1(c) and 1(d) corresponding to high and low exposure respectively. We see that when the exposure is high, a gaussian-shaped signal spectrum is clearly visible (Fig. 1(c)). However when a low exposure is used, we can only see the bright stripes along the axes and a peak at the origin. In order to see where the noise aggregates, two cross-sections of the spectrum are computed. One is the radial average of the spectrum on the f_x and f_y axes while the other is radially averaged over areas other than the axes. The cross-sections of the central band for both high and low signal cases are shown in Figs. 1(e) and 1(f), respectively. In the plots, the cutoff frequencies are marked by the vertical dashed lines and the spectrum beyond this cutoff frequency is purely due to noise. We can see that for the high signal case, both cross-sections within the cutoff frequency are much higher than the noise level. In the low signal case, the on-axis cross-section is significantly higher than the noise but most of the non-axis cross-section is only slightly higher than the noise level. This means that most of the spectrum in the overlap area is low for the central band except areas near the origin and axes. Thus when a low exposure is used, we can only see the bright stripes along the axes and a peak at the origin (Fig. 1(d)). Clearly the noise spectrum is not uniform and it aggregates more at the low frequencies and region along the axes(Fig. 1(d)). With this kind of central band, the estimated scaling factor tends to be much lower than in the high signal case. In this example, the amplitudes of the scaling factors for the 2^nd order bands of three directions are (0.37, 0.32, 0.40) for the high signal case while they are (0.07, 0.06, 0.09) for the low signal case. Using scaling factors this small in the reconstruction would enhance the noise and lower the contrast of the final reconstructed image.

Fig. 1 Central bands of beads data with high and low signal levels indicated by “H” or “L” inside the parenthesis. (Left column: high exposure with transmission filter set at 10%; Right: low exposure with transmission filter set at 0.1%. The exposure time for both cases is 5 ms). DL: diffraction limited image (a–b); 0^th band (the central band) is the Fourier transform of the average of 5 phase-shifted images, the amplitudes of which at the f_z = 0 plane are shown in (c–d). Cross-sections of 0^th bands (e–f): the vertical lines mark the lateral cutoff frequency of the objective; solid line (blue online) is averaged over the f_x and f_y axes while the dashed line (red online) is averaged over other areas at the f_z = 0 plane.

Download Full Size | PDF

In addition to estimation of the scaling factors, we must also estimate the spatial frequency of the illumination pattern, which we can determine using Eq. (2). If the frequency p_d corresponds to the exact experimental value, the recovered information from the overlap area, either from the central band (left side of Eq. (2)) or from the side band (right side of Eq. (2)), should be identical. The closer p_d is to the true value, the more the left and right side of Eq. (2) are correlated. Thus Eq. (4) describes the correlation between the recovered information from the overlap areas in both bands. By maximizing this correlation, we can find the spatial frequency for d^th direction of the illumination pattern. When the signal is low, the noise in the central band could lead to a wrong estimation of the frequency and direction of the illumination pattern. An Actin sample (FluoCells, Catalog Number: F36924, Life Technologies) was imaged by the same SIM microscope described above with high and low signals (T=30% and 1%, t=10ms). In the high signal case, the spatial frequency for the 3^rd direction is found to be about (190, 49) while in the low signal case, it was estimated as (217, 31). With this quite inconsistent estimate, the band is shifted to a wrong center in the Fourier space. Figure 2 shows the reconstructed images by the standard method where (a) is the high signal case while (b) is the low signal case. The final reconstructed image with low signal exhibits periodic artifacts compared to the image with high signal. For both signal levels, the same phase shifting algorithm is applied to separate the bands. The reconstruction is artifact free with high signal but with low signal, there are periodic artifacts. This means that it is the noise, not the phase shift between illumination patterns that causes the periodic artifacts in Fig. 2(b). In particular, it is the noise in the central band that leads to a wrong estimate of the spatial frequency of the illumination pattern. Although we could address this by limiting the search range, nevertheless the scaling factor would be too small and the noise will be further enhanced.

Fig. 2 Actin results with high and low signal levels indicated by “H” or “L” inside the parenthesis (Transmission filter is 30% for high signal and 1% for low signal; exposure times are 10 ms for both cases.) Wiener: reconstructed by the standard Wiener filtering where experimental parameters are estimated without pre-filtering; WP: Wiener method with pre-filtering. Scale bar: 5 μm.

Download Full Size | PDF

Thus the noise in low signal images significantly affects our estimation of experimental parameters. To address this, we propose to apply the following pre-filtering to the central band to eliminate some of this noise effect:

I_{d, 0} (k) = {\begin{array}{l} 0 & where | k, k_{x}, k_{y} | \leq Δ k, \\ I_{d, 0} (k) & elsewhere; \end{array}

With this filtering, the noise spectrum in the neighborhood of the origin and axes are blocked. Signal values in the unblocked regions are used to estimate the scaling factor and the spatial frequency of the illumination pattern. We apply this method to pre-process the same actin and beads data (shown in Figs. 2– 3) and reconstruct with a Wiener-type filter (Eq. (1)). The reconstructed results are shown in Figs. 2(c) and 2(d), and in Figs. 3(c) and 3(d). As a comparison, the results from the standard method are shown in Figs. 2(a) and 2(b), and in Figs. 3(a) and 3(b). In the actin case, our Wiener method actually gives a higher contrast image compared to the current method when the signal is high. This is because our estimate of the scaling factors, even at high signals are better than the standard Wiener method. The estimated scaling factors with the current standard method for 2^nd order bands of three directions are (0.17, 0.16, 0.18) at high signal levels, which are less than half of the values with the beads data even though both samples are thin and the data were taken with the same instrument. With our pre-filtering, these values are all closer to the scaling factor values for the beads data (∼ 0.4). With our estimated scaling factors, the various bands are properly OTF compensated in the Wiener method and the final image has a better contrast laterally and longitudinally (Fig. 2(c)). Note that the images shown here are plotted in the same scale and there is no subtraction or any other local processing. This means that our preprocessing of the data, namely filtering of the central band, helps in providing a better estimate of the scaling factors when the biological sample is complex. When the signal is low, the original estimate of the wave vector for the third direction of the illumination pattern is quite far from the value when the signal is high (high: (∼190, ∼49); low: (∼217, ∼31)). With our pre-processing, the estimate is much closer to the high signal value and the Wiener reconstructed image, while noisy due to the low signals, exhibits no periodic artifacts.

Fig. 3 Reconstructed beads images with two exposure levels as shown in Fig. 1. Wiener: reconstructed by standard Wiener method; WP: by Wiener method with pre-filtering; TV: by TV method with pre-filtering. The letters “H” and “L” in the parenthesis mean high and low signal respectively. Scale bar: 0.7 μm.

Download Full Size | PDF

With the beads data, our Wiener method gives a similar result to the standard Wiener method when the signal is high (comparing Figs. 3(a) and 3(c)). When the signal is 100 times lower, our Wiener method still gives reasonably good results (less noise) while the standard Wiener method enhances the noise (comparing Figs. 3(b) and 3(d)). This also verifies that our estimation method of the scaling factors in the low signal case is sufficient to give good quality reconstructions.

Note that the filtering of the central band through Eq. (5) is used only in the preprocessing of the data. For the subsequent reconstruction (Eq. (1)), the original central band is used instead to avoid artifacts in the final images. As to the value of Δk, we empirically set it to be 0.125μm⁻¹ in this paper.

3. Multi-channel TV constrained algorithm for SIM

In Section 2 we have discussed how to get a better estimation of the scaling factors and the wave vectors of the illumination when the signal is low. In this section we are going to describe how to filter and assemble the bands when the signal is low. The current method (Eq. (1)) is a Wiener-type filter that depends on strong signals. By contrast, using TV as a constraint has proven successful even with extremely low signals [9–11]. So far the TV constraint has been used in the fields of photography and conventional microscopy. In structured illumination microscopy, after the bands are separated the side bands have peaks quite far away from the zero-frequency. Adding a TV constraint directly to each side band will produce an overly smoothed reconstruction result. Therefore we must modify the traditional TV constrained reconstruction problem to account for this. In this paper we treat SIM as a multichannel imaging system where one object s is imaged by 15 channels. In each channel, the object is multiplied by a modulation function e^{j2πmp_d · r}, where r = (x, y, z), before passing through the channel. The channel-specific transfer function is a_d,mO_m. Thus our modified TV constrained optimization problem can be formulated as follows:

min_{s} \int ({| \nabla s |}_{1} + \frac{μ}{2} {\sum_{d, m} | s e^{j 2 π m p_{d} \cdot r} \otimes (a_{d, m} p_{d, m}) - i_{d, m} |}^{2}) d x d y d z,

where

{| \nabla s |}_{1} = \sqrt{{(\frac{\partial s}{\partial x})}^{2} + {(\frac{\partial s}{\partial y})}^{2} + {(\frac{\partial s}{\partial z})}^{2}}

defines the variation of the object at point (x, y, z) and an integral of it will be the total variation of the object. ⊗ is the convolution operator, p_d,m and i_d,m are the inverse Fourier transforms of the optical transfer function O_m and the band I_d,m, respectively (they can be thought as the point spread function and image of the channel corresponding to direction d and order m). The first term inside of the integral in Eq. (6) is the TV constraint and the second term is the data fidelity term. The parameter μ is the relative weight of these two terms during the optimization.

Because of the l1 term in the optimization function Eq. (6), the traditional gradient optimization method is very slow. A fast and convergence-guaranteed way to solve this optimization problem was developed recently through the Bregman variable splitting method [14]. In this method, new variables are introduced to approximate the derivatives of the sample function, i.e., W = (W_x, W_y, W_z) ⇀ ∇s. The difference between them, b = (W_x, W_y, W_z) − ∇s, is accounted for by reinserting it back into the optimization problem during each iteration. Thus we can solve the optimization problem of Eq. (6) by iterating

\begin{array}{l} (s^{n + 1}, W^{n + 1}) = \\ min_{s, W} \int ({| W |}_{1} + \frac{μ}{2} \sum_{d, m} {| s e^{j 2 π m p_{d} \cdot r} \otimes (a_{d, m} p_{d, m}) - i_{d, m} |}^{2} + \frac{β}{2} {| W - \nabla s - b^{n} |}^{2}) d x d y d z \end{array}

b^{n + 1} = b^{n} + (\nabla s^{n + 1} - W^{n + 1}) .

Here the parameter β is positive and does not need to be large to ensure W ⇀ ∇s [15]. For details of variable splitting in the Bregman method, we refer the readers to [14]. s and W can be further decoupled by solving them in two steps:

s^{n + 1} = min_{s} \int (\frac{μ}{2} \sum_{d, m} {| s e^{j 2 π m p_{d} \cdot r} \otimes (a_{d, m} p_{d, m}) - i_{d, m} |}^{2} + \frac{β}{2} {| W^{n} - \nabla s - b^{n} |}^{2}) d x d y d z

W^{n + 1} = min_{W} \int ({| W |}_{1} + \frac{β}{2} {| W - \nabla s^{n + 1} - b^{n} |}^{2}) d x d y d z .

Since Eq. (9) is quadratic, sⁿ⁺¹ can be directly solved as:

s^{n + 1} = ℱ^{- 1} (\frac{ℱ (\nabla^{T} (W^{n} - b^{n})) + \frac{μ}{β} \sum_{d, m} a_{d, m}^{*} O_{m}^{*} (k + m p_{d}) I_{d, m} (k + m p_{d})}{{| ℱ (\nabla) |}^{2} + \frac{μ}{β} \sum_{d, m} {| a_{d, m} O_{m} (k + m p_{d}) |}^{2}});

where the superscript T means “conjugate transpose”. And Eq. (10) can be solved efficiently by the shrinkage operation [9]:

W^{n + 1} = {\begin{array}{l} (| \nabla s + b | - \frac{1}{β}) \frac{\nabla s}{| \nabla s |}, & where | \nabla s + b | > \frac{1}{β}, | \nabla s | > 0 \\ 0, & elsewhere; \end{array}

Thus we formulate our reconstruction algorithm as follows:

Algorithm (TV-SIM)
Initialization: s⁰ = 0, W⁰ = 0, b⁰ = 0
For n = 1, 2,... repeat until a stopping criterion is reached
Update sⁿ according to Eq. (11)
Update Wⁿ according to Eq. (12)
Update bⁿ according to Eq. (8).
Output: sⁿ

Since in each step, sⁿ, Wⁿ and bⁿ can be solved analytically, each iteration is computationally cheap. At the same time, due to the updating of the difference between ∇s and W through application of Eq. (8), the convergence of the TV algorithm is achieved with only a few iterations [15]. The stopping criteria could be a maximum iteration number or based on the difference Δs = sⁿ⁺¹ − sⁿ. In this paper, we stop the algorithm after 20 iterations and usually the algorithm is close to convergence at this point. Regarding the parameters, smaller μ or larger β will give a smoother image with reduced resolution. We have used μ = 1 × 10⁶, β = 100 for the TV method throughout this paper. The Wiener parameter w = 10⁻³ is used for both Wiener methods. A larger w will give a more heavily smoothed image with reduced resolution.

4. Testing on the experimental beads data

We first test our algorithm on the beads data described in Section 2. The TV constrained reconstructed results for both high and low signals are shown in Figs. 3(e) and 3(j). We see that when the signal is high, the result from the TV method is similar to the Wiener methods (both the original Wiener method and our Wiener method with pre-filtering). When the signal is low, the TV result is still similar to the result with 100 times higher signal (compare Fig. 3(f) to Fig. 3(a), 3(c) or 3(e)). The result of our Wiener method at the low signal case, though much better than the result from the standard Wiener method, is still a bit noisy compared to the crisp image of the high signal or the TV result with low signal. Thus, together with the preprocessing step and the subsequent TV constrained reconstruction we are able to obtain an image of similar quality even with two orders of magnitude lower light exposure.

At the same time, we want to make sure that we did not sacrifice resolution when suppressing the noise. A resolution check is conducted by studying the intensity profiles of single beads along the x-,y- and z-axis. Figure 4 shows the profiles of 3 randomly chosen beads. The results from the standard method are plotted in the left side of each figure while the results of TV method are horizontally shifted and plotted on the right side. In the high signal case, we can see that the width of the beads reconstructed by both methods are similar to each other, about 110 nm in x or y and 300 nm in z. When the signal is low, beads reconstructed by the standard Wiener method start to become irregular while with the TV method they remain Gaussian-shaped and have similar widths as with the high signal. Note here we have used the same parameters for both high and low signal cases.

Fig. 4 Single beads profile along x,y,z. In each plot, the left side shows the result from the standard Wiener method and the right side shows the result from TV method.

Download Full Size | PDF

5. Testing on fixed samples

When it comes to biological samples, a rule of thumb for the commercial SIM in our lab to yield a good reconstructed image is that the maximum intensity of the acquired raw images should be 3 to 11K counts. In the next example we show results when actin data was taken with significantly lower exposure (using the same actin slide as described in Section 2 but with a different field-of-view). Figure 5 is the case when the maximum intensity is ∼ 100 counts. The reference diffraction limited image is shown in the top left panel. The result of the standard method is shown in the top right panel while the result from the pre-filtering and TV reconstruction is in the bottom left panel. A magnified region of interest corresponding to the white box of Fig. 5(a) is shown in the bottom right panel. Here we see once again that the low signal leads to noise enhancement with the standard method (Figs. 5(b) and 5(d)). However with our method, we are able to reconstruct the super-resolved image with high quality (Figs. 5(e) and 5(f)). The width of the actin filaments in our reconstructed image is about 100 nm. Note that the signal level in Fig. 5 is about 100 times lower than the current ”rule of thumb” level.

Fig. 5 Actin results where the exposure setting are: transmission filter 0.1% and exposure time 30 ms. The peak photon number detected is ∼ 100. (a) DL: diffraction limited image (shown with increased contrast to see structures); (b) Wiener: reconstructed by the standard Wiener method; TV: reconstructed by TV method. Zoom-in versions corresponding to the white box in (a) are shown in (c)–(e). Scale bar in (a): 10μm; in (f): 2μm.

Download Full Size | PDF

6. Test on the live cell data

In this section we will test our algorithm on live cells to see whether we can expand the number of time points imaged by structured illumination microscopy. The samples are HeLa cells. They were cultured in a 35 mm glass bottom imaging dish (Ibidi, Germany), stained with Mitotracker green (Life Technologies) and washed in Opti-MEM media. Samples were imaged at room temperature with time series taken every 20 seconds. First the samples were measured with different exposure settings and the imaging thickness is 3 μm (scanned over z with step size 125 nm). The raw data were submitted to the standard Wiener and our TV reconstruction algorithms to obtain the overall results of the structured illumination microscopy. For different exposure settings, different cells were imaged.

In Fig. 6 we show the results with a high exposure setting (transmission filter T=100%, exposure time t=2 ms) and in Fig. 7 the low exposure setting (T=10%, t=1ms). For the high setting, 10 time points of data were taken and for low setting, 50 time points. Here we can see clearly that due to photobleaching in the high exposure setting, the quality of the reconstructed images is decreasing with time for both methods even though the TV results are better than the Wiener results for all time points. With low exposures, the photobleaching is less and the noise is more. The Wiener method enhances the noise and fails to reconstruct the details of the mitochondria while the TV method suppresses the noise and the structure of mitochondria can be clearly seen in the reconstruction. For the time lapse images, please refer to the media files ( Media 1, Media 2, Media 3, and Media 4).

Fig. 6 SIM images reconstructed by Wiener filter and TV method where the data is taken with high exposures. Maximum intensity projection of the 3D images are shown here. Top row (a–c): results from the standard Wiener method (time lapse video shown in Media 1); Bottom Row (d–f): results from the TV method (time lapse video shown in Media 2). White scale bar in (a): 5 μm.

Download Full Size | PDF

Fig. 7 SIM images with low exposures. Maximum intensity projection of the 3D images are shown here. Top row (a–d): reconstructed with the standard Wiener method (time lapse video shown in Media 3); Bottom Row (e–h): reconstructed with our TV method (time lapse video shown in Media 4). White scale bar in (e): 5 μm.

Download Full Size | PDF

As the photobleaching is governed by a power law [16], $I (t) = I_{0} {(\frac{t_{A}}{t})}^{α}$ where I(t) and I₀ are the intensities at time point t and the beginning of the measurement respectively; t_A and α are parameters describing the intensity decay. In our experiments, we calculate the total intensity in the raw data and plot the decay of the intensity over time in Fig. 8 for both exposure levels. The curves were fitted with the power law and we found that t_A = 1.331 and 1.113 for high and low exposures respectively, α = 0.9805 and 0.1608 for high and low cases respectively. If we assume that the image quality at the 3^rd time point in the high signal case (Fig. 6(b)) is minimally acceptable, the image at the 50^th time point in the low signal case is of similar quality. Thus we can conclude that for this sample, we have increased the number of usable time points by at least 15 times. However, the exact improvement will depend on the sample and on the photobleaching property of the fluorophores. For instance, if there are both very bright and very dim structures in the sample, the quality of the reconstructed image will be lower and we might expect faster disappearance of those dim structures over time than the bright structures.

Fig. 8 Intensity decay with high and low exposure experiments whose results are shown in Fig. 6 and 7. Dashed line (blue online) is the power law curve for the high signal case and the solid line (blue online) is the power law curve for the low signal case.

Download Full Size | PDF

7. Discussion and summary

In this paper we have developed an image processing algorithm for structured illumination microscopy when the exposure is low. The noise hinders the current standard method and the reconstructed images are either quite noisy or have artifacts. We demonstrated the noise effect on estimating experimental parameters (such as the scaling factors and spatial frequency of the illumination pattern) and proposed pre-filtering on the central band to reduce the noise effect. When assembling all bands to form a super-resolved image, we also proposed a total-variation constraint to suppress the noise. With these two processes, we successfully reconstructed images with 100 times lower exposure on fixed samples (beads and actin). The final image maintains the same resolution as an image with high signals. With the reduced signal level requirement, photobleaching of fluorophores in the live cell should be slower and more time points of data can be obtained before the intensity drops too low. We tested our algorithm on live mitochondria and found that at least 15 times more time-points could be expected. All verifications of our method are on data taken with commercial structured illumination microscopes. Using a spatial light modulator to rapidly generate illumination patterns in SIM, Shao et al [17] reported 3D live cell images with 50-time points using the standard data processing [3]. With our proposed method, we believe many more time points are possible with SIM microscopes capable of live sample imaging. The current algorithm processes the images using MATLAB. For an image size of 512×512×375 with 50 time points, the total processing time was about 2 hours using a desktop PC equipped with a 3.4GHz,i7 processor and 12 GB of RAM. Although this time equates to only approximately 2 minutes to reconstruct a whole 3D volume per time point, we expect that the processing could be significantly sped up using a GPU to parallelize the processing.

Acknowledgments

We thank Fiona Houghton (The University of Melbourne) for preparing samples for the live cell imaging experiments. This work was funded by NSF Partnerships for Innovation: Accelerating Innovation Research (Award number: 1343479).

References and links

1. T. A. Klar and S. W. Hell, “Subdiffraction resolution in far-field fluorescence microscopy,” Opt. Lett. 24, 954–956 (1999). [CrossRef]

2. E. Betzig, G. H. Patterson, R. Sougrat, O. W. Lindwasser, S. Olenych, J. S. Bonifacino, M. W. Davidson, J. Lippincott-Schwartz, and H. F. Hess, “Imaging intracellular fluorescent proteins at naometer resolution,” Science 313, 1642–1645 (2006). [CrossRef] [PubMed]

3. M. G. Gustafsson, L. Shao, P. M. Carlton, C. 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]

4. B. Huang, W. Wang, M. Bates, and X. Zhuang, “Three-dimensional super-resolution imaging by stochastic optical reconstruction microscopy,” Science 319, 810–813 (2008). [CrossRef] [PubMed]

5. T. Dertinger, R. Colyer, R. Vogel, J. Enderlein, and S. Weiss, “Achieving increased resolution and more pixels with superresoltion optical fluctuation imaging (sofi),” Opt. Express 18, 18875–18885 (2010). [CrossRef] [PubMed]

6. F. C. Zanacchi, Z. Lavagnino, M. P. Donnorso, A. Del Bue, L. Furia, M. Faretta, and A. Diaspro, “Live-cell 3D super-resolution imaging in thick biological samples,” Nat. Methods 8, 1047–1049 (2011). [CrossRef]

7. L. M. Hirvonen, K. Wicker, O. Mandula, and R. Heintzmann, “Structured illumination microscopy of a living cell,” Eur. Biophys. J. 38, 807–812 (2009). [CrossRef] [PubMed]

8. L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D: Nonlinear Phenomena 60, 259–268 (1992). [CrossRef]

9. Y. Wang, J. Yang, W. Yin, and Y. Zhang, “A new alternating minimization algorithm for total variation image reconstruction,” SIAM J. Imaging Sci. 1, 248–272 (2008). [CrossRef]

10. S. Setzer, G. Steidl, and T. Teuber, “Deblurring poissonian images by split Bregman techniques,” J. Vis. Commun. Image Representation 21, 193–199 (2010). [CrossRef]

11. N. Dey, L. Blanc-Feraud, C. Zimmer, P. Roux, Z. Kam, J.-C. Olivo-Marin, and J. Zerubia, “Richardson-Lucy algorithm with total variation regularization for 3d confocal microscope deconvolution,” Microsc. Res. Tech. 69, 260–266 (2006). [CrossRef] [PubMed]

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

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

14. T. Goldstein and S. Osher, “The split Bregman method for L1-regularized problems,” SIAM J. Img. Sci. 2, 323–343 (2009). [CrossRef]

15. W. Yin and S. Osher, “Error forgetting of Bregman iteration,” J. Sci. Comput. 54, 684–695 (2013). [CrossRef]

16. F. Cichos, C. Von Borczyskowski, and M. Orrit, “Power-law intermittency of single emitters,” Curr. Opin. Colloid Interface Sci. 12, 272–284 (2007). [CrossRef]

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

Image reconstruction for structured-illumination microscopy with low signal level

Abstract

1. Introduction

2. Estimation of experimental parameters through pre-filtering

3. Multi-channel TV constrained algorithm for SIM

4. Testing on the experimental beads data

5. Testing on fixed samples

6. Test on the live cell data

7. Discussion and summary

Acknowledgments

References and links

Supplementary Material (4)

Cited By

Figures (8)

Equations (13)

Optics Express