Three-dimensional block-based restoration integrated with wide-field fluorescence microscopy for the investigation of thick specimens with spatially variant refractive index

Abstract. Development of a block-based restoration (BBR) method that addresses spatially variant (SV) imaging in wide-field fluorescence microscopy of thick samples is presented. The BBR method is based on a block-based imaging model, which approximates SV imaging using an efficient orthonormal basis decomposition of multiple SV point-spread functions computed at block vertices. The effect of reducing the number of blocks needed to account for SV imaging on the restoration accuracy was investigated with simulations using a numerical lung tissue phantom relevant to biological studies. Results show that reducing the number of blocks by 82% and 98% resulted in a 19% and 27% reduction in restoration accuracy, respectively, thereby establishing a reasonable tradeoff between computational resources and accuracy. Comparison of the BBR method to existing methods (deconvolution) that do not account for SV imaging demonstrates a 90% improvement in restoration accuracy. BBR results from synthetic and experimental images of a controlled test sample with SV refractive index (RI) show consistency, providing a validation of the BBR approach. In this study, information from DIC and fluorescence images was combined to identify regions with changing RI within the imaging volume. The BBR method provides a first step toward computationally tractable reconstruction of images from thick samples.

-Three-dimensional block-based restoration integrated with wide-field fluorescence microscopy for the investigation of thick specimens with spatially variant refractive index Sreya Ghosh and Chrysanthe Preza*

