White light diffraction tomography of unlabeled live cells

Abstract

Methods for obtaining a tomographic phase image of a specimen, either in transmission or in scatter. A specimen is illuminated by a temporally incoherent source and light collected in transmission or scattering is used to generate a scattered phase image of the specimen in multiple axial planes. The scattered field is solved for in wavevector space, and a derived instrument function is deconvolved to obtain specimen susceptibility in wavevector space. The specimen susceptibility is transformed to obtain a three-dimensional phase tomogram of the specimen.

Description

TECHNICAL FIELD

The present invention relates to methods of quantitative phase imaging, and, more particularly, to tomographic imaging using white light illumination of substantially transparent specimens.

DEFINITIONS

Objectives and techniques of quantitative phase imaging (QPI) are described in detail in Popescu, Quantitative Phase Imaging of Cells and Tissues, (McGraw-Hill, 2011), incorporated herein by reference. In QPI, the optical path length associated with substantially transparent specimens is measured and mapped, and translated into biomedically relevant information.

Off-axis phase imaging methods have provided the fastest acquisition rates by virtue of the fact that phase information, and thus optical path lengths, are extracted from a single recorded interferogram, as described, for example in Ikeda et al., Hilbert phase microscopy for investigating fast dynamics in transparent systems, Opt. Lett., vol. 30, pp. 1165-67 (2005), incorporated herein by reference. While diffraction-limited transverse resolution is intrinsically preserved in phase shifting methods, off-axis techniques may degrade transverse resolution. The highest temporal phase sensitivity (that is, the smallest frame-to-frame phase shift) is provided by common-path methods because they are the most stable. Finally, the highest spatial phase sensitivity (i.e., the smallest point-to-point phase change within the same frame) is obtained in the absence of speckles, which implies an incoherence in at least some respect—of the source of illumination.

Diffraction phase microscopy (DPM), as described in Popescu et al., Diffraction phase microscopy for quantifying cell structure and dynamics, Opt. Lett., vol. 31, pp. 775-77 (2006), incorporated herein by reference, is both off-axis and common-path, and thus combines the benefits of fast acquisition rates and high temporal sensitivity. These features enabled DPM to perform unprecedented biological studies, especially related to red blood cell membrane dynamics. Diffraction phase microscopy with white light is the subject of US Published Patent Application US 2014/0085715, which is incorporated herein by reference.

