## Abstract

We present a method for modeling image formation in optical projection tomographic microscopy (OPTM) using high numerical aperture (NA) condensers and objectives. Similar to techniques used in computed tomography, OPTM produces three-dimensional, reconstructed images of single cells from two-dimensional projections. The model is capable of simulating axial scanning of a microscope objective to produce projections, which are reconstructed using filtered backprojection. Simulation of optical scattering in transmission optical microscopy is designed to analyze all aspects of OPTM image formation, such as degree of specimen staining, refractive-index matching, and objective scanning. In this preliminary work, a set of simulations is performed to examine the effect of changing the condenser NA, objective scan range, and complex refractive index on the final reconstruction of a microshell with an outer radius of 1.5 μm and an inner radius of 0.9 μm. The model lays the groundwork for optimizing OPTM imaging parameters and triaging efforts to further improve the overall system design. As the model is expanded in the future, it will be used to simulate a more realistic cell, which could lead to even greater impact.

© 2012 Optical Society of America

## 1. INTRODUCTION

Optical projection tomographic microscopy (OPTM) is a novel three-dimensional (3D) imaging technology being developed at the University of Washington, Arizona State University, and Cornell University for the early detection and diagnosis of disease from cellular samples, such as lung cancer [1]. Motivation for developing OPTM is provided, in part, by a threefold reduction in false-negative rates for diagnosing cancer in a controlled comparison of 3D versus two-dimensional (2D) images of the same cells using automated classifier analysis [2]. For individual and isolated cell analysis, VisionGate Inc. (Phoenix, Arizona) is commercializing OPTM for early disease detection using sputum as the cellular sample. The instrument is referred to as the optical Cell-CT (trademark of VisionGate) because 3D cell images are reconstructed from 2D optical projections analogous to that of X-ray CT [3].

OPTM forms isotropic 3D images of single cells embedded in optical gel that flows through a rotating capillary tube. Cells are imaged in transmission using bright-field optical microscopy techniques allowing standard absorptive stains (e.g., hematoxylin) used in almost all clinical diagnostics to be used to form 3D images. Thus, OPTM loosely combines techniques from radiology (CT), analytical chemistry (flow cytometry), and cytopathology (bright-field microscopy) [3]. Projections are acquired as cells flow past the objective one by one in a rotating capillary tube. Each projection is formed by scanning the focal plane of a high numerical aperture (NA) objective through the object of interest, which for this paper is assumed to be a single cell.

While OPTM has already shown its diagnostic capabilities [2], it could have a greater impact if it were transformed into a quantitative device where pixel values are directly correlated with specimen properties. Clinicians would no longer qualitatively decide whether disease is present or not; instead a quantitative metric would diagnose disease, which should lead to higher sensitivity and specificity. Quantitative microscopy has been proposed previously [4], but the large number of variables associated with both the microscope and the sample make this a difficult proposition. For example, filament intensity variations, ambient light, and staining variations highlight just a few of the many variables that affect the final image. Rather than attacking this problem experimentally, computational modeling provides the opportunity to precisely vary individual parameters, which would otherwise be difficult to decouple from one another.

Recently, the role of simulations in developing quantitative microscopy has gained traction among different research groups [5] leading to many publications. For example, simulations played a fundamental role in correlating scattering coefficients with refractive index (RI) in partial-wave spectroscopy [6], analyzing spectral variations in microspheres [7], and analyzing scattering produced by a biological cell from a focused Gaussian beam [8]. In each case, the simulations provided increased understanding of the underlying results that would otherwise have been difficult to interpret.

This work aims to begin providing similar analysis by building a computational model of the Cell-CT where improved interpretation of light scattering, image reconstruction, and image-analysis techniques will eventually lead to quantitative measurements. Our simulation will be used to evaluate the effects of diffraction, refraction, and absorption on the reconstruction of 3D images of standards, such as nanoparticles and microparticles. The theoretical model will also provide a framework to quantify system errors and inconsistencies, including nonlinear scanning of the objective, hysteresis in focal-plane scanning, and cell motion. This is the first computational model of a high-NA optical microscope for the field of optical projection tomography.

#### A. Optical Projection Tomographic Microscopy

The 3D OPTM images are reconstructed from standard 2D bright-field microscope images, acquired over multiple projections by rotating a microcapillary around its central axis. A simplified diagram is shown in Fig. 1. Cells are fixed, stained, and immersed in a solution of optical gel, and pushed through a microcapillary that passes between a condenser and objective. Once a cell is within the field of view (FOV), flow is arrested so that projections can be acquired. At each step in rotation, the focal plane is scanned through the cell, leaving the camera (CCD) shutter open to obtain a single optical projection by optically integrating through the entire cell. The microcapillary (and cell) is rotated to the next position and another projection is acquired. After projections are acquired from all desired angles, a 3D image is reconstructed using the filtered backprojection method. A total of 500 projections are typically acquired over 360° equating to a rotational step of 0.72° between projection acquisitions. The objective focal-plane scanning is determined *a priori* in accordance with the size of the cells being imaged. For example, lung epithelial cells are approximately 10 μm in diameter, so the scan distance is set to 12 μm. The microcapillary has a similar RI to the optical gel and index-matching immersion fluid such that its cylindrical geometry is transparent and nonrefractive [3].

#### B. Modeling

A computational model of light propagation through a cell within the microcapillary is derived to improve imaging parameters in OPTM. Computational modeling of the imaging process associated with OPTM is complicated by the use of high-NA condensers and objectives, which preclude the use of most theoretical simplifications. For example, Fourier optics techniques, such as Fresnel and Fraunhoffer diffraction, are only valid for paraxial applications [9]; furthermore, such models ignore the effect of refraction, multiple scattering [10], and polarization. The most accurate method of modeling microscopes requires directly solving Maxwell’s curl equations. While analytical solutions are preferred, solutions are only available for spheres and cylinders using Mie theory [11,12], which cannot be extended to arbitrary objects, such as a cell.

