## Abstract

A depth-variant (DV) image restoration algorithm for wide field fluorescence microscopy, using an orthonormal basis decomposition of DV point-spread functions (PSFs), is investigated in this study. The efficient PSF representation is based on a previously developed principal component analysis (PCA), which is computationally intensive. We present an approach developed to reduce the number of DV PSFs required for the PCA computation, thereby making the PCA-based approach computationally tractable for thick samples. Restoration results from both synthetic and experimental images show consistency and that the proposed algorithm addresses efficiently depth-induced aberration using a small number of principal components. Comparison of the PCA-based algorithm with a previously-developed strata-based DV restoration algorithm demonstrates that the proposed method improves performance by 50% in terms of accuracy and simultaneously reduces the processing time by 64% using comparable computational resources.

© 2015 Optical Society of America

## 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:

*m*

^{t}^{h}stratum interpolation weights given by [8]:

### 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:

*n*base component or principal component (PC),$C(n,{z}_{o})$are the corresponding PCA coefficients for depth ${z}_{o}$, and

^{th}*B*is the number of base components. Theoretically, the larger the number of components the better the PSF representation is. In practice, a few components can reasonably represent all the PSFs in a large range of depth thereby reducing the dimensionality of the data set [11].

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,

where $\gamma $, the regularization parameter, is a constant that scales the penalty functional. In this study, we have used the Good’s Roughness penalty functional defined by [18]:## 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.).

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

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

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

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.

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^{3}. 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.

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.

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]:

*NA*is the objective lens’ numerical aperture;

*n*is the RI of the immersion medium and

*n*is the RI of the specimen embedding medium

_{s}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:

where${\sigma}_{s}$and${\sigma}_{N}$are the unbiased estimates of the standard deviation of the noisy image and the noise (difference between the noisy and noiseless image), respectively. The SNR of the measured image of the test sample was determined with the help of a simulated image similar to the experiment by adding noise to it. The amount of Gaussian noise was determined from the background of the experimental image and the Poisson noise was varied empirically to match the axial intensity profile of the noisy simulated image and the experimental image. The SNR was then determined by comparing the noisy simulated image (which resembles the experimental image) with the noiseless simulated image using Eq. (13).### 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^{−4},
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.

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

## Acknowledgments

This work is supported by the National Science Foundation (NSF CAREER award DBI-0844682 and NSF IDBR award DBI-0852847, PI: C. Preza) and the University of Memphis (UoM). N. Patwary acknowledges support through the Herff Fellowship of the Herff College of Engineering at UoM. The authors thank L. H. Schaefer (Advanced Imaging Methodology Consultation, Ontario, Canada) for discussions and Matlab code for PSF computation; S. Yuan and D. Dong (UoM) for implementation of the PCA-EM in COSMOS; M. M. Rahman (UoM) for help with noisy simulations; C. M. Waters (University of Tennessee Health Science Center) for preparation of the lung epithelial cells sample; O. Skalli (UoM) for preparation of the Glioblastoma cells sample and S. V. King (UoM) for data acquisition from this sample.

## References and links

**1. **D. A. Agard, “Optical sectioning microscopy: cellular
architecture in three dimensions,” Annu. Rev.
Biophys. Bioeng. **13**(1), 191–219
(1984). [CrossRef] [PubMed]

**2. **J. R. Swedlow, J. W. Sedat, and D. A. Agard, “Deconvolution in optical
microscopy,” Deconvolution of Images and
Spectra **285**, 284–309
(1997).

**3. **J. G. McNally, T. Karpova, J. Cooper, and J. A. Conchello, “Three-dimensional imaging by deconvolution
microscopy,” Methods **19**(3), 373–385
(1999). [CrossRef] [PubMed]

**4. **D. A. Agard, Y. Hiraoka, P. Shaw, and J. W. Sedat, “Fluorescence microscopy in three
dimensions,” Methods Cell Biol. **30**, 353–377
(1989). [CrossRef] [PubMed]

**5. **T. J. Holmes, “Expectation-maximization restoration of band-limited,
truncated point-process intensities with application in
microscopy,” J. Opt. Soc. Am. A **6**(7), 1006–1014
(1989). [CrossRef]