Introduction
Computational processing of images and the methodology of acquiring them using relevant optics are rarely seen as independent operations in modern microscopy. Computational optical sectioning microscopy (COSM) 1 refers to the process that results in three-dimensional (3-D) imaging after simultaneous processing of multiple two-dimensional (2-D) images, acquired while focusing the objective lens at different axial locations within a sample. Imaging optically thick samples [>5 μm thick with a spatially variant (SV) refractive index (RI) distribution within the imaging volume contributing to an optical path length change] is prone to spherical aberration (SA) 2 in addition to out-of-focus light 3 and noise, which can be addressed by COSM restoration algorithms, thereby improving resolution, contrast, and signal-to-noise ratio. Development of model-based COSM restoration algorithms requires solving two interconnected imaging problems: (1) modeling the image formation process to produce a forward imaging model 1 and (2) obtaining a solution to the inverse imaging problem using the forward imaging model to restore the true fluorescence intensities in the underlying sample. 4 In this paper, we report new developments that address both these challenges.
Aberrations in 3-D imaging of optically thick samples occur due to mismatch between imaging and design conditions. In this paper, we focus on ameliorating SA caused by changes in the spherical wave-front of the emitted light 5 due to an RI mismatch between the different layers that make up the imaging system, namely the immersion medium of the lens, coverslip, and specimen. 6 SA introduces asymmetry, as well as attenuation and spreading of intensities in the point spread function (PSF), which characterizes the system's response. 6,7 Biological samples are generally thicker than the depth of field of high NA lenses (<2 μm) and have a variable RI distribution consistent with their cellular components. For example, in lung tissue with alveolar edema, 8,9 there are air-filled (RI ¼ 1.00) and water-filled (RI ¼ 1.33) alveoli surrounded by cellular matter (RI ≈ 1.33). Images from such samples exhibit an SV amount of SA depending on the geometry and RI distribution of the sample. Ideally, in order to have the most accurate imaging model, a distinct PSF 6,7 is needed for every point in the object space; however, this is not computationally feasible for COSM algorithms, creating the need for an approximation.
Existing commercial COSM algorithms, for computational ease, assume that the microscope can be represented by a single PSF, 1,10 that is, that the system is space-invariant (SI), an assumption that is valid for thin samples only. This assumption has been relaxed in the development of depth-variant (DV) algorithms, [11][12][13][14] which address SA due to thick samples that can be modeled by a uniform RI distribution. However, depth variability alone is not a good approximation for many biological samples, such as lung tissue that exhibits a variable RI distribution. In this paper, we develop a 3-D restoration algorithm for COSM that addresses the need to account for SVaberrations due to thick samples with nonuniform RI. We present a brief background on SV imaging followed by a discussion of work accomplished in DV microscopy.
2-D SV restoration based on superposition was initially explored in astronomy to account for variance in images acquired from the Hubble telescope. The mathematical groundwork for this effort was laid down by Trussell and Hunt, 15 who divided a 2-D image into subsections, which were individually processed using a SI method. This approach can lead to restoration artifacts 16 due to a mosaic-like effect, which have been addressed using different interpolation techniques. [17][18][19] Nagy and O'Leary 20 compiled and categorized the concepts of SV superposition in 2-D with linear interpolation between PSFs to prevent restoration artifacts for astronomy. They formally defined two approaches of sub-block convolution, namely overlap-add and overlap-save. In the overlap-add approach, the division of blocks occurs in the object space and is accurate due to preservation of information, while in the overlap-save approach the block division occurs in the image space and can exploit parallelism, making it faster. These algorithms were studied and compared by Rahman et al., 12 establishing the accuracy of overlap-add and the speed of overlap-save.
DV imaging (i.e., imaging that varies along the axial direction only) is an acceptable approximation when imaging specimens with almost uniform RI, in which case the depth-induced SA is introduced due to the thickness of the sample and the RI mismatch between the immersion medium of the lens and an average RI for the sample. In the strata-based approximation model for DV imaging, 13 the object space is divided into a number of layers or strata. Within each stratum, the PSF is the weighted interpolation of two PSFs, one computed at the top of the stratum and one at the bottom. Preza and Conchello 13 introduced an expectation maximization algorithm based on their strata model that solves the 3-D DV imaging problem. 11 Maalouf 21 later introduced a mask-based algorithm (EMMA) based on an approach that divides the image space, as opposed to dividing the object space as in the case of the strata approach. 13 Hadj et al. 22 introduced a blind restoration technique to obtain improved DV images in confocal microscopy, which produces images with significantly reduced intensity spread compared to wide-field microscopy, allowing the user to obtain restoration with minimal information about the system parameters. The methods discussed so far for DV imaging as well as our first attempt for an SV imaging model 23 are computationally intensive.
The need for an efficient method (in terms of both accuracy and speed) practical for use in biological imaging led to the use of principal component analysis (PCA), a method widely used for dimensionality reduction 24 in COSM. PCA-based representation of 3-D DV-PSFs and their use in a 3-D DV imaging model were originally proposed by Arigovindan et al., 25 and Yuan and Preza 26 incorporated the PCA representation of PSFs into a DV expectation maximization algorithm (PCA-EM). The computational burden was further reduced by using fewer PSFs for the generation of the PCA coefficients without loss of accuracy in 3-D restoration. 27 Patwary and Preza 28 also showed that the regularized PCA-EM has better performance than the strata-based EM restoration. 13 Two additional DV algorithms based on an accelerated CG-iteration scheme by Schaefer et al. 29 have been recently developed: (1) the DV-CG algorithm based on the strata model, which was shown to have fast convergence, 12 and (2) the PCA-CG 30 algorithm based on the PCA representation of the DV-PSF, which requires fewer computational resources than the DV-CG algorithm. The block-based restoration (BBR) algorithm presented here extends the PCA-CG algorithm to SV imaging. Our first forward model for SV imaging 23 is based on a block-based (BB) approximation, which is a direct extension of the strata approach. 13 Specifically, the entire object space was divided into a collection of small 3-D blocks with PSFs computed at the faces of the blocks. An optimized combination of the overlap-save and overlap-add methods of interpolation was used to compute the final SV image of the object. The methodology was validated using the experimental and simulated data. 23 The primary drawback of the method is its computational inefficiency, which is ameliorated by the PCA BB model developed in this paper.
The development of the BBR method is motivated by studies of alveolar mechanics, 8,9,31 which use fluorescence microscopy to image the expansion and deflation of air-sacs (alveoli) in thick lung tissue (∼200 μm) with variability of alveolar contents. There are four integrated goals that need to be achieved in order to obtain the solution of this challenging imaging problem. These goals are the successful development of (1) a methodology to derive the sample RI distribution; (2) an accurate but practical PSF model that accounts for specimen RI variability, which we reported elsewhere; 32 (3) a computationally tractable SV forward imaging model; 23,33 and (4) a practical restoration algorithm based on the SV model. The determination of the sample RI map is a challenging problem that has been tackled by several groups [34][35][36][37][38] including ours, 39,40 and it is beyond the scope of this paper. In this paper, we focus on the mathematical formulation, simulation, and experimental validation in the development of the third and fourth goals, in which we integrate our new PSF model from the second goal. This paper is organized as follows: Sec. 2 discusses the mathematical formulation of the PCA BB forward imaging model and its use in the BBR method. Methodology and results from simulations with optically thick specimens and SV RI distribution, used to establish the usefulness of the BBR method, are discussed in Sec. 3. Fundamental experimental validation of the BBR method with images from test samples is presented in Sec. 4. Section 5 presents BBR simulation results obtained from images of a lung tissue phantom that approximate our motivating application and provide a controlled environment for testing the algorithm. A summary of the findings from this study and further developments are discussed in Sec. 6. Part of the research presented here was also presented elsewhere. 30,33 2 Block-Based Restoration Method In this section, the mathematical development of the BBR method is presented. We first describe the mathematical formulation of the BB approximation for the SV forward imaging model. A solution to the SV imaging problem is developed by incorporating a computationally efficient BB forward model as the kernel of a fast converging iterative algorithm.