A variety of methods are available to numerically solve Maxwell’s equations, including the finite-difference time-domain (FDTD) method, the method of moments, and finite-element modeling. Each of these methods is computationally limited, in that it is impossible to propagate light waves all the way from the condenser to an image detector due to the distance per number of wavelengths [13]. To compensate, several groups have designed computational models that break image formation into multiple steps whereby different techniques can be used over the course of propagation to reduce the computational load. Generally, these steps cover four main aspects: illumination [14], near-field scattering [15], resampling [16], and image formation [17]. This four-step process allows near-field scattering to be computed using a computationally intensive numerical method from which the final image is computed. Our model supplements the original illumination technique with one that provides a more straightforward method of performing transmission imaging [14].

Our model is also differentiated by a number of programmatic elements that provide greater computational expediency and accuracy. These include precomputing coefficients necessary for calculating near-field scattering for improved computational cost and reduced numerical dispersion [18]. Additionally, the resampling technique used for this work has produced more accurate results compared to previous methods [16]. The 2D model is also extended to produce simulated 3D images similar to the Cell-CT using microshells and idealized conditions for improved computational efficiency.

## 2. THEORY

The theory of OPTM image formation and its implementation are provided below, where the cell can be of arbitrary size, shape, and material properties and they can sit in media with any RI. Since the microcapillary and slide/coverslip provide near refraction-free imaging, they are ignored for this analysis [3]. Future research could be performed to include their minor effects in the model. A general diagram of the modeling components used in image formation is illustrated in Fig. 2. The illumination is introduced into a numerical method to compute the near-field scattering (scattering is used broadly to encompass the effects of refraction, reflection, and absorption) around the object, which is then resampled to determine the scattering that would be seen on an infinite plane in the direction of the condenser. The scattering calculated on this plane is mapped to the final image. Since each of these four steps is well described elsewhere [14,19], they are only summarized here where two of the steps are implemented in an alternative way to produce a more accurate and computationally effective algorithm. Each part is linked together in Section 2.E to fully understand the entire process required to simulate a projection. The projections are finally reconstructed to form a 3D representation of a microshell, which could represent an isolated nuclear membrane from a cell.

#### A. Illumination

Illumination in a high-NA bright-field microscope is normally partially coherent Köhler illumination where point sources on the back focal plane of the objective result in plane waves near the focal plane in the image space. In order to simulate partial coherence, scattering due to a large number of orthogonal sets of plane waves is computed using separate simulations for each plane wave and the intensity from each simulation is summed to produce a final image [7]. While this method would produce the best possible results, the computational cost of this technique is quite high and becomes impractical when the objective focal plane is scanned as is performed in OPTM. Instead, our model assumes a coherent, polarized focused pulse, which is scanned throughout the object space, and the intensity resulting from the scattered light at each scan position is integrated to form a final image. This assumption reduces the computational cost of one simulation, described below, by more than eight times to approximately 24 h using an eight-core workstation.

Beyond the computational cost considerations, modeling using coherent illumination is supported experimentally by the use of an incidence filter designed to emphasize absorption of hematoxylin bound to DNA. The illumination filter is a 60 nm bandpass filter centered at 588 nm, which boosts the coherence length to approximately five times that of standard white light, or 1.75 μm. While some of the objects that will be simulated are a little larger than 1.75 μm, the use of a coherent pulse will provide results with accuracy high enough for our application such that coherent waves will approximate partial coherence for the illumination bandpass.

Illumination is mathematically derived in the time domain using coherent, polarized light, which is focused through a thin aplanatic converging lens [14]. Light rays polarized in the $x$-direction travel parallel to the optical axis, interacting with the lens before being diffracted by the exit pupil. The direction of each ray after interacting with the lens is determined by geometric optics. The rays describe the direction and polarization of the plane waves propagating proximal to the focus. The resulting electric field around the focal plane is a superposition of plane waves,

The focused pulse represented by Eq. (1) can be approximated discretely using numerical integration methods. While a 2D integration would be preferred, it is often difficult to implement these techniques. Instead, this work follows the previously defined one-dimensional (1D) technique, in which integration is performed using extended midpoint. Capoglu *et al.* use spherical variables of integration to discretize the aperture [14], leading to a large number of plane waves propagating from the center of the aperture and an inadequate number of plane waves propagating from the outer circle of the aperture sampling. The only way to sample the outer dimensions of the aperture with enough rays in this scheme is to significantly oversample the inner portion of the aperture. Instead, an alternative procedure is performed where a Cartesian grid is placed on top of the circular aperture and the incident rays are introduced at equally spaced Cartesian intervals throughout the entire aperture. This procedure significantly reduces the oversampling performed using spherical decomposition.

The number of rays representing the focused pulse in the Cartesian grid is dictated by the sampling theorem such that there is no statistical correlation between two adjacent points on the sample plane. The so-called “mutual coherence function” [21] quantifies this correlation and provides a means to determine the appropriate anti-aliasing criteria. A closed-form solution provided in [22] is replicated here for completeness:

#### B. Numerical Method

The illumination is employed in the FDTD method as it provides a simple and elegant framework to solve Maxwell’s curl equations in the time domain. FDTD allows an entire broadband pulse to be simulated in a single run. A brief overview of FDTD is described below, while the reader is referred to Taflove and Hagness for a more detailed presentation of the technique [15]. Maxwell’s curl equations are often written

