Computational modeling of optical projection tomographic microscopy using the finite difference time domain method

Ryan L. Coe; Eric J. Seibel

doi:10.1364/JOSAA.29.002696

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

Fig. 1. Diagram of the image-formation process utilized in OPTM where projections are produced by scanning the objective to optically integrate through the cell. Projections are acquired from desired angles by rotating the microcapillary between the condenser and objective.

Download Full Size | PDF

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.

Fig. 2. Diagram of the image formation model for a high-NA bright-field, transmission microscope. The model consists of four main parts: illumination, numerical method, resampling, and image formation.

Download Full Size | PDF

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,

E (r, t) = \frac{f}{2 π c} \iint_{s} \frac{{\hat{e}}^{w} (θ^{'}, φ^{'})}{s_{z}^{'}} E (t^{'}) \cos^{1 / 2} (θ^{'}) d s_{x}^{'} d s_{y}^{'},

integrated over the domain

S = {(s_{x}, s_{y}) \in R^{2} | s_{x}^{2} + s_{y}^{2} \leq n \sin θ_{\max}}

, where

f

is the focal length of the condenser and

c

is the speed of light. The Cartesian coordinate variables of integration

s_{x}^{'}

,

s_{y}^{'}

, and

s_{z}^{'}

are directly matched with their spherical coordinate counterparts where

θ^{'} = \cos^{- 1} (s_{z}^{'})

and

φ^{'} = atan 2 (s_{y}^{'}, s_{x}^{'})

. The maximum

θ^{'}

is

θ_{\max}

, which is the maximum angle of intersection with the optical axis as defined by the NA (

NA = n \sin (θ_{\max})

). The unit vector

{\hat{e}}^{w} (θ^{'}, φ^{'})

is the polarization vector found using geometric optics that is directly related with the position vector of the ray interacting with condenser lens. The polarization angle is defined as the angle the electric vector makes with the meridonal plane. The

\cos^{1 / 2} (θ^{'})

factor is a result of the intensity law of geometric optics, dampening the intensity of waves coming from the edge of the lens as compared to those propagating a shorter distance from the center of the lens [20]. The dot above the incident electric field indicates the time derivative and the incident electric field,

E (t^{'})

, is a function used to describe the time variations in the electric field along the propagating plane wave. The plane waves can be introduced with any desired wavelength range using a Gaussian modulated sine wave.

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:

Δ s_{x} < \frac{2 π}{k \max {D, W_{c}}}, Δ s_{y} < \frac{2 π}{k \max {D, W_{c}}},

where

D

is the sample dimension,

W_{c}

is the value at which the mutual coherence function falls below a negligible value (

W_{c} = 1 / (k \sin θ_{\max}) = 1 / (k_{0} {NA}_{illum})

), and the smallest wavelength in the excitation dictates the wavenumber for the anti-aliasing condition.

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

- μ \frac{\partial H}{\partial t} = \nabla \times E,

σ E + ε \frac{\partial E}{\partial t} = \nabla \times H,

where

ε

is permittivity,

μ

is the permeability,

σ

is the conductivity,

E

is the electric-field vector, and

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{μ_{r} {\tilde{ε}}_{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{ε}}_{r}}

. The complex RI can also be written to explicitly show the absorption or extinction coefficient (

κ

), in

\tilde{n} = n + i κ

.

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:

E_{z}^{q + 1} [m, n, p + \frac{1}{2}] = \frac{1 - \frac{σ Δ_{t}}{2 ε}}{1 + \frac{σ Δ_{t}}{2 ε}} E_{z}^{q} [m, n, p + \frac{1}{2}] + \frac{1}{1 + \frac{σ Δ_{t}}{2 ε}} (\frac{Δ_{t}}{ε Δ_{x}} {H_{y}^{q + \frac{1}{2}} [m + \frac{1}{2}, n, p + \frac{1}{2}] - H_{y}^{q + \frac{1}{2}} [m - \frac{1}{2}, n, p + \frac{1}{2}]} - \frac{Δ_{t}}{ε Δ_{y}} {H_{x}^{q + \frac{1}{2}} [m, n + \frac{1}{2}, p + \frac{1}{2}] - H_{x}^{q + \frac{1}{2}} [m, n - \frac{1}{2}, p + \frac{1}{2}]}),

where

Δ_{t}

is the change in time per time step,

Δ_{x}

,

Δ_{y}

, and

Δ_{z}

are the spatial increments, and

m

,

n

, and

p

are the spatial indices. There are five similar equations to describe the other two electric-field components and all three magnetic-field components.

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 (Δ^{2})$ [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

E (r_{p}) = - ∯_{s} [j ω μ (\hat{n} \times H) G (r_{s}, r_{p}) - (\hat{n} \times E) \times \nabla G (r_{s}, r_{p}) - (\hat{n} \cdot H) \nabla G (r_{s}, r_{p})] d S^{'},

H (r_{p}) = ∯_{s} [j ω μ (\hat{n} \times E) G (r_{s}, r_{p}) + (\hat{n} \times H) \times \nabla G (r_{s}, r_{p}) + (\hat{n} \cdot E) \nabla G (r_{s}, r_{p})] d S^{'},

where

\hat{n}

is the outward surface normal, and

G (r_{s}, r_{p})

is the free-space Green’s function:

G (r_{s}, r_{p}) = \frac{e^{- j k | r_{s} - r_{p} |}}{4 π | r_{s} - r_{p} |},

| r_{s} - r_{p} | = \sqrt{{(x - x_{p})}^{2} + {(y - y_{p})}^{2} + {(z - z_{p})}^{2}},

with

r_{p} = (x_{p}, y_{p}, z_{p})

being the observation point and

r_{s} = (x, y, z)

being the coordinates of the fields on the surface of integration

S

. The spatial offsets among the electromagnetic fields make it impossible to obtain a single surface of integration directly from FDTD. Traditionally, a single surface is obtained by interpolating fields using the arithmetic mean. However, the accuracy of the arithmetic mean is strongly hindered by phase errors [27,28]. Though the geometric mean produces better results [28], an alternative mixed-surface method was implemented to dramatically reduce phase errors associated with the other two methods [16].

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_{img} (P) = \frac{- j k}{2 π} \iint_{S_{0}} (r \times (n \times E (Q))) \frac{\exp (- j k r)}{r^{2}} d S,

where

r = (x_{P} - x_{Q}, y_{P} - y_{Q}, z_{P} - z_{Q})

is a vector pointing from

Q

to

P

. This theory fits directly with the Debye–Wolf diffraction integral to determine the resulting image at the detector as a result of propagating through the objective. The effect of the objective is incorporated through expansion of generalized Jones matrices used to describe the transformation based on geometric optics simplifications. The image for each equivalent magnetic dipole,

n \times E (Q)

, is determined as though it sits on the optical axis and the lateral image position is subsequently altered using the actual dipole position and the magnification of the objective:

E_{img} (r_{d}) = \frac{- j k}{2 π} \iint_{S_{0}} E_{e q, o n -axis}^{n \times E^{i} (Q)} (r_{d} + β q) d S,

where

q

is a vector from the intersection of

S_{0}

and the

z

-axis to the point

Q

and

β

is the magnification. The total intensity or image at the detector is found from the electric field found at the detector:

I = {| E_{img} |}^{2} .

An image at each desired focal plane is a result of summing the contribution from each desired frequency. The process of computing images for each desired focal plane is simplified by altering the relative position of the resampled plane in relation to the focal plane of the objective. This means that only the final Debye–Wolf algorithm must be recomputed for each desired focal plane. Figure 3 illustrates the physical arrangement of the condenser and the objective where the condenser never changes position, meaning the scattered light calculated at the infinite plane does not change regardless of the position of the objective. The only difference is the distance between the infinite plane and the focus of the object. This idea is particularly important for modeling OPTM where projections are created by optically integrating through the cell.

Fig. 3. Representation of axially scanning the objective in comparison to the condenser illumination. Three different positions of the objective (on the right with red focal circles) are shown in relation to the condenser (on the left with black focal crosses) and the resampled infinite plane (vertical line). These diagrams show that the near-field scattered light produced by the object would be the same on the resampled infinite plane irrespective of the objective position. The only difference is the distance between the infinite plane and the focus of the objective (i.e., distance between red cross and infinite plane) meaning only the final image formation portion of the algorithm has to be repeated to axially scan through the object.

Download Full Size | PDF

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

P = \int_{z_{\min}}^{z_{\max}} \int_{λ_{\min}}^{λ_{\max}} \int_{0}^{2 π} \int_{0}^{r_{\max}} I (r, θ, λ) r d r d θ d λ d z,

where

(r, θ)

are the polar coordinates of the focused pulse scan position,

r_{\max}

is the maximum radial scan position dictated by the FOV and/or the cell size,

I

is intensity found due to each scan position using Eq. (12),

λ

is the wavelength interrogated between the minimum (

λ_{\min}

) and the maximum (

λ_{\max}

), and

z

defines the height of the objective focal plane between the minimum (

z_{\min}

) and the maximum (

z_{\max}

). The scan positions are implemented in Eq. (13) using a polar coordinate scheme due to the fact that our simulations employ radially symmetric objects, but they could be implemented using a Cartesian coordinate system.

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 ( $Δ r$ and $Δ θ$ ) are dictated by divisibles of the Rayleigh lateral resolution of an optical microscope, which is defined as

R_{lateral} = \frac{1.22 λ}{{NA}_{obj} + {NA}_{cond}} .

The actual sampling interval is found by altering the scan distance in divisibles of the Rayleigh criteria and comparing the correlation coefficient produced in the final image stack. The angular sampling interval (

Δ θ

) is also found using the lateral resolution by setting the radial scan interval as a chord on a circle whose radius is the maximum radial scan distance. The angle is found by relating the chord and the radius to find the maximum common denominator of 360° where noncommon denominators unequally weight rotations at the beginning/end of the rotations. Mathematically this relationship is described as

Δ θ = \frac{360}{⌈ \frac{360}{2 \sin^{- 1} (\frac{Δ r}{2 Δ r S_{N}})} ⌉} = \frac{360}{⌈ \frac{180}{\sin^{- 1} (\frac{1}{2 S_{N}})} ⌉},

where

S_{N}

is the maximum number of scans. Since the maximum radius is found by multiplying the scan interval and the number of scans, this relationship is simplified, eliminating the radial sampling interval.

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.

Fig. 4. Diagram of the process to simulate an OPTM projection using visual for-loops. The illumination is defined for the desired illumination source and NA, which is incorporated in a numerical method to determine the near-field scattering from a cell. The near-field scattering is resampled for a particular frequency and used to compute an image at every desired focal plane. The near-field scattering is then resampled for another desired frequency and images are computed for all of the same focal planes. This entire process is performed for every scanning position to illuminate the entire cell. Each of these images is finally integrated to produce a single projection. Projections are acquired for each rotation of the cell and these projections are reconstructed into a three-dimensional image in a final post-processing step (not shown).

Download Full Size | PDF

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

f (x, y, z) = \int_{0}^{2 π} \int_{0}^{\infty} H (ζ) G (ζ, θ, z) \exp (j 2 π ζ (x \cos θ + y \sin θ)) ζ d ζ d θ,

where the 1D Fourier transformed projection data is

G (ζ, θ)

and

H (ζ)

is the filter. The final reconstructed image is a 3D representation of the entire cell, which can be analyzed to understand the impact that each parameter has on the final reconstruction.

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

E^{i} (t) = \exp (- {((t - t_{0}) / τ)}^{2} / 2) \sin (2 π f_{0} (t - t_{0})),

where

τ = 2 fs

,

f_{0} = 4.283 \times 10^{14} Hz

, and

t_{0} = 9 τ

. 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 + j 0.001

. Complex refractive indices of

1.45 + j 0.001

and

1.50 + j 0.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 ( $Δ x = Δ y = Δ 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_{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.

Fig. 5. Axial cross section for varying NAs of scattered light with a constant microshell RI of $1.45 + j 0.001$ in media with a RI of 1.45. The idealized cross section is also provided for comparison purposes. Each cross section is individually normalized.

Download Full Size | PDF

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

Fig. 6. Axial cross sections produced using the 0.5 NA condenser (top row) and the 1.3 NA condenser (bottom row) for microshells with complex refractive indices of $1.45 + j 0.001$ (left column) and $1.5 + j 0.001$ (right column). The microshells are placed in media with a RI of 1.45. Each image is individually normalized, but the pure absorption microshells have approximately three orders of magnitude less contrast.

Download Full Size | PDF

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.

Fig. 7. Cross-sectional reconstruction plot (left) for varying NA for 0.5, 0.9, and 1.3 with a constant microshell RI of $1.5 + j 0.001$ . A 1.3 NA reconstruction (right) is also provided. Each NA intensity plot is individually normalized.

Download Full Size | PDF

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

Fig. 8. Cross-sectional reconstructions produced with a condenser whose NA is 1.3. The projection consists of axial scanning from $\pm 5$ or $\pm 2 μm$ for a microshell with a RI of 1.5. Each plot is individually normalized.

Download Full Size | PDF

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 first two points realize that the smaller depth of field of high-NA condensers seems to have both advantages and disadvantages, where it more accurately represents the object dimensions, but it also does not provide even illumination. This improved understanding appears to indicate that the small depth of field of a high-NA condenser requires that the condenser scan in conjunction with the objective in order to accurately depict larger objects. Without scanning the condenser, the object is not evenly illuminated, meaning the contribution to projections is disproportionately weighted toward those regions of the cell, which are closest to the condenser focal plane. Scanning the condenser in conjunction with the objective not only provides better illumination of the cell, but it might also be possible to reduce the scan range from 360° to 180°. The decision to scan over the entire 360 deg was made to average out differences between projections acquired with 180° of separation, which were believed to be a result of scattering. These results suggest that these differences are in fact more a result of disproportionate illumination rather than scattering, which has a smaller role in image reconstruction than initially believed. If these results continue to hold true it would play a vital role in improving projection acquisition time as only half the number of projections would be required. Future work must be done to optimize the amount and interval of rotation taking into account larger cell like objects, which could alter our scattering hypothesis.

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]

Computational modeling of optical projection tomographic microscopy using the finite difference time domain method

Abstract

1. INTRODUCTION

A. Optical Projection Tomographic Microscopy

B. Modeling

2. THEORY

A. Illumination

B. Numerical Method

C. Resampling

D. Image Formation

E. Simulating a Projection

F. Three-Dimensional Image Reconstruction

3. METHODS

4. RESULTS

5. DISCUSSION

6. CONCLUSION

ACKNOWLEDGMENTS

REFERENCES

Cited By

Figures (8)

Equations (17)

Journal of the Optical Society of America A