1. INTRODUCTION

It is challenging to image thick samples with a standard microscope. High-resolution objective lenses offer a shallow depth of field, which require one to axially scan through the sample to visualize its three-dimensional (3D) shape. Unfortunately, refocusing does not remove light from areas above and below the plane of interest. This longstanding problem has inspired a number of solutions, the most widespread being confocal designs, two-photon excitation methods, light sheet microscopy, and optical coherence tomography. These methods “gate out” light from sample areas away from the point of interest and offer excellent signal enhancement, especially for thick, fluorescent samples [1].

Such gating techniques also encounter several problems. First, they typically must scan out each image, which might require physical movement and can be time consuming. Second, the available signal (i.e., the number of ballistic photons) decreases exponentially with depth. To overcome this limit, one must use a high NA lens, which provides a proportionally smaller image field of view (FOV). Finally, little light is backscattered when imaging non-fluorescent samples that are primarily transparent, such as commonly seen in embryology, in model organisms such as zebrafish, and after the application of recent tissue clearing [2] and expansion [3] techniques.

Instead of capturing just the ballistic photons emerging from the sample, one might instead image the entire optical field, which includes light that has scattered. Several techniques have been proposed to enable depth selectivity without gating for ballistic light. One might perform optical sectioning through digital deconvolution of a focal stack [4]. Light-field imaging [5] and point-spread function engineering [6] are two other alternatives. All three of these methods primarily operate with incoherent light, e.g., from fluorescent samples. They are thus not ideal tools for obtaining the refractive index distribution of a primarily transparent and non-fluorescent medium.

To do so, it is useful to use coherent illumination. For example, the amplitude and phase of a digital hologram may be computationally propagated to different depths within a thick sample, much like refocusing a microscope. However, the field at out-of-focus planes still influences the final result. Several techniques also aim for depth selectivity by using quasi-coherent illumination or through acquiring multiple images [7–10].

A very useful framework to summarize how coherent light scatters through thick samples is diffraction tomography (DT), as first developed by Wolf [11]. In a typical DT experiment, one illuminates a sample of interest with a series of tilted plane waves and measures the resulting complex diffraction patterns in the far field. These measurements may then be combined with a suitable algorithm into a tomographic reconstruction. An early demonstration of DT by Lauer is a good example [12]. Typically, the reconstruction algorithm assumes the first Born [13,14] or the first Rytov [15] approximation. It is also possible to apply the projection approximation, which models light as a ray. As a synthetic aperture technique, DT comes with the additional benefit of improving the resolution of an imaging element beyond its traditional diffraction-limit cutoff [12].

However, as a technique that models both the amplitude and phase of a coherent field, most implementations of DT require a reference beam and holographic measurement, or some sort of phase-stable interference (including spatial light modulator coding strategies, e.g., as in Ref. [16]). Since it is critical to control for interferometric stability [13] and thus limit the motion and phase drift to sub-micrometer variations, DT has been primarily implemented in well-controlled, customized setups. Several prior works have considered solving DT from intensity-only measurements to possibly remove the need for a reference beam [17–26]. However, while some applied the first Born approximation, these works also required customized setups and typically imposed additional sample constraints (e.g., a known sample support). They did not operate within a standard microscope or connect their reconstruction algorithms to ptychography, from which improvements like computational aberration correction [27] and multiplexed reconstruction [28] may be easily adopted.

Here, we perform DT using standard intensity images captured under variable LED illumination from an array source. Our technique, termed Fourier ptychographic tomography (FPT), acquires a sequence of images while changing the light pattern displayed on the LED array. Then, it combines these images using a phase retrieval-based ptychographic reconstruction algorithm, which computationally segments a thick sample into multiple planes, as opposed to physically rejecting light from above and below one plane of interest. Similar to DT, FPT also improves the lateral image resolution beyond the standard cutoff of the imaging lens. The end result is an accurate three-dimensional map of the complex index of refraction of a volumetric sample obtained directly from a sequence of standard microscope images.

2. RELATED WORK

To begin, a number of techniques attempt 3D imaging without applying the first Born approximation. These include lensless on-chip devices [29], lensless setups that assume an appropriate linearization [30], and methods relying upon effects like defocusing (e.g., the transport of intensity equation [31]) or spectral variations [32]. These techniques do not necessarily fit within a standard microscope setup or offer the ability to simultaneously improve spatial resolution. Two related works for 3D imaging, which also do not use DT with the first Born approximation, are by Tian and Waller [33] and Li et al. [34]. These two setups are quite similar to ours, and we discuss their operations in more detail below.