where $\epsilon $ is permittivity, $\mu $ is the permeability, $\sigma $ is the conductivity, $\mathit{E}$ is the electric-field vector, and $\mathit{H}$ is the magnetic-field vector. These equations show that the magnetic fields vary temporally with the spatial variations in the electric fields and*vice versa*. In this scheme it is important to understand that the relative permittivity and permeability are directly linked to the complex RI ($\tilde{n}=\sqrt{{\mu}_{r}{\tilde{\epsilon}}_{r}}$), where the real portion represents pure refraction and the complex portion represents absorption. At optical frequencies, the relative permeability is very close to one, so the complex RI is simplified to $\tilde{n}=\sqrt{{\tilde{\epsilon}}_{r}}$. The complex RI can also be written to explicitly show the absorption or extinction coefficient ($\kappa $), in $\tilde{n}=n+i\kappa $.

Discretization of the electric field is performed using the Yee cube in which the electric and magnetic fields are spatially and temporally located at half-step differences [23]. These spatial and temporal differences allow the electromagnetic fields to be approximated using finite differences. As an example, the ${E}_{z}$ field is discretized using finite differences:

The spatial and temporal time increments are tightly linked together by the Courant number, which is required to maintain numerical accuracy in explicit time-marching schemes used to solve partial differential equations. It is generally accepted that 20 points per wavelength are required to maintain accuracy where the error of finite differences is $O\left({\mathrm{\Delta}}^{2}\right)$ [15]. Since discretization of light waves at such a high resolution is computationally costly, the computational grid should contain as few grid points as possible. For this reason, light traveling in an infinite space is modeled using absorbing boundary conditions [15]. This work relies on the convolution perfectly matched layer (CPML) (a variation of Berenger’s PML [24]), which absorbs waves at the boundary of the grid with much greater accuracy compared to previous boundary conditions [25].

Electromagnetic waves are introduced into the FDTD grid as either point sources or plane waves. Plane waves are more applicable to this work due to the plane-wave decomposition of the illumination described previously [14]. Plane waves could be introduced using the pure scatter-field formulation, which reduces the amount of energy in the grid thereby reducing potential errors, but at a high computational cost. Alternatively, they can be introduced with arbitrary direction and polarization using the total-field scatter-field (TFSF) formulation [26]. Using this method the grid is separated into two regions: the total-field and the scatter-field. The total field is in the inner portion of the grid enclosing the cell and includes both the incident and scatter fields. The scatter region, on the other hand, only contains the scattered electromagnetic radiation. The two regions are separated by altering the fields at the total-field/scatter-field boundary based on the time step, wave direction, and polarization. The cell interacts with the illumination in the total field and the waves that are scattered propagate to the scatter field and are mapped to the final image.

#### C. Resampling

The scattered fields collected in FDTD are resampled in preparation for the final image-formation component, which requires an infinite plane that is not capable of being acquired directly from the finite FDTD lattice. Additionally, the high-NA microscope theory requires the plane to be close to the focal plane. Together, these two requirements necessitate the use of a near-to-near field (NTNF) transformation [19]. The NTNF transformation is capable of transforming near-field scattering inside the FDTD grid to arbitrary observation points outside the grid. As waves pass through an artificial surface that encloses the object of interest, they are weighted by their relative position to the observation point. While there are time-domain versions of this transformation, our implementation utilizes a frequency-domain version owing to its reduced storage and computational costs for a large number of observation points [15]. Concurrent with FDTD time stepping, a discrete Fourier transform (DFT) is performed for each desired frequency, meaning only two values (real and imaginary parts) need to be saved for each point on the surface of integration. Following completion of the FDTD time stepping, the points on the surface of integration can be weighted to determine the electromagnetic fields at any point outside the surface. In this work, the NTNF transformation uses the Stratton–Chu vectorial diffraction integrals defined as

Borrowing ideas from the Huygens’s sources, the mixed-surface technique acquires electric fields from their regular positions, but assumes they come from the magnetic field’s position during the weighting process. The opposite is also true where the magnetic fields are acquired from their regular position but assumed to be located on the electric-field surface. This idea of switching calculated and perceived locations reduces the computational and storage costs required for two separate magnetic fields [27,29]. The mixed-surface Stratton–Chu diffraction integrals are implemented in FDTD as a box or rectangle and calculated for each desired frequency. The resulting electromagnetic fields on the infinite plane are transformed through the objective to the detector space for each frequency.

#### D. Image Formation

Image formation is treated as a coherent superposition of equivalent magnetic dipoles using the Debye–Wolf integral [19]. The image from each point is computed separately and subsequently summed to determine the overall electric field at the detector. The *m-theory* of diffraction determines the contribution from each equivalent magnetic dipole across the infinite plane to find the electromagnetic field at any observation point on the infinite half plane extending to $z=\infty $. The electric field at point $P=({x}_{p},{y}_{p},{z}_{p})$ is a result of the superposition of spherical waves from points $Q=({x}_{Q},{y}_{Q},{z}_{Q})$ situated on surface ${S}_{0}$ [19]:

#### E. Simulating a Projection

The image produced from the description above is from a single focused pulse whose point-spread function is potentially not wide enough to illuminate the entire cell. As such, the focused pulse is scanned throughout the imaging plane interacting with the cell. The intensity image resulting from each scan position and wavelength is subsequently integrated to produce the final image at each focal plane. A projection is the result of optically integrating different focal planes as the objective is axially scanned through the cell. The entire illumination and objective scanning process required to simulate a projection is mathematically encapsulated as

Integration of Eq. (13) can be approximated using various discrete integration techniques such as extended midpoint and Gauss–Legendre. This work utilizes extended midpoint due to its ease of implementation. The radial and angular sampling intervals ($\mathrm{\Delta}r$ and $\mathrm{\Delta}\theta $) are dictated by divisibles of the Rayleigh lateral resolution of an optical microscope, which is defined as