Block-Based Forward Imaging Model
The 3-D image, gð̱ x i Þ, formed by the microscope can be represented using a SV kernel, K, defined below as the integral of the product of the specimen function, fð̱ x o Þ, and the SV-PSF hð̱ x o ; ̱ x i Þ, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 6 3 ; 6 9 0 gð̱ x i Þ ¼ Kfð̱ (1) where ̱ x o ¼ ðx o ; y o ; z o Þ is a point in the object space and ̱ x i ¼ ðx i ; y i ; z i Þ is a point in the image space. In SV imaging, the image of every point in the specimen function is associated with a unique SV-PSF [ Fig. 1(a)], rendering evaluation of Eq. (1) computationally impractical.
In the BB approximation, the object space is conceptually sectioned into M sections along X; N sections along Y, and K sections along Z based on the object's RI distribution map, which results in the formation of blocks [ Fig. 1(b)] that can be associated with a single RI value (a local average of the RI distribution over the block volume). A block is not necessarily cubic, and the number of blocks needed is adjusted based on the variability of the specimen's RI; that is, a sample with a rapidly varying RI is approximated by more blocks than a sample with a slowly varying RI.
The minimum size of a block is defined by the smallest volume of uniform RI distribution determined based on the 3-D RI map of the sample. A block is mathematically defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 6 9 7 B m;n;k ¼ f̱ (2) The intensity in the sample can be expressed by the sum of all the nonoverlapping blocks as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 6 0 8 where E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 1 ; 3 2 6 ; 5 4 8 Eight SV-PSFs associated with each block (one for each block vertex location) are computed using the N-interface PSF model, 32 which models light propagation through N stratified layers within a block. SV-PSFs are calculated at discrete locations, ðx o ; y o ; z o Þ ¼ ðX m ; Y n ; Z k Þ, marking block vertices, using imaging conditions including thickness and RI of the sample at these unique locations. These SV-PSFs can be represented using a few principal components (PCs), thereby reducing not only the memory required but also the number of convolutions in the forward SV imaging model. 33 The PCA formulation used in the SV imaging model is an extension of the PCA approach developed previously to represent DV-PSFs, 25 in that it uses SV weighting coefficients, c j ð̱ x o Þ, instead of the DV ones, c j ðz o Þ. Each SV-PSF can be written in terms of the PCA as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 6 3 ; where P 0 ð̱ x i Þ is the mean of the SV-PSFs, P j ð̱ x i Þ is the j'th PC, and c j ð̱ where ⊗ is the 3-D convolution operator. Since multiple PCs are used to represent the SV image of the object, multiple (J) convolutions are required to generate the forward model. As evident in Eq. (5), every point ̱ x o in the object space is associated with its own unique PCA weighting coefficient c j ð̱ x o Þ. Obtaining a SV-PSF at every voxel presents a challenge; hence, fewer PSFs than the number of voxels in the object space are used, resulting in PCA coefficients at discrete locations within the object space. However, to compute the forward image, coefficients are required for all the points in the object space [Eq. (5)]. This has been addressed by taking advantage of the regular pattern of the coefficients as proposed by Patwary and Preza, 27 where a splinebased interpolation was used to predict the missing weighting coefficients, extending the weighting matrix ½c j ðx m ; y n ; z k Þ available at discrete locations in the object space to a new weighting matrix ½c j ðx o Þ for the entire object space. The BB forward model based on the PCA representation [Eq. (5)] requires (J) convolutions to compute the forward image while our first BB SV model, 23 which relied on the overlap-add and the overlap-save methods, requires (2R þ 1) convolutions, 41 where R > J as shown in Eq. (5) (in practice R ≫ J). Thus, the BB-PCA model, which forms the basis of the BBR method [summarized in Fig. 1(c) and discussed in what follows] is computationally more efficient than the BB model.

Block-Based Restoration Method
Iterative, statistical restoration algorithms are well suited for solving the ill-posed inverse imaging problem for 3-D fluorescence microscopy, 42 since they have the potential to restore missing frequencies, thereby achieving optical sectioning while avoiding direct inversion of the image formation kernel. The computation of the forward image [Eq. (5)] requires many resources (memory, processing time); hence a fast converging algorithm is a requirement for practical SV restoration. Therefore, for our BBR method, we chose an accelerated conjugate gradient (CG)-based iteration method, 29 which has been shown to provide (1) a twofold increase in processing speed while providing accurate restoration results in the solution of the SI inverse imaging problem for fluorescence microscopy, 29 and (2) faster convergence over the expectation maximization algorithm in the solution of the DV inverse imaging problem. 13 For the BBR method, the algorithm by Schaefer et al. 29 was adapted to include the PCA-based SV kernel [Eq. (5)] redefined on a finite discrete domain, K SV-PSF ∶R d → R d , where d is the number of voxels in the acquired data. Schaefer et al. 29 developed a regularized CG-based method that minimizes the generalized restoration functional E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 6 3 ; 3 8 1 ΦðsÞ ¼ L P ðs 2 Þ þ βRðs 2 Þ; (6) where L P ðs 2 Þ is the log-likelihood function for a Poisson distribution in terms of the vector s 2 ¼ f that contains samples from the object function, fð̱ x o Þ, and ensures positivity of the estimated intensities; β ≥ 0 is the regularization parameter that determines the tradeoff between the regularization and the data fit; and Rðs 2 Þ is a regularization functional that measures the smoothness of the object, f, by its distance to the SV image, which is represented by the vector g in a finite discrete domain. For a system with Poisson noise, the negative loglikelihood functional with the SV forward imaging kernel is given as 29 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; 2 2 4 L P ðs 2 Þ ¼ where ðK SV-PSF s 2 Þ n indicates the n'th component of the vector K SV-PSF s 2 . The gradient of the restoration functional [Eq. (6)] is given as 29 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 6 3 ; 1 4 9 ½∇φðsÞ T ¼ 2s where the sign "*" indicates pointwise multiplication of two vectors; ̱ 1 is a vector with all its components equal to 1; and K T SV-PSF is the transpose of the matrix. Equations used in the implementation of the accelerated CG iteration, and the step size computation can be found within Ref. 29.
A schematic diagram summarizing the steps in the BBR approach is shown in Fig. 1(c). Computation of the SV forward image [Eq. (5)] does not rely on sectioning the object or image space, but rather it uses interpolated PCA weighting coefficients to effectively and efficiently model the variance in SV-PSFs computed at the vertices of a small number of blocks defined in the object space (Sec. 2.1). Thus the BBR approach avoids edge artifacts due to sectioning, which have to be accounted for in other BBR techniques, 20 while reducing computational resources. These advantages make the BBR method a suitable approach for SV restoration. In the next section, the application of the BBR method to simulated and experimental data is discussed.

