1. INTRODUCTION

Quantitative phase imaging (QPI) is an emerging field of biomedical optics in which the refractive index (RI) of phase objects is indirectly imaged through interferometric analysis [1]. QPI is preferred over alternative forms of biomedical imaging in that it is label-free. Thus, live cells can be imaged in their natural state without issues associated with photo-toxicity or photo-bleaching as in fluorescence microscopy. Another important benefit is that it is quantitative, enabling measurement of intrinsic properties, as opposed to modalities providing qualitative (nonlinear) information such as phase and differential interference contrast microscopy. The data’s quantitative nature also lends itself toward image processing, thereby improving the extractability of various features and properties [2]. In addition to biomedical applications, QPI techniques have been useful for a variety of other applications including adaptive optics, semiconductor defect inspection [3], and optical fiber characterization [4].

QPI can refer either to two-dimensional QPI (2D QPI), in which the 2D phase image is interpreted as the integrated optical path length (OPL) through the sample, or three-dimensional QPI (3D QPI), in which the real part of the object’s complex RI is imaged in all three spatial dimensions [5]. In spite of the fact that most objects, including biological specimens, are essentially 3D phase objects, 2D QPI has found wide-spread biomedical applicability [1,5–9]. For example, it is known that cellular dry mass is linearly related to OPL [6]. Therefore, QPI can be used to monitor cell growth as a function of cell-cycle [7]. Another example would be the study of red blood cell membrane dynamics [8], for which 2D QPI is appropriate since RI is essentially homogeneous.

In a general sense, however, single phase projections are insufficient for characterizing heterogeneous objects since, without approximation or special measures, it is impossible to differentiate OPL variations owing to changes in thickness versus RI. Thus, for samples containing complex internal structure, such as eukaryotic cells, 3D QPI is necessary for a complete morphological characterization. Current trends in QPI methodology reflect this need as incorporating tomography and 3D microscopy is a major focus area for research [1]. Although most 3D QPI research has centered on methodology, applications areas, such as the biophysical characterization of malarial parasite exit from human erythrocytes [10] and the quantification of chromosomal dry mass values for human colon cancer cells [11], are being developed in parallel.

Conventionally, 3D QPI is realized via either tomographic [12] or deconvolution methods [13–16]. In tomography, the object is illuminated over a range of incident angles via either object rotation, in which the sample itself is rotated relative to the imaging system (usually along a principal Cartesian axis), or beam rotation, in which the angle of incidence of the illuminating beam is changed relative to the object and the optical axis of the imaging system [17]. If beam rotation within a nonmoving optical system is used, only a limited range of incident angles is possible due to the finite numerical apertures (NAs) of the illuminating and imaging optics, resulting in a missing conical region of the frequency domain support in which data is only recoverable using algorithmic approaches requiring a priori knowledge of the sample [18]. Without such recovery, spatial resolution will be degraded along the optical axis. Alternatively, although object rotation enables isotropic spatial resolution, it also introduces technical challenges associated with a moving object and limits acquisition speed [17].

In order to recover RI, phase is measured at each angle of incidence. Once the phase is measured, the choice of tomographic reconstruction algorithm depends on how the interaction with the object is modeled [12]. In projection tomography, for example, the measured phase is interpreted as in 2D QPI, in which the light propagates straight through the object in an undeviated manner so that phase is simply the RI of the object integrated along the optical axis multiplied by the incident wave vector magnitude. In this case, the object can be reconstructed using conventional algorithms such as filtered back-projection [12]. This model, however, is usually inappropriate at optical wavelengths since characteristic dimensions of the object are often of the same order as the illuminating wavelength, meaning that diffraction, in addition to refraction at object boundaries, contribute to image formation and degrade the line integral approximation.

For these reasons, optical diffraction tomography (ODT) is usually employed in 3D QPI studies [17]. ODT accounts for diffraction of the incident light by the object and thus provides a more accurate model for image formation. In ODT, one of two approximations, namely, the first-order Born or first-order Rytov approximations, are usually employed to linearize the relationship between the object’s complex scattering potential and the angularly resolved complex image data, for which amplitude absorption and phase of each pixel are measured [12]. The choice of approximation once again depends on the image formation model. The first Born approximation is known to be appropriate for “weakly scattering” objects in which the total phase delay through the object is less than around $\pi /2$ [12]. For biomedical applications, the first Rytov approximation is usually a better choice as it allows for a large total phase delay as long as the gradient of the complex scattered phase is not too large [19], as is usually the case for weak RI contrast [17]. In both cases, the complex image is related to the Fourier transform of the object along a semi-circular arc in the spatial frequency domain and reconstruction algorithms using either spatial or frequency domain interpolation [20] can be utilized.

Although a significant portion of current 3D QPI research is centered on developing and applying ODT methodologies such as tomographic phase microscopy [21], ODT has some negative features which encourage the development of alternatives in parallel. In general, ODT requires the illumination to be coherent, both temporally and spatially, resulting in difficulties associated with coherent noise sources such as phase jitter and speckle interference [22]. Another factor which may prohibit the wide-scale commercialization and adoption of ODT among biomedical users is the cost associated with such laser/optomechanical systems, as ODT must combine interferometric imaging with either object or beam rotation. Most often, ODT employs beam rotation using either single [21] or dual axis [18] galvanometer-controller mirrors to change the angle of incidence. A recent approach combines two modular units to attach on to conventional microscopes providing beam rotation and single-shot 2D QPI, respectively [23]. Object rotation has also been achieved on live cells via a hollow fiber capillary cell culture [24], patch-clamping with a micropipette [25], and via holographic tweezers [26].