As mentioned above, there are also several prior works applying DT under the first Born approximation that only use intensity measurements [17–26]. While some of these works examine phase retrieval as a reconstruction algorithm, they must either shift the focal plane [22,24] or source [25] axially between each measurement, or must assume constraints on the sample [19,20,23,26] to successfully recover the phase. These prior works do not connect DT to ptychographic phase retrieval (i.e., they do not recover the phase by using diversity between the variably illuminated DT images). Connections between phase retrieval and DT under the first Born approximation have also been explored within the context of volume hologram design [35].

The field of x ray ptychography also offers a number of methods to image 3D samples with intensity measurements [36–38]. However, none of these ptychography methods seem to directly modify DT under the first Born or Rytov approximation, to the best of our knowledge. A popular technique appears to use standard two-dimensional (2D) ptychographic solvers to determine the complex field for individual projections of a slowly rotated sample, which are subsequently combined using conventional DT techniques [39].

Fourier ptychography (FP) [40] uses a standard microscope and no moving parts to simultaneously improve image resolution and measure quantitative phase but is restricted to thin samples. FPT effectively extends FP into the third dimension. As noted above, Tian and Waller [33] and Li et al. [34] also examine the problem of 3D imaging from intensities in a standard microscope. These two examples adopted their reconstruction technique from a 3D ptychography method [37,38] that splits up the sample into a specified number of infinitesimally thin slices (each under the projection approximation) and applies the beam propagation method (i.e., assumes small-angle scattering) [41]. This “multi-slice” approach remains accurate within a different domain of optical scattering than the first Born approximation (i.e., for different types of samples and setups, see chapter 2 in Ref. [42]). For example, one may include both forward and backscattered light in a DT solver to accurately reconstruct a 3D sample under the Born approximation [12]. However, the multi-slice method does not directly account for backscattered light. Its projection approximation also assumes the lateral divergence of the optical field gradient at each slice is zero. Alternatively, the validity of first Born approximation breaks down when the total amount of absorption and phase shift from a sample is large [43], whereas this appears to impact the multi-slice approach less (e.g., it can image two highly absorbing layers separated by a finite amount of free space [33,34]).

Thus, while the experimental setup of FPT is similar to prior work [33,34,40], our reinterpretation of ptychography within the physical framework of DT (under the first Born approximation) allows us to accurately reconstruct new specimen types in 3D without measuring phase. For example, primarily transparent samples of continuously varying optical density, which are often encountered in biology, typically obey the first Born approximation. Accordingly, we have used FPT to compute some of the first quantitatively accurate 3D maps of clear samples that contain contiguous features (e.g., an unstained nematode parasite and starfish embryo) from standard microscope images.

In addition, FPT offers a clear picture of the location and amount of data it captures in 3D Fourier space. Such knowledge currently helps us to establish various solution guarantees for 2D image phase retrieval [44]. These guarantees may also extend to the current case of 3D tomographic phase retrieval (e.g., to support its potential clinical use). Furthermore, instead of specifying an arbitrary number of sample slices and their location in a 3D volume, FPT simply inserts measured data into its appropriate 3D Fourier space location and ensures phase consistency between each measurement. By solving for the first term in the Born expansion, we aim this approach as a general framework to eventually form quantitatively accurate tomographic maps of complex biological samples with sub-micrometer resolutions.

3. METHOD OF FPT

In this section, we develop a mathematical expression for our image measurements using the FPT framework and then summarize our reconstruction algorithm. We use the vector $\mathbf{r}=(r_x,r_y,r_z)$ to define the 3D sample coordinates and the vector $\mathbf{k}=(k_x,k_y,k_z)$ to define the corresponding k-space (wave vector) coordinates (see Fig. 1).

 

Fig. 1. Setup for Fourier ptychographic tomography (FPT). (a) Labeled diagram of the FPT microscope, including optical functions of interest. (b) FPT captures multiple images under varied LED illumination. (c) A ptychography-inspired algorithm combines these images in a 3D k-space representation of the complex sample. (d) FPT outputs a 3D tomographic map of the complex index of refraction of the sample. Included images are experimental measurements from a starfish embryo (real index component, threshold applied; see Fig. 7).

Download Full Size | PDF

A. Image Formation in FPT

It is helpful to begin our discussion by introducing a quantity termed the scattering potential, which contains the complex index of refraction of an arbitrarily thick volumetric sample,

Here, $n(\mathbf{r})$ is the spatially varying and complex refractive index profile of the sample, $n_b$ is the index of refraction of the background (which we assume is constant), and $k=2\pi /\lambda$ is the wavenumber in vacuum. We note that $n(\mathbf{r})=n_r(\mathbf{r})+1i·n_{im}(\mathbf{r})$, where $n_r$ is associated with the sample’s refractive index, $n_{im}$ is associated with its absorptivity, and we define $1i=\sqrt{-1}$ for notational clarity. We typically neglect the dependence of $n$ on $\lambda$ since we illuminate with quasi-monochromatic light. This dependence cannot be neglected when imaging with polychromatic light. Finally, we use the term “thick” for samples that do not obey the thin sample approximation, which requires the sample thickness to be much less than $2/k{\theta }_{\max }^2$, where ${\theta }_{\max }$ is the magnitude of the maximum scattering angle [45].