Evaluation of the Block-Based Restoration
Method with Simulated Images In this section, we use the BBR method to restore images from a simulated test sample. In what follows, we discuss the generation of the simulated phantom, computation of its SV image, and its restoration. We further compare the performance of the new BBR method with the performance of DV and SI algorithms applied to the same data.

Spatially Variant Imaging Using the Block-Based Forward Model
In this study, imaging a specimen with only laterally variant RI was investigated. This was done to ensure that the change along the axial direction in the final image could be attributed exclusively to the change in the sample thickness; on the other hand, the lateral change could be attributed to RI variability. This study was performed to investigate the impact of specimen RI variability on the 3-D fluorescence microscopy image of the specimen. The image of a 3-D numerical object [ Fig. 2(a)] with five spherical shell structures dispersed in a medium with variable RI over a grid size of 256 × 256 × 256 voxels was generated using the BB-PCA forward imaging model [Eq. (5)]. The spherical shells have an outer diameter of 2 μm and an inner diameter of 1 μm. The chosen structures are identical in shape to highlight the changes introduced due to SV imaging. Three of the spherical shells were aligned along the Z-axis to exhibit image variability introduced due to depth, and three were aligned along the X-axis to exhibit changes in the image due to the RI variability at the same axial location (i. e., at a constant Z). The specimen was assumed to have three distinct RIs, changing between the values of 1.47, 1.52, and 1.56 shown by the three shaded columns, representing the change in RI, in Fig. 2(b).
For this simulation, PSFs for a 63 × ∕1.4 NA oil lens (RI ¼ 1.515) at a 540-nm wavelength were calculated at different depths within the specimen (in the range of 20 to 40 μm below the coverslip and in depth steps of 5 μm). Due to the sample RI variability, the system experiences both negative (when the RI of immersion medium >RI of the sample) and positive (when the RI of immersion medium <RI of the sample) SA. For this example, the specimen was divided into three sections along the X and Y axes based on the RI change, and to four sections along the Z-axis to account for the change in the sample thickness; hence, the specimen was approximated by (3 × 3 × 4 ¼) 36 nonoverlapping blocks. All PSFs have a 0.1 μm × 0.1 μm × 0.1 μm cubic voxel. MATLAB ® code developed for the N-interface PSF model was used to compute the [ð3 þ 1Þ × ð3 þ 1Þ × ð4 þ 1Þ ¼] 80 SV-PSFs at the block vertices.
PCA was applied to the 80 SV-PSFs to obtain 25 PCs and the weighting coefficient at 80 discrete locations in the object space. The choice to generate 24 PCs was arbitrary, but we kept the number of PCs ∼30% of the total number of PSFs used in the computation. The PCA method sorts the PCs based on significance; that is, the first PC has the most information about the PSFs, while information content in subsequent PCs tapers off as the PC number increases. Comparing the values of the weighting coefficients for the 25 PCs generated, the value of the coefficient stabilized at 0 around the 15th PC for all 80 PSFs  5) and (8)] in this study.
To investigate the accuracy of the PCA BB model, using the imaging conditions described previously, two SV images were generated:   Fig. 2(g). Axial intensity profiles taken through the center of the XZ cross-section show that the intensity peaks in the object appear attenuated, spatially shifted, and blurred in the SV images computed with the two different models [ Fig. 2(g)]. Since the spread of intensities is anisotropic in wide-field microscopy, the axial intensity profiles show attenuation of intensity values and loss of structural integrity, while the lateral profiles preserve some of the structural details even though attenuation is still present. As evident from the profiles, there is a difference in the intensity of the two SV images [Figs. 2(c) and 2(f)] based on the BB and BB-PCA approximations; although overall the profiles demonstrate the same trend in intensity change. To further quantify this difference, we used the structural similarity index measure (SSIM), 43 which was found to be equal to 0.0349 and 0.0217 computed between the true object and the BB and BB-PCA SV images, respectively. In addition, the total integrated intensity (TI) of the numerical object (TI ¼ 1.8270 × 10 7 ) was found to be better preserved by the 3-D SV image predicted by the BB-PCA model (in which TI ¼ 1.7726 × 10 7 ) than the 3-D image predicted by the BB model (in which TI ¼ 2.4468 × 10 6 ). Thus, the BB-PCA model improves computational efficiency (the number of convolutions was reduced by 91% in the simulated study when the BB-PCA model was used over the BB model) without degrading achieved accuracy in the SV image compared to our previously developed approximate BB model.