Spatial light interference microscopy (SLIM), described, for example, in U.S. Pat. No. 8,184,298 (hereinafter, “Popescu '298”), entitled “Spatial Light Interference Microscopy,” and by Wang et al., Spatial Light Interference Microscopy (SLIM), Opt. Exp., vol. 19, pp. 1016-26 (2011), both of which are incorporated herein by reference, combines Zernike's phase contrast microscopy and Gabor's holography, has been shown to be of great value in measuring nanoscale structures and dynamics in live cells by generating endogenous contrast and mapping quantitative optical path-lengths at each point in the image.

SUMMARY OF EMBODIMENTS OF THE INVENTION

In accordance with embodiments of the invention, methods are provided for obtaining a tomographic phase image of a substantially transparent specimen. The methods have steps of:

• a. illuminating the specimen with substantially temporally incoherent light through an objective of high-numerical aperture, the temporally incoherent light characterized by a focus;
• b. stepping the focus of the temporally incoherent light to a sequence of successive focal planes within the substantially transparent specimen;
• c. forming a scattered phase image of the specimen in an image plane where the scattered phase image of the specimen corresponds to each of the successive focal planes;
• d. solving for the scattered field in wavevector space;
• e. deconvolving a derived instrument function to obtain specimen susceptibility in wavevector space; and
• f. transforming the specimen susceptibility to derive a three-dimensional phase tomogram of the specimen.

In accordance with alternate embodiments of the invention, the specimen may be substantially transparent in a visible portion of the electromagnetic spectrum, or, alternatively, the specimen may be opaque in the visible portion of the electromagnetic spectrum. The specimen may be substantially transparent in an infrared portion of the electromagnetic spectrum. More particularly, the specimen may be a biological cell.

In accordance with further embodiments of the present invention, forming a scattered phase image of the specimen in an image plane may correspond to measuring a temporal cross-correlation function between a scattered field and a reference plane wave, which, more particularly, may traverse the specimen.

In accordance with yet further embodiments of the present invention, forming a scattered phase image may include combining a plurality of interferograms obtained with distinct phase retardation rings. Forming a phase image may include obtaining a single interferogram from an off-axis type interferometer, such as a diffraction phase microscopy, in particular.

BRIEF DESCRIPTION OF THE DRAWINGS

The foregoing features of the invention will be more readily understood by reference to the following detailed description, taken with reference to the accompanying drawings, in which:

FIG. 1 is a general schematic depiction of a SLIM module used to perform a tomographic phase measurement in accordance with an embodiment of the present invention;

FIGS. 2A-2D show interferograms obtained with corresponding rings of phase retardation in FIGS. 2E-2H, processed, in accordance with SLIM, to obtain the phase image shown in FIG. 2I.

FIG. 3A depicts light scattering by a scattering potential under the first-order Born approximation; FIG. 3B is a three-dimensional rendering of an instrument transfer function using a WDT calculation in accordance with an embodiment of the present invention, while FIG. 3C is a cross-sectional slice at the k_y=0 plane;

FIG. 4 is a schematic depiction of optical phase sectioning of a specimen in accordance with embodiments of the present invention;

FIGS. 5A and 5B show the middle frame of a z-stack of raw data acquired in imaging an E. coli cell and a reconstructed z-stack, in accordance with an embodiment of the present invention. The 3D rendering of this E. coli cell is shown in FIG. 5C.

FIG. 6 shows a false-coloured 3D image of an HT-29 cell, employing WDT methods in accordance with an embodiment of the present invention.

DETAILED DESCRIPTION OF EMBODIMENTS OF THE INVENTION

Further details regarding the invention described herein may be obtained by reference to Kim et al., White-light diffraction tomography of unlabeled live cells, Nature Photonics, DOI: 10.1038/NPHOTON.2013.350 (Jan. 19, 2014) (hereinafter, “Kim (2014)”), and in the Supplementary Information appended thereto, both of which are incorporated herein by reference.

DEFINITIONS

The term “image” shall refer to any multidimensional representation, whether in tangible or otherwise perceptible form, or otherwise, whereby a value of some characteristic (amplitude, phase, etc.) is associated with each of a plurality of locations corresponding to dimensional coordinates of an object in physical space, though not necessarily mapped one-to-one thereonto. Thus, for example, the graphic display of the spatial distribution of some field, either scalar or vectorial, such as brightness or color, constitutes an image. So, also, does an array of numbers, such as a 3D holographic dataset, in a computer memory or holographic medium. Similarly, “imaging” refers to the rendering of a stated physical characteristic in terms of one or more images.

As used herein, “label-free” refers to a method of imaging a substantially transparent sample that does not require introduction of materials, such as fluorophores or contrast agents that are extrinsic to the imaged sample.

The terms “object,” “sample,” and “specimen” shall refer, interchangeably, to a tangible, non-transitory physical object capable of being rendered as an image.

The term “temporally incoherent” as applied to a source of irradiation shall refer to a source characterized by a coherence time τ_c that is no longer than 10 cycles of a central frequency of the light emitted by the source.

The term “substantially high spatial coherence” as applied to a beam of light shall denote a condition of coherence wherein if two slits were to be placed within the a coherence area exceeding the field of view in a plane transverse to the beam propagation axis, at least one interference null would result in the far field where the intensity is no greater than 10% of the peak intensity of the beam.

“High numerical aperture,” for purposes of the present description and claims, shall refer to NA>0.5, which for visible light insures a transverse resolution better than one micron.

Comment Regarding Notation:

For convenience of notation in the current description, the same symbol is used for a function and for its Fourier transform. In order to indicate the domain in which the function operates, all the arguments are carried explicitly, e.g., f(k_⊥,z; ω) is the Fourier transform of f(r_⊥, z; t) over r_⊥ and t.

Embodiments of the present invention, significantly, measure a complex field at an image plane. Complex field measurements in an image plane have been shown, as in Ding et al., Fourier transform light scattering of inhomogeneous and dynamic structures, Phys. Rev. Lett., vol. 101, 238102 (2008), incorporated herein by reference, to yield higher sensitivity than scatter measurements in the far zone. Once a complex field has been measured at successive depths to form a z-stack of complex images, the complex field stack is processed using white light diffraction tomography (WDT) in accordance with the following teachings of the present invention.

FIG. 1 shows a schematic of a SLIM module 100, designed as an add-on module to a commercial phase contrast microscope 102 (otherwise referred to herein as the “microscope”), such as an Axio Observer Z1, which may be obtained from Carl Zeiss Microscopy, LLC. The back focal plane of the objective lens (not shown) contains a thin metal annulus (the “phase ring”) which provides a π/2 phase shift between the scattered and unscattered light. In one embodiment of a SLIM system, the image plane of the phase contrast microscope 102 is first relayed and magnified through a 4f system formed by two focusing lenses, L1 and L2. The back focal plane of the objective lens is imaged onto a spatial light modulator (SLM) 104, such as a programmable liquid crystal spatial light modulator that may be obtained from Boulder Nonlinear Systems, Inc., Boulder, Colo., which is located in the Fourier plane of a second 4f system formed by L3 and L4. A desired ring pattern is projected onto the SLM, and is matched to the image of the phase ring on the SLM. In this manner the SLM is used to provide additional phase shifts of 0, π/2, π and 3π/2 (shown in FIGS. 2E-2H) between the scattered and un-scattered light. Thus, four phase shifted interferograms, shown in FIGS. 2A-2D, are recorded at the image plane (at the focus of the second Fourier lens L4) using a CCD detector 104. The SLM 104, CCD 106, and microscope 102 are synchronized through a processor 108 executing a software program. In a typical implementation, 8 frames are acquired per second, corresponding to 2 quantitative phase images (QPI) per second, although it is to be understood that other data acquisition rates may be practiced within the scope of the present invention. The acquisition rate is typically limited by the frame rate of the CCD and the refresh rate of SLM. At the rate discussed here by way of example, a z-stack of 100 slices can be obtained in 50 seconds. By using these four interferograms, the actual phase shift (shown in FIG. 2I) caused by the sample can be uniquely determined within a range of 2π. Any 2π ambiguities are corrected by Goldstein's unwrapping algorithm, which is described in Goldstein et al., Satellite radar interferometry: Two-dimensional phase unwrapping, Radio Science, vol. 23, pp. 713-20 (1988), which is incorporated herein by reference. Determination of the phase φ(r, t) of the imaged object based on the four phase-shifted interferograms is described, for example, in Popescu '298.

A sample (not shown) is transilluminated, in microscope 102, by a broadband source 103. In one example, the broadband source is a halogen lamp with a temporal coherence length of 1.2 μm. The broadband source 103 (also referred to, herein, as a “lamp”) assumes a significant role in increasing resolution since it does not suffer from speckle, and thus, improves the sensitivity.

Insofar as the SLIM module 100 is preferably an add-on module to a commercial microscope, it is possible to overlay SLIM with other microscopy modalities (e.g. epi-fluorescence, or DIC), within the scope of the present invention. In particular, within the scope of the present invention, a phase image may be obtained from a single interferogram from an off-axis type interferometer, such as a diffraction phase microscope, for example.

With these features, SLIM is capable of performing multimodal and functional studies. Furthermore, with the short coherence length and a high numerical aperture objective, SLIM may advantageously provide excellent depth sectioning (1.2 μm depth sectioning for a 63×/1.4 NA objective), which plays a large role in the reconstruction method described herein.

It is known that the Rytov approximation is more appropriate for reconstructing smooth objects with respect to the wavelength of light, i.e., for low-resolution imaging, while the Born approximation works better for imaging finer structures. Thus, for purposes of resolving features of imaged cells, for example, the latter is appropriate, and is used for purposes of the discussion herein. Under the first-order Born approximation, illustrated in FIG. 3A, with an incident plane wave, U_i=A(ω)e^iβ(ω)z, the forward scattered field U_s is solved in the wavevector space, as described above, instead of using the traditional Green's function and Weyl's formula approach.

For purposes of the present teachings, a temporal cross-correlation function, Γ₁₂ (r,τ), is usefully defined as follows:Γ₁₂(r,τ)=U_s(r,t)U_r*(r,t+τ),  (1)where U_s and U_r are the scattered and the reference (unscattered) field, respectively. Γ₁₂(r,τ) is a complex function that can be written asΓ₁₂(r,τ)=|Γ₁₂(r,τ)|e^iφ(r,τ),  (2)where φ(r,τ) is the phase associated with the image field, measured in SLIM using the four phase shifted interferograms, as discussed above. The generalized Wiener-Khintchine theorem (see Wiener, Generalized harmonic analysis, Acta Mathematica, vol. 55, pp. 117-258 (1930), incorporated herein by reference) allows Γ₁₂(r,τ) to be related to the cross-spectral density through Fourier transform of

$\begin{array}{cc}{\Gamma }_{12}\left(r,\tau \right)={\int }_0^{\infty }{W}_{12}\left(r,\omega \right)e^{\mathrm{i\omega \tau }}d\omega \mathrm{where},& \left(3\right)\\ W_{12}\left(r,\omega \right)=\langle U_s\left(r,\omega \right)U_r^{*}\left(r,\omega \right)\rangle .& \left(4\right)\end{array}$

During the phase shifting measurement, Γ₁₂(r,τ) is evaluated at τ=0, to yield:

$\begin{array}{cc}{\Gamma }_{12}\left(r,\tau =0\right)={\int }_0^{\infty }\langle U_s\left(r,\omega \right)U_r^{*}\left(r,\omega \right)\rangle d\omega ,& \left(5\right)\end{array}$ in which the unknown quantity is the scattered field, U_s. In the following, an analytical solution is derived for U_s which establishes the relationship between the measurable quantity, Γ₁₂, and the 3D object structure of interest.

The inhomogeneous Helmholtz equation, which describes the field U in a medium with index distribution of n(r) is∇²U(r,ω)+n²(r,ω)β₀²(ω)U(r,ω)=0,  (6)where β₀ (ω)=ω/c is the wavenumber in vacuum. Eq. (6) may be rearranged as∇²U(r,ω)+β²(ω)U(r,ω)=−χ(r,ω)β_o²(ω)U(r,ω),  (7)where β(ω)=nβ₀(ω) (we assume non-dispersive objects), n=n(r)_r is the spatially averaged refractive index, and χ(r,ω)=n²(r,ω)−n²(ω) is the scattering potential. The total field U(r,ω) can be written as U(r,ω)=U_i(r,ω)+U_s(r,ω), where U_i(r,ω) is the incident wave. U_i satisfies the homogeneous wave equation,∇²U_i(r,ω)+β²(ω)U_i(r,ω)=0.  (8)Eq. (8) has a plane-wave solution U_i(r,ω)=A(ω)e^iβ(ω)z where A(ω) is the spectral amplitude of the incident field. Subtracting Eq. (7) from Eq. (8) gives∇²U_s(r,ω)+β²(ω)U_s(r,ω)=−χ(r)β_o²(ω)U(r,ω).  (9)Eq. (9) makes clear that the driving term of the right hand side represents the interaction of the object scattering potential χ with the total field U. Under the first Born approximation, i.e., |U_s(r,ω)|<<|U_i(r,ω)|, U(r,ω) on the right hand side can be approximated as U_i(r,ω), allowing Eq. (9) to be re-written as,∇²U_s(r,ω)+β²(ω)U_s(r,ω)=−χ(r,ω)β_o²(ω)A(ω)e^iβ(ω)z.  (10)

Instead of employing the traditional Green's function approach and the angular spectrum representation (Weyl's formula), the scattered field directly may be solved for directly in wavevector space, using the 3D Fourier transformation. First, a 3D Fourier transform of Eq. (10) is performed, which yields[β²(ω)−k²]U_s(k,ω)=−β₀²(ω)A(ω)χ[k_⊥,k_z−β(ω)].  (11)