Next, to understand what happens to light when it passes through this volumetric sample, we define the complex field that results from illuminating the thick sample, $U(\mathbf{r})$, as a sum of two fields: $U(\mathbf{r})=U_i(\mathbf{r})+U_s(\mathbf{r})$. Here, $U_i(\mathbf{r})$ is the field “incident” upon the sample (i.e., from one LED) and $U_s(\mathbf{r})$ is the resulting field that “scatters” off of the sample. We may insert this decomposition into the scalar wave equation for light propagating through an inhomogeneous medium and use Green’s theorem to determine the scattered field as [11]

Here, $G(|{\mathbf{r}}^{\prime }-\mathbf{r}|)$ is the Green’s function connecting light scattered from various sample locations, denoted by $\mathbf{r}$, to an arbitrary location ${\mathbf{r}}^{\prime }$. $V(\mathbf{r})$ is the scattering potential from Eq. (1). Since $U(\mathbf{r})$ is unknown at all sample locations, it is challenging to solve Eq. (2). Instead, it is helpful to apply the first Born approximation, which replaces $U(\mathbf{r})$ in the integrand with $U_i(\mathbf{r})$. This approximation assumes that $U_i(\mathbf{r})\gg U_s(\mathbf{r})$. It is the first term in the Born expansion that describes the scattering response of an arbitrary sample [11]. It assumes a weakly scattering medium. Specifically, the first Born approximation remains valid when the relative index shift $\delta n=|n(\mathbf{r})-n_b|$ and sample thickness $t$ obey the relation, $kt\delta n/2\ll 1$ [43]. We expect FPT to remain quantitatively accurate with samples obeying this condition. By including higher-order terms, the above framework may in principle include samples with stronger scattering [46,47].

Our system sequentially illuminates the sample with an LED array, which contains $q=q_x\times q_y$ sources positioned a large distance $l$ from the sample (in a uniform grid, with inter-LED spacing $c$; see Fig. 1). It is helpful to label each LED with a 2D counter variable $(j_x,j_y)$, where $-q_x/2\le j_x\le q_x/2$ and $-q_y/2\le j_y\le q_y/2$, as well as a single counter variable $j$, where $1\le j\le q$. Assuming each LED acts as a spatially coherent and quasi-monochromatic source (central wavelength $\lambda$) placed at a large distance from the sample, the incident field takes the form of a plane wave traveling at a variable angle such that ${\theta }_{jx}={\tan }^{-1}(j_x·c/l)$ and ${\theta }_{jy}={\tan }^{-1}(j_y·c/l)$ with respect to the $x$- and $y$-axes, respectively. We may express the $j$th field incident upon the sample as

As ${\theta }_{jx}$ and ${\theta }_{jy}$ vary, ${\mathbf{k}}_j$ will always assume values along a spherical shell in 3D $(k_x,k_y,k_z)$ space (i.e., the Ewald sphere), since the value of $k_{jz}$ is a deterministic function of $k_{jx}$ and $k_{jy}$.

After replacing $U(\mathbf{r})$ in Eq. (2) with $U_i^{(j)}(\mathbf{r})$ from Eq. (3) and additionally approximating the Green’s function $G$ as a far-field response, the following relationship emerges between the scattering potential $V$ and the Fourier transform of the $j$th scattered field, ${\widehat{U}}_s^{(j)}(\mathbf{k})$, in the far field [11]:

We refer to $\widehat{V}(\mathbf{k})$ as the k-space scattering potential, which is the three-dimensional Fourier transform of $V(\mathbf{r})$, with $\mathbf{k}$ the scattered wave vector in the far field. For simplicity, we have left out a multiplicative pre-factor $(-1i\pi /k_z)$ on the right-hand side of Eq. (5), and instead assume it is included within the function $\widehat{V}$ for the remainder of this presentation. The field scattered by the sample and viewed at a large distance, ${\widehat{U}}_s^{(j)}(\mathbf{k})$, is given by the values along a specific manifold (or spherical “shell”) of the k-space scattering potential, here written as $\widehat{V}(\mathbf{k}-{\mathbf{k}}_j)$. We illustrate the geometric connection between $\widehat{V}(\mathbf{k}-{\mathbf{k}}_j)$ and ${\widehat{U}}_s^{(j)}(\mathbf{k})$ for a 2D optical geometry in Fig. 2(b). The center of the $j$th shell is defined by the incident wave vector, ${\mathbf{k}}_j$. For a given shell center, each value of $\widehat{V}(\mathbf{k}-{\mathbf{k}}_j)$ lies on a spherical surface at a radial distance of $|\mathbf{k}|=k$ [see colored arcs in Fig. 2(b)]. As ${\mathbf{k}}_j$ varies with the changing LED illumination, the shell center shifts along a second shell with the same radius [since ${\mathbf{k}}_j$ is itself constrained to lie on an Ewald sphere; see gray circle in Fig. 2(b)].

 

