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

Sreya Ghosh; Chrysanthe Preza

doi:10.1117/1.JBO.21.4.046010

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

Sreya Ghosh, Chrysanthe Preza

Author Affiliations +

Journal of Biomedical Optics, Vol. 21, Issue 4, 046010 (April 2016). https://doi.org/10.1117/1.JBO.21.4.046010

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.

1. 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)¹ 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)² in addition to out-of-focus light³ 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¹ 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.⁴ 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⁵ 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.⁶ SA introduces asymmetry, as well as attenuation and spreading of intensities in the point spread function (PSF), which characterizes the system’s response.⁶^,⁷ 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,⁸^,⁹ there are air-filled ( $RI = 1.00$ ) and water-filled ( $RI = 1.33$ ) alveoli surrounded by cellular matter ( $RI \approx 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⁶^,⁷ 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,¹^,¹⁰ 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,¹¹^–¹⁴ 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 SV aberrations 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,¹⁵ who divided a 2-D image into subsections, which were individually processed using a SI method. This approach can lead to restoration artifacts¹⁶ due to a mosaic-like effect, which have been addressed using different interpolation techniques.¹⁷^–¹⁹

Nagy and O’Leary²⁰ 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.,¹² 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,¹³ 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¹³ introduced an expectation maximization algorithm based on their strata model that solves the 3-D DV imaging problem.¹¹ Maalouf²¹ 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.¹³ Hadj et al.²² 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²³ 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²⁴ 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.,²⁵ and Yuan and Preza²⁶ 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.²⁷ Patwary and Preza²⁸ also showed that the regularized PCA-EM has better performance than the strata-based EM restoration.¹³ Two additional DV algorithms based on an accelerated CG-iteration scheme by Schaefer et al.²⁹ have been recently developed: (1) the DV-CG algorithm based on the strata model, which was shown to have fast convergence,¹² and (2) the PCA-CG³⁰ 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²³ is based on a block-based (BB) approximation, which is a direct extension of the strata approach.¹³ 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.²³ 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,⁸^,⁹^,³¹ which use fluorescence microscopy to image the expansion and deflation of air-sacs (alveoli) in thick lung tissue ( $\sim 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;³² (3) a computationally tractable SV forward imaging model;²³^,³³ 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³⁴^–³⁸ including ours,³⁹^,⁴⁰ 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.³⁰^,³³

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.

2.1.

Block-Based Forward Imaging Model

The 3-D image, $g ({\underset{̲}{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 ({\underset{̲}{x}}_{o})$ , and the SV-PSF $h ({\underset{̲}{x}}_{o}, {\underset{̲}{x}}_{i})$ ,

Eq. (1)

g ({\underset{̲}{x}}_{i}) = K f ({\underset{̲}{x}}_{o}) = \int_{R^{3}} f ({\underset{̲}{x}}_{o}) h ({\underset{̲}{x}}_{o}, {\underset{̲}{x}}_{i}) d {\underset{̲}{x}}_{o},

where

{\underset{̲}{x}}_{o} = (x_{o}, y_{o}, z_{o})

is a point in the object space and

{\underset{̲}{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.

Fig. 1

BB approximation for SV imaging. (a) Impact of imaging depth and RI variability on the 3-D SV-PSF. The solid-line box encloses a set of DV-PSFs. System parameters: $63 \times 1.4 NA$ , oil lens ( $RI = 1.515$ ); grid size: $256 \times 256 \times 256$ with axial and lateral spacing equal to $0.1 μ m$ ; emission wavelength at 540 nm. (b) Sectioning of the object enables computation of SV-PSFs at the eight vertices of each block (only 1 to 7 shown) within the 3-D volume, used by the SV model and the BBR method. (c) Schematic diagram of the BBR method.

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

Eq. (2)

B_{m, n, k} = {{\underset{̲}{x}}_{o} : X_{m} \leq x_{o} < X_{m + 1}; Y_{n} \leq y_{o} < Y_{n + 1}; Z_{k} \leq z_{o} < Z_{k + 1}} m = 1,2, 3, \dots, M; n = 1,2, 3, \dots, N; k = 1,2, 3, \dots, K .

The intensity in the sample can be expressed by the sum of all the nonoverlapping blocks as

Eq. (3)

f ({\underset{̲}{x}}_{o}) = \sum_{m = 1}^{M} \sum_{n = 1}^{N} \sum_{k = 1}^{K} f_{m, n, k} ({\underset{̲}{x}}_{o}),

where

f_{m, n, k} ({\underset{̲}{x}}_{o}) = {\begin{cases} f ({\underset{̲}{x}}_{o}) & for {\underset{̲}{x}}_{o} \in B_{m, n, k} \\ 0 & otherwise \end{cases}

Eight SV-PSFs associated with each block (one for each block vertex location) are computed using the $N$ -interface PSF model,³² 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.³³ The PCA formulation used in the SV imaging model is an extension of the PCA approach developed previously to represent DV-PSFs,²⁵ in that it uses SV weighting coefficients, $c_{j} (\underset{̲}{x_{o}})$ , instead of the DV ones, $c_{j} (z_{o})$ . Each SV-PSF can be written in terms of the PCA as follows:

Eq. (4)

h ({\underset{̲}{x}}_{o}, {\underset{̲}{x}}_{i}) = \sum_{j = 0}^{R} P_{j} ({\underset{̲}{x}}_{i}) c_{j} ({\underset{̲}{x}}_{o}) \approx \sum_{j = 0}^{J} P_{j} ({\underset{̲}{x}}_{i}) c_{j} ({\underset{̲}{x}}_{o}),

where

P_{0} (\underset{̲}{x_{i}})

is the mean of the SV-PSFs,

P_{j} (\underset{̲}{x_{i}})

is the

j

’th PC, and

c_{j} (\underset{̲}{x_{o}})

are the corresponding SV coefficients, with

c_{0} (\underline{x_{o}}) = 1

for all

{\underline{x}}_{o}

. The SV weighting coefficient is computed by the following inner product:

c_{j} (\underset{̲}{x_{o}}) = 〈 P_{j} ({\underset{̲}{x}}_{i}), [h ({\underset{̲}{x}}_{o}, {\underset{̲}{x}}_{i}) - P_{0} (\underset{̲}{x_{i}})] 〉

. Each location

(x_{o}, y_{o}, z_{o}) = (X_{m}, Y_{n}, Z_{k})

is associated with a unique

c_{j} (\underset{̲}{x_{o}})

for each PC. As

j

increases the value of

c_{j} (\underset{̲}{x_{o}})

becomes close to zero [see Fig. 2(h)], allowing the use of fewer components (

J < R

) in the approximation of Eq. (4).

Fig. 2

BB SV imaging of a simulated sample. $X Z$ section images taken through the center of 3-D volumes are shown in all cases. (a) Numerical object with five identical $2 - μ m$ in diameter spherical shells embedded in a medium with RI variance. (b) RI map of the object space showing 12 out of 36 blocks. Blocks with different intensity contrast levels denote RI values of 1.47, 1.52, and 1.56, from left to right. Image of the numerical object computed using (c) the BB SV imaging model²³ and superposition of 160 convolutions (not based on PCA); (d) the space-invariant model and a single nonaberrant PSF; (e) the DV strata-based model with five strata and superposition of 10 convolutions, for uniform sample ( $RI = 1.52$ ); and (f) the BB-PCA SV model [Eq. (5)] computed from the superposition of 15 convolutions. (g) Lateral and axial intensity profiles taken through the center of the object [along the dashed lines in (a)] and its SV images (c) and (f), denoted in the legend by “BB Img.” and “BB-PCA Img.,” respectively. (h) Graphical representation of the first 15 PCA coefficients used in (f) to represent PSF #1 (top) and PSF #66 (bottom). Imaging parameters: $63 \times 1.4 NA$ oil lens ( $RI = 1.515$ ); wavelength at 540 nm. Scale bar: $6.4 μ m$ .

The SV image, $g ({\underset{̲}{x}}_{i})$ , of an object, $f ({\underline{x}}_{o})$ , can be written using the PCA-represented SV-PSFs, by the kernel $K_{SV - PCA}$ ,

Eq. (5)

g ({\underset{̲}{x}}_{i}) = K_{SV - PCA} f ({\underset{̲}{x}}_{o}) = \sum_{j = 0}^{J} [f ({\underset{̲}{x}}_{o}) c_{j} ({\underset{̲}{x}}_{o})] \otimes P_{j} ({\underset{̲}{x}}_{i}),

where

\otimes

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

{\underset{̲}{x}}_{o}

in the object space is associated with its own unique PCA weighting coefficient

c_{j} ({\underset{̲}{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,²⁷ where a spline-based 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} ({\underline{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,²³ which relied on the overlap-add and the overlap-save methods, requires ( $2 R + 1$ ) convolutions,⁴¹ 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.

2.2.

Block-Based Restoration Method

Iterative, statistical restoration algorithms are well suited for solving the ill-posed inverse imaging problem for 3-D fluorescence microscopy,⁴² 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,²⁹ 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,²⁹ and (2) faster convergence over the expectation maximization algorithm in the solution of the DV inverse imaging problem.¹³

For the BBR method, the algorithm by Schaefer et al.²⁹ was adapted to include the PCA-based SV kernel [Eq. (5)] redefined on a finite discrete domain, $K_{SV - PSF} : R^{d} \to R^{d}$ , where $d$ is the number of voxels in the acquired data. Schaefer et al.²⁹ developed a regularized CG-based method that minimizes the generalized restoration functional

Eq. (6)

Φ (s) = L_{P} (s^{2}) + β R (s^{2}),

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 ({\underset{̲}{x}}_{o})

, and ensures positivity of the estimated intensities;

β \geq 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 log-likelihood functional with the SV forward imaging kernel is given as²⁹

Eq. (7)

L_{P} (s^{2}) = \sum_{n} {(K_{SV - PSF} s^{2})}_{n} - \sum_{n} g_{n} \ln {(K_{SV - PSF} s^{2})}_{n},

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²⁹

Eq. (8)

{[\nabla φ (s)]}^{T} = 2 s * K_{SV - PSF}^{T} (\underset{̲}{1} - \frac{g}{K_{SV - PSF} s^{2}}) + β {[\nabla R (s^{2})]}^{T},

where the sign “*” indicates pointwise multiplication of two vectors;

\underset{̲}{1}

is a vector with all its components equal to 1; and

K_{SV - PSF}^{T}

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,²⁰ 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.

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

3.1.

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 \times 256 \times 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 \times / 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 \times 3 \times 4 =$ ) 36 nonoverlapping blocks. All PSFs have a $0.1 μ m \times 0.1 μ m \times 0.1 μ m$ cubic voxel. MATLAB^® code developed for the $N$ -interface PSF model was used to compute the [ $(3 + 1) \times (3 + 1) \times (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 $\sim 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 [Fig. 2(h)]; that is,

c_{1, \dots, 14} (x_{m}, y_{n}, z_{k}) > c_{15} (x_{m}, y_{n}, z_{k}) \approx 0

for

m = 1,2, 3, \dots M; n = 1,2, 3, \dots, N; k = 1,2, 3, \dots, K

, indicating that the effect of the PCs beyond the 15th PC is negligible. This allowed us to eliminate PCs with negligible contribution, thereby reducing the number of convolutions required in the BBR approach [Eqs. (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:

1. One image using our first SV model (based on the overlap-save and overlap-add methods) described in Ref. 23. This method required 160 convolutions to generate the SV image [Fig. 2(c)] using 80 PSFs;
2. Another image using Eq. (5) and the PCA-represented PSFs [Fig. 2(f)]. This method required only 15 convolutions (using 14 PCs and the base component in the restoration) in addition to the PCA, which is precomputed once.

For comparison, we also generated an SI image [Fig. 2(d)] using a single nonaberrant PSF and a DV image [Fig. 2(e)] using a series of DV-PSFs extracted from the larger set of SV-PSFs (i.e., PSFs only affected by the change in the 20 to $40 μ m$ depth range) assuming that the sample has a uniform $RI = 1.52$ . This particular RI value was chosen since it is approximately equal to the average of all the RIs used in this study.

3.2.

Comparison of Simulated Images Using Different Forward Imaging Models

In the SV image [Fig. 2(c)], the rightmost shell experiences negative SA, while on the other hand, the leftmost shell experiences positive SA [Fig. 2(c)], indicated by the direction of the light spread. The SI image [Fig. 2(d)] and DV image [Fig. 2(e)] computed using earlier imaging approximations are qualitatively different from the SV images of this test object [Figs. 2(c) and 2(f)]. This comparison shows the effect of progressively incorporating more physics-based assumptions for the specimen and image formation in order to improve accuracy in approximation models used in COSM. To quantify the observable change in the SV images [Figs. 2(c) and 2(f)], intensity profiles are plotted and compared to the object intensities in Fig. 2(g). Axial intensity profiles taken through the center of the $X Z$ 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),⁴³ 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 \times 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 \times 10^{7}$ ) than the 3-D image predicted by the BB model (in which $TI = 2.4468 \times 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.

3.3.

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)²⁹ [Fig. 3(c)] and the DV-CG¹² [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.

Fig. 3

Evaluation of the BBR method in simulation. $X Z$ section images taken through the 3-D volumes are shown in all cases. (a) Numerical object as in Fig. 2(a). (b) BB-PCA SV image as in Fig. 2(f). The restoration obtained from (b) using (c) the SI-CG algorithm and a single nonaberrant PSF, (d) the DV-CG algorithm and five DV-PSFs in the depth range of 20 to $40 μ m$ (for uniform sample $RI = 1.52$ ), and (e) the BBR method with 14 PCs and coefficients computed from 80 SV-PSFs. Intensity profiles taken through the object at the red-dashed lines in (a) and the restored images (c, d, and e) along (f) the optical axis $z$ , and (g) the lateral axis $x$ . Imaging parameters are $63 \times 1.4 NA$ oil lens ( $RI = 1.515$ ); emission wavelength is 540 nm. Scale bar: $6.4 μ m$ .

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) because the DV-PSFs used by the DV-CG account for the DV SA. Quantitative comparison between the restored images and the object is demonstrated by axial [Fig. 3(f)] and lateral [Fig. 3(g)] intensity profiles taken through the center of the images. The intensity profile through Fig. 3(e) shows that the reconstruction of the shell structure has two distinct peaks per shell, while intensity profiles through the other two reconstructions [Figs. 3(c) and 3(d)] fail to do so. To further quantify the advantage of using BBR to restore SV images, the error between the true object and the restorations shown in Figs. 3(c)–3(e) is investigated using three discrepancy measures summarized in Table 1: the normalized mean square error (NMSE),²⁶ the SSIM, and the I-divergence (I-div).⁴⁴

Table 1

The advantage of using BBR over DV and SI approaches, quantified using NMSE, SSIM, and I-div.

Restoration	NMSE	SSIM	I-div
BBR	0.1180	0.8796	$3.83 \times 10^{4}$
DV	0.4322	0.7264	$2.90 \times 10^{6}$
SI	1.8410	0.0870	$6.17 \times 10^{6}$

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.

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

4.1.

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 \times / 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.

Fig. 4

Experimental evaluation of the BB-PCA SV imaging model. (a) Photo of glass slide showing the location of the different media in the experimental test sample. The red oval designates the location imaged. (b) Experimental DIC image (single focal plane of grid size: $2616 \times 2375 pixels$ ) showing the etched guide, the two regions with distinct RI (equal to 1.00 and 1.42), and the ROI (enclosed by the blue rectangle) imaged using different modes. Zoomed version of the ROI from (c) DIC, (d) bright-field, and (h) fluorescence modes. (e) Schematic of imaging layers. $X Y$ section image from (f) the center of the numerical object (used in simulation that approximates the experimental data), (g) the BB-PCA SV simulated image [Eq. (5)], and (h) the experimental fluorescence image. (i) Plot of lateral normalized intensity profiles taken through shells 1 and 2 of the numerical object (f), along the red-dashed lines in the simulated image (g), and the experimental image (h). Scale bar: $6.4 μ m$ (c, d, f–h). Imaging parameters: $20 \times / 0.5 NA$ air lens ( $RI = 1.00$ ); emission wavelength is 540 nm. All images are normalized and displayed on the same intensity scale.

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 \times 340 \times 300$ ) comprising seven shells, with a voxel size of $0.32 \times 0.32 \times 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.² 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.

4.2.

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³² 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(c)] is due to the change in RI. The boundary (dashed white line) is visible in the bright-field image [Fig. 4(d)] because of scattering, even though there is no change in the background intensity.

The simulated image [Fig. 4(g)] captures the main features in the experimental fluorescence image [Fig. 4(h)]. A cluster of three shells comes into focus $\sim 6.4 μ m$ away from the best focus of the other four shells in the field of view. This shift in best focus can be attributed to the different amounts of SA due to the different RIs in distinct volumes of the slide. The cluster of three shells comes into focus at a shallower depth than the other four; this indicates the presence of positive SA. Comparison of intensity profiles, from both the experimental and simulated images, taken through the center of two shells located in different media [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.

4.3.

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 $X Y$ 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. 4.1. The experimental image from this ROI [Fig. 5(e)] was predicted using our SV imaging model [Fig. 5(f)] from a numerical object [Fig. 5(c)] constructed based on an RI map inferred from the DIC image [Fig. 5(b)]. Shell 2 [Fig. 5(c)] is in a medium with $RI = 1.42$ and thus appears to be out of focus [Fig. 5(e)] due to the presence of SA caused by the RI mismatch. On the other hand, Shell 1 [Fig. 5(c)] is in an air vacuole (i.e., in a medium with $RI = 1.00$ ), which matches the immersion medium of the dry lens used, and thus it appears well focused in Fig. 5(e).

Fig. 5

Comparison of BBR results from the experimental and simulated images of a test sample. In all cases, only $X Y$ images taken through the center of the acquired volumes are shown. (a) DIC image of test sample shows the boundary between two media with RI equal to 1.00 and 1.42. (b) A magnified region identified by the red square in (a) shows two spherical shells, each located in one of the two media in the sample. (c) Numerical object constructed based on (b and e) used in simulation. (d) Fluorescent image of the same field of view as in (a). (e) Magnified subimage extracted from (d). (f) Simulated SV image of (c) computed with Eq. (5). BBR result from (g) the experimental image and (h) the simulated image. (i) Lateral intensity profiles through the center of the shells in (c), the experimental image (Img: Shell 1 and Img: Shell 2) before (e) and after BBR (g) (Rstrd: Shell 1 and Rstrd: Shell 2). Results shown after 25 iterations of the BBR algorithm. Grid size: $226 \times 226 \times 300$ for (b, c, e–h). Scale bar (white) for (b, c, e–h): $3.2 μ m$ . Scale bar (red) for (a and d): $6.4 μ m$ . Pixel size: $0.32 \times 0.32 \times 0.32 μ m$ . Imaging parameters: $20 \times / 0.5 NA$ air lens ( $RI = 1.00$ ); emission wavelength is 540 nm. All images are normalized and displayed on the same intensity scale.

A simulated numerical object [Fig. 5(c)], composed of two identical shells, was generated using the manufacturer’s information about the fluorescent shells. To generate a simulated image, the 3-D volume in the object space was divided into $2 \times 2 \times 10$ sections along $X, Y$ , and $Z$ , resulting in 40 blocks and the generation of 99 SV-PSFs. The BB forward model [Eq. (5)] was computed using 12 PCs (most significant) derived from the 99 PSFs and the numerical object to generate the SV image [Fig. 5(f)]. The simulated image [Fig. 5(f)] captures the main features in the experimental image [Fig. 5(e)], although some difference in intensity is evident. The agreement between the experimental and simulated images [Figs. 5(e) and 5(f), respectively] validates the SV forward imaging model used by the BBR method.

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 $X Y$ 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.

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

5.1.

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.⁴⁵ 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⁸ 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.⁴⁶ 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 \times 0.5 NA$ air lens with a sampling of $0.32 \times 0.32 \times 0.32 {μ m}^{3}$ .

5.2.

Application of the Block-Based Restoration Method to the Lung Phantom

$X Y$ and $X Z$ cross-sections from the numerical lung tissue phantom developed and used in our simulation are shown in Figs. 6(d) and 6(i), respectively. To generate the forward image, the object was divided into $7 \times 7 \times 7$ sections along $X, Y$ , and $Z$ , respectively; hence, the object space consists of 343 nonoverlapping blocks. An $X Y$ view of the arrangement of the 343 blocks is shown by the grid in Fig. 6(a). Because the phantom has regions filled with air, edema fluid, and cellular material, the RI map for this sample has values between 1.00, 1.33, and 1.336. 512 SV-PSFs were computed for a $20 \times 0.5 NA$ air lens and were represented using 10 PCs based on the PCA methodology [Eq. (4)] described in Sec. 2.1.

Fig. 6

Application of BBR to lung tissue phantom image. RI map of lung tissue sectioned into different number of blocks: (a) $7 \times 7 \times 7$ blocks (grid A); (b) $4 \times 4 \times 4$ blocks (grid B); and (c) $2 \times 2 \times 2$ blocks (grid C). (d–h) $X Y$ and (i–m) $X Z$ cross-section images taken through the center of 3-D volumes: (d and i) numerical lung tissue phantom with nodule 1 immersed in air and nodule 2 immersed in fluid; (e and j) the SV image produced with grid (A), that is, with 512 PSFs; BBR result computed from the grid (A) SV image using (f and k) grid (A), (g and l) grid (B), and (h and m) grid (C). Results shown after 100 iterations of the BBR algorithm. Axial profiles through the center of each nodule, its SV-image and restoration obtained with grid (A): (n) nodule 1 and (o) nodule 2. Scale bar: $10 μ m$ . Imaging parameters: $20 \times / 0.5 NA$ air lens ( $RI = 1.00$ ); emission wavelength is 635 nm. All images are normalized and displayed on the same intensity scale.

The BB forward model [Eq. (5)] was used to generate the SV image [Figs. 6(e) and 6(j)]. In the SV image [Figs. 6(e) and 6(j)], 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:

1. 343 ( $= 7 \times 7 \times 7$ ) blocks and 512 PSFs [grid A, Fig. 6(a)]. Grid (A) contains all the blocks used to compute the SV image of the lung tissue phantom.
2. 64 ( $= 4 \times 4 \times 4$ ) blocks and 125 PSFs [grid B, Fig. 6(b)]. Grid (B) contains 18% of the number of blocks used to compute the SV image of the phantom.
3. 8 ( $= 2 \times 2 \times 2$ ) blocks and 27 PSFs [grid C, Fig. 6(c)]. Grid (C) contains only 2% of the number of blocks used to compute the SV image of the phantom.

$X Y$ and $X Z$ 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 $X Z$ 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 $X Z$ 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.

Table 2

The effect of change in the information available to the BBR algorithm, quantified using NMSE, SSIM, and I-div.

Restoration	% of blocks used	NMSE	SSIM	I-div
Grid A	100	0.0130	0.9588	$1.47 \times 10^{7}$
Grid B	18	0.1747	0.7688	$2.22 \times 10^{13}$
Grid C	2	0.3245	0.6840	$6.72 \times 10^{19}$

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.

6. 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 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³⁰ 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.

Acknowledgments

We thankfully acknowledge our collaborator Lutz Schaefer (Advanced Imaging Methodology Consultation, Canada) for allowing us to extend and use in this study his MATLAB^® code for the SI-CG and DV-CG image restoration algorithms. We would also like to thank Sharon V. King (University of Memphis) for her useful comments and suggestions during this project. This work was supported by the National Science Foundation (CAREER award DBI-0844682, IDBR awards DBI-0852847, and DBI-1353904; PI: C. Preza) and the University of Memphis.

References

1.

D. A. Agard, “Optical sectioning microscopy: cellular architecture in three dimensions,” Annu. Rev. Biophys. Bioeng., 13 (1), 191 –219 (1984). http://dx.doi.org/10.1146/annurev.bb.13.060184.001203 ABPBBK 0084-6589 Google Scholar

2.

A. Egner, S. W. Hell, “Aberrations in confocal and multi-photon fluorescence microscopy induced by refractive index mismatch,” Handbook Of Biological Confocal Microscopy, 404 –413 Springer, US (2006). Google Scholar

3.

J. W. Goodman, Introduction to Fourier Optics, McGraw Hill, New York (1960). Google Scholar

4.

J. G. McNally et al., “Three-dimensional imaging by deconvolution microscopy,” Methods, 19 (3), 373 –385 (1999). http://dx.doi.org/10.1006/meth.1999.0873 MTHDE9 1046-2023 Google Scholar

5.

W. C. Chew, Waves and Fields in Inhomogeneous Media, IEEE Press, New York (1995). Google Scholar

6.

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,” JOSA A, 8 (10), 1601 –1613 (1991). http://dx.doi.org/10.1364/JOSAA.8.001601 Google Scholar

7.

P. Török and P. Varga, “Electromagnetic diffraction of light focused through a stratified medium,” Appl. Opt., 36 (11), 2305 –2312 (1997). http://dx.doi.org/10.1364/AO.36.002305 APOPAI 0003-6935 Google Scholar

8.

C. E. Perlman and J. Bhattacharya, “Alveolar expansion imaged by optical sectioning microscopy,” J. Appl. Physiol., 103 (3), 1037 –1044 (2007). http://dx.doi.org/10.1152/japplphysiol.00160.2007 Google Scholar

9.

C. E. Perlman, D. J. Lederer and J. Bhattacharya, “Micromechanics of alveolar edema,” Am. J. Respir. Cell Mol. Biol., 44 (1), 34 –39 (2011). http://dx.doi.org/10.1165/rcmb.2009-0005OC AJRBEL 1044-1549 Google Scholar

10.

D. S. Biggs, “3D deconvolution microscopy,” Curr. Protoc. Cytometry, 12 12.19. 1 –12.19. 20 (2010). http://dx.doi.org/10.1002/0471142956.cy1219s52 Google Scholar

11.

C. Preza and V. Myneni, “Quantitative depth-variant imaging for fluorescence microscopy using the COSMOS software package,” Proc. SPIE, 7570 757003 (2010). http://dx.doi.org/10.1117/12.840688 PSISDG 0277-786X Google Scholar

12.

M. M. Rahman et al., “Comparison of computational methods developed to address depth-variant imaging in fluorescence microscopy,” Proc. SPIE, 8589 85890A (2013). http://dx.doi.org/10.1117/12.2005227 Google Scholar

13.

C. Preza and J. Conchello, “Depth-variant maximum-likelihood restoration for three-dimensional fluorescence microscopy,” JOSA A, 21 (9), 1593 –1601 (2004). http://dx.doi.org/10.1364/JOSAA.21.001593 Google Scholar

14.

S. Ben Hadj et al., “Blind restoration of confocal microscopy images in presence of a depth-variant blur and poisson noise,” in IEEE Int. Conf. on Acoustics Speech and Signal Processing (ICASSP 2013), 915 –919 (2013). http://dx.doi.org/10.1109/ICASSP.2013.6637782 Google Scholar

15.

H. Trussell and B. Hunt, “Image restoration of space variant blurs by sectioned methods,” in IEEE Int. Conf. on Acoustics, Speech, and Signal Processing, (ICASSP’78), 196 –198 (1978). http://dx.doi.org/10.1109/ICASSP.1978.1170472 Google Scholar

16.

A. V. Oppenheim, R. W. Schafer and J. R. Buck, Discrete-Time Signal Processing, Prentice-hall, Englewood Cliffs (1989). Google Scholar

17.

T. P. Costello and W. B. Mikhael, “Efficient restoration of space-variant blurs from physical optics by sectioning with modified Wiener filtering,” Digital Signal Process., 13 (1), 1 –22 (2003). http://dx.doi.org/10.1016/S1051-2004(02)00004-0 DSPREJ 1051-2004 Google Scholar

18.

S. Weddell and R. Webb, “The restoration of extended astronomical images using the spatially-variant point spread function,” in 23rd Int. Conf. Image and Vision Computing New Zealand (IVCNZ 2008), 1 –6 (2008). Google Scholar

19.

H. Adorf, “Towards HST restoration with a space-variant PSF, cosmic rays and other missing data,” The Restoration of HST Images and Spectra-II, 72 Space Telescope Science Institute, Baltimore (1994). Google Scholar

20.

J. G. Nagy and D. P. O’Leary, “Restoring images degraded by spatially variant blur,” SIAM J. Sci. Comput., 19 (4), 1063 –1082 (1998). http://dx.doi.org/10.1137/S106482759528507X SJOCE3 1064-8275 Google Scholar

21.

E. Maalouf, “Contribution to fluorescence microscopy, 3D thick samples deconvolution and depth-variant PSF,” Thesis in Engineering Sciences [physics], Université de Haute Alsace—Mulhouse, (2010). Google Scholar

22.

S. Ben Hadj et al., “Blind restoration of confocal microscopy images in presence of a depth-variant blur and poisson noise,” in IEEE Int. Conf. on Acoustics, Speech, and Signal Processing (ICASSP 2013), 915 –919 (2013). http://dx.doi.org/10.1109/ICASSP.2013.6637782 Google Scholar

23.

S. Ghosh and C. Preza, “A block-based forward imaging model for improved sample volume representation in computational optical sectioning microscopy,” Proc. SPIE, 9330 93300T (2015). http://dx.doi.org/10.1117/12.2077001 Google Scholar

24.

J. B. Tenenbaum, V. de Silva and J. C. Langford, “A global geometric framework for nonlinear dimensionality reduction,” Science, 290 (5500), 2319 –2323 (2000). http://dx.doi.org/10.1126/science.290.5500.2319 SCIEAS 0036-8075 Google Scholar

25.

M. Arigovindan et al., “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). http://dx.doi.org/10.1364/OE.18.006461 OPEXFF 1094-4087 Google Scholar

26.

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). http://dx.doi.org/10.1117/12.909638 PSISDG 0277-786X Google Scholar

27.

N. Patwary and C. Preza, “Computationally tractable approach to PCA-based depth-variant PSF representation for 3D microscopy image restoration,” (2014). Google Scholar

28.

N. Patwary and C. Preza, “Image restoration for three-dimensional fluorescence microscopy using an orthonormal basis for efficient representation of depth-variant point-spread functions,” Biomed. Opt. Express, 6 (10), 3826 –3841 (2015). http://dx.doi.org/10.1364/BOE.6.003826 BOEICL 2156-7085 Google Scholar

29.

L. Schaefer, D. Schuster and H. Herz, “Generalized approach for accelerated maximum likelihood based image restoration applied to three dimensional fluorescence microscopy,” J. Microsc., 204 (2), 99 –107 (2001). http://dx.doi.org/10.1046/j.1365-2818.2001.00949.x JMICAR 0022-2720 Google Scholar

30.

S. Ghosh et al., “Further developments in addressing depth-variant 3D fluorescence microscopy imaging,” Proc. SPIE, 8949 89490Q (2014). http://dx.doi.org/10.1117/12.2045856 Google Scholar

31.

D. E. O’Donell, S. M. Revill and K. A. Webb, “Dynamic hyperinflation and exercise intolerance in chronic obstructive pulmonary disease,” Am. J. Respir. Crit. Care Med., 164 (5), 770 –777 (2001). http://dx.doi.org/10.1164/ajrccm.164.5.2012122 AJCMED 1073-449X Google Scholar

32.

S. Ghosh and C. Preza, “Fluorescence microscopy point-spread function model accounting for aberrations due to refractive index variability within a specimen,” J. Biomed. Opt., 20 (7), 075003 (2015). http://dx.doi.org/10.1117/1.JBO.20.7.075003 JBOPFO 1083-3668 Google Scholar

33.

S. Ghosh and C. Preza, “Space-variant image formation for 3D fluorescence microscopy using a computationally efficient block-based model,” in IEEE Int. Symp. on Biomedical Imaging (ISBI), (2015). http://dx.doi.org/10.1109/ISBI.2015.7163990 Google Scholar

34.

Y. Park et al., “Refractive index maps and membrane dynamics of human red blood cells parasitized by Plasmodium falciparum,” Proc. Natl. Acad. Sci. United States of America, 105 (37), 13730 –13735 (2008). http://dx.doi.org/10.1073/pnas.0806100105 PNASA6 0027-8424 Google Scholar

35.

C. L. Curl et al., “Refractive index measurement in viable cells using quantitative phase‐amplitude microscopy and confocal microscopy,” Cytometry Part A, 65 (1), 88 –92 (2005). http://dx.doi.org/10.1002/(ISSN)1552-4930 1552-4922 Google Scholar

36.

A. Golabchi et al., “Refractive errors and corrections for OCT images in an inflated lung phantom,” Biomed. Opt. Express, 3 (5), 1101 –1109 (2012). http://dx.doi.org/10.1364/BOE.3.001101 BOEICL 2156-7085 Google Scholar

37.

F. P. Bolin et al., “Refractive index of some mammalian tissues using a fiber optic cladding method,” Appl. Opt., 28 (12), 2297 –2303 (1989). http://dx.doi.org/10.1364/AO.28.002297 APOPAI 0003-6935 Google Scholar

38.

E. Barone‐Nugent, A. Barty and K. Nugent, “Quantitative phase‐amplitude microscopy I: optical microscopy,” J. Microsc., 206 (3), 194 –203 (2002). http://dx.doi.org/10.1046/j.1365-2818.2002.01027.x JMICAR 0022-2720 Google Scholar

39.

S. V. King, M. S. Hossain and C. Preza, “Dual acquisition of fluorescence and quantitative phase microscopy with high-speed switchable optics for DIC,” (2015). Google Scholar

40.

M. S. Hossain, S. King and C. Preza, “An integrated approach to determine prior information for improved wide-field imaging models from computational interference microscopy,” Proc. SPIE, 9330 933014 (2015). http://dx.doi.org/10.1117/12.2077322 Google Scholar

41.

J. G. Nagy et al., “Space-varying restoration of optical images,” JOSA A, 14 (12), 3162 –3174 (1997). http://dx.doi.org/10.1364/JOSAA.14.003162 Google Scholar

42.

M. Bertero and P. Boccacci, Introduction to Inverse Problems in Imaging, CRC Press(1998). Google Scholar

43.

Z. Wang et al., “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. Image Process., 13 (4), 600 –612 (2004). http://dx.doi.org/10.1109/TIP.2003.819861 IIPRE4 1057-7149 Google Scholar

44.

I. Csiszár, “I-divergence geometry of probability distributions and minimization problems,” Ann. Probab., 3 (1), 146 –158 (1975). Google Scholar

45.

D. Shier, J. Butler and R. Lewis, Hole’s Human Anatomy, McGraw-Hill, New York (1996). Google Scholar

46.

H. Kitaoka et al., “A 4-dimensional model of the alveolar structure,” J. Physiol. Sci., 57 (3), 175 –185 (2007). http://dx.doi.org/10.2170/physiolsci.RP000807 Google Scholar

Biography

Sreya Ghosh received her PhD in computer engineering from the University of Memphis in August 2015. She is a postdoctoral associate at the National Center for Microscopy and Imaging Research at the University of California San Diego, working on computational imaging for microscopy, since January 2016. Her research interests include signal processing, biophysical modeling, and inverse problems.

Chrysanthe Preza is an associate professor in the Department of Electrical and Computer Engineering at the University of Memphis, since 2006. She received her DSc degree in electrical engineering from Washington University in St. Louis in 1998. She leads the research in the Computational Imaging Research Laboratory at the University of Memphis. Her research interests are imaging science and estimation theory with applications in computational optical sensing and imaging applied to multidimensional multimodal light microscopy.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Sreya Ghosh and Chrysanthe Preza "Three-dimensional block-based restoration integrated with wide-field fluorescence microscopy for the investigation of thick specimens with spatially variant refractive index," Journal of Biomedical Optics 21(4), 046010 (28 April 2016). https://doi.org/10.1117/1.JBO.21.4.046010

Published: 28 April 2016

Access the abstract

JOURNAL ARTICLE
14 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 6 scholarly publications.

Explore citations on Lens.org

KEYWORDS

3D modeling

Point spread functions

3D image processing

Luminescence

Lung

Principal component analysis

Tissues

Show All Keywords

CHORUS Article. This article was made freely available starting 28 April 2017

1.

Introduction

2.

Block-Based Restoration Method

2.1.

Block-Based Forward Imaging Model

Eq. (1)

Fig. 1

Eq. (2)

Eq. (3)

Eq. (4)

Fig. 2

Eq. (5)

2.2.

Block-Based Restoration Method

Eq. (6)

Eq. (7)

Eq. (8)

3.

Evaluation of the Block-Based Restoration Method with Simulated Images

3.1.

Spatially Variant Imaging Using the Block-Based Forward Model

3.2.

Comparison of Simulated Images Using Different Forward Imaging Models

3.3.

Restoration of Simulated Spatially Variant Images

Fig. 3

Table 1

4.

Application of the Block-Based Restoration Method to Experimentally Acquired SV-Images

4.1.

Experimental Acquisition of Space-Variant Images

Fig. 4

4.2.

Comparison Between Block-Based Model Predictions and Experimental Images

4.3.

Restoration of Space Variant Images Using the Block-Based Restoration Method

Fig. 5

5.

Application of Block-Based Restoration to Space-Variant Imaging of a Simulated Lung Tissue Phantom

5.1.

Modeling the Lung Tissue Phantom

5.2.

Application of the Block-Based Restoration Method to the Lung Phantom

Fig. 6

Table 2

6.

Summary and Conclusions

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years