Figure 4 illustrates Eq. (13) and the entire process required to compute a single projection utilizing the theory presented in Sections 2.A–2.D. Rather than running each simulation step a single time to produce an image, a projection is the result of running each step several times to acquire all the required data. All the data is subsequently integrated to produce a single projection. For OPTM, projections need to be produced for every desired perspective prior to reconstructing a 3D representation of the cell.

#### F. Three-Dimensional Image Reconstruction

Experimentally acquired projections are reconstructed using the filtered backprojection method based on the projection-slice theorem [30]. This reconstruction technique mimics CT where X-rays are assumed to travel in a linear direction and absorption is the major contrast mechanism while refraction is the minor contrast mechanism. In OPTM, there are some obvious inconsistencies with the underlying linear, nonrefraction assumptions used in the filtered backprojection, and thus our model provides the opportunity assess the range of NAs and scattering coefficients that the filtered backprojection is adequate for reconstruction.

The filtered backprojection can be performed for each position along the optical axis such that 1D projections are transformed to 2D slices rather than transforming 2D data to 3D. In this way, 1D Fourier transforms are taken on projections along the dimension associated with the angle of rotation. These 1D transforms are placed at the angle of rotation in a 2D matrix and the final slice is computed by taking the inverse 2D transform. The backprojection by itself acts as a low-pass filter, producing blurry reconstructions. To reduce the effect of the transformation, a filter is often applied to the 1D transformed data to eliminate low-frequency information and emphasize high-frequency information. OPTM currently utilizes a variation of the ramp function where higher frequency is amplified and then removed beyond a point to reduce aliasing. The filtered backprojection is described mathematically as

## 3. METHODS

The OPTM model is implemented in this paper for the purpose of understanding the effect of condenser NA and object complex RI on image reconstruction of purely scattered light. This particular work reduces the computational complexity of the preceding model by assuming a microshell at the condenser focal plane, which reduces the number of projections to only one since every rotation produces the exact same projection. Further, laterally scanning the illumination is reduced as scans at a single radial distance for various angles are approximated as pure rotations of the final intensity image. While polarization has a larger impact at high NA, reducing the accuracy of this assumption, its overall contribution is actually fairly small and is neglected in the theory of many high-NA applications [31]. Without the preceding assumptions, the computational cost, especially for high NA, becomes time intensive and beyond the computational capabilities of a single PC workstation.

The incident plane waves described in Eq. (1) are modeled as Gaussian modulated sine waves of the form

where $\tau =2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{fs}$, ${f}_{0}=4.283\times {10}^{14}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{Hz}$, and ${t}_{0}=9\tau $. These parameters result in a pulse centered at 588 nm, mimicking OPTM incident light, which is bandpass filtered (see Fig. 1) to look specifically at the hematoxylin stain. The focused pulses, ranging in NA from 0.1 to 1.3, illuminate microshells with an outer radius of 1.5 μm and inner radius of 1.35 μm. The microshells are situated in media similar to optical gel whose RI is set at 1.45, similar to that of glycerol [32]. The RI inside the microshell matches that of the surrounding media so the only source of contrast is the shell. The RI of the shell is similar to experimentally tested microspheres [33] with a complex RI of $1.50+j0.001$. Complex refractive indices of $1.45+j0.001$ and $1.50+j0.0$ are also simulated to better understand the individual impact of absorption and refraction on the final reconstruction.The microshells are executed using the staircasing method, which is a standard technique in FDTD where curved boundaries are represented as a series of blocklike steps. The FDTD grid spacing is set at 15 nm ($\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{\Delta}z$) and the simulation is run at 95% of the 3D Courant limit ($0.95/\sqrt{3}$) for 2048 time steps. The FDTD grid is terminated using an 11-cell CPML with the total grid dimensions of $360\times 240\times 240$ cells. One dimension is extended such that scanning is performed in a single grid by moving the microshell instead of the focused pulse. For the first scan position, the microshell is placed at one end of the grid and the focused pulse is aimed at the center of the microshell. The remaining scan positions are the result of focusing the illumination at the same position and moving the microshell at predefined intervals toward the extended grid dimension. The scanning intervals, found using the method described in Section 2.E, are set to a maximum of $3/4{R}_{\text{lateral}}$, which produced correlation coefficients in the final image stack above 0.99. The scan intervals are further dictated by the grid spacing such that the actual sampling intervals are 0.36 μm for NA 0.1, 0.27 μm for NA 0.5, 0.225 μm for NA 0.9, and 0.195 μm for NA 1.3. Each scan is extended to 1.8 μm in radial distance from the center of the microshell. The FDTD algorithm is implemented in parallel utilizing message passage interface (MPI) to transfer magnetic fields between the subcubes in the 3D topology. MPI is also used for the remainder of the implementation to parallelize computation among the available computer resources.

The illumination TFSF coefficients are precomputed for the entire FDTD time sequence where the superposition of each plane wave is performed concurrently to reduce the overall computer-memory burden. Once the TFSF coefficients are computed for a specific grid size, the FDTD algorithm can be run in the time required to introduce a single plane wave. This has an obvious advantage for laterally scanning the microshell. Additionally, the TFSF coefficients are computed using the analytical field propagation (AFP) method, which offers significant advantages when dealing with numerical dispersion related to introducing plane waves at angles close to the Cartesian axes [18]. Numerical dispersion is an inherent effect in FDTD where waves travel at slower than desired velocities due to the number of cells per wavelength with which they are discretized and the angle that they travel with respect to the grid axes [15]. AFP requires a significant upfront computational demand, particularly when performed for a superposition of such a large number of plane waves, but the numerical dispersion advantages outweigh the computational cost disadvantage. TFSF coefficients are periodically loaded into memory throughout the FDTD algorithm to reduce the memory burden on the computer system. Most 3D implementations of the TFSF technique rely on either the single frequency compensation [15] or the matched numerical dispersion [34] techniques. These techniques are not as accurate as AFP [15], but they have reduced computational costs in situations where scanning relative to a focused pulse is not performed in a single grid.

