3D differential phase contrast microscopy

Michael Chen; Lei Tian; Laura Waller

doi:10.1364/BOE.7.003940

1. Introduction

Standard commercial microscopes use partially spatially coherent illumination for better spatial resolution, higher light throughput and reduction of speckle artifacts, as compared to coherent illumination. However, only absorption information is directly visible, not phase, which is a critical drawback for imaging unstained samples such as biological cells. Many qualitative and quantitative methods have been proposed to enable visualization of a transparent sample’s phase delay with partially coherent light. Zernike phase contrast (PhC) and differential interference contrast (DIC) render phase effects visible [1]; however, a single image is not quantitative. To measure phase quantitatively, PhC [2] and DIC [3–5] can be extended by phase-shifting approaches. Alternatively, other partially coherent 2D phase methods have been proposed, using defocus, illumination or interference [2,6–10].

When the sample is thicker than the microscope’s depth-of-field, we require 3D imaging methods. 2D phase measures integrated optical path delay through the sample, so refractive index and thickness are coupled. This ambiguity between shape and index is naturally removed in 3D phase imaging, which recovers the 3D refractive index distribution. Traditionally, 3D phase imaging is achieved tomographically - by capturing 2D projections at many angles [11–13], often employing priors to mitigate limited-angle artifacts [14,15]. In some cases, a ray-based model is sufficient (e.g. in X-ray [11,16,17]). However, when diffraction effects become prominent (e.g. in the visible regime), a diffraction tomography model [12,18] is needed. Generally, this assumes knowledge of the complex-field at each angle, requiring a two-step inverse problem: 2D phase retrieval, followed by tomography to reconstruct 3D. The 2D phase projection reconstructions may contain artifacts that propagate to the 3D reconstruction. Global reconstruction methods that relate all the measurements to the final estimate, without an intermediate 2D phase retrieval step, can be more robust to experimental errors [19,20].

Though one can digitally back-propagate a measured 2D complex-field to refocus, this does not provide the optical sectioning capabilities of true 3D imaging. Phase being a projected quantity means that one cannot simply measure 2D phase at different focus planes to reconstruct 3D phase with coherent illumination. Partially coherent imaging, however, provides a distinct focus plane and optical depth sectioning. Hence, with partially coherent light, one can use 2D phase retrieval at multiple focus planes to reconstruct 3D refractive index. Previous 3D phase imaging [21] used defocus with temporal partial coherence for sectioning. Here, we use defocus with spatial partial coherence for sectioning, then solve the corresponding 3D inverse problem.

Our method is an extension of differential phase contrast (DPC) microscopy [22–25], using asymmetric illumination for phase contrast. Four images are captured with rotated half-circle source patterns, for which 2D quantitative phase recovery models have been developed [24,25]. The dynamic source switching is achieved in a commercial microscope whose source has been replaced with a programmable LED array [25–29]. This flexible hardware platform has been used for gigapixel imaging [27,29,30], multi-contrast [26,28], 3D phase [20] and phase contrast [31], and aberration removal [30,32].

Here, we present a 3D DPC model that recovers 3D absorption and refractive index from intensity images taken at different focus planes with each of the 4 half-circle source patterns (Fig. 1). In order to properly account for out-of-focus contributions, we derive a full 3D model, rather than solving for 2D phase independently at each focus plane. Our algorithm is global in that it solves for 3D index directly from the measurements, without an intermediate phase retrieval step. In order to formulate a linear inverse problem, we consider only weak scattering (first Born/Rytov approximation) [33–36]. Finally, we explore the use of priors for mitigating both out-of-focus and halo artifacts. The resulting non-interferometric 3D quantitative phase method is simple to implement in a commercial microscope, achieves the incoherent resolution limit (2× the coherent resolution limit) and is accurate for most biological samples.

Fig. 1 3D Differential Phase Contrast (DPC) microscopy. The setup is a microscope equipped with LED array illumination and an axial motion stage. Through-focus intensity stacks are captured using 4 source patterns (top, bottom, right, and left half-circles). The intensity data is related to the 3D refractive index distribution by illumination-dependent transfer functions, according to the Born approximation. A deconvolution algorithm then recovers the 3D complex refractive index.

Download Full Size | PDF

2. Principles of 3D differential phase contrast (DPC) microscopy

2.1. Forward model

A 3D sample can be characterized by its scattering potential $V (\vec{r}) = k_{o}^{2} \cdot (n_{o}^{2} - n^{2} (\vec{r}))$ [33], where n_o is the refractive index of the surrounding media and n = n_Re + i · n_Im is the sample’s complex refractive index. The real part of the refractive index, n_Re, describes how phase delays accumulate as light passes through the sample, while the imaginary part, n_Im, describes absorption.

We start by describing the output field $U (\vec{r})$ of a 3D sample under illumination by an incident field $U_{I} (\vec{r})$ . The product of the two is convolved with the 3D Green’s function G, which characterizes the 3D field distribution from a point source,