To address these issues, 3D QPI solutions involving the deconvolution of 3D images have been proposed [13–16]. In such methods, a 3D image is constructed by collecting a through-focal series, after which RI is recovered via 3D deconvolution based on a linearized model. This approach is similar to 3D fluorescence deconvolution microscopy in which out-of-focus blur is removed numerically [27]. Partially coherent illumination is often employed, enabling compatibility with commercial microscopy, greatly reducing the anticipated cost of such systems. The optical sectioning capability of various methods is derived from differing mechanisms including coherence and high-NA gating [13] as well as extended optical transfer function (OTF) support using partially spatially coherent illumination [16,28]. Thus, it is possible to obtain similar spatial frequency domain support (ultimately limited by illuminating NA) to ODT under beam rotation using commercial microscopy hardware [28], as has been exploited in 3D fluorescence deconvolution microscopy [27,29]. In spite of these benefits, through-focal deconvolution methods, like tomography under beam rotation, result in degraded resolution along the optical axis, which may be limiting for samples possessing complex internal structure with rapidly varying features inconsistent with constraints imposed by iterative limited-angle tomographic recovery algorithms, such as known object support or piecewise constancy [18].

In what follows, we present a new numerical reconstruction method and approach for 3D QPI, called tomographic deconvolution phase microscopy (TDPM), which addresses the aforementioned issues by combining through-focal deconvolution with object rotation to enable isotropic resolution using commercial microscopy hardware. TDPM is analogous to methods used in 3D fluorescence microscopy [30–34] which are capable of 3D spatial resolution better than confocal microscopy. The extension to QPI was originally suggested by Cogswell et al. [32], although, to the authors’ knowledge, this concept was never realized. Although the recovery model will be based on 3D weak object transfer function (WOTF) theory [28,35], we show that TDPM recovery is possible for “nonweak” phase objects with large total phase delay. Altogether, TDPM is an attractive alternative to ODT for both biomedical and industrial applications due to its compatibility with commercial microscopy, experimental simplicity, isotropic spatial resolution, and tolerance of large phase objects.

2. PRINCIPLES OF TOMOGRAPHIC DECONVOLUTION PHASE MICROSCOPY

A. Relationship to First-Order Diffraction Tomography

First-order diffraction tomography is a scalar theory based on the inhomogeneous Helmholtz equation

In Eq. (1), $k(\mathit{r})=k_{0}n(\mathit{r})$, where $k_0=2\pi /\lambda$ is the free-space wave vector magnitude for the wavelength $\lambda$ and $n(\mathit{r})=n_0+\mathrm{\Delta }n(\mathit{r})$ in which $n_0$ is the average RI and $\mathrm{\Delta }n(\mathit{r})$ is the spatially varying component which defines the object, $u(\mathit{r})$ is the total complex field amplitude (single polarization component for electromagnetic fields), and $\mathit{r}=x\widehat{x}+y\widehat{y}+z\widehat{z}$ denotes spatial coordinates. We may rewrite Eq. (1) as Eq. (2) to isolate the driving terms

In Eq. (2), $V(\mathit{r})=k_0^2[n^2(\mathit{r})-n_0^2]$ is the complex scattering potential which is evidently zero outside the support of the object (given as $V^{\prime }$). Using the method of Green’s functions, we may write the solution for $u(\mathit{r})$ as [12]

It has been shown that we may also write the solution to Eq. (2) as

Let us now consider bright-field microscopy operating in transmission for which only forward propagating waves which fall within the system aperture exist in image space. Initially, the illumination is modeled as a spatially coherent quasi-monochromatic plane wave defined by $u_0(\mathit{r},{\mathit{\rho }}^{\prime })=\sqrt{S({\mathit{\rho }}^{\prime })}\exp (i2\pi {\mathit{\rho }}^{\prime }·\mathit{r})$, where $S({\mathit{\rho }}^{\prime })$ defines intensity over the illumination pupil. The quasi-monochromatic approximation, which implies that the illumination bandwidth is much smaller than the central wavelength, or $\mathrm{\Delta }\lambda \ll \lambda$, is easily obtained in microscopy through the use of interference filters and allows us to ignore partial temporal coherence from spectrally broadened sources [36]. Let us assume that the first Rytov approximation is valid, so that the scattered complex phase is well approximated by Eq. (6). Using Eq. (5a), we may write an expression for total intensity, in which the ${\mathit{\rho }}^{\prime }$ dependence has been made explicit,

Expanding the exponent of Eq. (8) in a Taylor series reveals that if

The Fourier transform of Eq. (11) is given by

In Eq. (13), $\mathit{\rho }={\rho }_x{\widehat{\rho }}_x+{\rho }_y{\widehat{\rho }}_y+{\rho }_z{\widehat{\rho }}_z$ denotes frequency coordinates conjugate to $\mathit{r}$ and $A(\mathit{\rho })$ as well as $P(\mathit{\rho })$ are the Fourier transforms of $A(\mathit{r})$ and $P(\mathit{r})$, respectively. Also used in Eq. (13), $G^{\prime }(\mathit{\rho })$ is the Fourier transform of the Green’s function filtered to transmit forward propagating waves within the system pupil as defined by $P({\mathit{\rho }}_{\perp })$ in

Assuming an extended incoherent source, we may incorporate partial spatial coherence by integrating over the illumination pupil [28,35], so that the final intensity spectrum may be written as

 

Fig. 1. Magnitudes of the on-axis coherent (a) AOTF and (b) POTF derived from simulating the scattered complex field amplitude due to a point scatterer. All figures are plotted as a function of normalized frequency coordinates $\mathit{\rho }\lambda /n_0$ and have rotational symmetry about ${\widehat{\rho }}_z$.

Download Full Size | PDF

From Fig. 1(a) we see that the AOTF has a strong low-pass characteristic since the Ewald sphere cap and its conjugate overlap additively due to the fact that the imaginary part of $g^{\prime }(\mathit{r},{\mathit{\rho }}^{\prime })$ is an even function of $\mathit{r}$. This low-pass characteristic implies that the total absorption through the object must be small, as is consistent with our interpretation regarding scattered intensity. For the POTF, the opposite is true since the real part of $g^{\prime }(\mathit{r},{\mathit{\rho }}^{\prime })$ is odd. Thus, in Fig. 1(b) we observe the cancellation of contrast near the origin of frequency space. This implies that large but “slowly varying” phase objects with most of their energy residing in lower spatial frequencies are well approximated by the 3D WOTF theory, which further cements the Rytov approximation used in its derivation. In Section 3, examples demonstrating the validity of this observation are provided using a split-step beam propagation method (BPM) validated against rigorous electromagnetic solutions to scattering from a homogenous cylinder.