In Eq. (11), we used the shift theorem of Fourier transforms, namely, χ(r,ω)e^iβ(ω)z→χ[k_⊥,k_z−β(ω)], where the arrow indicates Fourier transformation. Note that, throughout the manuscript, we use the same symbol for a function and its Fourier transform but carry all the arguments explicitly, which clearly identifies the domain in which the function operates. For example, χ(k) is the Fourier transform of χ(r). From Eq. (1), the scattered field U_s s obtained immediately as

$\begin{array}{cc}{U}_s\left(k,\omega \right)=-{\beta }_0^2\left(\omega \right)A\left(\omega \right)\frac{\chi \left[k_{\perp },k_z-\beta \left(\omega \right)\right]}{{\beta }^2\left(\omega \right)-k_{\perp }^2-k_z^2}.& \left(12\right)\end{array}$

Next, an expression is derived for the field as a function of axial distance z, i.e., we arrange the terms such that a 1D inverse Fourier transform with respect to k_z can be easily performed. Toward this end, a k_⊥-dependent propagation constant is defined, q=√{square root over (β²(ω)−k_⊥²)}, Eq. (12) is rewritten as

$\begin{array}{cc}{U}_s\left(k,\omega \right)=-{\beta }_0^2\left(\omega \right)A\left(\omega \right)\chi \left[k_{\perp },k_z-\beta \left(\omega \right)\right]\frac{1}{2q}\left(\frac{1}{k_z-q}-\frac{1}{k_z+q}\right).& \left(13\right)\end{array}$