Restoration of Simulated Spatially Variant Images
The BBR method described in Sec. 2 was applied to the SV image generated using the PCA-based model [ Fig. 2(b)], and restoration results are shown in Fig. 3. Restorations obtained with the space-invariant conjugate gradient (SI-CG) 29 [ Fig. 3(c)] and the DV-CG 12 [ Fig. 3(d)] algorithms, computed from the same SV image [ Fig. 3(b)], are compared to the BBR result [ Fig. 3(e)] to highlight the errors in restoring an SV image using an algorithm based on an inaccurate model assumption.
All algorithms were run for 50 iterations because the CG algorithm ensures fast convergence. The DV-CG and SI-CG algorithms use the same CG type iteration scheme as the BBR algorithm.
For the DV-CG restoration [ Fig. 3(d)], 4 strata and 5 DV-PSFs were used as two DV-PSFs are computed for each stratum, as opposed to eight SV-PSFs for each block in the case of the BBR approach. The 5 DV-PSFs used for the DV-CG restoration were exactly the same as those used to generate the DV image [ Fig. 2(e)]. With the SI-CG algorithm (aka deconvolution), a single nonaberrant PSF was used to restore the same SV image. Figure 3(e) qualitatively resembles the true object [ Fig. 3(a)] with significant structural details restored. The shells qualitatively appear identical in the restored image, highlighting the efficiency of the BBR algorithm to account for spatial variance. Figures 3(c) and 3(d) both show artifacts, although the DV-CG algorithm yields a slightly improved restoration compared to the SI-CG result (quantified by metrics in Table 1 Table 1: the normalized mean square error (NMSE), 26 the SSIM, and the I-divergence (I-div). 44 The BBR is the closest to the true object as judged by the lowest NMSE and I-div values, and the highest SSIM value compared to the DV and SI restoration results (Table 1). Table 1 reinforces the improvement achieved with the BBR method over existing methods currently used for COSM.
This simulation study demonstrated the differences in SV imaging, approximated by the BB forward model and existing DV and SI imaging models. The ability of the BBR method to restore an object from its simulated SV image was also shown to be superior to results from existing restoration methods based on DV and SI assumptions, which suffer from artifacts due to the effective mismatch between the forward imaging model and the actual imaging process.

Application of the Block-Based Restoration
Method to Experimentally Acquired SV-Images In this section, results from a controlled experimental test sample that highlights the effect of RI variance on COSM are discussed. Information derived from differential interference contrast (DIC) microscopy images was used to guide the application of the BBR approach to the experimental data and to generate simulated images that predict the data. Both experimental and simulated images were restored using the BBR method to further validate its usage. Details about data acquisition, comparison of experimental and simulated images, and their restoration are presented in what follows.

Experimental Acquisition of Space-Variant Images
Experimental images of 6 μm beads were acquired from a sample with SV RI using a Zeiss AxioImager (Carl Zeiss, GmbH, Germany) with an AxioCam MRm camera and different imaging modes: wide-field fluorescence, DIC, and bright-field microscopy. FocalCheck ™ microspheres, 6 μm in diameter and stained throughout with DAPI with a 1-μm fluorescent (Alexa Fluor ® ) outer shell, were imaged with a 20 × ∕0.5 NA air lens. The microspheres were air dried on a glass slide, on which guides had been etched with a diamond-head drill in order to facilitate determination of imaging depth from the slide, hence acting as guides. Half of the slide was sealed with UV cured optical cement (NOA 60) with RI ¼ 1.52, and the other half was sealed with ProLong ® Diamond antifade reagent with RI ¼ 1.42 (after being cured for 24 h). A layer of FluoSpheres ® dyed with DAPI was air dried on a coverslip to facilitate locating the exact location of the coverslip. When the coverslip was placed on the sample, a large air vacuole was introduced in the sample. This made the sample comprise of media with three different RIs. The sample was further cured to obtain a well-sealed arrangement. A photograph of the slide (acquired with a Samsung Galaxy S4) shows the boundaries between the different media with varying RI, and the approximate area imaged with the microscope indicated by the red oval [ Fig. 4(a)]. The location of the microscope stage recorded at the best focus of the various beads in the sample was used to compute the thickness of each layer in the sample. The DIC microscopy mode was used to bring into focus the etched guide on the slide. The location of the guide in the sample was identified by acquiring a wider field of view, encompassing the guide as well as the fluorescent shells, using the 3-D Panorama mode in the Zeiss ZEN software available on our microscope. A region-of-interest (ROI) [enclosed by the blue box Fig. 4(b)] was identified, and 3-D data were acquired from this region to evaluate the BB forward model while keeping the computations tractable. In the DIC image [ Fig. 4(b)], the guide is seen going out of focus within the field of view with a change in contrast, which indicates the change in RI. The ROI was imaged using the DIC [ Fig. 4(c)] and bright-field [ Fig. 4(d)] imaging modes, to determine the boundary where the RI transition occurs. The change in contrast at the boundary indicates the presence of a possible gradient of RIs between the two different RIs used in the sample preparation. The location of the boundary and information about the RI of the specimen were used to infer an approximate 3-D RI map by associating regions in the 3-D image with a unique RI. Images were acquired at the boundary between the air vacuole (RI ¼ 1.00) and the mountant (RI ¼ 1.42) to image the effect of a large RI difference similar to that encountered in alveolar imaging. A schematic of the imaging system, shown enclosed by the red circle in Fig. 4(a), is shown [ Fig. 4

(e)].
A 3-D image (grid size of 340 × 340 × 300) comprising seven shells, with a voxel size of 0.32 × 0.32 × 0.32 μm in the X; Y, and Z directions was acquired from the sample described above [ Fig. 4(g)]. The image was gathered using the AlexFluor illumination (emission wavelength at 535 nm). The thickness of the specimen layer was determined to be 143 μm from the imaging distance between the 6-μm beads and the subresolved marker beads (located on the coverslip). The imaging distance was first calibrated to account for the shift introduced due to SA. 2 In the next section, the experimental image described here is compared to a model prediction in order to validate the SV model proposed in this paper. Details about the simulation, comparison of simulated and measured images and restoration from these images are discussed in what follows.

Comparison Between Block-Based Model Predictions and Experimental Images
The information available from the experimental conditions including information about the spatial variability of RI derived from DIC images, described in Sec. 4.1, was used in the computation of the BB forward imaging model to generate the simulated images. A set of theoretical SV-PSFs were computed using the experimental parameters and the N-interface PSF generation model 32 and were used to generate a model prediction of the experimental data.
Based on the knowledge about the slide preparation, the right side of the slide is comprised of mountant while the left side is comprised of air vacuole [ Fig. 4(a)]. The boundary between the air vacuole and the mountant is highlighted using the white dashed line in Figs. 4(b) and 4(c). The transition in contrast can be attributed to the optical path difference (OPD), which can be imaged by DIC but not by bright-field microscopy [ Fig. 4(d)]. The lack of contrast change in Fig. 4(d) confirms that the contrast change in Fig. 4(c) can be attributed exclusively to OPD and not to the presence of a feature of interest, allowing us to conclude that the visible change in contrast  Fig. 4(i)], shows agreement with respect to the shape of the intensity distribution. BBR results from a measured 3-D image of a smaller ROI from this test sample, acquired using the same experimental conditions discussed here, are presented in the next section. Our current implementation of the BBR algorithm is not optimized, and thus the use of a smaller 3-D image allowed us to keep the computation time for the restoration tractable.