B. TDPM RI Recovery

Assuming that the validity conditions [Eq. (19) and first Rytov approximation] are met, the 3D WOTF [Eq. (17)] becomes the basis for TDPM. Shown in Fig. 2 are the AOTF and POTF (shown in the ${\rho }_x-{\rho }_z$ plane with rotational symmetry implied) that are used through the remainder of this paper which are calibrated to match the imaging properties of the microscope utilized (Olympus BX60). The OTFs in Figs. 2(a) and 2(b) were calculated by simulating the scattering due to a line absorber [$V(\mathit{r})=i\delta (x)\delta (z)$] and a line scatterer [$V(\mathit{r})=\delta (x)\delta (z)$], respectively, at $\lambda =546\mathrm{nm}$ with $n_0=1.46$.

 

Fig. 2. Partially coherent (a) AOTF and (b) POTF plotted as a function of normalized frequency coordinates $\mathit{\rho }\lambda /n_0$.

Download Full Size | PDF

In calculating the OTFs, the scattered intensity due to each coherent plane wave in the illumination pupil was summed incoherently as in Abbe’s method for partially coherent image formation [38]. In order to increase accuracy, primary spherical aberration owing to focusing through uncompensated media was modeled by adding an extra term in the exponent of the defocusing pupil [36], for which the microscope was assumed to be aberration-free at $z=0$. A circular illumination pupil with $N{A}_c=0.375$ was used to balance the trade-off between optical sectioning and image contrast. A Gaussian distribution for $S({\mathit{\rho }}^{\prime })$ was assumed based on 2D curve fitting [Fig. 3(b)] to an image [Fig. 3(a)] of the back-focal-plane (BFP) of the objective lens (Olympus UPlanFL $40\times /0.75$ $0.17\infty$) obtained by inserting a Bertrand lens into the optical train with no sample in place. Thus, each point source was weighted by the curve-fitted estimate [Fig. 3(b)].

 

Fig. 3. (a) Image of objective BFP with no sample in place. (b) Gaussian fit to (a) serving as input source distribution in the calculation of OTFs in Fig. 2. Both figures are plotted as a function of normalized frequency coordinates in the illumination pupil ${\mathit{\rho }}^{\prime }\lambda /n_0$.

Download Full Size | PDF

Since the linearization conditions depend only on the relative intensity contrast due to each point source, they are independent of $S({\mathit{\rho }}^{\prime })$ and depend only on the choice of system pupil $P({\mathit{\rho }}_{\perp })$. In order to remain compatible with commercial microscopy, the choices for $P({\mathit{\rho }}_{\perp })$ are limited. Although it may seem advantageous to utilize Zernike phase contrast, in this case the real and imaginary parts of $V(\mathit{r})$ effectively interchange roles in Eq. (19), thereby implying a weak phase condition. Differential interference contrast may be a much better option since contrast is related to lateral phase gradient rather than absolute value [39]. For bright-field microscopy, it may prove beneficial to optimize over the illumination pupil intensity distribution $S({\mathit{\rho }}^{\prime })$. Although no such optimization has been attempted, it has been observed that pupils with monotonically decreasing intensity with illumination angle provide more uniform contrast in the frequency domain and lead to increased stability.

It can be seen from Fig. 2 that spatial frequency coverage under partial spatial coherence is similar to the coverage obtained in ODT under beam rotation [28] as well as widefield deconvolution microscopy [40]. This is because each plane wave in the illumination pupil samples the same Ewald sphere cap as in ODT plus its complex conjugate. Thus, by measuring a through-focal series in a bright-field microscope one obtains similar information as in ODT with phase measurements over many angles. This is only true, however, in certain cases, such as imaging pure phase objects, in which a direct or regularized deconvolution between measured intensity and $P(\mathit{r})$ can be achieved based on $H_P(\mathit{\rho })$, as was first demonstrated by Noda et al. using annular illumination [16]. Another example would be when absorption is assumed to be proportional to phase, or $A(\mathit{r})=\epsilon P(\mathit{r})$, as in [41]. For a general object, with both weak absorption and phase, Streibl suggested that it should be possible to recover both components by measuring through-focal series under two different pupil functions [28].

In TDPM, this is realized via object rotation. If the sample is rotated by 180°, for example, the symmetries of $H_A(\mathit{\rho })$ and $H_P(\mathit{\rho })$ allow $P(\mathit{r})$ to be recovered uniquely by subtraction of their respective 3D image stacks relative to a single reference coordinate system because the absorption contrast is an even function about each scatterer. This is analogous to phase recovery using 2D WOTF theory based on subtraction of images on either side of the focus [42]. Another benefit of subtracting the through-focal series obtained from opposing perspectives is the ability to recover stronger pure phase objects because the second-order term, as well as all even-ordered terms, in the Taylor series expansion of Eq. (8) produce even contrast. In spite of the benefits of 360° coverage, for pure phase objects with weak RI contrast, complete object recovery is possible via object rotation over 180° and will be the basis of conventional TDPM RI recovery.

Assuming the addition of an experimental configuration for object rotation, which in practice will likely be a glass fiber or capillary coupled to a rotation stage/device [30,31,34], we are ready to devise a strategy for implementing TDPM. In order to sample the object spectrum in an isotropic fashion near the origin, the object must be rotated at least $N$ times, where $N\ge \lceil \pi /(2{\theta }_c)\rceil$, ${\theta }_c={\sin }^{-1}(N{A}_c/n_0)$ is the marginal illumination angle in radians, and $N{A}_c$ is the NA of the illuminating condenser lens adjustable via an aperture diaphragm in the front focal plane. In practice, however, it is often necessary to select $N\ge \lceil \pi /{\theta }_c\rceil$ in order to enable reasonable contrast across the entire spectrum. Assuming equiangular rotation, choosing $N$ results in the rotational increment $\mathrm{\Delta }\theta =\pi /N$. An optimal solution for $P(\mathit{\rho })$ is sought via the least squares formalism summarized by Eq. (20a), where $j$ is an index associated with object rotation angle ${\theta }_j=j\mathrm{\Delta }\theta$, $I_{{\theta }_j}(\mathit{\rho })/B$ are the zero-mean normalized 3D intensity spectra, and $\alpha$ is a regularization parameter. Equation (20) can be solved directly to yield $P(\mathit{\rho })$, as given by Eq. (20b),