Fig. 2. Mathematical summary of FPT. (a) The field from the $j$th LED scatters through the sample and exits its top surface as $U_j(x^{\prime })$ (in 2D). This field forms ${\widehat{U}}_j(k_x)$ at the microscope back focal plane, where it is bandlimited by the microscope aperture $a(k_x)$ before propagating to the image plane to form the $j$th sampled intensity image. (b) Under the first Born approximation, each detected image is the squared magnitude of the Fourier transform of one colored “shell” in $(k_x,k_z)$ space. (c) By filling in this space with a ptychographic phase-retrieval algorithm, FPT reconstructs the complex values within the finite bandpass volume ${\widehat{V}}_e(k_x,k_z)$ (color indicates expected bowl overlap for this example). The Fourier transform of this reconstruction yields our complex refractive index map with resolutions $\mathrm{\Delta }x$ and $\mathrm{\Delta }z$ along $x$ and $z$.

Download Full Size | PDF

The goal of DT is to determine all the complex values within the volumetric function $\widehat{V}$ from a set of $q$ scattered fields, ${\{{\widehat{U}}_s\}}_{j=1}^q$, which is often measured holographically [12,15]. Each 2D holographic measurement maps to the complex values of $\widehat{V}$ along one 2D shell. The values from multiple measurements [i.e., the multiple shells in Fig. 2(b)] can be combined to form a k-space scattering potential estimate, ${\widehat{V}}_e$. Nearly all stationary optical setups will yield only an estimate, since it is challenging to measure data from the entire k-space scattering potential without rotating the sample. Figures 1(c) and 2(c) display typical measurable volumes, also termed a bandpass, from a limited-angle illumination and detection setup. Once sampled, an inverse 3D Fourier transform of the band-limited ${\widehat{V}}_e(\mathbf{k})$ yields the desired complex scattering potential estimate, $V_e(\mathbf{r})$, which contains the quantitative index of refraction.

In FPT, we do not measure the scattered fields holographically. Instead, we use a standard microscope to detect image intensities and apply a ptychographic phase-retrieval algorithm to solve for the unknown complex potential. The scattered fields in Eq. (5) are defined at the microscope objective back focal plane (i.e., its Fourier plane), whose 2D coordinates ${\mathbf{k}}_{2\mathrm{D}}=(k_x,k_y)$ are Fourier conjugate to the microscope focal plane coordinates $(x,y)$. If we neglect the effect of the constant background plane wave term (i.e., $U_i$ in the sum $U=U_i+U_s$), we may now write the $j$th shifted field at our microscope back focal plane as ${\widehat{U}}^{(j)}({\mathbf{k}}_{2\mathrm{D}})=\widehat{V}({\mathbf{k}}_{2\mathrm{D}}-{\mathbf{k}}_{j2\mathrm{D}},k_z-k_{jz})$. These new coordinates highlight the 3D to 2D mapping from $\widehat{V}$ to $\widehat{U}$, where again $k_z=\sqrt{k-k_x^2-k_y^2}$ is a deterministic function of ${\mathbf{k}}_{2\mathrm{D}}$, and the same applies between $k_{jz}$ and ${\mathbf{k}}_{j2\mathrm{D}}$.

Each shifted, scattered field is then bandlimited by the microscope aperture function, $a({\mathbf{k}}_{2\mathrm{D}})$, before propagating to the image plane. The limited extent of $a({\mathbf{k}}_{2\mathrm{D}})$ (defined by the imaging system NA) sets the maximum extent of each shell along $k_x$ and $k_y$. The $j$th intensity image acquired by the detector is given by the squared Fourier transform of the bandlimited field at the microscope back focal plane:

Here, $F$ denotes a 2D Fourier transform with respect to ${\mathbf{k}}_{2\mathrm{D}}$, and we neglect the effects of magnification (for simplicity) by assuming the image plane coordinates match the sample plane coordinates, $(x,y)$. The goal of FPT is to determine the complex 3D function $\widehat{V}$ from the real, non-negative data matrix $g(x,y,j)$. A final 3D Fourier transform of $\widehat{V}$ yields the desired scattering potential, and subsequently the refractive index distribution, of the thick sample.

B. FPT Reconstruction Algorithm