**6. **T. J. Holmes, “Maximum-likelihood image restoration adapted for
noncoherent optical imaging,” J. Opt. Soc. Am.
A **5**(5), 666–673
(1988). [CrossRef]

**7. **S. F. Gibson, *Modeling the
Three-Dimensional Imaging Properties of the Fluorescence Light
Microscope* (Carnegie Mellon University, 1990).

**8. **C. Preza and J.-A. Conchello, “Depth-variant maximum-likelihood restoration
for three-dimensional fluorescence microscopy,” J.
Opt. Soc. Am. A **21**(9), 1593–1601
(2004). [CrossRef] [PubMed]

**9. **V. Myneni and C. Preza, “3-D
reconstruction of fluorescence microscopy image intensities using multiple
depth-variant point-spread functions,” Digital Image Processing and
Analysis (DIPA), Imaging and Applied Optics, OSA Optics and Photonics Congress,
DTuA2 (2010). [CrossRef]

**10. **C. Preza and V. Myneni, “Quantitative depth-variant imaging for fluorescence
microscopy using the COSMOS software package,”
Proc. SPIE **7570**, 757003
(2010). [CrossRef]

**11. **M. Arigovindan, J. Shaevitz, J. McGowan, J. W. Sedat, and D. A. Agard, “A parallel product-convolution approach for
representing the depth varying point spread functions in 3D widefield
microscopy based on principal component analysis,”
Opt. Express **18**(7), 6461–6476
(2010). [CrossRef] [PubMed]

**12. **S. Yuan and C. Preza, “Performance evaluation of an image estimation
method based on principal component analysis (PCA) developed for
quantitative depth-variant fluorescence microscopy imaging,”
Proc. SPIE **8227**, 82270H (2012). [CrossRef]

**13. **S. Yuan and C. Preza, “3D fluorescence microscopy imaging accounting
for depth-varying point-spread functions predicted by a strata interpolation
method and a principal component analysis method,”
Proc. SPIE **7904**, 79040M (2011). [CrossRef]

**14. **S. Yuan and C. Preza, “3D
computational microscopy with depth-varying point-spread functions using a
principal component analysis method,” in *Imaging and Applied Optics*, OSA Technical
Digest (online), paper IM3E.4 (2013). [CrossRef]

**15. **N. Patwary and C. Preza,
“Computationally tractable approach to PCA-based depth-variant PSF
representation for 3D microscopy image restoration,” in *Classical Optics 2014*, OSA Technical Digest (online), paper CW1C.5
(2014). [CrossRef]

**16. **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]

**17. **J.-A. Conchello and J. G. McNally,
“Fast regularization technique for expectation maximization algorithm for
optical sectioning microscopy,” Proc. SPIE 2655, Three-Dimensional
Microscopy: Image Acquisition and Processing III, 199 (1996).

**18. **I. J. Good, “Non-Parametric Roughness Penalty for
Probability Densities,” Nat. Phys. Sci
(Lond.) **229**(1), 29–30
(1971). [CrossRef]

**19. **Computational Imaging Research
Laboratory, Computational Optical Sectioning Microscopy Open Source (COSMOS)
software package; *http://cirl.memphis.edu/COSMOS*.

**20. **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]

**21. **O. Haeberlé, “Focusing of light through a stratified medium:
a practical approach for computing microscope point spread functions. Part
I: Conventional microscopy,” Opt. Commun. **216**(1), 55–63
(2003). [CrossRef]

**22. **D. L. Snyder, C. W. Helstrom, A. D. Lanterman, M. Faisal, and R. L. White, “Compensation for readout noise in CCD
images,” J. Opt. Soc. Am. A **12**(2), 272–283
(1995). [CrossRef]

**23. **N. Patwary, “Performance Aanalysis of
PCA-based Image Reconstruction in 3D Wide Field Fluorescence Microscopy (MS
Thesis),” University of Memphis, https://umwa.memphis.edu/etd/index.php (2014).

**24. **J.-B. Sibarita, *Deconvolution
Microscopy* (Springer Berlin / Heidelberg, 2005).

**25. **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]

**26. **M. Bertero and P. Boccacci,
*Introduction to Inverse Problems in Imaging* (CRC press,
2010).