In instances where the transmitted or forward-scattered field is acquired, the backscattered field (the second term in the right-most parentheses) may be ignored. Conversely, if the backscattered field is exclusively acquired, the forward scattering term may be ignored. However, implementation of the methods of the present invention to encompass the backscattered field is similarly within the capability of a person of ordinary skill in the art, once the present teachings have been appreciated. Solely for heuristic purposes, the teachings that follow will be expressed in terms of a forward scattered field, though their application to backscattered fields is within the scope of the present invention.

Performing an inverse Fourier transform on Eq. (13) with respect to k_z in order to obtain the forward-scattered field (for example) as a function of transverse wavevector, k_⊥, axial distance, z, and angular frequency, ω, one obtains:

$\begin{array}{cc}{U}_s\left(k_{\perp },z;\omega \right)=-\frac{{\beta }_0^2\left(\omega \right)A\left(\omega \right)}{2q}\left[\chi \left(k_{\perp },z\right)e^{\mathrm{i\beta }\left(\omega \right)z}\right]z_{}{e}^{i\mathrm{qz}}& \left(14\right)\end{array}$ In Eq. (14), {circle around (v)}_z indicates convolution along the z-dimension,

$f\left(z\right)z_{}{e}^{i\mathrm{qz}}={\int }_{-\infty }^{\infty }f\left(z^{\prime }\right)e^{iq\left(z-z^{\prime }\right)}d{z}^{\prime }.$ can be easily seen that the convolution of a function with a complex exponential yields the Fourier transform of that function multiplied the complex exponential, which yields the simple result [χ(k_⊥,z; ω)e^iβ(ω)z]{circle around (v)}_ze^iqz=e^iqzχ[k_⊥,q−β(ω)].

In order to insert the result of Eq. (14) into Eq. (3), Γ₁₂(r,τ) is Fourier transformed with respect to the transverse position vector, r_⊥=(x,y). Since U_r(r,ω) is a plane wave propagating in z-direction, Eq. (a4) in k_⊥-space isW₁₂(k_⊥,z,ω)=U_s(k_⊥,z,ω)U_r*(z,ω).  (15)Using the solution of U_s from Eq. (b9) and U_r(z,ω)=A(ω)e^iβ(ω)z, one has