Equation (6) closely resembles the data matrix measured by FP [40], but now the intensities are sampled from shells within a 3D space (i.e., the curves in Fig. 2). We use an iterative reconstruction procedure, mirroring that from FP [40], to “fill in” the k-space scattering potential with data from each recorded intensity image. Ptychography and FP require at least approximately 50%–60% data redundancy (i.e., overlapping measurements in k-space) to ensure the successful convergence of the phase retrieval process [48]. With such a similar problem structure, FPT will also require overlap between shell regions in 3D k-space. With one extra dimension, overlap is less frequent and more images are needed for an accurate reconstruction. Both a smaller LED array pitch and a larger array-sample distance along $z$ increase the amount of k-space overlap. An example cross section of FPT k-space overlap is shown in Fig. 2(c). As we demonstrate experimentally, several hundred images are sufficient for a complex reconstruction that offers a $4\times$ increase in resolution along $(x,y)$ and contains approximately 30 unique axial slices. Additional overlap (i.e., more images across the same angular range) will increase robustness to noise.

It is important to select the correct limits and discretization of 3D k-space (i.e., the FOV and resolution of the complex sample reconstruction). The maximum resolvable wave vector along $k_x$ and $k_y$ is proportional to $k({\mathrm{NA}}_o+{\mathrm{NA}}_i)$, where ${\mathrm{NA}}_o$ is the objective NA and ${\mathrm{NA}}_i$ is maximum NA of LED illumination. This lateral spatial resolution limit matches FP [49]. The maximum resolvable wave vector range along $k_z$ is also determined as a function of the objective and illumination NA as $k_z^{\max }=k(2-\sqrt{1-{\mathrm{NA}}_o^2}-\sqrt{1-{\mathrm{NA}}_i^2})$. As shown in Ref. [12], this relationship is easily derived from the geometry of the k-space bandpass volume in Fig. 2. We typically specify the maximum imaging range along the axial dimension, $z_{\max }$, to approximately match twice the expected sample thickness. This then sets the discretization level along $k_z$: $\mathrm{\Delta }{k}_z=2\pi /z_{\max }$. The total number of resolved slices along $z$ is set by the ratio $k_z^{\max }/\mathrm{\Delta }{k}_z$.

We now summarize the FPT reconstruction algorithm:

In practice, we also implement a pupil function recovery procedure [27] as we update each extracted shell from k-space, which helps remove possible microscope aberrations. As with other alternating projections-based ptychography solvers, the per iteration cost of the above FPT algorithm is $O(n\log n)$, using the big-$O$ notation. Additional details regarding algorithm robustness and convergence are in Supplement 1.

4. RESULTS AND DISCUSSION

We experimentally verify our reconstruction technique using a standard microscope outfitted with an LED array. The microscope uses an infinity corrected objective lens (${\mathrm{NA}}_o=0.4$, Olympus MPLN, $20\times$) and a digital detector containing 4.54 μm pixels (Prosilica GX 1920, $1936\times 1456$ pixel count). The LED array contains $31\times 31$ surface-mounted elements (model SMD3528, center wavelength $\lambda =632\mathrm{nm}$, 20 nm approximate bandwidth, 4 mm LED pitch, 150 μm active area diameter). We position the LED array 135 mm beneath the sample to create a maximum illumination NA of ${\mathrm{NA}}_i=0.41$. This leads to an effective lateral NA of ${\mathrm{NA}}_o+{\mathrm{NA}}_i=0.81$ and a lateral resolution gain along $(x,y)$ of slightly over a factor of 2 (from a 1.6 μm minimum resolved spatial period in the raw images to a 0.78 μm minimum resolved spatial period in the reconstruction). The associated axial Nyquist resolution is computed at 3.7 μm. We reconstruct samples across a total depth range of approximately $z_{\max }=110\mathrm{\mu m}$, which is approximately 20 times larger than the stated objective lens DOF of 5.8 μm.

For most of the reconstructions presented below, we capture and process $q=675$ images from the same fixed pattern of LEDs. Some LEDs from within the $31\times 31$ array produce images containing a shadow of the microscope back focal plane. We do not use these LEDs so as to avoid certain reconstruction artifacts. We typically use the following parameters for FPT reconstruction: each raw image is cropped to $1000\times 1000$ pixels, the reconstruction voxel size is $0.39\mathrm{\mu m}\times 0.39\mathrm{\mu m}\times 3.7\mathrm{\mu m}$ for Nyquist–Shannon rate sampling, the reconstruction array contains approximately $2100\times 2100\times 30$ voxels (110 μm total depth), and the algorithm runs for 5 iterations. The first Born approximation should remain valid across this total imaging depth. The primarily transparent and somewhat sparse samples that we test next justify the relatively large 110 μm depth that we typically use.

A. Quantitative Verification