The scattering is calculated for three different wavelengths spanning the 60 nm bandpass filter centered at 588 nm. There is no significant variation between the scattering produced by the upper and lower bounds of the bandpass filter, justifying the limited number of frequencies. Specifically in one test, the two wavelengths produce a Pearson correlation coefficient above 0.99. This high degree of similarity reduces the computational cost significantly. The previously defined simulation parameters and the DFT dictate the interrogated wavelengths are 555.3, 581.7, and 610.8 nm.

The infinite plane, used in the resampling, is placed at 2.5 μm along the optical axis from the focus of the illumination and the width is set to 20 μm in both lateral directions with a sampling of 200 nm. Both the width and the sampling interval are found by computing the correlation coefficient of final reconstructions using subsequently smaller sampling intervals and larger widths. The final width and sampling interval are determined where the correlation coefficient for the final reconstruction is above 0.99 produced using windows with adjacently computed parameters. Depending on the application, a correlation coefficient below 0.99 could be afforded to improve computational efficiency, but visual inaccuracies were noticeable for out-of-focus illumination in the axial cross section with lower correlation coefficient standards.

The light on the infinite plane is imaged through an objective with the same parameters across all simulations with a 1.3 NA and $100\times $ magnification. The detector space has an RI of 1.0 and the focal plane is scanned axially through the object from $-5$ to 5 μm at 0.2 μm increments. The axial scanning interval was determined to maintain a correlation coefficient of 0.99 in the reconstruction. In the magnified detector space, the images are calculated from $-400$ to 400 μm at 8 μm increments mimicking pixel spacing in standard CCD and CMOS cameras. All image units are demagnified in the figures presented below.

As previously mentioned, the final images are approximated by reducing the illumination to a single radial direction and all other angles are approximated as rotations of the intensity image for that radial distance. Rather than using Eq. (15) to determine the scan position of the focused pulse, it is instead used to determine the rotation interval for the final intensity images; all rotations are performed in MATLAB (Mathworks, Natick MA) using the *imrotate* function with linear interpolation. The projection is the result of integrating every frequency, focal plane, and scan position [Eq. (13) and Fig. 4]. Reconstructions using the projections are performed with the MATLAB *iradon* function using the Ram–Lak filter where the simulated projection is replicated 500 times to mimic OPTM acquisition every 0.72° over the entire 360° of rotation.

## 4. RESULTS

The best point to start analyzing OPTM image formation is prior to axially summing the 3D image stack to produce a projection. A set of microshell images produced with varying condenser NAs and an idealized representation are shown in Fig. 5. The microshells have an extinction coefficient of 0.001 while their real RI is matched to the surrounding media. The resulting image cross sections are individually normalized.

The image stack produced using the 0.1 NA condenser is unable to focus on the microshell where there are two “hot spots” symmetrically located at approximately $\pm 5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$, well outside the actual dimensions of the microshell. The actual location of the microshell only contains interference fringes. Since these results provide such an inadequate means to focus on the microshell, they are ignored for reconstruction and further analysis. The 0.5 NA condenser provides a more effective means to focus on the microshell as the object shape is starting to resemble the idealized cross section; however, the axial dimensions are extended to approximately $\pm 3.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$. This extension is due directly to the large depth of field, approximated as three times the lateral resolution [i.e., Eq. (14)]. As the condenser NA gets larger, the depth of field gets smaller and the image stacks begin to better resemble the idealized representation. The 0.9 NA condenser extends the maximum intensity to approximately $\pm 1.8\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$ while the 1.3 NA condenser does not extend the microshell axially by any appreciable amount.

The 0.9 NA condenser also slightly distorts the out-of-focus edges of the microshell creating a boxy-looking representation at either edge. This distortion is the result of constructive and destructive interference by out-of-focus waves, which is verified using a thinner microshell with an inner radius of 1.35 um. Using this particular sized microshell reduces the effect material properties have on imaging where the illumination and object shape have greater importance. While these results are not shown, the thinner microshell causes distinct interference fringes where the 0.9 NA condenser produces a large destructive interference pattern through the optical axis center of the microshell. The fringes produced by the 1.3 NA condenser are greater in number, but smaller in size. These constructive and destructive fringes become less noticeable as the microshell increases in thickness due to increased scattering.

While a smaller depth of field provides better representation of the physical dimensions, it does a worse job of illuminating out-of-focus edges of the microshell. The 0.9 NA condenser reduces the out-of-focus microshell edge contrast to only 75%–80% of the in-focus plane contrast, while the out-of-focus portion of the microshell illuminated using the 1.3 NA condenser is only 60%–65% of the in-focus plane. These variations in intensity at different focal positions impact the reconstructions, as the projections are unequally weighted toward the portion of the object closer to the condenser focal plane. The fact that depth of field gets smaller with increasing NA is verified by looking at the maximum value from each plane in the axial cross section. Ideally the maximum should be constant throughout the outer radius dimensions of the microshell and then be completely zero outside the microshell; however, it is apparent that the maximum intensity begins to fall immediately as the objective scans away from the condenser focal plane.

