1. Introduction

Computational optical sectioning microscopy (COSM) has facilitated 3D high resolution imaging from wide field (WF) fluorescence data [1]. Unlike confocal imaging, in COSM out of focus light intensity is accumulated, rather than discarded, and through processing it is reassigned to correct locations within the data volume resulting in images with improved resolution and contrast. COSM is thus preferable for low light, non-invasive in-vivo imaging. Traditionally, image restoration algorithms [2–6 ] developed for COSM were based on the assumption that the imaging system is space-invariant. This assumption makes the restoration process computationally feasible; however, it only holds when imaging thin (<5µm) samples. Imaging thicker samples requires focusing deeper into the sample, which often has an average refractive index (RI) that is different from the RI of the objective lens’ immersion medium. This RI mismatches causes depth induced spherical aberration (SA) [7] making the system depth-variant (DV).

An approximate DV forward model was developed by Preza and Conchello [8] where the object space is stratified into non-overlapping strata following the assumption that the PSF is invariant throughout each stratum. The strata-based maximum likelihood restoration problem was solved using a DV expectation-maximization (EM) algorithm (strata-EM) [8]. The strata-EM algorithm has been used successfully to restore the image of thick samples [9], but it can be computationally intensive if a large number of strata is used, and sometimes it can result in artifacts if a small number of strata is chosen [10]. The number of strata controls the tradeoff between computational complexity of the restoration algorithm and accuracy of the restoration, which both relate to the accuracy of the imaging model and the PSF representation. The development of a two-stage principal component analysis (PCA) method for the representation of 3D DV PSFs by orthonormal basis components [11] has improved the PSF accuracy. Arigovindan et al. [11] demonstrated that with a few number of principal components, PSFs over a desired depth can be computed with higher accuracy compared to the strata-based representation. In addition, they derived a PCA-based forward imaging model for 3D fluorescent microscopy; however, their study lacks experimental validation of the model. In this paper, the PCA DV forward image model is validated with a comparative study between experimental and simulated test-sample images. Yuan and Preza [12] proposed a 3D DV restoration algorithm based on the PCA imaging model, which showed promising preliminary simulation results presented in conference publications [12–14 ]. In this paper, PCA-based DV restoration is further investigated and applied to biological data.

The focus of this paper is an in-depth performance evaluation of a PCA-based expectation maximization algorithm (PCA-EM) for 3D image restoration. Toward this end, noisy simulated images (that predict experimental images) and experimentally acquired images, of a 3D test sample and of biological samples, were restored using the PCA-EM algorithm. Simulated DV-PSFs were used in the restoration of synthetic and experimental images. In addition, to facilitate application of the algorithm to experimental data from thick samples, we developed an approach that reduces the number of PSFs required for the computation of the PCA coefficients. Preliminary results of this approach have been presented in a conference publication [15]. Our implementation allows the computation of PCA components on a desktop computer and eliminates the need for a high performance computing facility for the computation of the PCA coefficients from a large set of 3D PSFs required in the case of thick specimen imaging.

The paper is organized as follows: Section 2 presents theoretical background for DV imaging using the strata- and PCA-based approaches. Section 3 describes methods used in the simulations, the experimental data acquisition, and the restoration process. Results from our investigations including restoration performance of the PCA-EM algorithm and a comparative study with the strata-EM algorithm, are presented in Section 4. Finally, Section 5 summarizes findings from this study.

2. Related work and background

In this Section, for completeness, we present the strata and PCA based forward imaging models [8, 11 ], which are the basis for the two DV restoration approaches presented in this paper. The proposed PCA-EM algorithm [12] is also described.

2.1 DV imaging model

In a discrete finite domain, the intensity in a 3D image formed by DV imaging can be described by:

2.1.1 Strata-based DV imaging model

The first approximation to Eq. (1) was based on the strata-based DV imaging model [8], in which the object is divided into M non-overlapping strata along the Z axis. The 3D DV PSF associated with a particular stratum is the weighted linear interpolation of two PSFs,$h_m(x,y,z)$and$h_{m+1}(x,y,z)$, computed at the depths, $Z_m$and$Z_{m+1}$, that define the top and the bottom of each stratum, respectively. The intensity in the 3D image can be approximated by:

2.1.2 PCA-based DV imaging model

In the PCA-based PSF representation model [11], a set of DV PSFs is expressed by the principal components (PCs) and coefficients computed from a two-stage tensor product-based PCA (TP-PCA) algorithm that facilitates the computations for a set of 3D images. In the TP-PCA algorithm, a 3D DV PSF due to a point source of light at depth$z_0$in the sample can be represented using the following equation:

Using the PCA-based PSF representation (Eq. (4)), Arigovindan et. al [11] developed a forward imaging model given by:

2.2 PCA-based solution to the inverse imaging problem

Using the PCA-based forward imaging model (Eq. (5)), Yuan and Preza proposed an image estimation algorithm [12], which estimates the object from the measured image using an iterative expectation-maximization algorithm (referred to as the PCA-EM algorithm). The PCA-EM algorithm is an extension of the well-known EM algorithm that solves the space-invariant maximum likelihood (ML) restoration problem for fluorescence microscopy [16]. In the space-invariant case, the estimated intensity at the (k + 1)^th iteration is computed using the current estimate$f^{(k)}(x,y,z)$of the sample and a back projection of an image ratio, computed from the measured image $d(x,y,z)$and the model prediction at the k^th iteration,$g^{(k)}(x,y,z)$, to the object space:

In the presence of noise, the EM algorithm without any regularization has been shown to produce isolated bright spots, which degrade the quality of the restored image [17]. The PCA-EM is regularized by adding a penalty functional,$R\left[f(x)\right]$, to the log-likelihood function,

3. Methods

This section describes the methodologies used in the investigation studies, which include: a method to compute PCA coefficients efficiently via an interpolation approach; DV PSF computation; computation of simulated objects and forward images; and PCA-EM restoration from simulated and measured data. Simulation studies were carried out to investigate and quantify the performance of the PCA-EM restoration method in various test cases.

3.1 Efficient implementation of the PCA computation

As mentioned above, the PCA coefficients$C(n,z)$used in the PCA-EM iteration (Eqs. (7) and (8) ), correspond to all discrete axial locations (planes) within the sample. In the case of a thick specimen, the number of axial planes in the 3D image of the specimen increases, which causes the number of coefficients needed for the restoration to increase. Consequently, because each $C(n,z_o)$ is associated with the DV PSF at depth $z_o$ (Eq. (4)), the number of DV PSFs needed for the PCA computation also increases. Taking advantage of the fact that the DV PSF does not vary significantly in a small depth interval, we have developed a PCA computation technique based on interpolated coefficients (referred to as PCA-IC) that reduces the number of PSFs required to compute the PCs and coefficients needed for DV restoration. PCA-IC is computationally efficient and suitable for the use on a desktop computer.

Our approach is based on the observation that the PCA coefficients vary smoothly as a function of depth (Fig. 1(a)-1(c) ). This allows computation of fewer coefficients from a smaller number of DV PSFs for 3D restoration of a thick sample. Then the coefficients corresponding to every axial location depth needed for the restoration (Eqs. (7) and (8) ) are computed using the cubic spline interpolation technique implemented in MATLAB (MathWorks, Inc.).

 

Fig. 1 The PCA coefficients of the first three components as a function of depth (axial location$z$) for a 20X/0.8 NA dry lens in the presence of spherical aberration. (a), (b), and (c) are the coefficients of the first, second, and third components, respectively. (d) Normalized mean square error (NMSE_3D) between the true PSFs and PSFs computed from the PCA-IC technique, and the full PCA computation as a function of significant PCs.

Download Full Size | PDF

To demonstrate the accuracy of the interpolation we compared the coefficients computed with the PCA-IC method from a set of 40 PSFs to PCA coefficients computed from a set of 200 PSFs without interpolation (i.e., computed with the full PCA computation). For this study, 200 PSFs for a 20X/0.8 NA dry lens were computed at an axial (depth) interval of 0.15 µm over a depth range of 30 µm in a watery sample. In the PCA-IC implementation, we used only 40 PSFs at a depth interval of 0.75 µm. Figures 1(a)-1(c) show the PCA coefficients corresponding to the first three PCs computed from the set of 40 PSFs (circles), and the interpolated coefficients (solid line). The mean square error (Eq. (14)) between the true PSFs and PCA-based PSFs represented using: (1) the interpolated coefficients and PCs computed with PCA-IC; and (2) the coefficients and PCs from the full PCA computation was computed as a function of the number of PCs used in the PCA representation of the PSF. Comparison of the normalized mean square error (NMSE) obtained in each case shows excellent agreement (Fig. 1(d)). The average difference in the NMSE computed for the case where only the three first significant PCs were used in the PCA representation of the PSF is only 0.05%, which is not significant and does not affect the PCA-EM restoration result. In Section 4.4 we show results that quantify the impact of this error on the accuracy of the restoration result.

3.2 Simulation studies

