We demonstrate 3D phase and absorption recovery from partially coherent intensity images captured with a programmable LED array source. Images are captured through-focus with four different illumination patterns. Using first Born and weak object approximations (WOA), a linear 3D differential phase contrast (DPC) model is derived. The partially coherent transfer functions relate the sample’s complex refractive index distribution to intensity measurements at varying defocus. Volumetric reconstruction is achieved by a global FFT-based method, without an intermediate 2D phase retrieval step. Because the illumination is spatially partially coherent, the transverse resolution of the reconstructed field achieves twice the NA of coherent systems and improved axial resolution.
© 2016 Optical Society of America
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 ; however, a single image is not quantitative. To measure phase quantitatively, PhC  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  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  and phase contrast , 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.
2. Principles of 3D differential phase contrast (DPC) microscopy
2.1. Forward model
A 3D sample can be characterized by its scattering potential , where no is the refractive index of the surrounding media and n = nRe + i · nIm is the sample’s complex refractive index. The real part of the refractive index, nRe, describes how phase delays accumulate as light passes through the sample, while the imaginary part, nIm, describes absorption.
We start by describing the output field of a 3D sample under illumination by an incident field . The product of the two is convolved with the 3D Green’s function G, which characterizes the 3D field distribution from a point source,33]),
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 described by , where represents 2D transverse frequencies, η′ is 1D axial frequency and . After passing through the sample, the field is low-pass filtered by the pupil function of the microscope, corresponding to the 2D coherent point spread function . The measured intensity is the incoherent sum of intensities from each point source. Under the first Born approximation and Eq. (2), the measured intensity isEq. (3). Taking a 3D Fourier transform of both sides of Eq. (3), we arrive at a Fourier space relation between the 3D measured intensity·spectrum and the scattering potential, Eq. (4) into the following three terms: Eq. (5) by dropping . 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, . Finally, we arrive at our forward model, relating the 3D intensity images under different illumination patterns to the 3D scattering potential:
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 . 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 . For DPC illumination, however, contrast comes mostly from phase  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 [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 , respectively.
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 (ℓ2) and TV regularized. The least squares algorithm minimizes the ℓ2-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. ℓ2 (Tikhonov) regularization is a standard way to alleviate this problem  by imposing a minimum total energy constraint,
Equation (9) has a closed-form solution that can be obtained by setting its derivative with respective to VRe and VIm to zero, giving
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 problem42].
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, . 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  or axial plane . 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 . 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 ℓ2 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 ℓ2 (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 ℓ2, as well as less noise. Both ℓ2 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 . 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 , 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.
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  and the transport of intensity equation (TIE) . 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 . Note that 2D DPC provides 2× better lateral resolution than TIE (see Fig. 5(a)) because it uses partially coherent illumination.
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.
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 (ℓ2) 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.
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]
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]
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]
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]
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]
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]
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]
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]