Figure 5 also shows symmetric scattering above and below the focal plane in the pure absorption cross sections. As real RI is incorporated into the microshell, the image stack is no longer symmetric along the zero focal plane and it becomes imbalanced to one side, depending on the mismatch. This particular feature is seen more explicitly in the 0.5 NA condenser results of Fig. 6 where the maximum intensity on the right side of the microshell is only 90% of that on the left side. If the microshell has real RI less than the surrounding media, then light is focused beyond the microshell in the positive direction. Alternatively, a microshell real RI greater than the surrounding media causes the hot spot to be reflected in the negative direction. This same lensing effect is described using microspheres in the optical trap literature [35–37]. As in the optical trap work, the effective focal length (i.e., the location of the hot spot) is reduced with increasing NA [35]. These results seem to indicate that changes in real RI have a smaller effect on the image stacks produced with higher NA condensers than those produced with smaller NA condensers. This conclusion has been verified using both microspheres and microshells with real refractive indices ranging from 1.4 to 1.55, which encompass the majority of biologically relevant refractive indices. It is also important to note that for these simulations the pure absorption images are three orders of magnitude smaller than those with the real refraction, meaning the real refraction is providing greater contrast than the absorption.

Finally, the axial stacks of images are summed to produce projections, which are reconstructed using the filtered backprojection and displayed with the idealized microshell reconstruction. Figure 7 shows superimposed plots of the cross section of the reconstructed images with varying NA where the cross section is the same in all directions. Both the axial and lateral cross-sectional reconstructions are the exact same due to the symmetry in the reconstruction; thus the axial cross section is arbitrarily displayed. The 1.3 NA reconstructed image is also provided in Fig. 7 for the reader’s curiosity. All cross sections reach their maximum at approximately the same position, but the peak widths are smaller for increasing NA where the peaks better represent the dimensions of the microshell.

Interestingly, the 0.9 NA cross section reaches a lower center value than the 1.3 NA cross section due to better illumination of out-of-focus edges with the worse depth of field. These parameters mean there is a greater distinction of shell versus nonshell, despite the dimensions not being represented as well. If a larger shell were used, then the 0.9 NA condenser would still illuminate more of the shell, but not be able to illuminate the out-of-focus edges as well as this example. These results could justify altering the illumination used in OPTM.

One apparent downside of increasing NA is the excess energy outside the dimensions of the microshell reconstruction as seen in Fig. 7. This excess energy is due to the “blooming effect” seen at out-of-focus points in Fig. 5. The image stacks resulting from lower NA condensers are cylindrically shaped patterns, while the higher NA condensers produce hourglass-shaped scattering as a result of incident rays propagating from larger angles. The hourglass shape introduces energy outside the actual dimensions and changes in the projections depending on the sweeping distance. In Figure 8, reconstructions with NA of 1.3 are displayed where projections produced by sweeping from $-2$ to 2 μm are compared to those produced by sweeping from $-5$ to 5 μm. These images show that by limiting the amount of axial information included in the projections, there is a major impact on the overall reconstruction where there is less energy outside the microshell dimensions while the reconstruction inside the dimensions of the microshell is approximately the same. In the inner diameter of the microshell, there is also less energy by limiting the sweeping distance. The negative value in the $\pm 2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$ sweep is due to the Matlab assuming the maximum value is at the center of the projection, which is not the case for microshells. The projections do not contain any negative values.

## 5. DISCUSSION

The model presented above provides the ability to analyze OPTM image formation in a much more systematic way than can be performed experimentally. For example, where analyzing different condenser NA configurations would each require their own condenser, costing thousands of dollars, the model can predict results by simply changing a few parameters. Also, the model provides the ability to isolate various RI components to determine their overall affect, which would be otherwise impossible to do experimentally. The preliminary data offered above is synthesized into four takeaway points that should be utilized in the future for improved image reconstructions using OPTM. These include the following:

- 1. High-NA condensers preserve axial shape information more effectively when compared to low NA condensers.
- 2. High-NA condensers should be scanned in conjunction with the objective to equally weight all axial information of the object in the projections.
- 3. Refractive-index differences do not significantly affect reconstructions for data acquired using high NA condensers.
- 4. Improved scanning based on
*a priori*information is necessary for high-NA condensers to eliminate blooming effect.

The chosen microshell material properties indicate for these particular simulations pure refraction provides greater contrast compared to pure absorption. However, experimental implementation of OPTM reduces the effect of RI contrast by closely matching the RI of the immersion oil, microcapillary, and surrounding media to the absorptive stained cell in order to highlight the absorption contrast mechanism. Obviously cells contain organelles whose refractive indices are not all the same and therefore cannot be matched with the surrounding media. Fortunately, the use of a high-NA condenser seems to mitigate distortions, caused by pure refraction, in the reconstruction. The greater number of incident angles offsets or averages out scattering resulting from pure refraction. Smaller NA condensers, however, are easily scattered in one direction or another, altering the reconstruction to a greater extent. Objects with refractive indices much different than the surrounding media may distort this claim, but this should not be a problem for biological cells.

Seemingly one of the few downsides of using high-NA condensers in OPTM is due to blooming of out-of-focus light, which is outside the lateral dimensions of the object itself. The easiest way to combat this blooming effect is to only scan the axial distance of the object and not beyond the object dimensions as simulated in Fig. 8. This method would eliminate a large percentage of the energy outside of the object dimensions. Alternatively, each focal plane could be acquired as a separate image, which would be individually segmented leaving only the object before summing the stack of images to produce a projection. This method could produce even more accurate reconstructions, but it would require alternative software and hardware to deal with this additional data.