The numerical implementation of Eq. (20) is nontrivial. Also, as can be seen in the upper-right inset of Fig. 2(b), contrast for low spatial frequencies near the origin is reduced resulting in ill-posed recovery, which is a well-known problem in 2D phase retrieval [43]. In fact, if $\mathrm{\Omega }$ is selected to be a rectangular prism which barely encompasses the object then the frequency domain is likely to be undersampled resulting in spatial aliasing artifacts. In this case, there are likely to be low spatial frequencies which are never sampled and are thus unrecoverable via 3D Fourier inversion. TDPM solves these issues by using a hybrid algorithm with different processing steps for low versus high spatial frequencies. The high-frequency algorithm is based on Eq. (20) and is summarized in Fig. 4. The low-frequency algorithm is the partially coherent analog to filtered back-propagation [44] and is summarized in Fig. 5. Although it has proven difficult to optimize the low-frequency algorithm as in Eq. (20), preliminary results indicate that such an optimization would only provide marginal improvement.

 

Fig. 4. Block diagram representation of TDPM RI recovery for high spatial frequencies.

Download Full Size | PDF

 

Fig. 5. Block diagram representation of TDPM RI recovery for low spatial frequencies.

Download Full Size | PDF

The high-frequency algorithm shown in Fig. 4 is divided into four stages. In the image capture and pre-processing stage, the through-focal series of an object are first acquired over the measurement domain $\mathrm{\Omega }$ as the object is rotated in equal increments of $\mathrm{\Delta }\theta$ about the y axis and then processed as inputs for RI recovery. The domain $\mathrm{\Omega }$ should at least encompass the object, which is usually either known or easily estimated. If possible, $\mathrm{\Omega }$ should also encompass the scattered intensity variations, thereby preventing errors associated with spatial aliasing and providing sufficient frequency domain resolution for a complete reconstruction via the optimized algorithm shown in Fig. 4. This will be difficult to achieve in practice, particularly for larger objects under increased spatial coherence. Thus far, however, simulation and experimentation suggest that errors due to spatial aliasing do not inhibit quantitative interpretation of the data.

Once the through-focal images are collected over all rotation angles, pre-processing steps include background intensity normalization and subtraction, $z$-slice registration, upsampling, and 3D image registration. Since background intensity $B$ is a conserved quantity [28], each $z$-slice image is first normalized by its average. Following normalization, background intensity variations which are characteristic of the system and not of the sample are removed via subtraction with a background through-focal series measured over the same domain $\mathrm{\Omega }$. The next pre-processing step entails the registration of each $z$-slice to its nearest neighbors. This is only necessary in the event of object movement during 3D image measurement and has been successfully implemented via normalized cross-correlation between neighboring slices. Once the through-focal series is aligned internally, it is then usually upsampled along the z axis. This is because the microscope depth of field is usually larger than the lateral resolution $\mathrm{\Delta }r=\mathrm{\Delta }x=\mathrm{\Delta }y$. An efficient sampling strategy, therefore, is to sample at the largest integer multiple of $\mathrm{\Delta }r$ which is less than or equal to the depth of field, so that $\mathrm{\Delta }z=\mathrm{\Delta }r$ after upsampling by the same integer. Note that care should be taken to ensure that defocus distances used in the measurement correspond to distances within object space, which may have a background RI $n_0$ greater than one, implying that the microscope should be defocused by $\mathrm{\Delta }z/n_0$ to realize a distance of $\mathrm{\Delta }z$ within object space.

The final pre-processing step is 3D image registration which is usually necessary due to radial runout of the rotating cylinder. Although phase correlation has previously been used for this step [30], herein another normalized cross-correlation between the measured through-focal series and the pre-simulated subimages which are characteristic of scattering from the cylindrical reference boundaries on a column by column basis along the axis of rotation (y axis) was found to be sufficient. Once positions of maximum correlation are identified, each column is circularly shifted to the center ensuring a single rotational reference.

The next stage consists of filtering the spectrum of each through-focal series with the POTF conjugate. Upon inverse Fourier transformation, the images are ready for inverse rotation via bilinear interpolation to compensate for their physical rotation angle. In the final stage, the high spatial frequencies of RI are synthesized by summing over the rotation angle, compensating for the frequency domain overlap between measurement angles via the denominator of Eq. (20b), filtering with a high-pass filter $H_{hp}(\mathit{\rho })$, offsetting the result so that the scattering potential $P(\mathit{r})$ is zero outside the object, and converting from scattering potential to RI. The denominator of Eq. (20b) is constructed by Fourier transforming (indicated by $\mathfrak{I}$) the sum of rotated phase PSF autocorrelations [Eq. (21)]. The rotations are made in real space via bilinear interpolation, which is found to yield less reconstruction error than frequency domain interpolation,

Thus far, a hard cutoff high-pass filter has been used [$H_{hp}(\rho )=U(\rho -\xi )$, where $\rho =|\mathit{\rho }|$] which retains all frequencies above a radial threshold $\xi =1/({\mathrm{\Omega }}_z\sin {\theta }_c)$(${\mathrm{\Omega }}_z$ is the extent of $\mathrm{\Omega }$ along to optical axis) which defines the boundary between sampled frequencies and frequencies which are never sampled due to insufficient frequency resolution.