Restoration of Space Variant Images Using the Block-Based Restoration Method
Results from applying the BBR method to a 3-D experimental image of the test sample (Sec. 4.1) and a corresponding simulated image are summarized in Fig. 5, where XY section images through the volume are shown. The ROI [red square in Fig. 5(a)] shows two 6-μm spherical shells known to be physically located at the same depth by construction of the sample as described in Sec In the simulation study, the BBR method was found to converge after 25 iterations based on the I-div value computed between the restored image and the numerical object at each iteration. Thus, results reported here from both experimental and simulated images are after 25 BBR iterations. Regularization [β ¼ 0.0001 in Eq. (6)] was used during the restoration of the experimental image to address instability in the solution due to the presence of noise introduced during data acquisition. The XY cross-section of the restoration obtained from applying the BBR method to the 3-D experimental image is shown in Fig. 5(g), while the corresponding section of the restoration obtained from the 3-D simulated image is shown in Fig. 5(h). Both restored images show that the spherical shell structures are restored well as evident by the sharp high intensity ring. The BBR method was effective in eliminating the SA-induced axial shift, which caused the image of Shell 2 to appear blurred in Figs. 5(e) and 5(f), thereby bringing both spherical shells into focus at the same axial plane [Figs. 5(g) and 5(h)], as expected (because both spherical shells are located at the same depth within the sample by construction). The lateral intensity profiles through the center of the experimental image and the restoration [Figs. 5(e) and 5(f)] are compared in Fig. 5(i) to the true intensity of the numerical spherical shells [ Fig. 5(c)]. As shown by the profiles, the intensity in both restored shells has two peaks and follows the intensity of the true object.
To quantify the ability of the BBR to restore the experimental and simulated image, the SSIM over a ROI around bead 1 and 2 was calculated and averaged. The SSIM between the numerical object and the restored simulated image is 0.6083, while the restored experimental image yielded an SSIM of 0.5755. The SSIM values are comparable, indicating consistency between the experimental and simulated restoration results.
In this section, the experimental images acquired from a controlled test sample show SV imaging due to change in RI of the media. This section highlights the ability of the BB forward model to predict SV images from the test sample, using a simple 3-D RI map inferred from a DIC image. The BBR results obtained from the experimental and simulated SV images of a test sample were found to be consistent, validating the implementation of the approach for the conditions tested. In the next section, results obtained from applying the BBR method to a lung tissue phantom are presented to demonstrate the effect of SV imaging and BBR in a biological application of interest.

Application of Block-Based Restoration to
Space-Variant Imaging of a Simulated Lung Tissue Phantom In this section, we present results from applying the BBR method to simulated data relevant to biological studies that motivate our work. Toward this end, we combined available information from physiological and mechanical studies of lung alveoli to generate a simulated lung tissue phantom. In what follows, we discuss how the phantom and its RI map were used with the BB forward model to simulate its SV image, and we report results from processing the SV image using the BBR approach. A study of the effect of the block size used in the BBR approach on the restoration accuracy is also presented.

