The performance of structured illumination microscopy (SIM) is hampered in many biological applications due to the inability to modulate the light when imaging deep into the sample. This is in part because sample-induced aberration reduces the modulation contrast of the structured pattern. In this paper, we present an image restoration approach suitable for processing raw incoherent-grid-projection SIM data with a low fringe contrast. Restoration results from simulated and experimental ApoTome SIM data show results with improved signal-to-noise ratio (SNR) and optical sectioning compared to the results obtained from existing methods, such as 2D demodulation and 3D SIM deconvolution. Our proposed method provides satisfactory results (quantified by the achieved SNR and normalized mean square error) even when the modulation contrast of the illumination pattern is as low as 7%.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
Structured illumination microscopy (SIM) provides improved optical sectioning (OS)  and super-resolution [2–4] compared to widefield microscopy; however, SIM performance is restricted in imaging optically-thick biological samples . This is because of light scattering and depth-induced spherical aberration (SA), which significantly reduce the modulation contrast of the structured pattern when imaging deep into the sample [5,6]. Adaptive optics systems have been used in SIM for aberration correction by which imaging of an approximately 40-µm thick sample has been demonstrated . Other approaches are based on either the combination of structured illumination with light-sheet microscopy [7–9], or the use of multiple focal spots as the excitation plane in a fluorescence microscope . In this paper, we focus on a computational approach to restore images from incoherent-grid-projection SIM raw data affected by depth-induced SA.
The SIM forward model presented in Ref  describes the intermediate raw image captured from a grid-projection SIM instrument . Later, the model was extended to include the impact of SA and a non-iterative deconvolution method based on this model was developed to recover the 3D image from the SIM raw data in the presence of SA [11–13]. In Ref  an expectation maximization (EM) iterative algorithm based on a model that did not account for SA was presented, in which an improvement in the resolution was demonstrated; however, this study was limited to noiseless simulated data. Other model-based SIM reconstruction methods that do not account for SA and do not require knowledge of the illumination patterns to obtain an artifact-free restoration have also been developed [15–18]. Although these methods work well in the case when the sample is sufficiently modulated, they do not specifically address the reduced modulation contrast that occurs while imaging thick biological tissue.
In this study, we concentrate on the effect of the reduced modulation contrast on the performance of SIM, and specifically its impact in imaging optically thick biological specimen. In our investigation, we have used SIM raw data obtained from a commercially available SIM device (Zeiss ApoTome.2), in which the structured illumination patterns are generated by the incoherent projection of a sinusoidal grating . Since light emits incoherently from a fluorescent sample, the modulated light in SIM is always subject to the incoherent optical transfer function (OTF)  of the imaging system. The OTF attenuates the fringe contrast (especially at high modulation frequencies), which in turn reduces the signal-to-noise ratio (SNR) in the final restored images. In the ApoTome.2 SIM, the structured pattern is doubly attenuated (by the illumination and the detection OTFs); thereby, the modulation contrast of the illumination pattern is significantly reduced. Because of this, in the ApoTome.2, the recommended modulation frequency, fM, is very low (for example, fM is equal to only 9% of the cut-off frequency of the imaging system, fc, for a 63 × /1.4 numerical aperture (NA) oil immersion lens), and as a consequence this type of SIM is primarily used to improve the OS capability in the widefield imaging system. Although the low modulation frequency recommended by the manufacturer for ApoTome SIM could be increased further (e.g. it could be doubled), in practice its value is limited because of the reduction of the contrast through the illumination system, which results in images with poor SNR.
In this paper, we present an image restoration approach for ApoTome.2 SIM data that suffers from reduced modulation contrast. The novelty of our proposed approach is the cascaded application of existing methods that are sequentially applied to the SIM raw data to first increase the modulation fringe contrast and then, restore the 3D final image using a 3D iterative EM method  and a synthetic SIM point spread function (PSF) that accounts for spherical aberration derived in Section 2.2. In the proposed method, we approximate the OS image by a linear model using a 3D synthetic PSF, which we define after the demodulation step. In the ApoTome.2 instrument, the OS image obtained after the demodulation step contains only the in-focus information; therefore, for processing, we only need to capture data from the region of interest (ROI), where the object is located. Our proposed approach improves the OS capability of the system and increases the SNR and the contrast in the restored images, thereby revealing more details within the object’s structure compared to the 2D demodulation (2D-DMOD) [1,21] and the non-iterative 3D deconvolution (3D-DCV) [12,13] methods. We demonstrate performance of the proposed approach by restoring simulated and experimental images obtained using the Zeiss ApoTome.2 SIM instrument. We also compare our results with the results obtained from the ZEN.2 software (Carl Zeiss, Jena, Germany). Initial findings of this study were reported in a conference publication .
Despite the fact that we have carried out this study with the raw image obtained from the ApoTome SIM, this computational approach can be extended for other SIM implementations in general. We have recently extended and applied our method to a new system integrating wavefront encoding (WFE), developed to reduce the impact of SA, with ApoTome SIM . Because WFE causes the PSF to spread further , this results in a reduction of the SNR in the recorded image and in a lower modulation fringe contrast.
The paper is organized as follows: Section 2 describes background and theory of the proposed restoration approach. Section 3 describes the methodologies used in the investigation, which include simulation, sample preparation, and image restoration procedures. Results of this investigation study are reported in Section 4, while Section 5 summarizes the findings.
2. Background and theory
2.1 Forward imaging model
As stated earlier, in the case of ApoTome SIM, the 3D illumination pattern is produced by the incoherent projection of a grating onto the object plane, and is mathematically represented by:21] was later extended to account the effect of depth-variant (DV) SA [11, 12]. The DV forward SIM model uses multiple DV PSFs computed at different depth locations and computes the raw SIM image using a strata-based model instead of a convolution . In this paper, we use an approximation of the DV forward imaging model, which is accurate if the object has a small axial extent at a particular depth (as is the case of the numerical objects and biological samples used in the studies presented here). This approximation was motivated by an earlier study  in which SA was accounted using a single aberrant PSF. Specifically, we use the model in  and a single aberrant PSF at a particular depth. Thus, the corresponding N forward (raw) SIM images at an imaging depth d below the cover slip are described by:Eq. (2) reflects the case where the lateral magnification of the illumination and the recording system are equal to one.
As evident by Eq. (2) the sinusoidal pattern of Eq. (1) is attenuated by both the illumination and the detection processes. Because in epi-illumination fluorescence microscopy, the same objective lens is used for the illumination and the detection, consequently the corresponding PSFs at the best focus (bf) plane can be assumed to be equal . Note that the illumination and the detection PSFs depend on the corresponding excitation and emission wavelengths; however, we have ignored the difference for the sake of simplicity. Therefore, the modulation contrast (mc) of the illumination pattern at the bf plane of the objective lens, can be expressed in terms of the corresponding 2D detection OTF by the following:Eq. (3) is the theoretical maximum limit, which is valid for a thin object; however, in practice, it may further reduce depending on the sample’s optical thickness. In our previous work , we demonstrated that for a 63 × /1.4 NA oil immersion lens, the theoretical fringes’ visibility drops to mc < 10% at depths d > 20 µm for fM 0.3fC (see Fig. 1(a) in ).Eq. (4) it is evident that when the value of mc is low, the modulated components, and , are attenuated, which will result in a poor SNR in the demodulated image.
2.2 Inverse imaging problem solution
Because the synthetic SIM-PSF, , is axially-confined by the axial modulation function , the contribution of out-of-focus light at each axial location in the modulated sectioned components of the raw SIM data is significantly less compared to the out-of-focus light in the widefield component. Therefore, the intensity at the jth axial plane, of the modulated sectioned component imaged at the best focal plane of the objective lens can be approximated by:Equation (7) shows that each transverse section is diffraction-limited due to the aberrated 2D PSF at the best focal plane,, and thus, by applying a 2D Wiener filter , the diffraction effects can be removed, thereby, improving the fringe contrast in the modulated component. As we will show later, this Wiener preprocessing step improves the overall SNR of the demodulated (i.e. optical-sectioned) image obtained using the magnitude demodulation technique [21, 26]:
It should be mentioned that residual fringes (RF) may be evident in the OS image [Eq. (8)] due to experimental mismatches such as, the difference between the desired and the implemented fringe phases , and the unequal intensity in each phase-shifted images . Furthermore, RF artifacts can also occur because in commercial instruments rectangular grating patterns are used instead of ideal sinusoidal patterns [5,21]. To remove the RF artifact, different methods have been proposed such as phase estimation from the raw data [27–29], and a parameter optimization method to minimize the mismatch between the theoretical and the experimental illumination patterns . Since the RF have frequencies at the harmonics of the sinusoidal illumination pattern [5,21], in this study, the RF artifact is removed by applying a 2D inverted Gaussian filter (IGF) to each 2D optical-sectioned plane, as shown below:Eq. (8); ur = n fM with n being a positive integer; and σ is the standard deviation of the Gaussian function. The value of σ is chosen carefully (i.e. σ ≤ 0.05) so that the object spectrum is minimally affected. If the structured pattern is not purely sinusoidal, multiple harmonics may need to be removed to obtain the RF free image. A similar filter was also used in Ref , while restoring SIM data in order to emphasize the frequency components that fill the missing cone of the OTF to improve the OS capability.
The data processing due to Eqs. (8-9) consists of 2D plane-wise operations; therefore, in order to account for the axial response of the system in the restoration process, we determine a modified synthetic 3D SIM-PSF, , that characterizes the system after the demodulation process [Eq. (8)]. Using the assumption that the OS image intensity is linear to the object intensity (an evaluation of this assumption is given below), we approximate the OS image by the 3D convolution between the object and as:Eq. (11)], the final estimate of the object, is obtained using the iterative EM algorithm  regularized with the Good’s penalty functional .Fig. 1(g) for a 63 × /1.4NA oil immersion objective lens, at d = 30 µm and fM = 0.2 fC.
To investigate the accuracy of the approximation in Eq. (11), we have conducted a study, where we compared in simulation the result obtained with Eq. (11) and the OS image of a test object (a 6-µm in diameter spherical shell with shell thickness equal to 1 µm) obtained with Eq. (9) from SIM data with Poisson noise (SNR = 30 dB, using the methods described in Section 3.1.2). The structure similarity (S-SIM)  and the normalized mean square error (NMSE)  between the approximated OS image [right-hand side of Eq. (11)] and the filtered OS image [Eq. (9)] were found to be equal to 0.94 and 0.10, respectively, when the SNR in the raw data is equal to 30 dB and the modulation contrast, mc = 0.07. Based on this simulated result, we conclude that Eq. (11) approximates the OS image with reasonable accuracy.
In summary, the image restoration approach presented in this paper consists of the following steps: at first, each 2D axial plane of the captured raw image is deconvolved by an aberrated 2D detection PSF using a Wiener filter. In the second step, the OS image is obtained by demodulating the Wiener-filtered images using the magnitude demodulation technique. The third step is the removal of any RF (if needed) from the OS image using an IGF. And, finally, the object is restored by a regularized iterative EM algorithm using a synthetic SIM-PSF that characterizes the system after the initial three steps.
This section describes the methodologies used in the investigation studies, which include: computation of the widefield PSF (WF-PSF), , and corresponding SIM-PSF, ; SIM raw image simulation, experimental data acquisition techniques, and implementation of the proposed restoration approach method as well as other restoration methods used for the comparison studies.
All simulations were carried out using Matlab (MathWorks, Inc., Natick, Massachusetts).
3.1.1 PSF computation
To compute the WF-PSFs we used the Gibson & Lanni optical path distance model  on a 128 × 128 × 256 grid with a voxel dimension of 0.05 µm3 and assuming a sample embedding medium refractive index (RI) equal to 1.47. PSFs for a 63 × /1.4 NA oil (RI = 1.518) immersion objective lens were computed using an emission wavelength equal to 515 nm. The synthetic SIM-PSFs were computed using Eq. (5) in which the function was calculated using Eq. (6). PSFs were computed for different depths below the coverslip.
3.1.2 Simulation of SIM images
Simulated images were computed using Eq. (4), which requires both the WF-PSF and the SIM-PSF. To demonstrate the improved OS capability achieved with the proposed restoration approach under ideal conditions, i.e., without any noise, we simulated objects (Fig. 1) with one and two fluorescence layers separated by 0.5 µm, a distance smaller than the axial resolution of a 1.4 NA oil-immersion lens at 515 nm, µm. In both cases, three phase-shifted raw SIM images were simulated for φi = 0, 120 and 240 deg using an aberrated PSF computed at a 30-µm depth (the SA is quantified by w40 = −1.08, which is computed using Eq. (4) in Ref ).
The third simulated object is a 3-µm in diameter spherical shell with the shell thickness equal to 200 nm. Inside this spherical shell, there are three 250-nm in diameter beads separated by 300 nm [Fig. 2 (a)]. The forward SIM images were simulated using Eq. (4) discretized on a 128 × 128 × 256 grid at three different modulation frequencies: 0.5, 1.0, and 1.5 µm−1, which correspond to 0.1 ×, 0.2 × and 0.3 × the cutoff frequency of the objective lens (fc), respectively. We simulated the object at three different depth locations below the coverslip: 0 µm, 30 µm, and 60 µm, to show the effect of depth-induced SA on the SIM system (for these different depths the SA coefficient w40 is, respectively, 0, −1.08, and −2.16). Noisy images were simulated as the realizations from a Poisson random process  for a SNR = 29 dB. The noise was incorporated in the simulated data using the poissrnd built-in function of Matlab and the SNR was calculated by comparing the energy of the noiseless image to the energy of the noise using a standard equation reported in an earlier publication . In the case of experimental images, we computed the peak SNR by comparing the maximum energy of the image to the average background intensity computed from the mean of a cubic dark area of 10 × 10 × 10 pixels. The peak SNR of restored images was computed to compare the performance of different restoration methods.
3.2 Experimental data acquisition
For experimental evaluation of our method, we have used images from a test sample, and two biological specimens. The test sample consisted of FocalCheck microspheres, 6-µm in diameter fluorescent green ring stain/blue throughout with shell thickness equal to 1 µm (Invitrogen, Molecular Probes, Carlsbad, CA). A solution of the beads was diluted in deionized water and dried both on a slide and on a ZEISS high precision cover slip (RI = 1.522). To determine the location of the spherical shells, 170-nm in diameter spheres were also dried on the cover slip, which were used as a reference for the coverslip location (considered as the 0-μm imaging depth). The microspheres were mounted in ProLong Diamond Antifade Mountant (Invitrogen, Molecular Probes) with a RI = 1.47, which was cured for 24 hours. During the curing process, the beads localize themselves at random locations. While imaging the beads, their axial location was determined with respect to the cover slip location (marked by the nanospheres) and they were found to be centered at a 3-µm, a 24-µm and a 50-µm depth below the cover slip, respectively. The interface between the cover slip and the mounting medium was also confirmed using differential interference contrast microscopy.
The first biological sample is the Vimentin protein of U-373 MG human Glioblastoma cells mounted in ProLong Diamond (Invitrogen, Molecular Probes). The second biological sample is a commercially-available microscopy slide of multi-pollen grains (Item # 304264, Carolina, Burlington, North Carolina), in which pollens are mounted in Practamount, which is a resin (RI ~1.47). All experimental images were recorded using a 63 × /1.4 NA oil-immersion objective lens using the ApoTome.2 slider in the illumination path that is implemented in the ZEISS AxioImager.Z1. The images of the test sample at different depths and the Glioblastoma cell were captured using the L-illumination grating of the ApoTome.2, whereas the image of a pollen grain was captured using the M-illumination grating. The modulation frequency of the M- and L-illumination grating is, respectively, fM = 0.436 µm−1 (0.09fc) and fM = 0.848 µm−1 (0.17fc) at the object plane. Although the M-illumination grating is recommended for the 63 × /1.4 NA oil immersion lens, we used the L-illumination grating in the case of the test object as the sample is sparse, and in the case of the Glioblastoma cell as it is a thin specimen (<5 µm); thus it provided sufficient modulation contrast at the highest grid frequency. In all the experimental data, three phase-shifted images were captured at phase shifts 0, 120 and 240 deg of the illumination pattern. The samples were scanned axially at an interval of 0.1 µm, and the images were recorded digitally by a CCD camera (AxioCam MRm), which has a pixel dimension equal to 6.45 × 6.45 µm2 (equivalent to 0.1 × 0.1 µm2 in the object space).
3.3 Image restoration
The first three steps of the proposed image restoration approach, (2D Wiener filtering, magnitude demodulation, and the RF removal) were implemented in Matlab. The Wiener filter was implemented using the deconvwnr function of Matlab with the Wiener parameter set to 0.05, which was chosen to reduce the impact of noise without losing significant details. The OS images were obtained by using the Eq. (8). As the ApoTome uses a periodic illumination grating, the RF occurs at the harmonics of the illumination grating. We removed the first harmonic in the case of the Glioblastoma cells, and the first three harmonics in the case of pollen grain using the IGF. Note that the RF of higher order harmonics appeared in the case of the pollen image, which may occur due to the thickness of the pollen. In the IGF [Eq. (10)], the standard deviation of the Gaussian filter, σ was chosen to be equal to 0.05 (a smaller value of σ may not remove the RF completely, and a higher value may distort the object’s spectrum). Finally, the open source software package COSMOS  was used to restore the image using the EM algorithm. The appropriate regularization parameters were empirically determined to reduce the impact of noise without losing fine details in different cases and found them to be within the range of 10−2 to 10−4, in the case of simulated images, and 5 × 10−1 to 5 × 10−2 in the case of experimental images. Note that, although we simulated the same noise level in the raw SIM images, the SNR in the demodulated result reduced with increasing fM and increasing depth below the coverslip, thereby a larger amount of regularization was needed in the EM restoration. Because the OS image obtained from Eq. (9) has a reduced amount of out-of-focus blur, the EM algorithm applied to the image converges fast and the object is restored successfully within 100 iterations. Thus, in this paper, we display and analyze the results after 100 iterations of the EM algorithm. To have an idea about the execution time of the proposed method, we estimated the restoration time required for an image with a 128 × 128 × 256 grid on the Intel(R) Core(TM)i5 CPU @3.40 GHz processor. We observed that the total wall clock time required was less than two minutes (the first two steps of the algorithm took 1.3 seconds, step 3 took 0.33 seconds, and the iterative-EM step took one minute and 40 seconds to complete 100 iterations).
To compare the performance of the proposed approach with other available methods, the images were also restored using the 2D magnitude demodulation technique (2D-DMOD) , which demodulates the image plane by plane, and the non-iterative 3D deconvolution (3D-DCV)  method, which uses the 3D SIM-PSF and a regularized inverse filter in the reconstruction process. The simulated images were restored using the 2D-DMOD and the 3D-DCV methods implemented in Matlab, while the experimental images were restored using the commercial implementation of these methods in the ZEN.2 software, which we refer to as the ZEN-OS and ZEN-DCV methods, respectively. To restore the simulated objects using the 3D-DCV and the proposed method, SIM-PSFs computed at different depths (depending on the location of the test objects) below the coverslip were used. In the case of experimental images, an aberrated PSF computed at a 30-µm depth was used in the restoration, which accounts for a pre-aberration effect of our microscope system . In the ZEN-DCV method, the restoration can be regularized in a scale of 1 to 10, where 1 provides the highest amount of regularization. We determined the appropriate regularization level empirically, and found it to be equal to 2 in the case of experimental data.
In this section, we present restored images obtained from the simulated and the experimental raw SIM images. Results from the proposed method, and from the two other methods discussed above (2D-DMOD/ZEN-OS and 3D-DCV/ZEN-DCV) are shown in a comprehensive comparison to demonstrate the improved performance of the proposed restoration approach.
4.1 Demonstration of achieved OS capability in noiseless simulation
The OS capability of the proposed restoration method is investigated in Fig. 1, where we show a single simulated phase-shifted image [Fig. 1(b)] of a thin fluorescent layer [Fig. 1(a)], and corresponding restored images obtained with different methods [Fig. 1(c-e)]. Only the XZ-sectional images are shown as the object does not vary laterally, and thus restoring the object in the axial direction is challenging. The axial normalized intensity profiles from the different restored images are compared in Fig. 1(h), which demonstrate the improvement in the axial confinement (i.e. achieved OS) obtained with the proposed method. This OS improvement is also evident in Fig. 1(c)-(e), where the XZ cross-sectional view images at different restoration methods are shown. The axial confinement due to the demodulation process is shown in Fig. 1(f), where the function for fM = 0.2fc and d = 30 μm (w40 = −1.08) has one strong lobe, whereas the corresponding function has four dominant lobes. Comparison of intensity profiles taken along the red dashed horizontal line shown in Fig. 1(b), before and after application of the Wiener filter to the raw data (Fig. 1(g)) shows a 2.1 × improvement in the contrast of the illumination fringe pattern due to the application of the Wiener filter. To quantify the OS capability further, we determined the full width at half maximum (FWHM) of the axial profiles plotted along the red dashed vertical line in Fig. 1(c) [see Fig. 1(h)]. The FWHM value of the profile from the true object is 0.1 µm, whereas the FWHM values of the profiles from the restored images are 1 µm, 0.8 µm, and 0.2 µm for 2D-DMOD, 3D-DCV, and the proposed method, respectively. Note that, the peak of the intensity profile of the restored images in the case of 2D-DMOD is shifted due to the presence of depth-induced SA [Fig. 1(h)]. In the case of both the 3D-DCV and the proposed method, the original location of the object is properly restored as these two methods use an aberrated 3D PSF computed at 30 µm depth below the cover slip in the restoration process.
To further evaluate the OS capability achieved with the proposed method, we restored images of an object with two fluorescent layers separated by 0.5 µm, a distance lower than the axial resolution of the system [Fig. 1(i)]. As expected, the layers cannot be distinguished in the 2D-DMOD result while their separation is hinted in the 3D-DCV result. The improved performance obtained by the proposed method in resolving the two layers is evident in the restored image and further quantified by the axial intensity profile taken through the center of the XZ view image [Fig. 1(j)]. From the profiles in Fig. 1(j), we measure the drop in intensity between the maximum value and the value at the midpoint between the two restored fluorescent layers. For the 3D-DCV and the proposed method results, the drop in intensity is 11% and 70%, respectively. Following the Rayleigh criterion, the two layers are considered resolvable if the drop in intensity is higher than 25%. This result demonstrates the improved OS capability of the proposed method.
4.2 Impact of SA and modulation frequency on SIM restored images
The efficacy of the proposed method in the case when the fringe contrast is compromised is investigated through noisy simulation of SIM image of a test object [Fig. 2(a)] described in Section 3.1. Because the fringe contrast can be compromised by both the depth location of the object and the modulation frequency, restored images are shown when the test object is placed at different depths, d and modulated with different frequencies, fM. The XZ cross-sectional image of the true object, and the corresponding simulated forward raw SIM image for a single phase-shift (i. e. φ = 0 deg) at d = 30 µm and fM = 0.5 µm−1 are shown in Fig. 2(a). The XZ cross-sections from restored images obtained using different methods are shown in Fig. 2(b-d), where in each case the data was modulated using a different modulation frequency: (a) fM = 0.1fc; (b) fM = 0.2fc; and (c) fM = 0.3fc. The first, second and the third columns in each case of Fig. 2(b)-(d) show results obtained using a different method: (i) the 2D-DMOD; (ii) 3D-DCV; and (iii) the proposed method, respectively. All the methods can recover the object at the 0 µm depth (i.e. in the absence of SA), but the images restored using the 2D-DMOD method become noisy with increasing fM. In the presence of SA (at 30-µm and 60-µm depths), the SNR is reduced significantly in the images restored using the 2D-DMOD method for fM ≥ 1.0 µm−1 [Fig. 2(c & d)]. Although the 3D-DCV algorithm produces improved results compared to the 2D-DMOD method, the images restored using the proposed method show qualitatively the best match with the true object in the respective cases, where there is presence of depth-induced aberration (at 30-µm and 60-µm imaging depths).
For a quantitative comparison, we consider the results obtained at depth 30 µm and fM = 1.5 μm−1, where the theoretical modulation contrast is 0.07. The peak SNR in the restored image, the SNR reduction in the OS image with respect to the raw data, the NMSE  and the S-SIM  computed between the true object and the restored images are shown in Table 1. These metrics demonstrate the improved performance achieved with the proposed restoration method.
4.3 Experimental image results from a test sample
The proposed method is validated by restoring experimentally acquired images from the test objects described in Section 3.2. The XZ cross-sectional image from the raw SIM image of a 6-µm spherical shell located at a depth 3 µm below the coverslip is shown in Fig. 3(a). Figure 3(b) shows restored images using the 2D-DMOD, the 3D-DCV and the proposed method from the spherical shells located at different depths: 3 µm (top row), 24 µm (middle row), and 50 µm (bottom row). The modulation contrast in the raw data at increasing depth was found to be equal to 0.12, 0.10 and 0.09, which was increased to 0.20, 0.18, and 0.15, respectively, after the Wiener filter step. These observations suggest an average 1.71 × improvement in the modulation contrast due to the Wiener filter step, which consequently improved the SNR in the restored images. Images restored by the 3D-DCV method are similar to those obtained by the 2D-DMOD technique except that the SNR is improved by the 3D-DCV method, which is further improved by the proposed method. A quantitative comparison using the axial intensity profiles [Fig. 3(c)] taken through the center of the restored images from the spherical shells located at a 50-µm depth below the coverslip, shows the SNR improvement in the result obtained with the proposed method. The peak SNR values were found to be equal to 19.1 dB, 21.1 dB and 24.3 dB in the results obtained by the ZEN-OS, the ZEN-DCV, and the proposed method; respectively, the differences in these values is with the SNR values achieved in simulation (Table 1).
As mentioned earlier, the outer diameter of the spherical shell is 6 µm, and the shell thickness is 1 µm based on the manufacturer specification. We computed the FWHM distance of the outer shell along the axial and the lateral directions and the FWHM distance of the shell itself [from the intensity profile shown in Fig. 3(c)] for all the restored beads. The corresponding FWHM distances shown in Table 2 are within a 95% confidence interval, where we observe that the axial and lateral diameter are restored similarly by all the methods (with an error < 8%). However, the shell thickness in the axial direction is closer to the truth (with an error ~4%) in the case of the proposed method compared to the results obtained from the other two techniques (which show an error > 20%). These results demonstrate that all three methods perform well in restoring object dimensions that are much larger than the diffraction limit, while smaller dimensions (< 1 µm) are restored better by the proposed method. To quantify how isotropic our results are, we also computed the ratio between the axial and lateral values obtained for the shell diameter and thickness (Table 2). From these results, we can also conclude that our proposed method does not introduce any artifact while restoring the image.
4.4 Experimental results: biological specimen
To evaluate the performance of the proposed approach using biological data, we reconstructed the 3D image from a Glioblastoma cell and a pollen grain. Because the Glioblastoma cell has two axial planes close to each other (within 1 µm), and the pollen grain is a thick object (around 30 μm of thickness), these two samples are suitable to evaluate the OS performance and the effect of sample thickness, respectively. XY cross-sectional demodulated images of the Glioblastoma cell, and a zoomed region-of-interest (ROI) of the pollen grain image are shown in Fig. 4 (a-b), and (c-d), respectively. From these results, it is evident that application of the Wiener filtering prior to the demodulation [Eq. (8)] provides signal improvement. This signal increase (found to be 50% in the case of the pollen grain) is due to the increase of the fringe contrast by the Wiener filter [Fig. 4 (e) & (f)]. We observe that the modulation contrast increases within the range of 1.49 × to 1.77 × , in the Glioblastoma cell case, and in the range 1.6 × to 2.0 × in the pollen grain case, respectively.
4.4.1 Experimental validation of the achieved OS capability
The OS capability achieved with the proposed method in the case of experimental images is demonstrated in Fig. 5, where the restored image of the vimentin protein of a Glioblastoma cell (described in Section 3.2) is shown. As the vimentin surrounds the nucleus, this layer has a hole inside. The restored images at three different axial locations separated by 0.3 µm, and the corresponding XZ cross-sectional images are shown in Fig. 5(a-c). It is seen that in the case of the ZEN-OS result, each axial plane is affected by the out-of-focus light from other planes, which is reduced in the ZEN-DCV result and the result from the proposed method. In the zoomed ROI of the restored images at d0 + 0.3 μm [Fig. 5(i-iii)] it is clearly seen that the background intensity is minimized resulting in a higher contrast image in the case of the proposed method (panel iii). To visualize the axial confinement clearly, the zoomed versions of the rectangular-bound area of XZ cross-sectional images of the corresponding ROI [4th row of panels (a-c)] are shown in Fig. 5(iv-vi), which show the improvement in the axial intensity of the restored image obtained with the proposed method. The axial confinement of the two layers is quantified by the intensity profile [Fig. 5(d)] drawn through the dashed vertical line shown in Fig. 5(vi). It is observed that in the case of the ZEN-OS result, the two lobes are not separated; whereas, in the case of the ZEN-DCV result and the proposed method’s result, both lobes are visibly discriminated. The OS performance of the ZEN-DCV method and the proposed method is compared using the FWHM widths of the two corresponding lobes. The FWHM widths of the first and the second lobes in the case of ZEN-DCV result are found to be equal to 520 nm and 470 nm whereas, the corresponding lobes in the case of the proposed method’s result are found to be equal to 390 nm and 420 nm. These values demonstrate the improved axial confinement obtained with the proposed method. From Fig. 5(iv-vi), and its corresponding axial intensity profiles shown in Fig. 5(d), it is clearly seen that, the restored image from the proposed method appears to have fewer artifacts compared to the one from the ZEN-DCV method. The peak SNR computed from the entire 3D image is 19.9 dB, 22.3 dB, and 28.7 dB in the results from the ZEN-OS, ZEN-DCV, and the proposed method, respectively.
4.4.2 Experimental evaluation of thick sample restoration
The performance of the proposed method in restoring a thick sample is demonstrated in Fig. 6, where we show the restored images of a pollen grain in which the imaging thickness is approximately 30 µm. Figure 6 (a-c) shows the XY and XZ view images restored using the ZEN-OS, ZEN-DCV, and the proposed method, respectively. From a qualitative comparison, we observe that the signal in the restored images from the proposed method is improved compared to the results from the other two methods. As before, the signal improvement is quantified by the peak SNR, found to be equal to 15.3 dB, 17.6 dB, and 20.3 dB, in the 3D image from ZEN-OS, ZEN-DCV, and the proposed method, respectively. For better visualization of the comparison between the restored images from different methods, zoomed images [Fig. 6(i-iii)] of the ROI marked in the top row of Fig. 6(a-c), show the improved contrast achieved by the proposed method compared to the other two methods. Further zoom of a 5 × 5 µm2 ROI from Fig. 6(i-iii) shown in Fig. 6(iv-vi), respectively, clearly demonstrates the improved SNR and contrast achieved with the proposed method. The improvement in the image contrast is quantified by measuring the maximum and minimum values of the axial intensity profile [Fig. 6(d)] drawn along the arrow shown in Fig. 6(iii), where three point structures separated by approximately 0.5 µm are evident. The average contrast is found to be more than 3 × higher in the case of the proposed method’s result compared to the result from the other two methods.
5. Summary and conclusion
In this paper, an image restoration approach for ApoTome-based SIM is presented for recovering the image when the modulation contrast of the illumination fringes is poor at the image plane. In the first step of our proposed restoration method, we apply a 2D Wiener filter to the raw SIM data, which improves the modulation contrast [Fig. 1(g) and Fig. 4(e) & (f)] by removing the diffraction effect due to the 2D detection PSF, and then we demodulate the Wiener-filtered phase-shifted images. As evident from Fig. 4, the application of the Wiener filter provides a signal enhancement in the restored demodulated image. Finally, we use an EM algorithm (previously developed for widefield microscopy) using a synthetic SIM-PSF [Eq. (12)] to iteratively restore the image intensity.
For a better understanding of the impact of each step in the proposed restoration method, Fig. 7 compares a zoomed ROI from restored images of the pollen grain obtained using different steps of the proposed method. Comparison of Fig. 7(a) and (b) shows that the presence of the residual fringes (marked by the red arrows) is reduced by the IGF step, which can be used in any sinusoidal-based SIM restoration approach. The impact of the Wiener filter step is evident in the comparison of Figs. 7(a) and (c), where a signal improvement is obtained with the use of this step (the maximum intensity is doubled). As expected, the EM step provides signal and contrast improvement in the result as evident by the comparison between Figs. 7(a) and (d). Specimen features are less evident in Fig. 7(d) than they are in Fig. 7(a). The contrast improvement due to the EM step is also evident by the improved visibility of the fringes in Fig. 7(b) compared to Fig. 7(d). It is important to mention that the biggest benefit of the EM step is the correction of the SA through the use of an aberrant PSF. This effect was shown in Fig. 2 where the performance of the proposed method was evaluated in the presence of SA. It is noted that other iterative or non-iterative restoration algorithms can also be used in this step using the synthetic SIM-PSF.
The theoretical improvement in the modulation contrast and the OS capability of the system achieved with the use of our proposed method is demonstrated through the simulation of a thin fluorescent layer. Our results show a 2.1 × improvement in the modulation contrast after application of the Wiener filter step to the SIM raw data. The improvement in the OS performance of the proposed method compared to the 2D-DMOD and the non-iterative 3D-DCV methods is evident from the reduced axial FWHM distance of the restored thin layer and in resolving two thin layers separated by 0.5 µm. Restored noisy simulated images of a test object show that the proposed restoration approach can recover the object with improved accuracy even when the modulation contrast is as low as 0.07, which is demonstrated with an SNR improvement of 7.2 dB and 3.8 dB compared to the results obtained with the 2D-DMOD and the 3D-DCV methods, respectively. Similar SNR improvement was also observed in the experimental studies presented here.
The proposed method was experimentally validated by restoring images from a test sample with 6-µm in diameter spherical shells, and from two biological samples (one with Glioblastoma cells and the other with the pollen grains). Restored images of the test sample were consistent with simulations and validated the method. Results from the proposed method show approximately 16% less error in restoring the 1-µm shell thickness compared to the results obtained with the ZEN-OS and the ZEN-DCV methods.
In our experimental biological images from a Glioblastoma cell and a pollen grain we obtained a 1.63 × and a 1.8 × average improvement in the modulation contrast, respectively, after applying the Wiener filter step. Restored images of a Glioblastoma cell demonstrate improvement in both the axial confinement and the SNR in the result from our proposed method compared to the results obtained using the ZEN.2 software (both ZEN-OS and ZEN-DCV methods). Furthermore, the restored image of a pollen grain validates the ability of our proposed method to reveal fine details within a large structure with improved contrast compared to results from the other two methods. The restored images of the pollen grain are consistent with published images of similar pollen grains . As mentioned earlier, the approximation of Eq. (7) is based on the axial confinement of the illumination pattern. This approximation is especially true for a high NA ( = 1.4) objective lens, where the depth-of-field (DoF) is small (less than < 0.5 µm). For a low NA objective lens with a higher DoF, and for the case where there is no confinement in the illumination pattern, the 2D approximation may not be valid. In such cases, one should consider performing a 3D deconvolution with the 3D synthetic SIM-PSF .
In summary, we have demonstrated a SIM image restoration approach, which is suitable for imaging conditions in which the modulation contrast is reduced due to different experimental conditions. Our simulated and experimental results show that the proposed approach improves the achieved OS capability and SNR in the final restored image, compared to results obtained from existing 2D demodulation and non-iterative 3D deconvolution methods. Although the study was carried out with raw SIM images obtained from the ApoTome SIM system, the proposed framework could be used for other 3D-SIM imaging systems, in which the illumination is confined axially . However, our method would not probably be easily extended to 3D-SIM based on 3-wave interference  because the 3D illumination pattern is not axially confined.
National Science Foundation (IDBR award DBI-1353904, PI: CP); University of Memphis.
Authors thank O. Skalli and L. Thompson (Integrated Microscopy Center, The University of Memphis, TN) for providing the Glioblastoma cell sample and a biological interpretation of our results; G. Saavedra (Dept. of Optics, Univ. of Valencia, Spain) for useful technical discussion and L. H. Schaefer (Advanced Imaging Methodology Consultation, Ontario, Canada) for Matlab implementation of the 2D-DMOD and the 3D-DCV SIM methods.
The authors declare that there are no conflicts of interest related to this article.
References and links
1. 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(24), 1905–1907 (1997). [CrossRef] [PubMed]
4. M. G. L. Gustafsson, L. Shao, P. M. Carlton, C. J. 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(12), 4957–4970 (2008). [CrossRef] [PubMed]
7. B.-C. Chen, W. R. Legant, K. Wang, L. Shao, D. E. Milkie, M. W. Davidson, C. Janetopoulos, X. S. Wu, J. A. Hammer 3rd, Z. Liu, B. P. English, Y. Mimori-Kiyosue, D. P. Romero, A. T. Ritter, J. Lippincott-Schwartz, L. Fritz-Laylin, R. D. Mullins, D. M. Mitchell, J. N. Bembenek, A. C. Reymann, R. Böhme, S. W. Grill, J. T. Wang, G. Seydoux, U. S. Tulu, D. P. Kiehart, and E. Betzig, “Lattice light-sheet microscopy: imaging molecules to embryos at high spatiotemporal resolution,” Science 346(6208), 1257998 (2014). [CrossRef] [PubMed]
9. L. Gao, L. Shao, C. D. Higgins, J. S. Poulton, M. Peifer, M. W. Davidson, X. Wu, B. Goldstein, and E. Betzig, “Noninvasive imaging beyond the diffraction limit of 3D dynamics in thickly fluorescent specimens,” Cell 151(6), 1370–1385 (2012). [CrossRef] [PubMed]
10. A. G. York, S. H. Parekh, D. Dalle Nogare, R. S. Fischer, K. Temprine, M. Mione, A. B. Chitnis, C. A. Combs, and H. Shroff, “Resolution doubling in live, multicellular organisms via multifocal structured illumination microscopy,” Nat. Methods 9(7), 749–754 (2012). [CrossRef] [PubMed]
11. C. Preza, “Simulating structured-illumination microscopy in the presence of spherical aberrations,” Proc. SPIE 7904, 79040D (2011). [CrossRef]
12. C. Preza, L. Schaefer, D. Schuster, A. Ghaffar, S. Yuan, and G. Lobo, “Impact of spherical aberration on structured-illumination microscopy,” In Focus on Microscopy, Singapore, April 1–4, (2012).
13. L. Schäfer and D. Schuster, [Method and device for reconstructing images] US Patents 20060285632, (2011).
14. J. Conchello, “Image estimation for structured-illumination microscopy,” Proc. SPIE 5701, 1–8 (2005). [CrossRef]
16. A. Jost, E. Tolstik, P. Feldmann, K. Wicker, A. Sentenac, and R. Heintzmann, “Optical sectioning and high resolution in single-slice structured illumination microscopy by thick slice blind-SIM reconstruction,” PLoS One 10(7), e0132174 (2015). [CrossRef] [PubMed]
17. F. Ströhl and C. F. Kaminski, “A joint Richardson-Lucy deconvolution algorithm for the reconstruction of multifocal structured illumination microscopy data,” Methods Appl. Fluoresc. 3(1), 014002 (2015). [CrossRef] [PubMed]
18. S. Dong, P. Nanda, R. Shiradkar, K. Guo, and G. Zheng, “High-resolution fluorescence imaging via pattern-illuminated Fourier ptychography,” Opt. Express 22(17), 20856–20870 (2014). [CrossRef] [PubMed]
19. J. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1996).
20. J. A. Conchello, “Superresolution and convergence properties of the expectation-maximization algorithm for maximum-likelihood deconvolution of incoherent images,” J. Opt. Soc. Am. A 15(10), 2609–2619 (1998). [CrossRef] [PubMed]
21. L. H. Schaefer, D. Schuster, and J. Schaffer, “Structured illumination microscopy: artefact analysis and reduction utilizing a parameter optimization approach,” J. Microsc. 216(2), 165–174 (2004). [CrossRef] [PubMed]
22. N. Patwary, A. Doblas, and C. Preza, “Computational approach to address reduced modulation contrast in structured-illumination microscopy,” in Imaging and Applied Optics, OSA Technical Digest Series (Optical Society of America, 2017), paper JTu5A.9. [CrossRef]
23. N. Patwary, J. Sola-Pikabea, A. Doblas, G. Saavedra, and C. Preza, “Evaluation of the use of wavefront encoding to reduce depth-induced aberration in structured illumination microscopy,” Proc. SPIE 10499, 10499 (2018).
24. N. Patwary, H. Shabani, A. Doblas, G. Saavedra, and C. Preza, “Experimental validation of a customized phase mask designed to enable efficient computational optical sectioning microscopy through wavefront encoding,” Appl. Opt. 56(9), D14–D23 (2017). [CrossRef] [PubMed]
25. N. Wiener, Extrapolation, Interpolation, and Smoothing of Stationary Time Series (MIT Press, 1949).
26. T. Tkaczyk, M. Rahman, V. Mack, K. Sokolov, J. Rogers, R. Richards-Kortum, and M. Descour, “High resolution, molecular-specific, reflectance imaging in optically dense tissue phantoms with structured-illumination,” Opt. Express 12(16), 3745–3758 (2004). [CrossRef] [PubMed]
27. Q. Yang, L. Cao, H. Zhang, H. Zhang, and G. Jin, “Method of lateral image reconstruction in structured illumination microscopy with super resolution,” J. Innov. Opt. Health Sci. 9(3), 1630002 (2016). [CrossRef]
28. F. Orieux, E. Sepulveda, V. Loriette, B. Dubertret, and J. C. Olivo-Marin, “Bayesian estimation for optimized structured illumination microscopy,” IEEE Trans. Image Process. 21(2), 601–614 (2012). [CrossRef] [PubMed]
29. 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(2), 413–424 (2009). [CrossRef] [PubMed]
31. J. A. Conchello and J. G. McNally, “Fast regularization technique for expectation maximization algorithm for optical sectioning microscopy,” Proc. SPIE 2655, 199–208 (1996).
32. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. Image Process. 13(4), 600–612 (2004). [CrossRef] [PubMed]
33. S. Yuan and C. Preza, “Point-spread function engineering to reduce the impact of spherical aberration on 3D computational fluorescence microscopy imaging,” Opt. Express 19(23), 23298–23314 (2011). [CrossRef] [PubMed]
34. S. F. Gibson and F. Lanni, “Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy,” J. Opt. Soc. Am. A 9(1), 154–166 (1992). [CrossRef] [PubMed]
35. G. Saavedra, I. Escobar, R. Martínez-Cuenca, E. Sánchez-Ortiga, and M. Martínez-Corral, “Reduction of spherical-aberration impact in microscopy by wavefront coding,” Opt. Express 17(16), 13810–13818 (2009). [CrossRef] [PubMed]
36. M. Bertero, P. Boccacci, G. Desiderà, and G. Vicidomini, “Image deblurring with Poisson data: from cells to galaxies,” Inverse Probl. 25(12), 123006 (2009). [CrossRef]
37. [Computational Imaging Research Laboratory, Computational Optical Sectioning Microscopy Open Source (COSMOS) software package; http://cirl.memphis.edu/COSMOS].
38. N. Patwary, S. V. King, G. Saavedra, and C. Preza, “Reducing effects of aberration in 3D fluorescence imaging using wavefront coding with a radially symmetric phase mask,” Opt. Express 24(12), 12905–12921 (2016). [CrossRef] [PubMed]
39. M. Sivaguru, L. Mander, G. Fried, and S. W. Punyasena, “Capturing the surface texture and shape of pollen: a comparison of microscopy techniques,” PLoS One 7(6), e39129 (2012). [CrossRef] [PubMed]
40. H. Shabani, A. Doblas, G. Saavedra, and C. Preza, “3D structured illumination microscopy using an incoherent illumination system based on a Fresnel biprism,” Proc. SPIE 10499, 10499 (2018).