Because frequencies below $\xi$ are unrecoverable when $\mathrm{\Omega }$ does not encompass the scattered intensity, the algorithm shown in Fig. 5 is necessary. The first stage of Fig. 5 is the same as the first stage of Fig. 4. The first stage shown in Fig. 5, which is the second stage in the low-frequency algorithm, consists of applying 2D phase retrieval to solve for phase at each $z$-slice through the object. Since it needs to operate on partially coherent intensity data, an algorithm is selected which can easily model partial coherence, such as TIE phase recovery [45] or methods based on inversion of the 2D WOTF [46]. In this work, a recent phase reconstruction method referred to as POTF recovery [42] was utilized. This method is the 2D analog of TDPM and results in phase recovery from multiple defocused planes which is optimal in the sense of minimized noise in the final phase image.

Actually, the general use of depth-resolved phase recovery is unjustified in the case of partially coherent illumination of a 3D phase object because phase has no meaning in the out-of-focus planes from each slice. Below the previously defined threshold $\rho <\xi$, however, the projection approximation is actually rather good, justifying the use of depth-resolved phase recovery for these frequencies since the phase contributed from each slice is roughly independent of the plane in which it is measured. Small phase variations do occur over the same length scale as ${\mathrm{\Omega }}_z$, which is why the algorithm detailed in Fig. 5 is used, as opposed to conventional filtered back-projection [12].

In implementing depth-resolved phase retrieval, a selected number of defocused intensity images on either side of the focus are used as inputs to the POTF recovery algorithm. This number should be chosen so that phase can be reconstructed over the entire extent of the object without needing to use defocused images estimated via circular padding. Another good reason to select $\mathrm{\Omega }$ to be as large as possible is that the eventual phase SNR is roughly proportional to the defocus range used. Planes near the top and bottom edges which are unrecoverable in this manner are estimated by extension of their nearest recoverable neighbor. The remaining stages shown in Fig. 5 correspond to the conventional filtered back-propagation algorithm [44]. In order to compensate for the increased sampling density near the spatial frequency origin, depth resolved phases are filtered using normalized Ram–Lak filters with cylindrical symmetry ($\mathrm{\Delta }r{\rho }_{xz}/N$ in which ${\rho }_{xz}=\sqrt{{\rho }_x^2+{\rho }_z^2}$) after which rotation is achieved via bilinear interpolation in the spatial domain [44]. In the synthesis stage, the filtered phases are summed over rotation angle, low-pass filtered using $H_{lp}(\rho )=1-U(\rho -\xi )$, and level shifted to compensate for the lack of absolute phase information. The final step in TDPM RI recovery is to add the results from Figs. 4 and 5 to obtain the overall RI.

3. SIMULATION RESULTS

A. Modified Split-Step BPM

In order to model the imaging of 3D phase objects using bright-field microscopy, a modified split-step wide-angle beam propagation method (BPM) is used. In this method, contributions arising from each point source in the illumination pupil are added incoherently to form the final 3D bright-field image. Each coherent simulation consists of implementing an off-axis wide-angle BPM [47] which incorporates an obliquity factor (OF) given by Eq. (21) and which is associated with the local phase gradient magnitude [48], which compensates for additional phase delay associated with propagating through an effective thickness of $\mathrm{OF}(\mathit{r})\mathrm{\Delta }z$, where $\mathrm{\Delta }z$ is the longitudinal resolution of the simulation. In Eq. (21), ${\nabla }_{\perp }$ is a gradient operating on lateral coordinates only and $\varphi (\mathit{r})$ is the phase of the total field $u(\mathit{r})$:

The BPM is initialized with $u_z(\mathit{r})=u_0(\mathit{r},{\mathit{\rho }}^{\prime })$, as given in Section 2.A. At each $z$-slice, the wave is first propagated by a half-step using the angular spectrum method in accordance with Eq. (22) [49], after which additional phase delay is added according to Eq. (23), and then Eq. (22) is applied once again:

For simplicity, the OFs are capped at $\sqrt{2}$ to avoid modeling the propagation between neighboring pixels in a given $z$-slice. At each $z$-slice, the phase gradient is estimated via a central difference approximation acting on the unwrapped phase values of the previous slice. Once the aforementioned algorithm has been applied through the entire structure, the result is then filtered by $N{A}_o$ and back-propagated through the simulation space. The squared magnitude is the intensity contribution associated with the source point ${\mathit{\rho }}^{\prime }$. Integrating over $S({\mathit{\rho }}^{\prime })$ completes the partially coherent 3D image.

In order to validate this model, its coherent outputs are compared with rigorous vectorial solutions to Maxwell’s equations for plane waves polarized along the axis of cylinders with real homogeneous RI, as found in [50]. Simulations were conducted to compare plane wave scattering at both normal and marginal incidence for a range of RI values. In order to enable a direct comparison with forward-scattered waves detectable in transmission, the complex electric field amplitude from rigorous solutions which incorporate multiple scattering due to both forward and backward propagating waves was calculated at positions coincident with the last $z$-slice. Since the problem is now constrained to 2D, this 1D wave is then filtered by $N{A}_o$ and back-propagated through the same simulation space, resulting in a solution for forward propagating waves only.

The BPM simulation is then compared against this solution using the normalized mean square error of scattered intensity [NMSE, given by Eq. (24)] as a metric with results shown in Fig. 6. In Eq. (24), $I_S(\mathit{r})$ is the simulated forward-scattered intensity and $I_A(\mathit{r})$ is the analytic solution with mean value ${\overline{I}}_A$:

 

Fig. 6. (a) NMSEs of various BPM simulations with and without OF correction [Eq. (21)] for normal (${\mathit{\rho }}^{\prime }=\mathbf{0}$) and marginal (${\mathit{\rho }}^{\prime }={\rho }_c{\widehat{\rho }}_x=N{A}_c{\widehat{\rho }}_x/\lambda$) incidence. (b) Simulated (with OF correction) and (c) analytic intensities with $\mathrm{\Delta }n=0.025$. Simulation parameters: $\lambda =546\mathrm{nm}$, $n_0=1$, $\mathrm{\Delta }x=\mathrm{\Delta }z=245\mathrm{nm}$, $N{A}_o=0.75$, and $N{A}_c=0.375$.

Download Full Size | PDF

