Image restoration approach to address reduced modulation contrast in structured illumination microscopy

: 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%.


Introduction
Structured illumination microscopy (SIM) provides improved optical sectioning (OS) [1] and super-resolution [2][3][4] compared to widefield microscopy; however, SIM performance is restricted in imaging optically-thick biological samples [5]. 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 [6]. Other approaches are based on either the combination of structured illumination with light-sheet microscopy [7][8][9], or the use of multiple focal spots as the excitation plane in a fluorescence microscope [10]. 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 [21] describes the intermediate raw image captured from a grid-projection SIM instrument [1]. 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][12][13]. In Ref [14] 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][16][17][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 [1]. Since light emits incoherently from a fluorescent sample, the modulated light in SIM is always subject to the incoherent optical transfer function (OTF) [19] 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, f M is equal to only 9% of the cut-off frequency of the imaging system, f c , 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 [20] 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 [22].
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 [23]. Because WFE causes the PSF to spread further [24], 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.

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: is the Dirac delta function; M f and i ϕ are the modulation frequency, and the phase of the illumination fringes, respectively, and 3 N ≥ is the number of phase-shifted patterns. The forward imaging model developed for ApoTome SIM assuming space-invariant imaging [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 [11]. 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 [12] in which SA was accounted using a single aberrant PSF. Specifically, we use the model in [11] 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: where 0 x and 1 x are the coordinates in the grating and sample space, respectively; ( ) o x is the 3D fluorescence intensity distribution of the object; ( ; ) h d x are the detection, and the illumination PSFs, respectively at a depth d below the cover slip. Without any loss of generality, 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 [17]. 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: where ( , ) u v are the lateral frequency coordinates. Note that the value of m c obtained from 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 [22], we demonstrated that for a 63 × /1.4 NA oil immersion lens, the theoretical fringes' visibility drops to m c < 10% at depths d > 20 µm for f M ≥ 0.3f C (see Fig. 1(a) in [22]).
As previously shown, Eq. (2) can be reformulated as the sum of the widefield and the modulated components, x , respectively [11,12]: where 3 ⊗ is the 3D convolution operator, is the synthetic SIM-PSF at a depth d below the coverslip, and

Inverse imaging problem solution
Because the synthetic SIM-PSF, , is axially-confined by the axial modulation function ( ; , ) M m z f d , 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 j th axial plane, , j z z = of the modulated sectioned component imaged at the best focal plane of the objective lens can be approximated by: 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 [21], and the unequal intensity in each phase-shifted images [21]. 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][28][29], and a parameter optimization method to minimize the mismatch between the theoretical and the experimental illumination patterns [21]. 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: where, ( ) is the 2D IGF, 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 [30], 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 ( ; , ) and An example of ( ; , ) M m z f d  is shown in Fig. 1(g) for a 63 × /1.4NA oil immersion objective lens, at d = 30 µm and f M = 0.2 f C .
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) [32] and the normalized mean square error (NMSE) [33] 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.

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

PSF computation
To compute the WF-PSFs we used the Gibson & Lanni optical path distance model [34] on a 128 × 128 × 256 grid with a voxel dimension of 0.05 µm 3 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 ( ; , ) M m z f d was calculated using Eq. (6). PSFs were computed for different depths below the coverslip.

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, 2 2 / (NA) 0.8 z λη ∆ = = µ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 w 40 = −1.08, which is computed using Eq. (4) in Ref [35]).
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 (f c ), 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 w 40 is, respectively, 0, −1.08, and −2.16). Noisy images were simulated as the realizations from a Poisson random process [36] 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 [24]. 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.

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 oilimmersion 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 f M = 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 µm 2 (equivalent to 0.1 × 0.1 µm 2 in the object space).

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 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 [37] 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 f M 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) [21], which demodulates the image plane by plane, and the non-iterative 3D deconvolution (3D-DCV) [13] 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 [38]. 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.

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

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

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 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) f M = 0.1f c ; (b) f M = 0.2f c ; and (c) f M = 0.3f c . 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 f M . 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 f M ≥ 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 f M = 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 [33] and the S-SIM [32] 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. Table 1. Comparison between the restored images shown in Fig. 2(b)-(d)

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.

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.

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

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

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. Fig. 7. Impact of different steps used in the proposed restoration method. XY cross-sections of a zoomed ROI from restored images of the pollen grain obtained with different steps of the proposed method: (a) using all the steps; (b) IGF step skipped; (c) Wiener filter step skipped; and (d) EM restoration and IGF step skipped. All panels have been displayed using the same color scale for direct comparison of signal strength. Residual fringes shown by a red arrow in (b) are not evident in (a) due to the IGF step. The dashed red rectangles highlight signal enhancement due to the Wiener filter and further signal and contrast improvement due to the EM restoration step.
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 [39]. 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-offield (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 [13].
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 [40]. However, our method would not probably be easily extended to 3D-SIM based on 3-wave interference [4] because the 3D illumination pattern is not axially confined.