For the simulation studies, PSFs were generated using MATLAB (MathWorks, Inc.) code. The PCA components and coefficients, as well as the forward simulated images, were computed using the CosmTools module of the COSMOS open source software package [19]. The PCA-EM restoration algorithm investigated in this paper and the strata-EM algorithm used in the comparative study (Section 4.2) are both implemented in the CosmEstimation module of COSMOS [19], which was used for all the restorations presented in this study. The Linux version of the COSMOS software, available on the high computing facility (HPC) at the University of Memphis, facilitated completion of time-consuming simulations (Fig. 2 ).

 

Fig. 2 Determination of the OPR for: (a) SNR = 11.2 dB; (b) SNR = 10 dB; (c) SNR = 6.8 dB. The ORP is the value of γ for which the NMSE_3D between the true object and the restored image becomes the minimum. The ORPs for these three cases were found to be equal to 0.00004, 0.0001, and 0.003, respectively.

Download Full Size | PDF

3.2.1 Validation of PCA-based forward model

For the PCA model evaluation, a simulated image of a physical test sample was generated using a numerical object and the imaging conditions that matched those used in the data acquisition of a 6-µm in diameter fluorescent bead sample described in Section 3.3.1. The numerical spherical shell (referred to as object 1, shown in Fig. 3(a) ) was generated on a 256 × 256 × 512 grid, and centered at the grid point (128,128,256), where each voxel was mapped to a physical dimension of 0.1 µm, assuming that its center is at a depth equal to 65 µm below the coverslip. The simulated 3D image was generated using the PCA depth-variant imaging model (Eq. (5)) implemented in the variant tab of the CosmTools using 60 DV PSFs. The DV PSFs (based on the Gibson and Lanni PSF model [20] and the vectorial field approximation theory [21]) were computed at an interval of 0.1 µm within the depth range of 62.5 µm to 67.5 µm, which includes the location of the physical object below the coverslip. The PSFs for a 63X/1.4 NA oil-immersion objective lens were computed on a 128 × 128 × 128 grid with a voxel dimension equal to 0.1 µm, assuming that the specimen embedding medium is glycerol (RI = 1.47). Comparisons between the model prediction of the numerical object and the experimental data are discussed in Section 4.1.

 

Fig. 3 Evaluation of the PCA forward imaging model. XZ medial section of: (a) the numerical spherical shell (object 1 as described in Section 3.2.1); (b) the simulated image of (a) obtained by the PCA model (described in Section 3.2.1); (c) experimentally acquired image of a 6-μm bead test sample (described in Section 3.3.1). (d) Axial normalized intensity profile through the center of the images shown in (b) and (c). To plot the intensity profile, both simulated and experimental images were normalized, and were registered properly to match the peak intensities. Lens: 63X/oil immersion; Emission Wavelength: 515 nm.

Download Full Size | PDF

3.2.2 Comparison between two DV-restoration methods: the PCA-EM and the strata-EM