From Fig. 6(a) it is observed that correcting for the OF through each slice reduces error significantly, which is useful since split-step BPMs assume $\mathrm{\Delta }n/n_0\ll 1$ [47], yet the RI contrast for dehydrated cells in water is about $(1.55-1.33)/1.33\approx 0.17$, where $n=1.55$ for dehydrated cells [51] has been assumed. Also from Fig. 6, it is observed that NMSEs are less than 20% for $\mathrm{\Delta }n/n_0<0.15$ with OF correction, validating the model and indicating its usefulness for simulating bright-field imagery from 3D phase objects. For the purposes of modeling TDPM reconstruction, this model is especially well suited as the intensity contrast must satisfy Eq. (9). Therefore, $\mathrm{\Delta }n$ must be weak as in ODT under the first Rytov approximation [12], resulting in improved model accuracy. This can be seen qualitatively in Figs. 6(b) and 6(c) for which Eq. (9) begins to break down.

B. Simulated TDPM Reconstruction

Due to memory and time constraints imposed by modeling 3D intensity distributions under partial coherence, the tomographic reconstructions presented in this section are, without loss of generality in 3D, based on a 2D cylindrical phantom (shown in Fig. 7). Even though the object is 2D, off-axis waves emanating from the entire illumination aperture [Fig. 3(b)] have been incorporated in the partially coherent image calculation, so that the modeled intensities coincide with cylindrical scattering under Köhler illumination. The parameters used in the simulation are the same as in Fig. 6 except that $n_0=1.46$ as opposed to $n_0=1$.

 

Fig. 7. RI contrast $\mathrm{\Delta }n(\mathit{r})$ of modified Shepp–Logan phantom.

Download Full Size | PDF

Figure 7 shows a modified version of the Shepp–Logan phantom in which the outer skull material has been replaced by the surrounding head material. For such a phantom the surrounding head material may represent cytoplasm with internal ellipses representing organelles.

In Fig. 8, multiple tomographic reconstructions, including reconstructions obtained via ODT using filtered back-propagation under the first Born [Figs. 8(a), 8(d), and 8(g)] and Rytov [Figs. 8(b), 8(e), and 8(h)] approximations as well as the TDPM method [Figs. 8(c), 8(f), and 8(i)], are compared directly. The simulations are also differentiated by row according to maximum RI contrast $\mathrm{\Delta }{n}_{\max }$. In the results obtained via filtered back-propagation, the object was rotated 825 times about the y axis so that the rotational increment corresponded roughly with the angular resolution (along the y axis) used in the partially coherent image simulation.

 

Fig. 8. Reconstructions obtained using filtered back-propagation under the [(a), (d), (g)] Born and [(b), (e), (h)] Rytov approximations as well as [(c), (f), (i)] TDPM for maximum RI contrast values of [(a), (b), (c)] $\mathrm{\Delta }{n}_{\max }=0.004$, [(d), (e), (f)] $\mathrm{\Delta }{n}_{\max }=0.02$, and [(g), (h), (i)] $\mathrm{\Delta }{n}_{\max }=0.1$. The resulting RMSEs are (a) 0.00042, (b) 0.00015, (c) 0.00017, (d) 0.00580, (e) 0.00088, (f) 0.00106, (g) 0.02681, (h) 0.01124, and (i) 0.01293.

Download Full Size | PDF

For each rotation angle, the aforementioned modified spilled-step BPM was used to calculate the total field, and the scattered field was used in the Born approximation and the scattered phase was used in the Rytov approximation. In calculating the scattered phase, MATLAB’s unwrap function was utilized. The TDPM reconstructions utilized the methods outlined in Figs. 4 and 5 using the POTF shown in Fig. 2. The TDPM reconstruction parameters were set to $N=15$ and $\alpha =0$.

Immediately it is seen that TDPM is not as restrictive as ODT in the first Born approximation in the size of objects which may be reconstructed, verifying that the object need not be “weak” as in the conventional interpretation of the first Born approximation [12], but must be “slowly varying” as given by Eq. (19). This is evident by comparing the reconstructions made under the first Born approximation [Figs. 8(a), 8(d), and 8(g)] with the TDPM reconstructions shown in Figs. 8(c), 8(f), and 8(i). Even when RI contrast is relatively weak, as in Fig. 8(a) in which the total phase delay through the object is approximately 1.1 rad propagating left to right, the first Born approximation results in rotations in the complex plane which degrade the result [12]. Likewise, similar behavior is observed and is more pronounced in Figs. 8(d) and 8(g), in which the total phase delays are 5.6 and 28.2 rad, respectively. By contrast, TDPM reconstructions shown in Figs. 8(c) and 8(f) display no such rotation effect, and the associated reconstruction RMSEs are comparable to reconstructions in the first Rytov approximation [Figs. 8(b) and 8(e)] without any visible degradation of resolution or image quality.

For the case of $\mathrm{\Delta }{n}_{\max }=0.1$ the scattered phase is highly wrapped, leading to phase unwrapping errors and associated artifacts in the Rytov reconstruction shown in Fig. 8(h). In addition to errors due to phase unwrapping, deformation errors exist in Fig. 8(h), such as enlargement of the phantom’s left eye ellipse, which are associated with the breakdown of the Rytov approximation. In TDPM, no phase unwrapping is necessary and thus the reconstruction shown in Fig. 8(i) does not contain the same artifacts as in Fig. 8(h). The artifacts associated with TDPM appear to be predominant near object boundaries, resulting in a degradation of spatial resolution when RI contrast is too high. This results from an asymmetrical intensity distribution induced via multiple scattering events through the object which is especially pronounced near edges. It may be possible to counterbalance this asymmetry by illuminating from opposing angles over 360° as alluded to in Section 2.B. The application of this characteristic, however, will be the subject of future work. In spite of these errors, it is clear that TDPM, like ODT under the Rytov approximation, will be useful over an appreciable range of RI contrast and will thus be relevant and applicable to biomedical studies.