$\begin{array}{cc}{W}_{12}\left(k_{\perp },z,\omega \right)=\frac{-{\beta }_o^2\left(\omega \right)S\left(\omega \right)e^{iz\left[q-\beta \left(\omega \right)\right]}}{2q}\chi \left[k_{\perp },q-\beta \left(\omega \right)\right],& \left(16\right)\end{array}$ where S(ω)=|A(ω)|² is the power spectrum of the illumination field. Using Eq. (3), one obtains the temporal cross correlation at zero-delay as a function of the frequency integral

$\begin{array}{cc}{\Gamma }_{12}\left(k_{\perp },z,0\right)=-{\int }_0^{\infty }\frac{{\beta }_o^2\left(\omega \right)S\left(\omega \right)e^{iz\left[q-\beta \left(\omega \right)\right]}}{2q}\chi \left[k_{\perp },q-\beta \left(\omega \right)\right]d\omega .& \left(17\right)\end{array}$ With the relation β(ω)=nβ₀=nω/c, the integral from dω to dβ, S(ω) can be changed to S(βc/n), such that Eq. (17) becomes

$\begin{array}{cc}{\Gamma }_{12}\left(k_{\perp },z,0\right)=-\frac{c}{2{\stackrel{\_}{n}}^3}{\int }_0^{\infty }\frac{{\beta }^{2}S\left(\beta c/\stackrel{\_}{n}\right)e^{iz\left(q-\beta \right)}}{q}\chi \left(k_{\perp },q-\beta \right)d\beta .& \left(18\right)\end{array}$

In practicing the present invention, what is measured is S(λ). So, to obtain the spectrum distribution for S(βc/n) from S(λ), one considers the Jacobian transformations, S(βc/n)nS(β)/c and S(β)λ²S(λ)/2π. In order to evaluate the integral in Eq. (18), a new variable Q=q−β=√{square root over (β²−k_⊥²)}−β is defined, from which are derived

$\beta =-\frac{Q^2+k_{\perp }^2}{2Q}$ $\mathrm{and}$ $q=\left(-\frac{Q}{2}+\frac{k_{\perp }^2}{2Q}\right).$ Substituting dβ for dQ, the Jacobian becomes

$\frac{d\beta }{dQ}=\left(-\frac{1}{2}+\frac{k_{\perp }^2}{2Q^2}\right),$ and Eq. (18) becomes:

$\begin{array}{cc}\begin{array}{c}{\Gamma }_{12}\left(k_{\perp },z;0\right)=\frac{1}{8{\stackrel{\_}{n}}^2}{\int }_0^{\infty }\frac{{\left(Q^2+k_{\perp }^2\right)}^2}{Q^3}S\left(-\frac{Q^2+k_{\perp }^2}{2Q}\right)\chi \left(k_{\perp },Q\right)e^{i\mathrm{zQ}}dQ\\ =\frac{1}{8{\tilde{n}}^2}{\mathrm{FT}}_Q^{-1}\left[\frac{{\left(Q^2+k_{\perp }^2\right)}^2}{Q^3}S\left(-\frac{Q^2+k_{\perp }^2}{2Q}\right)\right]z_{}\chi \left(k_{\perp },-z\right)\\ =\Sigma \left(k_{\perp },-z\right)z_{}\chi \left(k_{\perp },-z\right).\end{array}& \left(19\right)\end{array}$

In Eq. (19), Σ is the function that incorporates all the details of the instrument response. Note that the k_⊥ coverage of Σ is limited to a maximum value k_⊥^max=β₀NA, where NA is the numerical aperture of the objective. By measuring z-stacks in SLIM, it becomes possible to reconstruct the object's 3D distribution through deconvolution of Eq. (19). Alternatively, Eq. (19) may be written in the spatial frequency domain as a product, namely,

$\begin{array}{cc}{\Gamma }_{12}\left(k_{\perp },Q;0\right)=\Sigma \left(k_{\perp },Q\right)\chi \left(k_{\perp },Q\right),\mathrm{where},& \left(20\right)\\ \Sigma \left(k_{\perp },Q\right)=\frac{1}{8{\stackrel{\_}{n}}^2}\frac{{\left(Q^2+k_{\perp }^2\right)}^2}{Q^3}S\left(-\frac{Q^2+k_{\perp }^2}{2Q}\right).& \left(21\right)\end{array}$ Calculating the Point Spread Function (PSF)

Since Eq. (21) gives the coherent transfer function, i.e., the Fourier transform of the point spread function (PSF) of the imaging system, therefore, by measuring the source spectrum and filtering according to the NA of the objective, the coherent transfer function of the system can be calculated directly. The PSF of the imaging system may also be referred to herein equivalently as an “instrument function.”