To demonstrate the advantages of the PCA-EM over the strata-EM restoration algorithm, a second simulated object (referred to as object 2, shown in Fig. 4(a) ) was used. Object 2 is also a spherical shell with outer diameter equal to 6 µm and shell thickness equal to 1 µm, embedded in water (RI = 1.33) and with its center located at 5 µm below the coverslip, and simulated in a 25.6µm x 25.6 µm x 10 µm volume. It was simulated at a shallower depth than object 1 in order to investigate the effect of the rapid PSF changes due to increasing SA in this region. The image of the object was computed using the depth-variant imaging model (Eq. (2) and 100 DV PSFs computed at an interval of 0.1 µm within the depth range of 0 µm to 10 µm.

 

Fig. 4 Algorithm performance comparison of the PCA-EM and strata-EM algorithms using a simulated image from a 6 µm in diameter spherical shell (object 2) shown in (a). (b)-(d) the PCA-EM restored images using 2, 4, and 6 PCs, respectively; (e)-(g) the strata-EM restored images using 2, 4, and 6 strata, respectively. (h) Algorithm performance quantified by the NMSE_3D (Eq. (14)) as a function of the number of PCs/strata used in the restoration. (f) Time needed to compute the restoration result after 1000 iterations of each algorithm as a function of the number of PCs/strata used in each case. Lens: 63X/1.4 oil immersion; Emission Wavelength: 515 nm.

Download Full Size | PDF

The simulated image was restored using two, four, and six principal components (PCs) in the case of the PCA-EM algorithm, and a comparable number of strata (determined using the relation between PCs and strata in [13]) in the case of the strata-EM algorithm. The restored images obtained after 1000 iterations of both algorithms are displayed in Fig. 4 and discussed in Section 4.2.

3.2.3 Noisy image restoration

The performance of the regularized PCA-EM was evaluated in simulation using a numerical test object with two spherical shells (referred to as object 3, shown in Fig. 5(a) ) generated using the same imaging condition used for object 2. The physical dimension of object 3 is 12.8 µm × 12.8 µm × 25.6 µm, with the centers of the two shells located at the (x,y,z) points (4.8 µm, 6.4 µm, 11 µm) and (4.8 µm, 6.4 µm, 14 µm), respectively. The outer diameter of the shells is 4 µm, and the shell thickness is 1 µm. The DV image of object 3 was computed using 256 DV PSFs defined at depths within the range of 0 µm to 25.6 µm separated by a 0.1 µm interval. The simulated DV image (Fig. 5(b)) was corrupted with Poisson and Gaussian noise to simulate the actual imaging and acquisition system [22]. Three different noise levels with overall SNR (Eq. (13)) equal to 11.2 dB, 10 dB and 6.8 dB were used in this algorithm performance investigation.

 

Fig. 5 Performance evaluation of the regularized PCA-EM algorithm demonstrated by restoring noisy simulated images using PCA-EM (with 10 PCs) and a roughness regularization penalty (Eq. (11)). XZ medial section of: (a) the synthetic object 3; and (b) the noisy simulated image with SNR 6.8 db. PCA-EM restored images computed from noisy images using an ORP: (c) γ = 0.0004 for SNR = 11.2 db; (d) γ = 0.0001 for SNR = 10 db; and (e) γ = 0.003 for SNR = 6.8 db. Specimen embedding medium is water. Lens: 63X/1.4 NA oil immersion. Emission wavelength: 515 nm.

Download Full Size | PDF

The noisy simulated images were restored using 10 PCs and the corresponding optimal regularization parameter (ORP). More PCs were used in the case of noisy image restorations compared to the noiseless studies, in order to study restoration accuracy in terms of the SNR of the noisy images. To determine the ORP, PCA-EM restorations were obtained using different values of the regularization parameter and compared in terms of the NMSE_3D between the true object and its restorations. The regularization parameter (γ) corresponding to the best restored image judged by the minimum NMSE_3D value was considered as the ORP (Fig. 2). The ORP was found to be equal to 0.00004, 0.0001, and 0.003 for SNR equal to 11.2 dB, 10 dB and 6.8 dB, respectively [23]. Restored images after 10,000 iterations of the regularized PCA-EM algorithm are displayed in Fig. 5 and discussed in Section 4.3.

3.2.4 DV restoration using a reduced set of DV-PSFs

In order to test the performance of the efficient implementation of the PCA computation using the proposed interpolated coefficient technique (discussed in Section 3.1), a thick test object was chosen (referred to as object 4, shown in Fig. 6(a) ), thereby increasing the requirements on the number of PSFs needed. Object 4 is 30 µm thick and has five spheres (3 µm in diameter) centered at depths 5 µm, 10 µm, 15 µm, 20 µm, and 25 µm, respectively. The whole volume was generated on a 336 × 336 × 336 grid. We simulated a 20X/0.8 NA dry objective lens assuming that the micro spheres are embedded in water. 200 DV PSFs were calculated in the depth range of 0 µm to 30 µm, using the same methods described above except that the voxel dimension was 0.15 µm³. These PSFs were used to compute the forward simulated image using Eq. (2). We used a lower NA lens for this study so that the number of PSFs required for adequate axial sampling was small enough so that computation of the PCA without interpolation was possible.

 

Fig. 6 Performance comparison between the PCA-EM restorations (after 600 iterations) from simulated data using 3 PCs computed with and without interpolated coefficients (IC). XZ sectional view from: (a) the numerical object; (b) the simulated image of the numerical object; (c) the PCA-EM restored image using PCA data from a full set of PSFs; and (d) the PCA-EM restored image using PCA data from a reduced set of PSFs computed with PCA-IC. A 20X/0.8 NA air objective lens was simulated assuming that the specimen is embedded in water.

Download Full Size | PDF

The simulated image was restored using the PCs obtained from the full PCA computation (where all the PSFs were used) and the PCA-IC (where a reduced set of 40 PSFs was used). For the full PCA computation, PCs were calculated from all the 200 PSFs computed at a depth interval of 0.15 µm. This PCA computation was performed on the HPC facility to accommodate the large memory requirement. On the other hand, in the case of the PCA-IC, PCs were calculated from a set of 40 PSFs computed at an interval of 0.75 µm using a desktop computer. Results from this study are discussed in Section 4.4.

3.3 Application of the PCA-EM to measured data

3.3.1 Sample preparation and data acquisition

Images from a test sample and two different biological samples were captured for this study. The test object was a slide with 6 µm in diameter fluorescent microspheres (FocalCheck, Invitrogen, Molecular Probes) with dual staining: blue throughout label, and green label on a 1 µm thick outer shell of the microsphere. The microspheres were dried on the microscope slide and glycerol (RI = 1.47) was added on top. A layer of 170 nm green fluorescent spheres (FluoSphere, Invitrogen, Molecular Probes) were dried on the coverslip (to mark its location) which was sealed onto the slide. The specimen was illuminated using a mercury-arc source and an excitation filter with a 450-460 nm pass band. The 3D image was captured using a ZEISS Imager.Z1 with a 63X/1.4 NA oil (RI = 1.518) immersion objective lens using the green emission fluorescent filter (with a 515-565 nm pass band), thus the microspheres appear as spherical shells. The sample was scanned axially using an interval equal to 0.1 µm, and the image was recorded digitally by a CCD camera (AxioCam MRm) with a pixel size equal to 6.45 µm × 6.45 µm (equivalent to 0.102 µm × 0.102 µm/pixel in the object space). The distance between the layer of the 175-nm microspheres at the coverslip and a spherical shell of interest allowed us to obtain the depth of the shell below the coverslip. This depth of the center of the shell was found to be equal to 65 µm.

Two biological samples were used in this study: (1) a slide with lung epithelial cells, and (2) a slide with cultured Glioblastoma cells. The image acquired from the first sample was used to demonstrate the optical sectioning capability of the PCA-EM. Images acquired from the second sample were used in a quantitative comparison of the resolution achieved with two different 3D imaging techniques: COSM and structured illumination microscopy (SIM). SIM was also used to evaluate the optical sectioning capability. To obtain SIM images, the ZEISS ApoTome.1 device was inserted in the illumination path and images were captured at three different illumination angles (0°, 120°, and 240°). SIM optical sectioned images were computed using the standard method available with the ApoTome.1 and deconvolved using regularization parameter γ = 0.0001. For the 20X/0.8 NA dry lens we used L1 grating with grating frequency equal to17.5 lines/mm while for the 63X/1.4 NA oil lens we used the phase H1 grating with frequency equal to 7.5 lines/mm.

Immortalized mouse lung epithelial cells (MLE-12) were stained using the actin cytoskeleton and Focal Adhesion Staining Kit (EMD Millipore, catalog no. FAK100) and the SlowFade (Life Sciences) mounting protocols. The image was captured using the same microscope described above with a 20X/0.8 NA dry lens, and an axial scanning interval equal to 0.32 µm.

In the cultured Glioblastoma cells, the actin was labeled with green phalloidin and the permeabilized membrane was infiltrated with Alexa Fluor. The cells were cultured with green fluorescent nano-spheres (170 nm in diameter), which were eventually engulfed by the cells. The cells with the nano-spheres were finally fixed and mounted in Prolong Gold (RI = 1.46). The WF image of the sample was recorded with a 63X/1.4 NA oil immersion objective lens using an axial scanning interval of 0.1 µm. Another image of the same field of view of this sample was captured using SIM.

3.3.2 Image restoration from measured data

The measured image acquired from the test sample (Fig. 7(a) ) was restored using the PCA-EM algorithm and PCs computed from: 1) a set of 200 DV PSFs at a depth interval of 0.1 µm; and 2) from a reduced set of 40 DV PSFs at a depth interval of 1 µm using the PCA-IC approach. The PSFs were computed for a 63X/1.4 NA oil immersion objective lens using a RI = 1.47 for the specimen embedding medium. The images from a mouse lung epithelial cell (Fig. 8 ) and a cultured Glioblastoma cell (Fig. 9 ) were restored using PCs computed from the reduced set of 40 DV PSFs at a depth interval of 1 µm for the same lens, but in these cases the RI for the specimen embedding was equal to 1.33 and 1.46, respectively.

 