Beyond the application of OPTM, the model presented could be used to better understand many different microscopy applications, especially those based on coherent illumination where some of the assumptions used for this work would not be necessary. For example, the model could be altered to simulate transmission-mode optical coherence tomography microscopy or phase-contrast techniques, such as holography, for improved device implementation [38]. While OPTM tries to reduce the RI contrast mechanism, other techniques characterize cells by their organelle scattering. The model could be applied to these techniques to better understand the organelle originated light-scattering signature of cancerous and noncancerous cells under focused illumination. Previous work individually evaluated the impact of each organelle on scattering utilizing plane-wave illumination [39,40], which produces a much different scattering pattern compared to focused illumination. Other work has utilized focused illumination to show there are different scattering processes associated with cancerous and noncancerous cells [41], but it does not pinpoint each organelle’s direct impact on the scattering process. Application of this model could bridge the gap between these manuscripts and provide a greater understanding of the light-scattering properties associated with different organelles, both cancerous and noncancerous, for focused illumination. The results of this type of study could be applicable outside of microscopy in areas, such as flow cytometry, where cells are analyzed based on their scattering signature. The model presented above would require little alterations to simulate collection of scattered light produced by coherent, laser light incident on cells. Computational modeling has many applications throughout microscopy and general imaging where it provides the ability to better understand the impact of light absorption and scattering on image formation for improved image analysis and interpretation.

## 6. CONCLUSION

An algorithm is presented to examine OPTM, which is designed with the use of high-NA condensers and objectives. The algorithm is used to examine the effect of condenser NA and microshell complex RI as a first example. The results of such analysis showed the benefits of using high-NA condensers when it comes to representing the dimensions of the cell, but that high NA requires scanning the condenser in conjunction with the objective to better illuminate the entire cell. Also, the use of a high-NA condenser reduces the effect of scattering on the reconstruction. The one downside is the introduction of blooming, which can lead to misrepresentations of the cell and could require alternative implementations. In the future, the algorithm will be experimentally validated using microspheres with known refractive indices and radii. Also, future algorithm development will incorporate other effects, such as nonlinear scanning and incomplete microcapillary rotation, to understand how each one affects the entire image-formation process. From these results, we will be able to triage efforts to improve the image-formation process and enhance the diagnostic capabilities of the Cell-CT.

## ACKNOWLEDGMENTS

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-0718124. This work is also funded by the National Science Foundation Grant No. CBET-1014976. The authors would like to thank Anthony Reeves, Ph.D., for providing computational facilities at Cornell University, Michael Meyer at VisionGate Inc. (Phoenix, Arizona) for his valuable insight, and Paul Wiggins, Ph.D., for helpful comments. Cell-CT is a trademark of Visiongate Inc. (www.visiongate3d.com).

## REFERENCES

**1. **Q. Miao, A. P. Reeves, F. W. Patten, and E. J. Seibel, “Multimodal 3D imaging of cells and tissue: bridging the gap between clinical and research microscopy,” Ann. Biomed. Eng. **40**, 263–276 (2012). [CrossRef]

**2. **M. G. Meyer, M. Fauver, J. R. Rahn, T. Neumann, and F. W. Patten, “Automated cell analysis in 2D and 3D: a comparative study,” Pattern Recogn. **42**, 141–146 (2009). [CrossRef]

**3. **M. Fauver, E. J. Seibel, J. R. Rahn, M. G. Meyer, F. W. Patten, T. Neumann, and A. C. Nelson, “Three-dimensional imaging of single isolated cell nuclei using optical projection tomography,” Opt. Express **13**, 4210–4223 (2005). [CrossRef]

**4. **I. T. Young, “Quantitative microscopy,” IEEE Eng. Med. Biol. **15**, 59–66 (1996). [CrossRef]

**5. **N. N. Boustany, S. A. Boppart, and V. Backman, “Microscopic imaging and spectroscopy with scattered light,” Annu. Rev. Biomed. Eng. **12**, 285–314 (2010). [CrossRef]

**6. **H. Subramanian, H. K. Roy, P. Pradhan, M. J. Goldberg, J. Muldoon, R. E. Brand, C. Sturgis, T. Hensing, D. Ray, A. Bogojevic, J. Mohammed, J. Chang, and V. Backman, “Nanoscale cellular changes in field carcinogenesis detected by partial wave spectroscopy,” Cancer Res. **69**, 5357–5363 (2009). [CrossRef]

**7. **I. R. Çapoglu, C. A. White, J. D. Rogers, H. Subramanian, A. Taflove, and V. Backman, “Numerical simulation of partially coherent broadband optical imaging using the finite-difference time-domain method,” Opt. Lett. **36**, 1596–1598 (2011). [CrossRef]

**8. **C. Smithpeter, A. K. Dunn, R. Drezek, T. Collier, and R. Richards-Kortum, “Near real time confocal microscopy of cultured amelanotic cells: sources of signal, contrast agents and limits of contrast,” J. Biomed. Opt. **3**, 429–436 (1998). [CrossRef]

**9. **N. Nakajima, “Phase retrieval from a high-numerical-aperture intensity distribution by use of an aperture-array filter,” J. Opt. Soc. Am. A **26**, 2172–2180 (2009). [CrossRef]

**10. **M. Totzeck, “Numerical simulation of high-NA quantitative polarization microscopy and corresponding near-fields,” Optik **112**, 399–406 (2001). [CrossRef]

**11. **F. Slimani, G. Grehan, G. Gouesbet, and D. Allano, “Near-field Lorenz–Mie theory and its application to microholography,” Appl. Opt. **23**, 4140–4148 (1984). [CrossRef]

**12. **G. Gouesbet, “Generalized Lorenz–Mie theory and applications,” Part. Part. Syst. Charact. **11**, 22–34 (1994). [CrossRef]

**13. **P. Török, P. R. T. Munro, and E. E. Kriezis, “Rigorous near- to far-field transformation for vectorial diffraction calculations and its numerical implementation,” J. Opt. Soc. Am. A **23**, 713–722 (2006). [CrossRef]

**14. **I. R. Çapoglu, A. Taflove, and V. Backman, “Generation of an incident focused light pulse in FDTD,” Opt. Express **16**, 19208–19220 (2008). [CrossRef]