First, the optical spectrum of the lamp 103 is measured at a specified temperature, such as 3200K, using a fiber optic spectrometer, for example. Insofar as the spectrum SW is measured as in air, we perform the Jacobian transformation introduced in the previous section to obtain the spectrum in terms of the variable β=−(Q²+k_⊥²)/2Q. Next through numerical calculations, this spectrum is resampled onto a 3D grid in spatial-frequency space (k_x, k_y, Q). For a 63×/1.4 NA objective, each pixel in space corresponds to 45 nm and each z-slice is separated by 200 nm. Because of the quadratic relationship between Q and β, the resampling yields two duplicates of the spectrum. Therefore, the second of the two duplicates which appears at high Q is removed by applying a spatial low-pass filter with cutoff Q=k_⊥. Further filtering in Q is performed to incorporate the physical minimum and maximum value of Q and the maximum value, which is determined by the NA of the objective, k_⊥^max=β₀NA. In order to avoid well known numerical artifacts due to sharp cutoffs in the frequency domain, we used a simple apodization procedure to smooth the edges of the filter function. The filter function may be convolved with a narrow Gaussian function, of a width that is 5% of the system's maximum transverse frequency. Finally, the coherent transfer function is Fourier transformed to obtain the PSF. The complex PSF is then used for the 3D reconstruction, as detailed below.

Deconvolution

As discussed above, SLIM measures the 3D complex field distribution, i.e. phase and amplitude, which can be expressed as a convolution between the point spread function (PSF) of the imaging system, Σ(r), and the structure of the object, χ(r)=n²(r)−n². The real field may be represented as a complex analytic signal, Ū(r)=χ(r){circle around (v)}Σ(r)+ζ(r), which also contains a signal independent noise term, ζ(r). Considering a weakly scattering phase object with very small absorption, where most of the useful information is contained in the phase image, it is reasonable to assume that the amplitude of the field is constant over the entire space. Therefore, the deconvolution can be performed on the phase term only, exp[iΦ(r)].

The function that is minimized, in accordance with a preferred embodiment of the invention, is

$\begin{array}{cc}\hat{\Phi }\left(r\right)=\underset{\Phi \left(r\right)}{\mathrm{argmin}}\frac{1}{2{\sigma }^2}{\exp \left[i\stackrel{\_}{\Phi }\left(r\right)\right]-\Sigma \left(r\right)\exp \left[\mathrm{i\Phi }\left(r\right)\right]}^2+\rho R\left[\Phi \left(r\right)\right]& \left(22\right)\end{array}$ where σ² is the noise variance, ρ is the regularization parameter, and R is the regularization functional. For simplification, these purely phase-dependent fields can be expressed in vector forms:

$\begin{array}{cc}\hat{f}=\underset{f}{\mathrm{argmin}}\frac{1}{2{\sigma }^2}{g-\Sigma f}^2+\rho R\left(f\right)& \left(23\right)\end{array}$ where g represents the measured field, exp[iΦ(r)], f represents the unknown field from the structure, exp[iΦ(r)], and Σ represents the convolution matrix corresponding to the PSF.

Phase contrast (PC) images are highly sensitive to the sharp object boundaries, but not to slow variations in the background region. For small scale objects, such as a biological cell, these characteristics of PC permit the use the sparse representation, described, for example, in Shechtman et al., Sparsity-based sub-wavelength imaging with partially incoherent light via quadratic compressed sensing, Opt. Express, vol. 19, pp. 14807-22 (2011), incorporated herein by reference. When an appropriate transform is applied to a phase image, only a few of the transform coefficients contain most of the signal energy while all the other coefficients become very small. This situation is known as transform sparsity. In a preferred embodiment of the present invention, the first- and second-order directional difference operators, [−1 1] and [−1 2 −1], along with 45° and −45° first-order derivative filters,

$\left[\begin{array}{cc}-1& 0\\ 0& 1\end{array}\right]\mathrm{and}\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right],$ are used as transforms. For each plane in the image, x-y, y-z and z-x, these 2D transforms are applied and total of 12 transforms are generated to hold the spatial variations. Therefore, by applying the sparse deconvolution principle with these filters, the problem in Eq. (23) can be expressed as

$\begin{array}{cc}\hat{f},{\hat{\alpha }}_{\mathrm{ki}}=\underset{f,{\alpha }_{\mathrm{ki}}}{\mathrm{argmin}}\frac{1}{{\sigma }^2}{g-\Sigma f}^2+\sum _k{\left(D_{k}f\right)}^T{A}_k\left(D_{k}f\right)& \left(24\right)\end{array}$ where α_ki are the weighting coefficients for each plane, i, D_k are the transform matrices, and A_k are diagonal matrices with α_ki in the diagonal. This problem is solved by an alternating iterative minimization scheme where only one unknown is estimated at a time while others are fixed. As a result, the complex image f is estimated by taking derivative of Eq. (24) and setting it to zero. Therefore, the solution is of the structure is,