Fig. 7 Restoration performance evaluation of the PCA-EM algorithm using 3 PCs. (a) XZ medial section of the measured image of a 6 µm in diameter spherical shell with a shell thickness equal to 1 µm. Restoration using PCs computed from the set of: (b) 200PSFs; and (c) 40 PSFs. The restorations were carried out with a regularization parameter equal to 0.0001. The restored images after 1000 iterations are displayed. Images were resampled in the axial direction to compensate for the axial shift due to SA (Eq. (12)).

Download Full Size | PDF

 

Fig. 8 Evaluation of the optical sectioning capability of the PCA-EM algorithm. (a)-(d) are four XY sections (planes at different axial locations) from the WF image of mouse lung epithelial cells captured using a 20X/0.8 NA dry lens; (e)-(h) are corresponding planes captured using the Zeiss ApoTome attachment (SIM); (i)-(l) are corresponding planes from the PCA-EM restored image computed with 3 PCs and a regularization parameter equal to 0.0001 after 1000 iterations. Image intensities are displayed on a separate scale for the best visualization of features in each image.

Download Full Size | PDF

 

Fig. 9 Comparison of PCA-EM restoration with results obtained from other optical sectioning methods. XY sectional image from a Glioblastoma cell acquired with two different modalities and restored using different methods: (a) WF unprocessed image with three marker beads enclosed by a rectangle, a circle, and a triangle; (b) Deconvolved SIM image captured using ApoTome; (c) Restoration with ZEN iterative deconvolution algorithm; (d) Restoration with PCA-EM algorithm with 3 PCs. In the case of iterative restoration, results were obtained after 1000 iterations using a regularization parameter equal to 0.0001. Image intensities are displayed on a separate scale to facilitate visualization of features in each image. Lens: 63x/1.4 NA oil; Wavelength: 509 nm.