In order to quantify the reconstruction error associated with noise in the bright-field imagery, white Gaussian noise with a normalized standard deviation set to 1% of the background intensity level ($\sigma =0.01$) was added to the simulated intensity values. This value was selected to match the noise produced by the camera (QImaging Retiga 1300R) used in the experimental results shown in Section 4 and is typical of many scientific imagers used for microscopy. The TDPM reconstruction results obtained with noise added are shown in Fig. 9(a) for $\mathrm{\Delta }{n}_{\max }=0.004$ and Fig. 9(b) for $\mathrm{\Delta }{n}_{\max }=0.04$.

 

Fig. 9. TDPM reconstructions obtained with additive noise with a normalized standard deviation of $\sigma =0.01$ and (a) $\mathrm{\Delta }{n}_{\max }=0.004$ and (b) $\mathrm{\Delta }{n}_{\max }=0.04$. The resulting RMSEs are (a) 0.00026 and (b) 0.00294.

Download Full Size | PDF

The reconstruction noise is visible primarily in Fig. 9(a) and dominated by the signal in the case of Fig. 9(b). In addition to demonstrating potential SNR, Fig. 9(b), for which the total phase delay through the object is approximately 11.3 rad, further supports the fact that TDPM does not require small total phase delay and can result in good reconstruction quality for objects with wrapped phase. Theoretically, the RI error associated with image noise will, at least for the high-frequency algorithm shown in Fig. 4, be spatial frequency dependent and inversely proportional to the square root of the denominator in Eq. (20b). Here the RI error associated with image noise is estimated by subtracting the reconstruction obtained without noise [Fig. 8(c)] from the reconstruction including noise [Fig. 9(a)] to yield a single-valued RI sensitivity of $\sim 2\times {10}^{-4}$ RI units.

4. EXPERIMENTAL RESULTS

In order to demonstrate TDPM experimentally, optical fibers are used as test 2D phase objects and a single exposure from a ${\mathrm{CO}}_2$-laser-induced azimuthally asymmetric long-period fiber grating (LPFG) is used as a test 3D phase object. In order to implement TDPM, all that is required, in addition to a commercial microscope with automated defocus control, is an external stage for object rotation. The configuration utilized in the present work is designed to implement tomography on optical fibers using an upright microscope (Olympus BX60) and has been described elsewhere [52,53]. In addition to this configuration, TDPM can immediately benefit from the groundwork laid for object rotation in similar fluorescent [30,31,34] and phase [24–26] methods. For the measurements presented here, the illuminating source was a mercury-arc lamp with a $\lambda =546\mathrm{nm}$ green interference filter with a full width at half-maximum bandwidth of $\mathrm{\Delta }\lambda =10\mathrm{nm}$. The imaging parameters and components are as outlined in Section 2.B in order to match the calculated POTF. The microscope defocusing was automated using a piezoelectric microscope objective scanner (Physik Instrumente P-721.SL2 with E-709.SR controller).

Figure 10 shows the results of implementing 2D TDPM, using the same reconstruction parameters outlined in Section 3 ($N=15$ and $\alpha =0$), with a single-mode fiber [Corning SMF-28, Fig. 10(a)], polarization-maintaining fiber [Thorlabs HB980T, Fig. 10(b)], and photonic-crystal fiber [Blaze Photonics ESM-12-01, Fig. 10(c)] used as test phase objects. For all cases, the defocused imagery were sampled at $\mathrm{\Delta }z=4\mathrm{\Delta }x=980\mathrm{nm}$, with 147 images per rotation angle for a total of 2,205 images. Since the objects are cylindrical, averaging along the fiber axis is possible to improve SNR and for the tomograms shown in Fig. 10, 51 columns were averaged.

 

Fig. 10. 2D TDPM cross-sectional tomograms obtained on (a) single-mode, (b) polarization-maintaining, and (c) photonic-crystal fibers. All colorbars indicate RI units.

Download Full Size | PDF

In all cases, the fibers were surrounded by RI matching oil ($n_0=1.46$ at $\lambda =589\mathrm{nm}$, Cargille Labs Series A) to match the RI of the fused silica cladding ($n=1.4601$ at $\lambda =546\mathrm{nm}$). For normal glass dispersion results in a slightly higher oil index $n_0>1.46$, however, for the results presented here, RI values are offset to $n_0=1.46$, which was the background index used in calculating the POTF. In Fig. 10(a), the expected step profile between fiber core and cladding is observed. Also seen is the well-known “center dip” in RI in the fiber’s core associated with dopant burnoff effect. This effect is also clearly seen in Fig. 10(b), in which the RI of the polarization-maintaining (bow-tie) fiber’s stress members is also clearly visible and well resolved. Last, the reconstructed hexagonal lattice structure of the photonic-crystal fiber shown in Fig. 10(c) highlights TDPM’s capability, and the results may be compared directly with a recently published state-of-the-art optical fiber tomographic algorithm which is based on ODT in the Rytov approximation [54]. In addition to the lattice structure, RI features resulting from the modification residual stresses in the fiber, such as the ring surrounding the air-hole lattice, are visible in Fig. 10(c).

In order to demonstrate the 3D capability of TDPM, we have implemented the full 3D reconstruction procedure over a field of view of $293\times 651\times 293$ cubic voxels, each with a volume of $\mathrm{\Delta }{r}^3={0.49}^3{\mathrm{\mu m}}^3$, corresponding to physical dimensions of $\sim 143\mathrm{\mu m}\times 318.5\mathrm{\mu m}\times 143\mathrm{\mu m}$. The resolution $\mathrm{\Delta }r=490\mathrm{nm}$ is a factor of two larger than previous cases to prevent excessive memory usage for 3D arrays. As in the 2D results, $N=15$ angles were used with 147 images ($\mathrm{\Delta }z=2\mathrm{\Delta }r=980\mathrm{nm}$) taken at each angle. The reconstruction procedure utilized the full 3D POTF as opposed to a column-by-column implementation of the 2D procedure used before.