First, we verify the ability of FPT to improve the lateral image resolution. The sample consists of 800 nm-diameter microspheres (index of refraction $n_s=1.59$) immersed in oil (index of refraction $n_o=1.515$). We highlight a small group of these microspheres in Fig. 3. The single raw image in Fig. 3(a) (generated from the center LED) cannot resolve the individual spheres gathered in small clusters. Based upon the coherent Sparrow limit for resolving two points ($0.68\lambda /{\mathrm{NA}}_o$), this raw image cannot resolve points that are closer than 1.1 μm. After FPT reconstruction, we obtain the complex index of refraction in Fig. 3(b), where we show the real component of the recovered index. The FPT reconstruction along the $\mathrm{\Delta }z=0$ slice clearly resolves the spheres within each cluster. This 800 nm distance is close to the expected Sparrow limit for the FPT reconstruction: $0.68\lambda /({\mathrm{NA}}_o+{\mathrm{NA}}_i)=540\mathrm{nm}$. The ringing features around each sphere indicate a jinc-like point-spread function, as expected theoretically from the circular shape of the finite FPT bandpass along $k_x,k_y$ [in Fig. 1(c)] for all reconstructions. The resulting constructive interference forms the undesired dip feature at the center of each cluster.

 

Fig. 3. Improved lateral resolution with FPT. (a) Single raw image of 0.8 μm microspheres. Beads within each cluster are not resolved. (b) Refractive index (real) from $\mathrm{\Delta }z=0$ slice (1 of 30) of the FPT reconstruction, resolving each microsphere.

Download Full Size | PDF

Second, we check the quantitative accuracy of FPT by imaging microspheres that extend across more than just a few reconstruction voxels. Figure 4 displays a reconstruction of 12 μm diameter microspheres (index of refraction $n_s=1.59$) immersed in oil (index of refraction $n_o=1.58$). We use the same data capture and post-processing steps as in Fig. 3. Here, we display a cropped section ($200\times 200\times 15$ voxels) of the full 3D reconstruction, which required 221 seconds of computation time on a standard laptop. We again display the real (non-absorptive) component of the recovered index across both a lateral slice (along the $\mathrm{\Delta }z=0$ plane) and a vertical slice (along the $\mathrm{\Delta }y=25\mathrm{\mu m}$ plane). We also include detailed 1D traces along the center of the vertical slice.

 

Fig. 4. FPT quantitatively measures refractive index in 3D. (a) Tomographic reconstruction of 12 μm microspheres in oil with lateral ($\mathrm{\Delta }z=0$) slice on left, axial ($\mathrm{\Delta }y=25\mathrm{\mu m}$) slice on right, and one-dimensional plots of index shift along both $x$ and $z$. (b) Digitally propagated FP reconstruction (middle) and refocused light field (right) created from the same data. FPT (left) best matches the expected spherical bead profile.

Download Full Size | PDF

Three observations here are noteworthy. First, the measured index shift approximately matches the expected shift of $\mathrm{\Delta }n=n_s-n_o=0.01$ across the entire bead, thus demonstrating quantitatively accurate performance across this limited volume. Currently, we only expect quantitative recovery for samples that meet the first Born approximation condition ($kt\delta n/2\ll 1$). This 12 μm microsphere sample approximately satisfies the condition ($kt\delta n/2=0.59$). Second, for any given one-dimensional trace through a microsphere center, we would ideally expect a perfect rect function (between $\mathrm{\Delta }n=0$ to $\mathrm{\Delta }n=0.01$). This is unlike 2D FP, which reconstructs the phase delay though each sphere and forms a parabolically shaped phase measurement (due to the varying thickness of each sphere along the optical axis). While FPT resolves an approximate step function through the center of the sphere along the lateral ($x$) dimension, it does not along the axial ($z$) direction. This is caused by the limited volume of 3D k-space that FPT measures (i.e., the limited bandpass or “missing cone” of information surrounding the $k_z$ axis). While our stationary sample/detector setup cannot avoid this missing cone, various methods are available to computationally fill it in [51].

Finally, we compare FPT with two alternative techniques for 3D imaging in Fig. 4(c). First, we use the same dataset to perform 2D FP and then holographically refocus its reconstructed optical field. We obtain this FP reconstruction using the same number of images ($q=675$) and follow the procedure in Ref. [40] after focusing the objective lens at the axial center of the 12 μm microspheres. The “out-of-focus noise” above and below the plane of the microsphere, created by digital propagation of the complex field via the angular spectrum method, noticeably hides its spherical shape. Second, we interpret the same raw image set as a light field and perform light-field refocusing [5]. While the refocused light field approximately resolves the outline of microsphere along $z$, it does not offer a quantitative picture of the sample interior, nor a measure of its complex index of refraction. The areas above the microsphere are very bright due to its lensing effect (i.e., the light field displays the optical intensity at each plane and thus displays high energy where the microsphere focuses light).