Download Full Size | PDF

The experimental data were also resampled to compensate for the axial shift due to the RI mismatch between the objective immersion medium and the specimen embedding medium according to the following equation [24]:

Restorations from the experimental data were obtained using 3 PCs in all cases. The signal to noise ratio (SNR) of the experimental data in the case of spherical shell object was found to be approximately equal to 10 dB using Eq. (13). Thus, in the restoration the regularization parameter was set equal to 0.0001 (determined as in Section 3.2.3). Restored images after 1000 iterations are shown and discussed in Section 4.5.

PCA-EM restoration obtained from the 3D Glioblastoma cell image and the PCA-IC result for the PCs and coefficients, were compared to optical sectioned images obtained with other restoration methods. Towards this end, the Glioblastoma cell was also processed using a constrained iterative algorithm (with the regularization parameter equal to 0.0001) available in the ZEN 2012, blue edition commercial image restoration software (Carl Zeiss Microscopy GmbH).

3.4 Restoration metrics

3.4.1 Signal to noise ratio (SNR)

To quantify the level of noise in the images used in noisy simulations and in restoration of experimental data, we computed the SNR using the following equation:

3.4.2 3D normalized mean square error (NMSE_3D)

To quantify the difference between two images$a(x,y,z)$and$b(x,y,z)$we used a 3D normalized mean square error (NMSE) which is calculated according to the following equation [25]:

4. Results and discussion

Results from several studies are presented in this section. Validation of the PCA-based forward model and results obtained from a comparative study between the PCA-EM and the strata-EM algorithms are presented in Section 4.1 and Section 4.2, respectively. The robustness of the PCA-EM algorithm to noise is investigated in Section 4.3. Section 4.4 shows results that investigate the impact of the interpolated coefficients on the accuracy of the PCA-EM. Image restorations from experimentally acquired images are shown in Section 4.5.

For qualitative comparison of images, only XZ medial sections are shown, unless otherwise specified, because XY sections do not provide any additional information while using circularly symmetric test objects. For the case of biological samples, XY sections are shown to help visualize useful structures in the cells. Quantitative comparisons are shown using either intensity profiles taken through the XZ medial sections or a cumulative normalized mean square error (NMSE_3D) computed over the 3D image (Eq. (14)). Displayed images are cropped to show the region of interest rather than the restoration grid size, which is often larger in order to minimize the effects due to the discrete finite representation of the data. Scale bars shown in figures apply to all the images in the corresponding figure.

4.1 Experimental validation of PCA forward imaging model

Evaluation of the PCA forward model is shown in Fig. 3 where a model prediction for the image of object 1 (Fig. 3(b)), computed as described in Section 3.2.1 is compared to the experimentally acquired image (Fig. 3(c)). The XZ sectional images show a qualitative agreement while the intensity profiles (Fig. 3(d)) quantify the difference between the simulation and observation. The NMSE_3D between the normalized experimental image intensity and the simulated image intensity is 0.0623, which sufficiently demonstrates the accuracy of the PCA-based forward imaging model using three PCs.