U (\vec{r}) = U_{I} (\vec{r}) + \int \int \int U ({\vec{r}}^{'}) \cdot V ({\vec{r}}^{'}) \cdot G (\vec{r} - {\vec{r}}^{'}) d^{3} {\vec{r}}^{'},

where

\vec{r}

and

{\vec{r}}^{'}

denote 3D spatial coordinates. The convolution term can be thought of as an equilibrium solution from the multiple scattering interactions between the output field U and the 3D sample V. For weakly scattering objects, multiple scattering is negligible and the convolution term can be approximated by the single-scattered field (first Born approximation [33]),

U (\vec{r}) \approx U_{I} (\vec{r}) + \int \int \int U_{I} ({\vec{r}}^{'}) \cdot V ({\vec{r}}^{'}) \cdot G (\vec{r} - {\vec{r}}^{'}) d^{3} {\vec{r}}^{'} .

Next, we consider quasi-monochromatic (with wavelength λ) spatially partially coherent illumination from an extended source in Köhler geometry. Each point on the source generates a distinct plane wave at the sample, having 3D spatial frequency $({\vec{u}}^{'}, η^{'})$ described by $U_{I} ({\vec{u}}^{'}, η^{'}) = e^{- i 2 π ({\vec{u}}^{'} \cdot \vec{x} + η^{'} z)}$ , where ${\vec{u}}^{'}$ represents 2D transverse frequencies, η′ is 1D axial frequency and ${| {\vec{u}}^{'} |}^{2} + {η^{'}}^{2} = λ^{- 2}$ . After passing through the sample, the field is low-pass filtered by the pupil function $P (\vec{u})$ of the microscope, corresponding to the 2D coherent point spread function $h (\vec{x})$ . The measured intensity $I_{img} (\vec{x}, z)$ is the incoherent sum of intensities from each point source. Under the first Born approximation and Eq. (2), the measured intensity is

I_{img} (\vec{x}, z) = \int \int S ({\vec{u}}^{'}) \cdot {| [e^{- i 2 π ({\vec{u}}^{'} \cdot \vec{x} + η^{'} z)} + (e^{- i 2 π ({\vec{u}}^{'} \cdot \vec{x} + η^{'} z)} \cdot V (\vec{x}, z)) \otimes G (\vec{x}, z)] \otimes h (\vec{x}) |}^{2} d^{2} {\vec{u}}^{'},

where ⊗ is the 3D convolution operator,

\vec{x}

denotes transverse coordinates, and z the axial coordinate, with its origin at the focal plane of the microscope. The 2D intensity distribution of the source is described by

S ({\vec{u}}^{'})

. In practice, we can only measure 2D intensity in a single image. Thus, we scan the sample through z and capture images at many focus planes in order to build up the 3D intensity measurement on the left side of Eq. (3). Taking a 3D Fourier transform

\tilde{\cdot}

of both sides of Eq. (3), we arrive at a Fourier space relation between the 3D measured intensity·spectrum and the scattering potential,

\begin{array}{r} {\tilde{I}}_{img} (\vec{u}, η) = \int \int S ({\vec{u}}^{'}) \cdot [P (\vec{u}) \cdot (δ (\vec{u} + {\vec{u}}^{'}, η + η^{'}) + \tilde{V} (\vec{u} + {\vec{u}}^{'}, η + η^{'}) \cdot \tilde{G} (\vec{u}, η))] \otimes \\ [P * (- \vec{u}) \cdot (δ (\vec{u} - {\vec{u}}^{'}, η - η^{'}) + \tilde{V} * ({\vec{u}}^{'} - \vec{u}, η^{'} - η) \cdot \tilde{G} * (- \vec{u}, - η))] d^{2} {\vec{u}}^{'}, \end{array}

where

\vec{u}

denotes transverse spatial frequency variables and η represents the axial frequency variable. We can further expand Eq. (4) into the following three terms:

\begin{matrix} {\tilde{I}}_{img} (\vec{u}, η) = {\tilde{I}}_{o} + {\tilde{I}}_{s s} + \int \int S ({\vec{u}}^{'}) \cdot [P (- {\vec{u}}^{'}) \cdot \tilde{V} * (- \vec{u}, - η) \cdot \tilde{G} * (- {\vec{u}}^{'} - \vec{u}, - η - η^{'}) \cdot P * (- {\vec{u}}^{'} - \vec{u}) \\ + P * (- {\vec{u}}^{'}) \cdot \tilde{V} (\vec{u}, η) \cdot \tilde{G} (- {\vec{u}}^{'} + \vec{u}, - η^{'} + η) \cdot P (- {\vec{u}}^{'} + \vec{u})] d^{2} {\vec{u}}^{'}, \end{matrix}

where

{\tilde{I}}_{o} = \int \int S ({\vec{u}}^{'}) {| P (- {\vec{u}}^{'}) |}^{2} d^{2} {\vec{u}}^{'} \cdot δ (\vec{u}, η)

is the DC term representing the background intensity.

{\tilde{I}}_{s s}

is a nonlinear term describing 2^nd order scattering interactions, which involves a convolution between

\tilde{V}

and

\tilde{V} *

. If the sample and surrounding medium have small index contrast and weak absorption, this 2^nd order term becomes negligible. Hence, under the weak object approximation (WOA), we can linearize Eq. (5) by dropping

{\tilde{I}}_{s s}

. The implication is that the linear relation only holds when a majority of the scattered energy is contributed from the interaction between the scattered field and the illumination [25,34,37], as opposed to the interaction of the scattered field with itself. This tends to be a valid approximation for transparent biological cells, which provide a strong DC term and minimal scatter-scatter interactions.

To evaluate the phase and absorption effects separately, we consider the real and imaginary parts of the scattering potential, $\tilde{V} = {\tilde{V}}_{Re} + i \cdot {\tilde{V}}_{Im}$ . Finally, we arrive at our forward model, relating the 3D intensity images under different illumination patterns to the 3D scattering potential:

{\tilde{I}}_{img} = {\tilde{I}}_{o} + H_{Re} \cdot {\tilde{V}}_{Re} + H_{Im} \cdot {\tilde{V}}_{Im},

where the 3D phase transfer function H_Re is

\begin{matrix} H_{Re} (\vec{u}, η) = \int \int S ({\vec{u}}^{'}) \cdot [P * (- {\vec{u}}^{'}) \cdot \tilde{G} (- {\vec{u}}^{'} + \vec{u}, - η^{'} + η) \cdot P (- {\vec{u}}^{'} + \vec{u}) \\ + P (- {\vec{u}}^{'}) \cdot \tilde{G} * (- {\vec{u}}^{'} - \vec{u}, - η^{'} - η) \cdot P * (- {\vec{u}}^{'} - \vec{u})] d^{2} {\vec{u}}^{'} \end{matrix}

and the 3D absorption transfer function H_Im is

\begin{matrix} H_{Im} (\vec{u}, η) = i \int \int S ({\vec{u}}^{'}) \cdot [P * (- {\vec{u}}^{'}) \cdot \tilde{G} (- {\vec{u}}^{'} + \vec{u}, - η^{'} + η) \cdot P (- {\vec{u}}^{'} + \vec{u}) \\ - P (- {\vec{u}}^{'}) \cdot \tilde{G} * (- {\vec{u}}^{'} - \vec{u}, - η^{'} - η) \cdot P * (- {\vec{u}}^{'} - \vec{u})] d^{2} {\vec{u}}^{'} . \end{matrix}

Since phase and absorption transfer functions only depend on the source S and pupil P, which are known, they can be pre-computed. We use here the first Born approximation, however the Rytov approximation could instead be used to arrive at the same result [38]. They are equivalent here because a small phase assumption (i.e. |ϕ| ≪ 1) is made in both cases to build a linear relationship between the scattering potential of the sample and the measured intensity. Note that WOTFs have been shown to be valid beyond the small-phase regime, for weakly scattering large-phase samples that have slowly varying phase/refractive index (i.e. |∇ϕ| ≪ 1) [38–40]. Here, we use non-paraxial WOTFs to also go beyond the paraxial regime [34, 38].

3D renderings and orthogonal slices of the phase and absorption transfer functions in Eq. (7) and (8) are shown in Fig. 2(a). For brightfield images, only the absorption transfer function contains non-zero values (i.e. no phase information is encoded). This is the familiar 3D incoherent transfer function in a wide-field microscope [34]. For DPC illumination, however, contrast comes mostly from phase [22] since the asymmetric illumination converts information about the real part of the scattering potential into measured intensity.

The spatial resolution of our method is determined by wavelength λ and NA. The transverse Fourier coverage of the DPC transfer function spans a bandwidth 2× larger than the NA of the objective (N A_⊥,_max = 2 × N A), so phase and absorption may be recovered with resolution that is 2× the coherent resolution limit. Axial bandwidth is also determined by NA, but is not directly proportional. It improves non-linearly as the illumination and objective NAs increase $(N A_{∥, m a x} = 1 - \sqrt{1 - N A^{2}})$ [18,33,34], as shown in Fig. 2(b). For instance, if an objective lens with 0.65 NA and a light source with λ = 500 nm are used, the Abbe diffraction-limited lateral and axial resolution are λ/(2 × N A) = 385 nm and $λ / (1 - \sqrt{1 - N A^{2}}) = 2083 nm$ , respectively.

Fig. 2 (a) Absorption (H_Im) and phase (H_Re) 3D transfer functions for the brightfield and DPC stacks. The NA of both illumination and detection is 0.65. (b) Lateral and axial resolution (measured by Fourier bandwidth normalized to units of NA) improve as the NA increases. Gray circles indicate the case shown in (a).

Download Full Size | PDF

2.2. Inverse problem

Our inverse problem estimates the sample’s 3D complex refractive index from the captured intensity measurements. In the following, we consider two reconstruction algorithms: least squares (ℓ₂) and TV regularized. The least squares algorithm minimizes the ℓ₂-norm of the difference between the actual measurements and the predicted measurements based on the current estimate and our forward model. In our problem, since both the illumination and imaging optics only cover a limited range of angles, there is a ‘missing cone’ along the axial dimension in both the phase and absorption transfer functions (see Fig. 2). Direct inversion will therefore result in high-frequency artifacts due to the missing information. ℓ₂ (Tikhonov) regularization is a standard way to alleviate this problem [25] by imposing a minimum total energy constraint,

\min_{{\tilde{V}}_{Re}, {\tilde{V}}_{Im}} \sum_{l} ∥ {\tilde{I}}_{l}^{'} - {H^{'}}_{Re, l} \cdot {\tilde{V}}_{Re} - {H^{'}}_{Im, l} \cdot {\tilde{V}}_{Im} ∥_{2}^{2} + α ∥ {\tilde{V}}_{Re} ∥_{2}^{2} + β ∥ {\tilde{V}}_{Im} ∥_{2}^{2},

where l indexes the illumination patterns corresponding to brightfield and DPC, α and β are regularization parameters, Ĩ′ denotes the Fourier transform of the normalized intensity data stack and

{H^{'}}_{Re} = H_{Re} / | {\tilde{I}}_{o} |

and

{H^{'}}_{Im} = H_{Im} / | {\tilde{I}}_{o} |

are normalized transfer functions. The brightfield image stack is computed by averaging the intensities from all DPC illumination patterns (top, bottom, left, right). We normalize each brightfield image by first subtracting a background image, and then dividing it by the (scalar) DC term. The DPC images are computed by taking the difference of the pairs of images from complimentary illumination patterns, then normalized by the illumination intensity. Mathematically, these procedures are written as:

I_{BF} = \frac{I_{t o p} + I_{b o t t o m} + I_{r i g h t} + I_{l e f t}}{2}, {I^{'}}_{BF} = \frac{I_{BF} - | {\tilde{I}}_{o} |}{| {\tilde{I}}_{o} |}, {I^{'}}_{DPC} = \frac{I_{t o p / r i g h t} - I_{b o t t o m / l e f t}}{| {\tilde{I}}_{o} |} .

where I_top, I_bottom, I_right and I_left represent measured intensities using top, bottom, right, and left half-circle illumination, respectively, I_BF is brightfield intensity and

{I^{'}}_{DPC}

is the normalized DPC intensity. One can acquire the DC term by measuring the background brightfield intensity without the sample. However, it can be approximated by the average intensity of the brightfield images when the WOA holds, i.e. |Ĩ_o | = 〈I_BF〉, which we adopt here.

Equation (9) has a closed-form solution that can be obtained by setting its derivative with respective to V_Re and V_Im to zero, giving

\begin{array}{l} V_{Re} = F^{- 1} {\frac{(\sum_{l} | {H^{'}}_{Im, l} |^{2} + β) \cdot \sum_{l} {H^{'}}_{Re, l}^{*} \cdot {\tilde{I}}_{l}^{'} - \sum_{l} {H^{'}}_{Re, l}^{*} \cdot {H^{'}}_{Im, l} \cdot \sum_{l} {H^{'}}_{Im, l} \cdot {\tilde{I}}^{'}}{(\sum_{l} | {H^{'}}_{Re, l} |^{2} + α) \cdot (\sum_{l} | {H^{'}}_{Im, l} |^{2} + β) - \sum_{l} {H^{'}}_{Re, l} \cdot {H^{'}}_{Im, l}^{*} \cdot \sum_{l} {H^{'}}_{Re, l}^{*} \cdot {H^{'}}_{Im, l}}} \\ V_{Im} = F^{- 1} {\frac{(\sum_{l} | {H^{'}}_{Re, l} |^{2} + α) \cdot \sum_{l} {H^{'}}_{Im, l}^{*} \cdot {\tilde{I}}_{l}^{'} - \sum_{l} {H^{'}}_{Re, l} \cdot {H^{'}}_{Im, l}^{*} \cdot \sum_{l} {H^{'}}_{Re, l}^{*} \cdot {\tilde{I}}^{'}}{(\sum_{l} | {H^{'}}_{Re, l} |^{2} + α) \cdot (\sum_{l} | {H^{'}}_{Im, l} |^{2} + β) - \sum_{l} {H^{'}}_{Re, l} \cdot {H^{'}}_{Im, l}^{*} \cdot \sum_{l} {H^{'}}_{Re, l}^{*} \cdot {H^{'}}_{Im, l}}}, \end{array}

where F⁻¹{·} represents inverse Fourier transform. Hence, the scattering potential can be recovered using deconvolution with fast Fourier transforms (FFTs). The two regularization parameters, α and β, should scale with the magnitude of V_Re and V_Im, respectively. Overly large regularization parameters result in reduction of absorption and phase values, while overly small regularization parameters result in high frequency noise.

Due to the missing cone problem, phase reconstructions suffer from halo artifacts (see Fig. 3(c)). To mitigate this, we exploit two different priors. First, a piece-wise constant approximation of the 3D scattering potential (i.e. the sample can be represented as multiple regions of constant refractive index) is made by using a total variation (TV) regularizer [17, 41]. Second, in our experiments, the sample is immersed in a medium whose refractive index is always smaller than the sample’s. Thus, we can further impose a positivity constraint on the refractive index. Including both priors, the 3D phase reconstruction solves the following minimization problem

\min_{{\tilde{V}}_{Re}, {\tilde{V}}_{Im}} \sum_{l} ∥ {\tilde{I}}_{l}^{'} - {H^{'}}_{Re, l} \cdot {\tilde{V}}_{Re} - {H^{'}}_{Im, l} \cdot {\tilde{V}}_{Im} ∥_{2}^{2} + γ_{Re} ∥ V_{Re} ∥_{TV} + γ_{Im} ∥ V_{Im} ∥_{TV}, s . t . V_{Re}, V_{Im} \leq 0,

where γ_Re and γ_Im are the regularization coefficients. We adopt a generalized 3D TV operator ‖·‖_TV= ∑ Φ_τ (∇·) in which Φ_τ (·) is the Huber loss function that allows a differentiable TV functional and ∇ is the 3D gradient operator, all implemented in TVReg toolbox [42].

Fig. 3 3D refractive index reconstructions of a 10µm diameter polystyrene bead immersed in oil of refractive index 1.58 (40× 0.65 NA objective). (a) Cross-sectional views of 3D brightfield and DPC measurements. (b) Absorption and (c) phase reconstructions with ℓ₂ and TV regularization. (d) 1D cross-sections of recovered refractive index.

Download Full Size | PDF

3. Experimental results

Our experimental setup is a commercial microscope (Nikon TE300) with source patterning achieved by a programmable LED array (32×32, 4 mm spacing, central λ = 513 nm). The LED array is placed at 69 mm above the focal plane in order to provide illumination NA 0.65. An automated piezo-stage (Thorlabs, MZS500-E) implements axial scanning while capturing 4 images at each position using the 4 half-circle illumination patterns shown in Fig. 1. In order to avoid aliasing, the axial spacing between focus planes within the 3D intensity stack should be smaller than the depth-of-field, $λ / (2 - 2 \sqrt{1 - N A^{2}})$ . We also need to capture the full axial range of the sample’s intensity variations. Hence, for 10–20 µm thick samples, we collect data at 100 defocus planes equally spaced by 1 µm. Note that our method requires similar amounts of data as compared to other diffraction tomography schemes, which use 4 phase-shifted interferograms (or an off-axis hologram with 4× more pixels) at each angle [43] or axial plane [21]. Both the LED array and the piezo-stage are synchronized to our camera (PCO.Edge 5.5), so acquisition speed is limited by the axial translation speed (∼ 90 seconds for 400 images).

During the reconstruction process, the 4 images at each z are first used to synthesize a brightfield and two DPC images along orthogonal directions [25]. These brightfield and DPC image stacks provide complementary information. Direct solutions from our algorithms are complex scattering potentials, which are converted to refractive index distributions. Using FFT-based implementations, reconstruction of 512×512×100 voxels using Eq. (9) on a desktop computer (Intel i7 CPU) with Matlab takes ∼10 seconds, whereas implementation of Eq. (12) typically requires about 30 iterations to converge (∼ 18 minutes). Computation times could be further reduced through parallel processing on Graphics Processing Units (GPUs).

3.1. Comparison between ℓ₂ and TV regularization

To compare reconstruction algorithms for a known test sample, we use a single polystyrene bead (Sigma-Aldrich) immersed in index-matching fluid (Cargille, index 1.58). This sample is a pure phase object with piecewise-constant refractive index, so satisfies the TV constraint. Two cross-sectional views of the 3D brightfield and DPC intensity stacks from Eq. (10) are shown in Fig. 3(a). Brightfield images appear to have absorption at the edges of the bead, where strong diffraction causes light to scatter outside the the objective’s collection aperture. Since this loss is small, the recovered absorption map using ℓ₂ (Tikhonov) regularization is noisy (Fig. 3(b)).

Looking at the reconstructed 3D refractive index (Fig. 3(c,d)), we see significantly less halo artifacts in the TV regularization case, as compared to ℓ₂, as well as less noise. Both ℓ₂ and TV reconstruction suffer from elongation along the axial direction - the axial width is ∼1.6× larger than the bead. This is due to the missing frequencies along the axial direction (see Fig. 2(a)). In order to reduce the elongation, one can use higher NA or rotate the object [38]. Alternatively, a nonlinear forward model [20,41] that accounts for multiple scattering can be used to go beyond the first Born approximation and the WOA, potentially allowing recovery of higher frequency components as well as missing cone information. However, such algorithms are not guaranteed to converge and require significantly more computation power and time. Despite using a weak object approximation, our reconstructed refractive index (∼ 1.61) matches well with previous experimental measurements [44], providing accurate quantitative phase results. In the rest of our experimental results, we will use the TV regularization algorithm.

3.2. Effect of NA

To experimentally study the effects of NA on our 3D phase reconstructions, we compare results from two objective lenses (20× 0.4 NA and 40× 0.65 NA). The test sample is composed of several polystyrene beads immersed in oil of refractive index 1.54. As expected, higher NA gives better resolution in both lateral and axial dimensions, as shown in Fig. 4. Poor depth sectioning from the lower NA also causes slowly varying background artifacts, which are mitigated in the high NA reconstructions. The improved axial resolution both reduces the axial elongation and enables observation of finer structures within the large beads (indicated by the white arrows). The diffraction-limited lateral and axial resolutions are 0.64 µm and 6.14 µm for 20× objective, and 0.39 µm and 2.14 µm for 40× objective, respectively. Figure 4 also illustrates a caveat of 3D DPC. As the index difference between the surrounding media and the beads becomes larger (Δn > 0.07 in this case), the WOA or slowly-varying phase approximation does not hold but our model attempts to fit the measurements with those assumptions. Hence, although the shape remains the same, the retrieved refractive index is lower than the actual value.

Fig. 4 Comparison of recovered 3D refractive index (6 μm polystyrene beads) with 20 × 0.4 NA and 40 × 0.65 NA objectives. Larger NA provides better lateral and axial resolution.

Download Full Size | PDF

3.3. Comparison between 2D and 3D phase reconstructions

Next, we illustrate the differences between 2D and 3D phase imaging, by comparing with two different 2D phase methods: 2D DPC [25] and the transport of intensity equation (TIE) [45]. Both use through-focus intensity stacks to recover the on-axis projected phase, ignoring sample thickness. We use a 40× 0.65 NA objective and fixed saline-immersed human mammary epithelial MCF10A cell samples. To implement TIE, we take 15 images with exponential axial spacing from 1 µm to 64 µm, using coherent (single-LED) illumination. The phase at the focal plane is recovered using a modified TIE algorithm [45]. Note that 2D DPC provides 2× better lateral resolution than TIE (see Fig. 5(a)) because it uses partially coherent illumination.

Fig. 5 Comparison between 2D phase reconstructions (TIE and 2D DPC) and 2D slices of our 3D DPC refractive index reconstruction for a human mammary epithelial MCF10A cell. (a) TIE and 2D DPC, next to three slices of the recovered refractive index from 3D DPC at three different depths. (b) 3D view of the recovered index’s 3D Fourier spectrum. (c) 3D rendering of the recovered refractive index distribution.

Download Full Size | PDF

Since the 2D phase reconstructions are proportional to the total projected optical path length of the sample, this quantity alone cannot distinguish features at different axial positions. Our 3D DPC method, on the other hand, can clearly distinguish sub-cellular features such as cytosol, cytoplasm and nucleus, which are located at different axial planes (see Fig. 5). The sharper appearance of the reconstructed slices is a result of the axial sectioning that removes out-of-focus features. 3D DPC provides the same lateral resolution as 2D DPC and the reconstructed 3D Fourier spectrum (Fig. 5(b)) covers the same volume in Fourier space as the transfer functions in Fig. 2(a). In Fig. 5(c), we show the 3D rendered refractive index reconstruction with different refractive index values mapping to different colors as in [13,21].

We also implement our method to image embryo cells (Fig. 6 and Visualization 1). These cells have diameter ∼100 µm and are surrounded by zona pellucida (bottom arrow). The top arrow points to the polar body of the oocyte. The nucleus as well as cytoplasm distribution along the axial direction are clearly observed within the region in the orange box. The blue dashed box shows a 4-cell stage embryo at three distinct focus planes, where two sets of two cells sit on top of each other, with cleavage planes on the diagonal and off-diagonal directions. Recovering the 3D information enables us to figure out the location of individual cells, providing a non-invasive way to isolate sub-cellular features with high resolution and accuracy.

Fig. 6 Reconstructed 3D refractive index of embryo cells. (Left) Full field-of-view at a single focus plane. (Right) 3 axial slices for regions in the blue dashed and orange boxes.

Download Full Size | PDF

4. Conclusion

We proposed a new quantitative 3D phase imaging technique that captures images with asymmetric partially coherent illumination at different focus planes. Under the Born approximation, we derived a linear 3D model that relates illumination-dependent intensity to 3D refractive index and absorption. The reconstruction algorithm is fast and efficient, using either least squares (ℓ₂) or TV regularization with a positivity constraint. Our experimental setup is simple and inexpensive (an LED array add-on), enabling label-free and stain-free single-cell imaging with sub-cellular feature specificity.

Acknowledgments

We thank Zachary Phillips, Li-hao Yeh, Jingshan Zhong, Olivia Scheideler and Colin Sheppard for assistance and the David and Lucile Packard Foundation (Grant 88685) for funding.

References and links

1. F. Zernike, “Phase contrast, a new method for the microscopic observation of transparent objects,” Physica 9, 686–698 (1942). [CrossRef]

2. Z. Wang, L. Millet, M. Mir, H. Ding, S. Unarunotai, J. Rogers, M. U. Gillette, and G. Popescu, “Spatial light interference microscopy (SLIM),” Opt. Express 19, 1016–1026 (2011). [CrossRef] [PubMed]

3. D. Murphy, Fundamentals of Light Microscopy and Electronic Imaging (Wiley-Liss, New York, NY, USA, 2001).

4. W. Lang, Nomarski Differential Interference-Contrast Microscopy (Oberkochen, Carl Zeiss, 1982).

5. M. R. Arnison, K. G. Larkin, C. J. R. Sheppard, N. I. Smith, and C. J. Cogswell, “Linear phase imaging using differential interference contrast microscopy,” J. Microsc. 214, 7–12 (2004). [CrossRef] [PubMed]

6. D. Paganin and K. A. Nugent, “Noninterferometric phase imaging with partially coherent light,” Phys. Rev. Lett. 80, 2586–2589 (1998). [CrossRef]

7. C. J. R. Sheppard, “Defocused transfer function for a partially coherent microscope and application to phase retrieval,” J. Opt. Soc. Am. A 21, 828–831 (2004). [CrossRef]

8. J. C. Petruccelli, L. Tian, and G. Barbastathis, “The transport of intensity equation for optical path length recovery using partially coherent illumination,” Opt. Express 21, 14430–14441 (2013). [CrossRef] [PubMed]

9. J. Zhong, L. Tian, J. Dauwels, and L. Waller, “Partially coherent phase imaging with simultaneous source recovery,” Opt. Express 6, 257–265 (2015). [CrossRef]

10. J. A. Rodrigo and T. Alieva, “Rapid quantitative phase imaging for partially coherent light microscopy,” Opt. Express 22, 13472–13483 (2014). [CrossRef] [PubMed]

11. P. Cloetens, W. Ludwig, J. Baruchel, D. V. Dyck, J. V. Landuyt, J. P. Guigay, and M. Schlenker, “Holotomography: Quantitative phase tomography with micrometer resolution using hard synchrotron radiation X rays,” Appl. Phys. Lett. 75, 2912–2914 (1999). [CrossRef]

12. Y. Sung, W. Choi, C. Fang-Yen, K. Badizadegan, R. R. Dasari, and M. S. Feld, “Optical diffraction tomography for high resolution live cell imaging,” Opt. Express 17, 266–277 (2009). [CrossRef] [PubMed]

13. Y. Cotte, F. Toy, P. Jourdain, N. Pavillon, D. Boss, P. Magistretti, P. Marquet, and C. Depeursinge, “Marker-free phase nanoscopy,” Nature Photon. 7, 113–117 (2013). [CrossRef]

14. D. Brady, K. Choi, D. Marks, and R. Horisaki, “Compressive holography,” Opt. Express 17, 13040–13049 (2009). [CrossRef] [PubMed]

15. Y. Sung and R. R. Dasari, “Deterministic regularization of three-dimensional optical diffraction tomography,” J. Opt. Soc. Am. A 28, 1554–1561 (2011). [CrossRef]

16. A. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” J. Opt. Soc. Am. A 19, 472–480 (2002). [CrossRef]

17. L. Tian, J. C. Petruccelli, Q. Miao, H. Kudrolli, V. Nagarkar, and G. Barbastathis, “Compressive X-ray phase tomography based on the transport of intensity equation,” Opt. Lett. 38, 3418–3421 (2013). [CrossRef] [PubMed]

18. E. Wolf, “Three-dimensional structure determination of semi-transparent objects from holographic data,” Opt. Commun. 1, 153–156 (1969). [CrossRef]

19. G. Gbur, M. A. Anastasio, Y. Huang, and D. Shi, “Spherical-wave intensity diffraction tomography,” J. Opt. Soc. Am. A 22, 230–238 (2005). [CrossRef]

20. L. Tian and L. Waller, “3D intensity and phase imaging from light field measurements in an LED array microscope,” Optica 2, 104–111 (2015). [CrossRef]

21. T. Kim, R. Zhou, M. Mir, S. D. Babacan, P. S. Carney, L. L. Goddard, and G. Popescu, “White-light diffraction tomography of unlabelled live cells,” Nature Photon. 8, 256–263 (2014). [CrossRef]

22. D. Hamilton and C. Sheppard, “Differential phase contrast in scanning optical microscopy,” J. Microsc. 133, 27–39 (1984). [CrossRef]

23. T. N. Ford, K. K. Chu, and J. Mertz, “Phase-gradient microscopy in thick tissue with oblique back-illumination,” Nat. Methods 9, 1195–1197 (2012). [CrossRef] [PubMed]

24. S. B. Mehta and C. J. Sheppard, “Quantitative phase-gradient imaging at high resolution with asymmetric illumination-based differential phase contrast,” Opt. Lett. 34, 1924–1926 (2009). [CrossRef] [PubMed]

25. L. Tian and L. Waller, “Quantitative differential phase contrast imaging in an LED array microscope,” Opt. Express 23, 11394–11403 (2015). [CrossRef] [PubMed]

26. G. Zheng, C. Kolner, and C. Yang, “Microscopy refocusing and dark-field imaging by using a simple LED array,” Opt. Lett. 36, 3987–3989 (2011). [CrossRef] [PubMed]

27. G. Zheng, R. Horstmeyer, and C. Yang, “Wide-field, high-resolution Fourier ptychographic microscopy,” Nature Photon. 7, 739–745 (2013). [CrossRef]

28. Z. Liu, L. Tian, S. Liu, and L. Waller, “Real-time brightfield, darkfield, and phase contrast imaging in a light-emitting diode array microscope,” J. Biomed. Opt. 19, 106002 (2014). [CrossRef] [PubMed]

29. L. Tian, Z. Liu, L.-H. Yeh, M. Chen, J. Zhong, and L. Waller, “Computational illumination for high-speed in vitro Fourier ptychographic microscopy,” Optica 2, 904–911 (2015). [CrossRef]

30. L. Tian, X. Li, K. Ramchandran, and L. Waller, “Multiplexed coded illumination for Fourier ptychography with an LED array microscope,” Biomed. Opt. Express 5, 2376–2389 (2014). [CrossRef] [PubMed]

31. L. Tian, J. Wang, and L. Waller, “3D differential phase-contrast microscopy with computational illumination using an LED array,” Opt. Lett. 39, 1326–1329 (2014). [CrossRef] [PubMed]

32. X. Ou, G. Zheng, and C. Yang, “Embedded pupil function recovery for Fourier ptychographic microscopy,” Opt. Express 22, 4960–4972 (2014). [CrossRef] [PubMed]

33. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light, 7th ed. (Cambridge University Press, 1999). [CrossRef]

34. N. Streibl, “Three-dimensional imaging by a microscope,” J. Opt. Soc. Am. A 2, 121–127 (1985). [CrossRef]

35. Y. Sung and C. J. R. Sheppard, “Three-dimensional imaging by partially coherent light under non-paraxial condition,” J. Opt. Soc. Am. A 28, 554–559 (2011). [CrossRef]

36. T. H. Nguyen, C. Edwards, L. L. Goddard, and G. Popescu, “Quantitative phase imaging of weakly scattering objects using partially coherent illumination,” Opt. Express 24, 11683–11693 (2016). [CrossRef] [PubMed]

37. R. A. Claus, P. P. Naulleau, A. R. Neureuther, and L. Waller, “Quantitative phase retrieval with arbitrary pupil and illumination,” Opt. Express 23, 26672–26682 (2015). [CrossRef] [PubMed]

38. M. H. Jenkins and T. K. Gaylord, “Three-dimensional quantitative phase imaging via tomographic deconvolution phase microscopy,” Appl. Opt. 54, 9213–9227 (2015). [CrossRef] [PubMed]

39. J. Guigay, “Fourier transform analysis of Fresnel diffraction patterns and in-line holograms,” Optik 49, 121–125 (1977).

40. Y. I. Nesterests and T. E. Gureyev, “Partially coherent contrast-transfer-function approximation,” J. Opt. Soc. Am. A 33, 464–474 (2016). [CrossRef]

41. U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “Learning approach to optical tomography,” Optica 2, 517–522 (2015). [CrossRef]

42. T. L. Jensen, J. H. Joergensen, P. C. Hansen, and S. H. Jensen, “Implementation of an Optimal First-Order Method for Strongly Convex Total Variation Regularization,” BIT Numer. Math. 52, 329–356 (2012). http://www.imm.dtu.dk/~pcha/TVReg/ [CrossRef]

43. 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] [PubMed]

44. S. Jones, M. King, and A. Ward, “Determining the unique refractive index properties of solid polystyrene aerosol using broadband Mie scattering from optically trapped beads,” Phys. Chem. 15, 20735–20741 (2013).

45. Z. Jingshan, R. A. Claus, J. Dauwels, L. Tian, and L. Waller, “Transport of intensity phase imaging by intensity spectrum fitting of exponentially spaced defocus planes,” Opt. Express 22, 10661–10674 (2014). [CrossRef] [PubMed]

3D differential phase contrast microscopy

Abstract

1. Introduction