We verify the axial resolution of FPT in Fig. 5 using a sample containing two closely separated layers of 2 μm microspheres ($n_s=1.59$) distributed across the surface of a glass slide with oil in between ($n_o=1.515$). The axial separation between the two microsphere layers, measured from the center of each sphere along $z$, is 3.9 μm [i.e., the separation between the microscope slide surfaces is 5.9 μm; see Fig. 5(a)]. This almost matches the expected axial resolution limit of 3.7 μm for the FPT microscope.

 

Fig. 5. Testing the axial resolution of FPT. (a) The sample contains two layers of microspheres separated by a thin layer of oil. Raw images (b) focused at the center of the two layers and (c) on the top layer do not clearly resolve overlapping microspheres. (d)–(f) Slices of the FPT tomographic reconstruction, showing $|\mathrm{\Delta }n|$, clearly resolve each sphere within the two individual sphere layers.

Download Full Size | PDF

Conventional microscope images of the sample, using the center LED for illumination, are in Figs. 5(b) and 5(c). Here, we focus on the center of the two layers ($\mathrm{\Delta }z=0$) as well as the top microsphere layer ($\mathrm{\Delta }z=1.9\mathrm{\mu m}$) in an attempt to distinguish the two separate layers. At the top of each image (where microspheres in the two layers overlap), it is especially hard to resolve each sphere or determine which sphere is in a particular layer. These challenges are due in part to the limited amount of information contained within the optical intensity at each plane, as opposed to the sample’s complex refractive index.

Next, we return the focus to the $\mathrm{\Delta }z=0$ plane and implement FPT. We display three slices of our 3D scattering potential reconstruction in Figs. 5(d)–5(f). Here, we show the absolute values of the potential near the plane of the top layer, at the center, and near the plane of the bottom layer. The originally indistinguishable spheres within the top and bottom layers are now clearly resolved in each $z$-plane. Due to the system’s limited axial resolution, the reconstruction at the middle plane ($\mathrm{\Delta }z=0$) still shows the presence of spheres from both layers. Comparing Figs. 5(b) and 5(c) with Figs. 5(e) and 5(f), it is clear that the axial resolution of FPT is sharper than manual refocusing. Not only is each sphere layer distinguishable (as predicted theoretically), but we now also have quantitative information about the sample’s complex refractive index.

B. Biological Experiments

For our first biological demonstration, we reconstruct a Trichinella spiralis parasite in 3D (see Fig. 6 and Visualization 1 for the complete tomogram). Since the worm extended along a larger distance than the width of our detector, we performed FPT twice and shifted the FOV between to capture the left and right sides of the worm. We then merged each tomographic reconstruction together with a simple averaging operation (matching that from FP [40], 10% overlap). The total captured volume here is $0.8\mathrm{mm}\times 0.4\mathrm{mm}\times 110\mathrm{\mu m}$. If our setup included a digital detector that occupied the entire microscope FOV, the fixed imaging volume would be $1.1\mathrm{mm}\times 1.1\mathrm{mm}\times 110\mathrm{\mu m}$, and no movement would be needed for this example.

 

Fig. 6. Tomographic reconstruction of a Trichinella spiralis parasite. (a) The worm’s curved trajectory resolved within various $z$ planes. (b) Refocusing the same distance to each respective plane does not clearly distinguish each in-focus worm segment (marked by white arrows). Since the worm is primarily transparent, in-focus worm sections exhibit minimal intensity contrast, presenting significant challenges for segmentation (see intensity along each black dash in inset plots, where black dash location is in-focus in left image). FPT, on the other hand, exhibits maximum contrast at each worm voxel. See Visualization 1.

Download Full Size | PDF

A thresholded 3D reconstruction of the parasite index is at the top of Fig. 6 (real component, threshold applied at $\mathrm{Re}[\mathrm{\Delta }n]>0.7$ after $|\mathrm{\Delta }n|$ normalized to 1, under-sampled for clarity). The maximum real index variation across the tomogram before normalization is approximately 0.06, which can be seen without thresholding in Visualization 1. Its 3D curved trajectory is especially clear in the three separate $z$-slices of the reconstructed tomogram in Fig. 6(a). The two downward bends in the parasite body are lower than the upward bend in the middle, as well as at its front and back ends. It is very challenging to resolve these depth-dependent sample features by simply refocusing a standard microscope. Figure 6(b) displays such an attempt, where the same three $z$ planes are brought into focus manually. Since the sample is primarily transparent, in-focus areas in each standard image actually exhibit minimal contrast, as marked by arrows in Fig. 6(b). We plot the intensity through a fixed worm section (black dash) in each of the three insets. The intensity contrast drops by over a factor of 2 at in-focus locations, which will pose a significant challenge to any depth segmentation technique (e.g., focal stack deconvolution [4]). Since FPT effectively offers 3D phase contrast, points along the parasite within its reconstruction voxels instead show maximum contrast, which enables direct segmentation via thresholding, as shown in the plots in Fig. 6(a).