Modeling the Lung Tissue Phantom
The BBR method as elucidated before has been developed for imaging specimens with SV RI, such as lung tissue, which exhibits RI values in the range of 1.00 to 1.342. The lung is composed of clusters of small air-sacs (alveoli) divided by thin, elastic walls, that is, membranes. Alveoli normally have a thin wall that allows for air exchange to occur, and fluids are usually kept out of the alveoli unless these walls lose their integrity. Pulmonary edema occurs when the alveoli fill up with excess fluid seeping out of the blood vessels in the lung instead of air. 45 For this study, a single alveolus with edema surrounded by air-filled alveoli was modeled.
Reported information about alveoli depends on parameters specific to the animal tissue, sample preparation, and the imaging tool used. To guide the generation of our numerical phantom, results from Perlman and Bhattacharya 8 were used because they studied alveoli with and without edema using fluorescence microscopy. In their study, 3-D images were acquired using a confocal microscope at a subplueral depth of 80 to 100 μm in lung tissue from male Sprague-Dawley rats weighing 300 to 600 g. The alveolar walls were stained with calcein red-orange AM for visibility in fluorescence. To establish edema a single alveolus was filled with an albumin solution (RI ¼ 1.336).
A cluster of alveoli mechanically modeled as a 3-D arrangement of dodecahedrons by Kitaoka et al. 46 forms the basis of our phantom. Because a slice from a dodecahedral structure appears as a honeycomb, we modeled an 80-μm thick slice of rat alveoli by stacking 2-D honeycomb arrangements to represent air-filled and fluid-filled alveoli surrounded by alveolar epithelial cells. The 3-D RI map of the phantom was generated using existing information from the literature about different media. The value used for the RI of the rat alveoli was obtained from Ref. 37. In the phantom, a 30-μm layer of a mountant, such as the ProLong ® Diamond antifade reagent (RI ¼ 1.46) often used in samples to prevent fluorescence bleaching, was also introduced, thereby introducing an additional layer of a medium with different RI. To mimic images from real samples and present restoration challenges, small nodular structures similar to cancerous nodules were added in alveoli with edema and without edema. The sample was imaged with a 20 × 0.5 NA air lens with a sampling of 0.32 × 0.32 × 0.32 μm 3 .

Application of the Block-Based Restoration
Method to the Lung Phantom  the alveoli suffering from edema are affected by a larger amount of SA (characterized by loss of structure due to spreading of intensities and axial focal shift) than the other two alveoli. Two nodules [labeled in Figs. 6(d) and 6(i)] were selected to assess the effect of SV imaging and BBR on their images. The two selected nodules, although initially located at the same axial location, appear at different imaging depths (6.4 μm apart) in the SV image [ Fig. 6(j)] due to SV SA. This phenomenon is quantified by the axial intensity profiles taken through the center of nodules 1 and 2, shown in Figs. 6(n) and 6(o), respectively.
The application of the BBR algorithm to the SV image of the lung tissue phantom was investigated using three different block sizes used in the restoration, in order to study the effect of this choice on the restoration accuracy. For a fixed image size, using a smaller number of blocks (or larger block size) reduces the number of SV-PSFs that need to be computed, and it also reduces the computation time of the PCs used in the PCA representation of the PSFs. The number of blocks and PSFs used in applying the BBR algorithm to the SV image of the lung tissue phantom in this study is as follows: XY and XZ cross-sections of the restored images after 100 iterations of the BBR algorithm using the three different block sizes are shown in Figs. 6(f)-6(h) and 6(k)-6(m). A minimal change in the I-div computed between the true object and the restored object at every iteration suggested that the algorithm converged after 100 iterations. Restorations obtained with grid (A) provide a reasonable approximation of the true object with the nodules clearly separated from the alveolar walls [Figs. 6(f) and 6(k)]. The overlap of the normalized axial intensity profiles plotted in Figs. 6(n) and 6(o) from the center of nodules 1 and 2, respectively, quantifies the accuracy achieved by the BBR method.
Some deterioration in the restoration obtained using grid (B) is evident in Figs. 6(g) and 6(l). In this case, there is some loss in structural integrity of the alveolar walls, but the nodules are restored with separation from the walls as they appear in the true object. In Fig. 6(g), the nodules are clearly restored to their circular shape. However, nodule 1 is restored better qualitatively than nodule 2, particularly when viewed in the XZ images [ Fig. 6(l)]. Increasing the restoration block size further, as in grid (C), results in more restoration artifacts, particularly in the region of the alveolar walls in the presence of edema [Figs. 6(h) and 6(m)]. In this case, the intensity of nodule 2 is not restored to the true object intensity. As evident in the XZ section image, nodules immersed in fluid are not restored to a spherical shape [ Fig. 6(m)].
To highlight the ability of the BBR to reconstruct the SV image using different block sizes, three error measures computed between the restorations and the true object are reported in Table 2. The change in the block size was introduced to evaluate the ability of the BBR to restore intensity using an approximate model based on a reduced number of PSFs. Results show that the BBR method is able to reconstruct the shape of the nodules and their distance from the alveolar wall with a SSIM ¼ 0.7688 when grid (B) is used, which contains only 18% of the original number of blocks used to compute the simulated image. In the case of grid (C), even though the BBR method uses only 5% of the original number of blocks, it is still able to reconstruct the object with a SSIM ¼ 0.6840, which represents a 27% drop in the restoration accuracy [as quantified by comparing it to the SSIM obtained in the grid (A) restoration]. As evidenced in Fig. 6, nodules in fluid experiencing larger amount of SA are harder to restore when the block size is increased because SV SA is not adequately accounted for. Overall, the achieved restoration of the primary structures of the lung tissue phantom, in the three grid cases tested, is promising, and it highlights the ability of the BBR method and the effect of the block size (or number of blocks) used on the restoration.
In this section, we proposed a lung tissue phantom generated based on available information in the literature in order to study and address imaging challenges of such a sample with fluorescence microscopy. Noiseless simulation results shown in this section suggest that application of the BBR method to SV images from lung tissue can provide successful restoration if an adequate block size (or number of blocks) is used. The numerical phantom developed for this study can be used in future simulation studies to investigate further the sensitivity of the BBR method to the block grid size as well as to the presence of noise.