The 3D sample consisted of a ${\mathrm{CO}}_2$-laser-induced azimuthally asymmetric LPFG period ($\mathrm{\Lambda }=335\mathrm{\mu m}$) [55,56]. In spite of the success of ${\mathrm{CO}}_2$-laser-induced LPFGs [57] since their introduction by Davis et al. [55], the mechanisms for grating formation in these [57], as well as arc-induced LPFGs [58], have not yet been fully characterized. Different investigations suggest differing mechanisms, such as residual stress modification [59], cladding densification [60], or geometric deformation [61]. A recent report by Hutsel and Gaylord presents 3D QPI data on a single-mode fiber exposed to focused ${\mathrm{CO}}_2$-laser radiation from exposures of successive durations [56]. The LPFG period measured here was fabricated using the same experimental configuration with two pulses of 200 and 100 ms which may be directly compared with the results in [56] for one 300 ms pulse. The results of TDPM applied to this sample are summarized in Figs. 11 and 12.

 

Fig. 11. RI modification induced via ${\mathrm{CO}}_2$-laser exposure. (a) Unperturbed single-mode fiber reconstructed $\sim 150\mathrm{\mu m}$ away from the center of exposure. (b) Reconstruction near the center of the exposed region showing azimuthal variation. (c) Multiple slices showing the 3D nature of TDPM data plotted with a reduced colormap range to highlight both azimuthal and axial changes in the fiber cladding facing the exposure. All colorbars indicate RI units.

Download Full Size | PDF

 

Fig. 12. Line profiles showing the axial variation of RI in selected regions of ${\mathrm{CO}}_2$ exposed single-mode fiber.

Download Full Size | PDF

In Fig. 11(a), a cross-sectional RI tomogram reconstructed $\sim 150\mathrm{\mu m}$ from the laser exposure is displayed and is similar to the result shown in Fig. 10(a) except that the RI difference between the matching oil and the fiber cladding is smaller. This most likely results from the temperature dependence of the oil ($\mathrm{\Delta }n=-0.00038/°\mathrm{C}$) since the observed difference corresponds to a temperature difference of $\sim 2°\mathrm{C}$ and the measurements were performed on different days with different pre-stabilization periods (the microscope is normally turned on for around three hours prior to imaging to stabilize the temperature of the oil). In Fig. 11(b), another tomogram, reconstructed near the center of ${\mathrm{CO}}_2$-laser exposure, is shown and clearly demonstrates the expected azimuthal asymmetry (see the lower-left inset in which the exposure is clearly incident from the upper-right).

The axial extent of this cladding index change can be easily visualized using Figs. 11(c) and  12. In Fig. 11(c), multiple tomograms throughout the volume are represented simultaneously and clearly show the asymmetric cladding perturbation and its axial variation. Selected line profiles, corresponding to the index matching oil, cladding on the opposite side of exposure [lower-left quadrant of Fig. 11(b)], cladding facing the exposure [upper-right quadrant of Fig. 11(b)], and core are overlaid on Fig. 11(c) and shown in Fig. 12. As expected, the RI of the oil remains constant throughout the extent of the sample. Similarly, the cladding side opposite the exposure appears to remain relatively unaffected. In contrast, the cladding side facing the exposure is clearly modulated and the magnitude and extent of RI change shown here ($\sim 5\times {10}^{-4}$ RI units and $\sim 100\mathrm{\mu m}$) is consistent with the results in [56], which adds validity to the proposed method since the results in [56] were based on a quantitative phase tomography technique [52], which has been used in a variety of fiber investigations [4,62]. The source of this cladding index modulation has been predicted to be glass densification caused by the relaxation of viscoelasticity frozen-in to the fiber during draw [63].

Also shown in Fig. 12 is the core RI values which indicate a slight increase due to laser exposure, which is in contradiction to the lowering predicted by residual mechanical stress relaxation [59]. Such an increase may indicate some form of densification in the core as well, however, a complete characterization and study of grating formation mechanisms is a subject of future work. Overall, the experiments presented here for both 2D and 3D reconstructions of optical fibers demonstrate the ability of TDPM to achieve 3D QPI in samples possessing RI variation in all three spatial dimensions.

To assess the spatial and RI resolution of the aforementioned measurements, the cladding oil step response of the single-mode fiber shown in Fig. 10(a) as well as the standard deviation of RI in the oil region of the LPFG period shown in Fig. 11(c). The 10%–90% rise distance determined by the cladding oil step response was $\sim 735\mathrm{nm}$ and provides a practical measure for spatial resolution in all three spatial dimensions. The standard deviation of RI in the oil region yielded a value of $\sim 7.7\times {10}^{-5}$, which provides an estimate of the RI sensitivity. Improved estimates of these specifications could be obtained by using polystyrene nanospheres to measure the system point spread function and corresponding POTF, which is another subject of future work.

5. SUMMARY

In summary, a new method, called tomographic deconvolution phase microscopy or TDPM, is described which enables 3D QPI using commercial microscopy with minimal hardware modification. The linearization conditions for TDPM, comprising both weak absorption and slowly varying phase, have been elucidated and indicate applicability with large phase objects in which reconstructions based on a first Born approximation are known to fail. The spatial frequency domain support of TDPM is roughly isotropic and requires no a priori knowledge of the sample or phase unwrapping as in limited-angle ODT in the Rytov approximation.

Due to its compatibility with commercial microscopy, TDPM is particularly well suited to wide-scale application among biomedical users. Techniques associated with 3D cell fixation and culture using glass capillaries and other cylindrical housings are still a matter of research [24–26,34]. However, it is anticipated that such techniques will be enabling in the application of TDPM to both fixed and living cells. Another important application area is in optical fiber characterization, such as the study of grating mechanisms in LPFGs, which was only briefly mentioned here. For applications which preclude the use of object rotation, such as high-speed 3D QPI for imaging live cell dynamics [64], the modification of TDPM to incorporate algorithmic recovery of the “missing cone” of spatial frequencies from a single through-focal series should be possible and will be the subject of future work. In such an approach no hardware modification would be necessary, permitting 3D QPI with similar frequency domain coverage as limited-angle ODT with only a fraction of the complexity and associated costs.