For our second 3D biological example, we tomographically reconstruct a starfish embryo at its larval stage [see Fig. 7(a) and Visualization 2]. Here, we again show three different closely spaced $z$-slices of the reconstructed scattering potential ($\mathrm{Re}[\mathrm{\Delta }n]$, no thresholding applied). Each $z$-slice contains sample features that are not present in the adjacent $z$-slices. For example, the large oval structure in the upper left of the $\mathrm{\Delta }z=0$ plane, which is a developing stomach, nearly completely disappears in the $\mathrm{\Delta }z=-3.7\mathrm{\mu m}$ plane. Now at this $z$-slice, however, small structures, which we expect to be developing mesenchyme cells [52] and various epithelial cells [53], clearly appear in the lower right. We confirm the presence of these structures with a differential interference contrast (DIC) confocal microscope in Fig. 7(c) (Zeiss LSM 510, 0.8 NA, $\lambda =633\mathrm{nm}$, 0.2 μm scan step). A DIC confocal scan is one of the few possible imaging options for this thick and primarily transparent sample, but suffers from a much smaller FOV (approximately 4% of the total effective FPT FOV). It also does not quantitatively measure the refractive index and requires mechanical scanning. Finally, we attempt to refocus through the embryo using a standard microscope ($\mathrm{NA}=0.4$) in Fig. 7(b). Both the particular plane of the developing stomach and even the presence of the mesenchyme cells are completely missing from the refocused images. This is due to the inability of the standard microscope to segment each particular plane of interest, the inability to accurately reconstruct transparent structures without a phase contrast mechanism, and an inferior lateral resolution with respect to FPT.

 

Fig. 7. 3D reconstruction of a starfish embryo at larval stage. (a) Three different axial planes of the FPT tomogram show significant feature variation (e.g., protocol is completely missing from $\mathrm{\Delta }z=-3.7\mathrm{\mu m}$ plane, expected developing mesenchyme cells are only visible in $\mathrm{\Delta }z=-3.7\mathrm{\mu m}$ plane). (b) Such axial information, and even certain structures (e.g., mesenchyme cells and various epithelia cells, marked in (a)) are completely missing from standard microscope images after manual refocusing. (c) A high-resolution DIC confocal scan of the $\mathrm{\Delta }z=-3.7\mathrm{\mu m}$ plane confirms presence of structures of interest. See Visualization 2.

Download Full Size | PDF

5. CONCLUSIONS

We have performed diffraction tomography using intensity measurements captured with a standard microscope and an LED illuminator. The current system offers a lateral resolution of approximately 400 nm at the Nyquist–Shannon sampling limit (550 nm at the Sparrow limit and 800 nm full period limit) and an axial resolution of 3.7 μm at the sampling limit. The maximum axial extent attempted thus far was 110 μm along $z$, and we demonstrated quantitative measurement of the complex index of refraction through several types of thick specimen with contiguous features.

To improve the experimental setup, an alternative LED array geometry that enables a higher angle of illumination will increase the resolution. Also, we set the number of captured images here to match the data redundancy required by ptychography [48]. However, we have observed that reconstructions are successful with much fewer images than otherwise expected. Along with using a multiplexed illumination strategy [28], this may help speed up the tomogram capture time. In addition, we did not explicitly account for the finite LED spectral bandwidth or attempt polychromatic capture, which can potentially provide additional information about volumetric samples [32,35]. Finally, we set our reconstruction range along the $z$-axis somewhat arbitrarily at 110 μm. We expect the ability to further extend this axial range in the future.

Subsequent work should also examine connections between FPT, multi-slice-based techniques for ptychography [33,34], and machine learning for 3D reconstruction [54]. While prior work already specifies sample conditions under which the first Born [43] and multi-slice [42] approximations remain accurate, it is not yet clear if this translates directly to reconstructions that require phase retrieval. By merging these approaches, it may be possible to increase the domain of sample validity beyond what is currently achieved by each technique independently.

Finally, FPT may also adopt alternative computational tools to help improve ptychographic DT under the first Born approximation. We used the well-known alternating projections phase retrieval update. Other solvers based upon convex optimization [55] or alternative gradient descent techniques [56,57] may perform better in the presence of noise. Alternative approximations besides first Born approximation (e.g., Rytov [15]) are also available to simplify the Born series. In addition, the resolution is currently impacted by the missing cone in 3D k-space, and various methods are available to fill this cone in by assuming the sample is positive only, sparse, or of a finite spatial support [51]. Finally, methods exist to solve for the full Born series by taking into account the effects of multiple scattering [46,47]. Connecting this type of multiple scattering solver to FPT may aid with the reconstruction of increasingly turbid biological samples.