**15. **A. Taflove and S. C. Hagness,*Computational Electrodynamics: The Finite-Difference Time-Domain Method* (Artech House, 2005).

**16. **R. L. Coe and E. J. Seibel, “Improved near-field calculations using vectorial diffraction integrals in the finite-difference time-domain method,” J. Opt. Soc. Am. A **28**, 1776–1783 (2011). [CrossRef]

**17. **P. R. T. Munro and P. Török, “Calculation of the image of an arbitrary vectorial electromagnetic field,” Opt. Express **15**, 9293–9307 (2007). [CrossRef]

**18. **J. B. Schneider and K. Abdijalilov, “Analytic field propagation TFSF boundary for FDTD problems involving planar interfaces: PECs, TE, and TM,” IEEE Trans. Antennas. Propag. **54**, 2531–2542 (2006). [CrossRef]

**19. **P. Török, P. R. T. Munro, and E. E. Kriezis, “High numerical aperture vectorial imaging in coherent optical microscopes,” Opt. Express **16**, 507–523 (2008). [CrossRef]

**20. **B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems. II. Structure of the image field in an aplanatic system,” Proc. R. Soc. A **253**, 358–379 (1959). [CrossRef]

**21. **M. Born and E. Wolf, *Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light* (Pergamon, 1980).

**22. **I. R. Capoglu, J. D. Rogers, A. Taflove, and V. Backman, “The microscope in a computer: Image synthesis from three-dimensional full-vector solutions of Maxwell’s equations at the nanometer scale,” Prog. Opt.57, 1–91 (2012). [CrossRef]

**23. **K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. **14**, 302–307 (1966). [CrossRef]

**24. **J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. **114**, 185–200 (1994). [CrossRef]

**25. **J. A. Roden and S. D. Gedney, “Convolution PML (CPML): an efficient FDTD implementation of the CFS–PML for arbitrary media,” Microw. Opt. Technol. Lett. **27**, 334–339 (2000). [CrossRef]

**26. **D. Merewether, R. Fisher, and F. W. Smith, “On implementing a numeric Huygen’s source scheme in a finite difference program to illuminate scattering bodies,” IEEE Trans. Nucl. Sci. **27**, 1829–1833 (1980). [CrossRef]

**27. **T. Martin, “An improved near- to far-zone transformation for the finite-difference time-domain method,” IEEE Trans. Antennas Propag. **46**, 1263–1271 (1998). [CrossRef]

**28. **D. J. Robinson and J. B. Schneider, “On the use of the geometric mean in FDTD near-to-far-field transformations,” IEEE Trans. Antennas Propag. **55**, 3204–3211 (2007). [CrossRef]

**29. **T. Martin, “On the FDTD near-to-far-field transformations for weakly scattering objects,” IEEE Trans. Antennas Propag. **58**, 2794–2795 (2010). [CrossRef]

**30. **A. C. Kak and M. Slaney, *Principles of Computerized Tomographic Imaging* (IEEE, 1988).

**31. **M. G. L. Gustafsson, “Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy,” J. Microsc. **198**, 82–87 (2000). [CrossRef]

**32. **N. Martini, J. Bewersdorf, and S. W. Hell, “A new high-aperture glycerol immersion objective lens and its application to 3D-fluorescence microscopy,” J. Microsc. **206**, 146–151 (2002). [CrossRef]

**33. **X. Ma, J. Q. Lu, R. S. Brock, K. M. Jacobs, P. Yang, and X.-H. Hu, “Determination of complex refractive index of polystyrene microspheres from 370 to 1610 nm,” Phys. Med. Biol. **48**, 4165–4172 (2003). [CrossRef]

**34. **C. Guiffaut and K. Mahdjoubi, “Perfect wideband plane wave injector for FDTD method,” in *Proceedings of IEEE Antennas Propagation Society International Symposium* (IEEE, 2000), pp. 236–239.

**35. **T. H. Chow, W. M. Lee, K. M. Tan, B. K. Ng, and C. J. R. Sheppard, “Resolving interparticle position and optical forces along the axial direction using optical coherence gating,” Appl. Phys. Lett. **97**, 231113 (2010). [CrossRef]

**36. **J. P. Brody and S. R. Quake, “A self-assembled microlensing rotational probe,” Appl. Phys. Lett. **74**, 144–146 (1999). [CrossRef]

**37. **J. J. Schwartz, S. Stavrakis, and S. R. Quake, “Colloidal lenses allow high-temperature single-molecule imaging and improve fluorophore photostability,” Nat. Nanotechnol. **5**, 127–132 (2009). [CrossRef]

**38. **W. Choi, C. Fang-Yen, K. Badizadegan, S. Oh, N. Lue, R. R. Dasari, and M. S. Feld, “Tomographic phase microscopy,” Nat. Methods **4**, 717–719 (2007). [CrossRef]

**39. **R. Drezek, A. Dunn, and R. Richards-Kortum, “Light scattering from cells: finite-difference time-domain simulations and goniometric measurements,” Appl. Opt. **38**, 3651–3661 (1999). [CrossRef]

**40. **R. Drezek, A. Dunn, and R. Richards-Kortum, “A pulsed finite-difference time-domain (FDTD) method for calculating light scattering from biological cells over broad wavelength ranges,” Opt. Express **6**, 147–157 (2000). [CrossRef]

**41. **H. Subramanian, P. Pradhan, Y. Liu, I. R. Capoglu, J. D. Rogers, H. K. Roy, R. E. Brand, and V. Backman, “Partial-wave microscopic spectroscopy detects subwavelength refractive index fluctuations: an application to cancer diagnosis,” Opt. Lett. **34**, 518–520 (2009). [CrossRef]