$\begin{array}{cc}\hat{f}={\left({\Sigma }^T\Sigma +{\sigma }^2\sum _k{D}_k^T{A}_k{D}_k\right)}^{-1}{\Sigma }^{T}g& \left(25\right)\end{array}$ and the weighting coefficients are estimated to be,

$\begin{array}{cc}{\hat{\alpha }}_{\mathrm{ki}}=\frac{1}{{\left(D_k\hat{f}\right)}_i^2+\varepsilon }& \left(26\right)\end{array}$ where ε is a small number 10⁻⁹, in a preferred embodiment) used to avoid the trivial solution.Application to Tomographic Reconstruction

It is known that It is known that the Rytov approximation is more appropriate for reconstructing smooth objects with respect to the wavelength of light, i.e., for low-resolution imaging, and the Born approximation works better for imaging finer structures. Thus, for purposes of resolving features of imaged cells, for example, the latter is appropriate, and is used for purposes of the discussion herein. Under the first-order Born approximation, with an incident plane wave, U_i=A(ω)e^iβ(ω)z, the forward scattered field U_s is solved in the wavevector space, as described above, instead of using the traditional Green's function and Weyl's formula approach. In the transverse wavevector domain, k_⊥, U_s can be expressed, per Eq. (13), as

$\begin{array}{cc}{U}_s\left(k_{\perp },z;\omega \right)=-\frac{{\beta }_0^2\left(\omega \right)A\left(\omega \right)e^{i\mathrm{qz}}}{2q}\chi \left[k_{\perp },q-\beta \left(\omega \right)\right],& \left(27\right)\end{array}$ with symbols as defined above. The dispersion in the object may be neglected, because most biological samples of interest here are weakly absorbing. This is true even for single red blood cells. Even though the haemoglobin absorbs strongly in blue, the overall absorption of visible light through a single red blood cell is very small. This is so because the absorption length of haemoglobin in a normal red blood cell is around 10 μm in blue (averaged over 400-500 nm wavelength) and around 3 mm in red (averaged over 600-750 nm wavelength), while the thickness of the cell is only 2-3 μm.

In conventional phase shifting interferometry, the cross-correlation of the scattered and reference fields is measured as Γ₁₂(r_⊥,z,τ)=U_s(r_⊥,z,t)U_r*(z,t+τ) at τ=0 which is equivalent to integrating the spectrum over ω. Knowledge of the spatial frequency response of our instrument, or coherent transfer function (see Supplementary Information), Σ(k_x, k_y, k_z), allows us to write the main result of our calculation, i.e., the solution to the inverse scattering problem, in terms of the measured data, Γ₁₂, and the instrument function (or the coherent transfer function), Σ, in the wave-vector domain as,

$\begin{array}{cc}\chi \left(k\right)=\frac{{\Gamma }_{12}\left(k;0\right)}{\Sigma \left(k\right)}.& \left(28\right)\end{array}$

In practice, the operation in Eq. (28) requires regularization, as discussed above. Σ(k) is given by

$\begin{array}{cc}\Sigma \left(k\right)=\frac{1}{8{\stackrel{\_}{n}}^2}\frac{{\left(Q^2+k_{\perp }^2\right)}^2}{Q^3}S\left(-\frac{Q^2+k_{\perp }^2}{2Q}\right).& \left(29\right)\end{array}$

In Eq. (29), S is the optical spectrum of the imaging field as a function of the wavenumber and Q=√{square root over (β²−k_⊥²)}−β. The 3D point spread function can be obtained, as shown above, through an inverse Fourier transform of Eq. (29).

Qualitatively, Σ has a physically intuitive behaviour, of which salient features are depicted in FIGS. 3B and 3C. Specifically, its dependence on z is related to the optical spectrum, S(ω), via a Fourier transform, meaning that a broader optical spectrum gives a narrower function, Σ(z). This relationship explains the inherent optical sectioning capabilities of the instrument. This type of optical gating, in which axial resolution is determined by the coherence properties of the illumination light, has been successfully employed in optical coherence tomography (OCT) of deep tissues, for example. However, there are significant differences between WDT and OCT. In OCT, the cross-correlation Γ₁₂ is resolved over a broad delay range, which provides the depth dimension of the object. In WDT, the z-information is collected by scanning the object through focus.

Most importantly, in methods in accordance with the present invention, the coherence gating works in synergy with the high numerical aperture (NA) optics and, thus, allows for high resolution tomography. In other words, in WDT, coherence gating by itself would not work at zero NA and, conversely, high NA gating would not work with monochromatic light.

Use of broadband light from source 103 and high numerical aperture objectives results in optical sectioning capabilities highly advantageous for high resolution tomography. The function Σ(k_x, k_y, k_z) for an imaging system in accordance with embodiments of the present invention is illustrated in FIGS. 3B and 3C: a three-dimensional rendering in FIG. 3B, and a cross-section at the k_y=0 plane, in FIG. 3C. As expected, the width of the k coverage increases with k_x, which indicates that the sectioning is stronger for finer structures or, equivalently, higher scattering angles. The structure of the object is recovered through a sparse deconvolution algorithm, as described above. Acquisition of phase images, in each plane of a z-stack of images, proceeds in accordance with the paradigm discussed above and depicted in FIG. 4. Experimentally, Σ was measured by imaging a microsphere much smaller than the resolution of the system, in accord with theoretical expectations.