Comparison between PCA-EM and strata-EM restoration in noiseless simulation Fig. 4. It is observed qualitatively that in the case of the PCA-EM algorithm, two PCs provide a restoration of the spherical shell (Fig. 4(b)) that may be reasonable for some applications. In the case of the strata-EM algorithm, at least six strata (Fig. 4(g)) are required to obtain similar result. Figure 4(h) compares the NMSE_3D values with respect to the number of PCs (for the PCA-EM algorithm) or strata (for the strata-EM algorithm). As it can be observed (Fig. 4(h)), a higher restoration accuracy is obtained with the PCA-EM algorithm compared to the strata-EM algorithm, at a comparable computational cost defined in terms of PCs and strata. In the case of 4 PCs or strata, the NMSE_3D is found to be equal to 0.022 and 0.044, respectively, which represents a 50% improvement in the restoration in terms of the NMSE_3D. The number of PCs required for an accurate restoration with the PCA-EM is clearly evident from the NMSE_3D curve.

The performance comparison between the PCA-EM and the strata-EM algorithm is also shown in terms of the execution time in hours required for the completion of 1000 iterations (Fig. 4(i)). It is observed that the PCA-EM algorithm with two PCs requires ~5 hours to perform 1000 iterations (for a 256x256x400 grid size on a 2.3GHz Quad Core processor with 4 GB RAM). On the other hand, in the case of strata-EM algorithm with two strata, ~15 hours are needed for the restoration using the same computer system as for the PCA-EM algorithm. The execution time increases linearly at a rate of 5.61 hours/stratum in the case of the strata-EM algorithm (Fig. 4(i)); whereas, in the case of PCA-EM the rate of increase is 2 hours/PC, which represents a 64% reduction in the execution time.

4.2 PCA-EM performance evaluation using noisy simulated data

Restorations from the noisy simulated images obtained as described in Section 3.2.3 are shown in Fig. 5 in order to evaluate the performance of the regularized PCA-EM algorithm. It is evident from Fig. 5 that appropriate levels of regularization result in desirable restorations. The effect of noise on the restored image at low SNR (6.8 dB) can be mitigated further using a higher value of the regularization parameter, but at the cost of reduced image resolution as predicted by Fig. 2(c) and [26]. The restoration accuracy is quantified in terms of the NMSE_3D (computed between the restored and true object), which was found to be equal to 0.0139, 0.0175 and 0.0289 in the cases of SNR values equal to 11.2 dB, 10 dB and 6.8 dB, respectively.

4.3 DV restoration using a reduced number of DV-PSFs

Results from the study evaluating the impact of computing PCA coefficients from a reduced number of PSFs are shown in Fig. 6. Restored images obtained by processing the DV simulated image (Fig. 6(b)) using the PCA-EM and PCA coefficients and components from the full and reduced set of PSFs (as described in Section 3.2.4) are shown in Fig. 6(c) and Fig. 6(d), respectively. The NMSE_3D between the two restorations is 1.03 × 10⁻⁴, which suggests that the approximation in the values of the PCs and coefficients due to the use of the PCA-IC approach has an insignificant effect on the PCA-EM image restoration. In addition, it is demonstrated here that for a 20X/0.8 NA dry lens used to image a 30 µm thick specimen in water, three PCs are sufficient for adequate restoration.

4.4 PCA-EM performance evaluation using experimental data

4.4.1 Restoration from test-sample data

The performance of the PCA-EM algorithm was evaluated using experimental data by restoring the image of a test sample (described in Section 3.3.1). The sectional view of the measured image along the XZ plane is shown in Fig. 7(a). PCA-EM restorations using PCs and coefficients computed from a set of 200 and 40 DV PSFs (with PCA-IC) are shown in Fig. 7(b) and Fig. 7(c), respectively. From these restored images it can be observed that PCA-EM restores the spherical shell considerably well. The FWHM diameters along the axial and lateral directions of the restored images (measured from the intensity profiles in Fig. 7(d)) are 5.8 µm and 5.9 µm, respectively, whereas the true diameter of the spherical shell is 5.9 µm (per vendor specification).The axial intensity profile shown in Fig. 7(d) shows quantitatively that computing PCA components from a reduced PSF set does not affect appreciably the restoration result, which is consistent with the results observed in simulation (Section 4.4).

4.4.2 Restoration from biological-sample data

Results from applying the PCA-EM algorithm to biological data are discussed in this section. The absence of ground truth in this case makes algorithm performance evaluation more challenging. The use of other optical sectioning methods facilitated the evaluation of the PCA-EM restorations. In addition, our samples offered unique assessment opportunities, which we present here.