Summary and Conclusions
The goal of the work presented in this paper is the development of a computationally tractable algorithm for efficient restoration of 3D images suffering from SV aberration from optically thick samples, such as lung tissue. The paper details the development and evaluation of an iterative BBR method that uses an accelerated conjugate gradient type minimization scheme to solve the SV inverse imaging problem. The BBR method is based on a PCA BB forward imaging model, which approximates SV imaging of thick specimens with spatially varying RI. In this BB approach, the object space is divided into blocks, thereby providing a grid of points where SV-PSFs are computed using our previously developed N-layered PSF model. Therefore, the BB approach models SV imaging more accurately than previously developed models based on a SI or DV assumption, which can only account for imaging of thin specimens or thick specimens with uniform RI, respectively. The use of PCA for the SV-PSF representation efficiently reduces the information (and computational resources) needed to compute the BB model and restoration while preserving accuracy. Simulated images computed using the BB-PCA approximation required only 11% of the total number of convolutions required by the BB approximation Table 2 The effect of change in the information available to the BBR algorithm, quantified using NMSE, SSIM, and I-div. without the PCA representation for the same imaging conditions, without appreciable disagreement between the simulated images, quantified using the total intensity of the 3-D images and the SSIM. The effect of SV imaging was modeled using different imaging conditions, and the ability of the BBR algorithm to restore SV images in these cases was investigated. Simulated images of a simple 3-D object (with uniform structures and nonuniform RI) were effectively restored using the BBR with a 90% improvement in accuracy over using SI restoration (aka deconvolution) and 57% improvement over using DV restoration (quantified using the NMSE measure). In this study, the RI values were intentionally chosen to exaggerate the space variability for testing purposes.
The ability of the BBR method to restore intensity values in the experimental images of microspheres was evaluated quantitatively and qualitatively using a controlled test sample created with an RI space-variance due to an air-mountant interface, which mimics the air-water interface in lung tissue in the presence of edema. An approximate 3-D RI map of the sample was created using available information about the RI of the sample media as well as images of the same field of view of the sample obtained with fluorescence and DIC microscopy. A simulated image, generated based on the experimental condition during data acquisition, was shown to capture the main features in the experimentally acquired image, thereby validating the BB SV imaging model. SSIM values of 0.5755 and 0.6083, computed between the numerical object and the restorations from the experimentally acquired image and the synthetic image, respectively, were found to be comparable, demonstrating the ability of the BBR to account for SV SA in the experimental data consistent with the simulated study.
Finally, application of the BBR method to a simulated SV image of a numerical phantom for lung tissue (generated using parameters acquired from the literature and the proposed BB approach) provided promising results with respect to future application of the BBR method to a biological application of interest. This study was also used to study the effect of the block size (a user defined parameter), critical in the BB approach, on the accuracy of the restoration. Results from this study showed that as the block size was increased (thereby reducing the number of blocks in a fixed volume size) the restoration accuracy dropped by 19% and 27% when the number of blocks was reduced to only 18% and 2% of the original number of blocks used to approximate the SV image of the phantom (synthetic data), respectively. Restored images demonstrated the ability of the BBR to effectively restore the primary features of the lung tissue phantom with reduced computational resources (quantified by the reduced number of SV-PSFs needed for the BBR).
Based on our current investigations, the BBR algorithm demonstrates a similar tradeoff between speed and accuracy observed in DV restoration with our previously developed PCA-CG 30 algorithm. In the future, we plan to use the BB forward model in other types of iterative restoration algorithms to investigate possible improvements in accuracy. The BBR method is computationally efficient in terms of the number of convolutions it computes per iteration. In the future, we plan to optimize and parallelize the method to provide further improvements in the processing speed.
The BBR method uses SV-PSFs that capture the SA due to the sample. In the approach presented here, the specimen RI and depth within the sample were used directly to compute the SV-PSFs. However, the quantitative phase image of the sample (which is proportional to the integral of RI over sample thickness) could be used instead in the PSF computation. In future work, we plan to investigate methods for sample segmentation based on its phase image, thereby eliminating the need to compute the RI distribution first from the phase image and the manual determination of the block size. Additionally, noisy simulation studies will be used to evaluate the performance of the BBR method at different noise and regularization levels in order to facilitate its application to biological data.
This study presents a computationally efficient SV forward imaging model and restoration method and provides a first step toward accurate reconstruction of 3-D images from thick biological samples with varying RI. This study contributes to the larger goal of 3-D imaging using COSM.