Example I

Tomography of Escherichia coli

The approach described herein in accordance with the present invention was applied to obtain three-dimensional images of Escherichia coli (E. coli) cells. A z-stack consisting of 17 slices with a step size of 280 nm was obtained. FIGS. 5A and 5B show the middle frame of raw data and the reconstructed z-stack. Upon deconvolution, the previously invisible protein helical subcellular structure of E. coli is resolved. The 3D rendering of this E. coli cell is shown in FIG. 5C. Only the bottom half of the cell is shown in the figure to emphasize the helical subcellular structure.

WDT was also employed to obtain 3D images of speculated red blood cells, shown in Kim (2014), to which the interested reader is directed.

Example II

Tomography of HT29 Cells

In order to study more complex subcellular structures, WDT was used to image human colon adenocarcinoma cells (HT29), using a 63×/1.4 NA oil immersion objective. The z-stack consisting of 140 frames was acquired in 150 nm z-steps. The results obtained on a cell that has recently divided are summarized in FIG. 6 which shows a false-coloured 3D image of an HT-29 cell, where sub-cellular structures are clearly observed: cell membrane in blue, nuclei in green, and nucleoli in red. These areas are first chosen by thresholding based on the phase values (0.6 rad for nucleoli and 0.1 rad for membrane) and then detailed based on morphology.

In preferred embodiments of the present invention, the disclosed methods for tomographic quantitative phase imaging of transparent structures using white light diffraction phase microscopy may be implemented as a computer program product for use with a computer system. Such implementations may include a series of computer instructions fixed either on a tangible medium, such as a computer readable medium (e.g., a diskette, CD-ROM, ROM, or fixed disk) or transmittable to a computer system, via a modem or other interface device, such as a communications adapter connected to a network over a medium. The medium may be either a tangible medium (e.g., optical or analog communications lines) or a medium implemented with wireless techniques (e.g., microwave, infrared or other transmission techniques). The series of computer instructions embodies all or part of the functionality previously described herein with respect to the system. Those skilled in the art should appreciate that such computer instructions can be written in a number of programming languages for use with many computer architectures or operating systems. Furthermore, such instructions may be stored in any memory device, such as semiconductor, magnetic, optical or other memory devices, and may be transmitted using any communications technology, such as optical, infrared, microwave, or other transmission technologies. It is expected that such a computer program product may be distributed as a removable medium with accompanying printed or electronic documentation (e.g., shrink wrapped software), preloaded with a computer system (e.g., on system ROM or fixed disk), or distributed from a server or electronic bulletin board over the network (e.g., the Internet or World Wide Web). Of course, some embodiments of the invention may be implemented as a combination of both software (e.g., a computer program product) and hardware. Still other embodiments of the invention are implemented as entirely hardware, or entirely software (e.g., a computer program product).

The embodiments of the invention described above are intended to be merely exemplary; numerous variations and modifications will be apparent to those skilled in the art. All such variations and modifications are intended to be within the scope of the present invention as defined in any appended claims.

Claims

1. A method for obtaining a tomographic phase image of a specimen, the method comprising: a. illuminating the specimen with substantially temporally incoherent light through an objective of high-numerical aperture, the temporally incoherent light characterized by a focus;b. stepping the focus of the temporally incoherent light to a sequence of successive focal planes within the specimen;c. forming a scattered phase image of the specimen in an image plane where the scattered phase image of the specimen corresponds to each of the successive focal planes;d. solving for the scattered field in wavevector space;e. deconvolving a derived instrument function to obtain specimen susceptibility in wavevector space; andf. transforming the specimen susceptibility to derive a three-dimensional phase tomogram of the specimen.

2. A method in accordance with claim 1, wherein the specimen is substantially transparent in a visible portion of the electromagnetic spectrum.

3. A method in accordance with claim 1, wherein the specimen is opaque in a visible portion of the electromagnetic spectrum.

4. A method in accordance with claim 1, wherein the specimen is substantially transparent in an infrared portion of the electromagnetic spectrum.

5. A method in accordance with claim 1, wherein forming a scattered phase image of the specimen in an image plane corresponds to measuring a temporal cross-correlation function between a scattered field and a reference plane wave.

6. A method in accordance with claim 5, wherein the reference plane wave traverses the specimen.

7. A method in accordance with claim 1, wherein the specimen is a biological cell.

8. A method in accordance with claim 1, wherein forming a scattered phase image includes combining a plurality of interferograms obtained with distinct phase retardation rings.

9. A method in accordance with claim 1, wherein forming a phase image includes obtaining a single interferogram from an off-axis type interferometer.

10. A method in accordance with claim 9, wherein the off-axis type interferometer is a diffraction phase microscope.