To demonstrate the optical sectioning capability of the PCA-EM algorithm restored images from mouse lung epithelial cells (described in Section 3.3.1) are shown in Fig. 8. Four wide field (WF) images from consecutive axial planes separated by 1.92 µm are shown in Fig. 8(a)-8(d). The corresponding planes from SIM and PCA-EM restored images are shown in Fig. 8(e)-8(h) and Fig. 8(i)-8(l), respectively. As expected, the SIM image set exhibits better optical sectioning compared to the set of WF images. The PCA-EM result shows consistency with the SIM result, which is used here to validate the PCA-EM result. Two different types of cells are observed in these XY section images, which are better distinguished in the PCA-EM restoration. The observed differences between the images shown in Fig. 8(i) and Fig. 8(j) suggest that PCA-EM can remove out-of-focus light from adjacent planes thereby achieving the desired optical sectioning.

To quantify achieved resolution in restorations from biological data, the PCA-EM algorithm was applied to images of a Glioblastoma cell cultured with 170 nm in diameter fluorescent nanoparticles (as described in Section 3.3.1). Results from this study are summarized in Fig. 9. A single XY section or axial plane (where most of the structures are visible) taken from the WF unprocessed image (Fig. 9(a)) and from restorations obtained using the commercial deconvolution software package ZEN (Fig. 9(c)) and the PCA-EM algorithm (Fig. 9(d)) are compared to the corresponding deconvolved SIM image (Fig. 9(b)).

Restoration accuracy is quantified using the images of three beads identified inside the cell in all the images. In Fig. 9(a), bead #1, bead #2, and bead #3 are denoted by a rectangle, a circle, and a triangle, respectively. Normalized axial intensity profiles of the 3 beads from all the restored images are shown in Fig. 10(a)-10(c) , and they are compared with the intensity profile from an ideal PSF (computed at 0-µm depth below the coverslip without aberration). It is observed that in all 3 bead cases, the intensity profile from the PCA-EM restoration has the highest peak intensity and it is the narrowest compared to the other profiles, indicating that a higher axial resolution is achieved than in the images from the other methods. A lateral intensity profile through bead #2 is shown in Fig. 10(d). The FWHM diameter of bead #2 along the lateral direction is 0.17 µm, which agrees with the actual diameter of the nanoparticle. However, the FWHM diameter along the axial direction is 0.24 µm, which exhibits that the achieved resolution in the PCA-EM result is not isotropic, albeit improved.

 

Fig. 10 Quantitative comparison between the processed images of 175-nm in diameter microspheres (PSF beads) cultured in the Glioblastoma cell. (a-c) Axial profile through the three beads marked in Fig. 9(a), which demonstrate that the PCA-EM restored images have better contrast and axial resolution compared to SIM images and ZEN (ZEISS deconvolution software) restored images; (d) Lateral intensity profile through bead#2 for all the cases. Legend shown in (a) corresponds to all the figures. Lens: 63x/1.4 NA oil; Wavelength: 509 nm.

Download Full Size | PDF

5. Discussion and conclusion

The PCA-EM restoration algorithm (Eqs. (7)-(9) ) presented in this paper is based on a depth-variant (DV) imaging model (Eq. (5)) that uses an efficient PCA representation of DV PSFs defined over the imaging depth within the underlying sample. The required principal components and corresponding coefficients are pre-computed from a reduced set of DV PSFs and then a subset of these are used in the restoration. We have shown here that using a small number (3-6) of principal components the PCA-EM restoration algorithm renders faster and more accurate image estimation in a DV imaging framework than the previously developed strata-EM restoration algorithm. When the PCA-EM was executed using comparable resources with the strata-EM algorithm (i.e., number of principal components vs. number of strata), an improvement in accuracy by 50% along with a 64% reduction in the execution time was found.

Noisy image restoration with proper regularization based on a smoothness penalty demonstrated the robustness of the algorithm when applied to low SNR data. An efficient method of computing PCA coefficients from a reduced set of DV PSFs has made the PCA approach computationally tractable and useful for processing images from thick samples without loss of accuracy.

Restoration of an experimentally acquired image of a test sample (6 µm in diameter spherical shell embedded in a mounting medium) established the applicability of the algorithm to experimental data. In the case of this test sample, we showed through comparison of a simulated image with the experimental image of the spherical shell that the PCA-based forward imaging model used by the PCA-EM algorithm is accurate sufficiently as it captures the main features in the data with a 94% match computed based on a normalized mean square metric. Restored images from mouse lung epithelial cell and Glioblastoma cell with engulfed nanoparticles demonstrated that COSM with PCA-EM restoration can outperform other imaging methods compared in this paper, both in terms of optical sectioning capability (Fig. 8) and image resolution (Fig. 10).

The DV restoration algorithm presented in this paper provides computationally feasible and efficient fluorescence imaging in the presence of depth-varying aberration making it suitable for practical use in the investigation of thick samples of interest.