Diffraction tomography with Fourier ptychography

This article presents a method to perform diffraction tomography in a standard microscope that includes an LED array for illumination. After acquiring a sequence of intensity-only images of a thick sample, a ptychography-based reconstruction algorithm solves for its unknown complex index of refraction across three dimensions. The experimental microscope demonstrates a spatial resolution of 0.39 $\mu$m and an axial resolution of 3.7 $\mu$m at the Nyquist-Shannon sampling limit (0.54 $\mu$m and 5.0 $\mu$m at the Sparrow limit, respectively), across a total imaging volume of 2.2 mm $\times$ 2.2 mm $\times$ 110 $\mu$m. Unlike competing methods, the 3D tomograms presented in this article are continuous, quantitative, and formed without the need for interferometry or any moving parts. Wide field-of-view reconstructions of thick biological specimens demonstrate potential applications in pathology and developmental biology.


I. INTRODUCTION
iment, 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. As a synthetic aperture technique, DT comes with the additional benefit of improving the limited resolution of an imaging element beyond its traditional diffraction-limit cutoff [14]. Thus, it appears like a well-suited method to study thick, transparent samples at high resolution.
However, as a technique that models both the amplitude and phase of a coherent field, nearly all prior implementations of DT required a reference beam and holographic measurement, or some sort of phase-stable interference (including SLM coding strategies, e.g. as in [22]). Reference fields require sub-micrometer stability in terms of both motion and phase drift, which has thus far limited DT to well-controlled, customized setups. While several prior works have considered solving the DT problem from intensity-only measurements from a theoretical perspective [23][24][25][26][27], none have implemented a DT system within a standard microscope, or connected their reconstruction attempts to ptychography.
Here, we perform DT based upon intensity images captured under variable LED illumination from an array source. Our technique, termed Fourier ptychographic tomography (FPT), captures 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 (as opposed to physically) rejects light from all areas above and below each plane of interest. FPT also improves the lateral image resolution beyond the standard cutoff of the objective lens used for imaging. The end result is a quantitatively accurate three-dimensional map of the complex index of refraction of a volumetric sample, obtained directly from a sequence of standard microscope images.

II. RELATED WORK
The theoretical foundations of DT were first developed by Wolf [13]. A number of implementations based upon holography have followed. An early demonstration by Lauer is a good example [14]. Prior methods have also implemented tomography within a microscope-like setup, but required the addition of a phase-stable reference beam. Early results applied the projection approximation, which models light as a ray [15]. Subsequent work has taken the effects of diffraction into account [16][17][18].
Instead of relying upon holography, this work measures intensity images and computationally recovers the missing phase of each field. As mentioned above, a few prior works consider the reconstruction problem from detected optical intensities, but must either move the focal plane between measurements, or must assume a sample support constraint. They do not attempt ptychographic phase retrieval. One related strategy worth mentioning is lifting-based phase retrieval for tomography [19]. The connection between the first Born approximation and phase retrieval has also been explored within the context of volume hologram design [20].
Related efforts to reconstruct volumetric samples from wide-field intensity-only measurements outside of the realm of DT include lensless on-chip setups [28,29], lensless techniques that assume an appropriate linearization [30], and methods relying upon effects like defocus (e.g., the transport of intensity equation [31]) or spectral variations [21]. None of these techniques fit within a standard microscope setup, nor offer the ability to simultaneously improve spatial resolution.
Based upon a standard microscope, Fourier ptychography (FP) [32] can simultaneously improve image resolution and measure quantitative phase [33]. However, it is restricted to thin samples.
FPT effectively extends prior developments of FP into the third dimension. Two recent works also examined the problem of 3D imaging from intensities in a standard microscope [34,35]. These two examples adopted their 3D reconstruction technique from the related field of 3D ptychography [36,37], where the sample under examination is split up into a specified number of infinitesimally thin slices, and the beam propagation method (i.e., multi-slice approximation) is applied [38]. Unlike the multi-slice approach, which works well with distinctly separated absorbing layers, FPT is best suited for continuous, primarily transparent samples. A number of related methods to perform 3D X-ray ptychography have also been proposed [39][40][41]. However, none 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 2D ptychographic solvers to determine the complex field for individual projections of a slowly rotated sample, which are subsequently combined using conventional DT techniques, as shown with both crystallographic [42] and unordered specimens [43].
Here, we first outline a solid foundation for the application of ptychographic phase retrieval to DT. Unlike approaching the problem from a projection-based or multi-slice perspective, the framework of DT (under the first Born approximation) follows directly from the scalar wave equation.
It offers a clear picture of achievable resolution in 3D, spells out sampling and data redundancy requirements for an accurate reconstruction, and presents a clear path forward for future extensions to account for multiple scattering [44]. Furthermore, our method does not require the arbitrary assignment of the number slices in the 3D volume, or their location, or for us to select a particular order in which to address each slice as iterations proceed. Instead, it simply inserts the measured data into its appropriate location in 3D Fourier space and ensures phase consistency between each measurement, given a sufficient amount of data redundancy (just like ptychography). From the initial starting point of solving for the first term in the Born expansion, we aim this approach as a general framework to eventually solve the challenging problem of forming tomographic maps of volumetric samples, at sub-micrometer resolution, in the presence of significant scattering.

III. METHODS
In this section, we develop a mathematical expression for our image measurements using the FPT framework, and then summarize our reconstruction algorithm. We describe our setup and reconstruction in 3D with the vector r = (r x , r y , r z ) defining the sample coordinates and the vector k = (k x , k y , k z ) defining the k-space (wavevector) coordinates (see Fig. 1).

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(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π/λ is the wavenumber in vaccuum. It is informative to point out that n(r) = n r (r) + in i (r), where n r is associated with the sample's refractive index and n i is associated with its absorptivity. We typically neglect the dependence of n on λ since we illuminate with quasi-monochromatic light. This dependence cannot be neglected when imaging with polychromatic light.
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 (r), as a sum of two fields: U (r) = U i (r) + U s (r). Here, U i (r) is the field incident upon the sample (i.e., from one LED) and U s (r) is the resulting field that scatters off of the sample. As detailed in [13], we may insert this decomposition into the scalar wave equation for light propagating through an inhomogeneous medium. Using Green's theorem, we may determine the scattered field as, Here, G(|r − r|) is the Green's function connecting light scattered from various sample locations, denoted by r, to an arbitrary location r . V (r) is the scattering potential from Eq. 1. Since U (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 (r) in the integrand with U i (r). This approximation assumes that U i (r) U s (r). It is the first term in the Born expansion that describes the scattering response of an arbitrary sample [13].
Our system sequentially illuminates the sample with an LED array, which contains q = q x × 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 ≤ j x ≤ q x /2 and −q y /2 ≤ j y ≤ q y /2, as well as a single counter variable j, where 1 ≤ j ≤ q. Assuming each LED acts as a spatially coherent and quasi-monochromatic source (central wavelength λ) 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 θ jx = tan −1 (j x · c/l) and θ jy = tan −1 (j y · c/l) with respect to the x and y axes, respectively. We may express the jth field incident upon the sample as, where k j = (k jx , k jy , k jz ) = k · sin θ jx , sin θ jy , 1 − sin 2 θ jx − sin 2 θ jy is the wavevector of the jth LED plane wave. As θ jx and θ jy vary, 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 (r) in Eq. 2 with U (j) i (r) in 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 jth scattered field,Û (j) s (k), in the far field [13]: Here,V (k), which we refer to as the k-space scattering potential, is the three-dimensional Fourier transform of V (r), k denotes the scattered wavevector in the far field, and we have neglected constant multiplication factors for simplicity. The field scattered by the sample and viewed at a large distance,Û (j) s (k), is given by the values along a specific manifold (or spherical "shell") of the k-space scattering potential, here written asV (k − k j ). We illustrate the geometric connection s (k) for a 2D optical geometry in Fig. 2(b). The center of the jth shell is defined by the incident wavevector, k j . For a given shell center, each value of interest lies on a spherical surface at a radial distance set by |k| = k (see colored arcs in Fig. 2(b)). As k j varies with the changing LED illumination, the shell center shifts along a second shell with the same radius (since k j is itself constrained to lie on an Ewald sphere, see gray circle in Fig. 2(b)).
The goal of DT is to determine all complex values within the volumeV , from a set of q scattered fields, {Û s } q j=1 . It is common to measure these scattered fields holographically [14,18]. Each 2D holographic measurement maps to the complex values ofV that lie along one 2D shell. Values from multiple measurements (i.e., the multiple shells in Fig. 2(b)) can be combined to form a k-space scattering potential estimate,V e [45]. 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. Fig. 1(c) and Fig. 2(c) each display a typical measurable volume, also termed a bandpass, from a limited-angle illumination and detection setup. Once sampled, an inverse 3D Fourier transform of the band-limitedV e (k) yields the desired complex scattering potential estimate, V e (r), from which the quantitative index of refraction is directly obtained.
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. 4 are defined at the microscope objective back focal plane (i.e., its Fourier plane), whose 2D coordinates k 2D = (k x , k y ) are 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 jth shifted field at our microscope back focal plane as,Û (j) (k 2D ) =V (k 2D − k j2D , k z − k jz ). These new coordinates highlight the 3D to 2D mapping fromV toÛ , where again k z = k − k 2 x − k 2 y is a deterministic function of k 2D , and the same applies between k jz and k j2D .
Each shifted scattered field is then bandlimited by the microscope aperture function, a(k 2D ), before propagating to the image plane. The limited extent of a(k 2D ) (defined by the imaging system NA) sets the maximum extent of each shell along k x and k y . The jth 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 k 2D 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 functionV from the real, positive data matrix, g(x, y, j). A final 3D Fourier transform ofV yields the desired scattering potential, and subsequently the refractive index distribution, of the thick sample.

B. FPT reconstruction algorithm
Eq. 5 closely resembles the data matrix measured by Fourier ptychography (FP, see [32]). Now, however, intensities are sampled from a volumetric function along shells in a 3D space (i.e., the curves in Fig. 2). We use an iterative reconstruction procedure, mirroring that from FP [32], to "fill in" the k-space scattering potential with data from each recorded intensity image. Just like FP requires a certain amount of data redundancy (i.e., overlap in k-space) to accurately recover the unknown optical phase, FPT also requires overlap between shell regions in 3D k-space. Since our discretized k-space now has an extra dimension, overlap is less frequent and more images are required for successful algorithm convergence. A coarser k-space discretization, a smaller LED array pitch and/or a larger array-sample distance along z will help increase overlap. As we demonstrate experimentally, several hundred images are sufficient for a complex reconstruction that contains approximately 30 unique axial slices.
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 wavevector along k x and k y is proportional to k(NA o +NA i ), where NA o is the objective NA and NA i is maximum NA of LED illumination. This lateral spatial resolution limit matches FP [46]. The maximum resolvable wavevector range along k z is also determined as a function of the objective and illumination NA as, k max This relationship is easily derived from the geometry of the k-space bandpass volume in Fig. 2, as shown in [14]. We typically specify the maximum imaging range along the axial dimension, z max , to equal approximately twice the expected sample thickness. This then sets the discretization level along k z : ∆k z = 2π/z max . The total number of resolved slices along z is set by the ratio k max z /∆k z .
We now summarize the FPT reconstruction algorithm in the following 5 steps: 1. Initialize a discrete estimate of the unknown k-space scattering potential,V e (k), selecting the appropriate 3D array size following the discussion above. Either a single raw image may be padded along all three dimensions and then Fourier transformed for this initialization, or the raw intensity image set may be used to form a refocused light field [34]. A constant matrix is also often an adequate initialization.
2. For j = 1 to q images, compute the center coordinate, k j , and select its associated shell (radius k, maximum width 2k · NA o ). This selection process samples a discrete 2D function, d j (k x , k y ), from the 3D k-space volume. The selected voxels must partially overlap with voxels from adjacent shells. Currently, no interpolation is used to map voxels from the discrete shell to pixels withind j (k x , k y ).
3. Fourier transformd j (k x , k y ) to the image plane to create d j (x, y), and constrain its amplitudes to match the measured amplitudes from the jth image. For example, the update may take the simple form, d j (x, y) = g(x, y, j) · d j (x, y)/|d j (x, y)|. More advanced alternating projection-based updates are also available [47]. Continue for a fixed number of iterations, or until satisfying some error metric. At the end, 3D inverse Fourier transformV e (k) to recover the complex scattering potential, V e (r).
In practice, we also implement a pupil function recovery procedure [48] as we update each extracted shell from k-space. This allows us to simultaneously estimate and remove possible aberrations present in the microscope back focal plane.

IV. RESULTS
We experimentally verify our reconstruction technique using a standard microscope outfitted  For most of the reconstructions presented below, we capture and process q = 675 images from the same fixed pattern of LEDs. We typically use the following parameters for FPT reconstruction: each raw image is cropped to 1000 × 1000 pixels, the reconstruction voxel size is set at 0.39 × 0.39 × 3.7 µm 3 for sampling at the Nyquist-Shannon rate, the reconstruction array contains approximately 2100 × 2100 × 30 voxels, and the algorithm runs for 5 iterations.

A. Quantitative verification
We include three different quantitative verifications of FPT performance using polystyrene microspheres as reference targets. First, we verify the ability of FPT to improve lateral image resolution. This matches the goal of FP for thin 2D samples. 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 (a) (generated from the center LED) cannot resolve the individual spheres gathered in small clusters.  resolve points that are closer than 1.1 µm. After FPT reconstruction, we obtain the complex index of refraction in Fig. 3(b), where here we show the real component of the recovered index. The FPT reconstruction along the ∆z = 0 slice clearly resolves the spheres within each cluster. This 800 nm distance is close to the expected Sparrow limit for FPT of 0.68λ/ (NA o + NA i ) = 540 nm.
The ringing features around each sphere indicate a sinc-like point-spread function, as we expect theoretically, and this ringing constructively interferes to form the undesired dip feature at the center of each cluster.
Second, we check the quantitative accuracy of FPT by imaging microspheres that extend across more than just a few reconstruction voxels. Fig. 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, and display a cropped section (200×200×15 voxels) of the full reconstruction. We again display the real (non-absorptive) component of the recovered index across both a lateral slice (along the ∆z = 0 plane) and a vertical slice (along the ∆y = 25 µm plane). We also include detailed 1D traces along the center of the vertical slice.
Three observations are noteworthy regarding this experiment. First, the measured index shift approximately matches the expected shift of ∆n = n s − n o = 0.01 across the entire bead, thus demonstrating quantitatively accurate performance. Second, for each 1D trace through the center of each microsphere, we would ideally expect a perfect rect function (from ∆n = 0 to ∆n = 0.01 and then back down). This is unlike 2D FP, which reconstructs the phase delay though each sphere, leading to a parabolic function (due to their varying thickness along the optical axis).
While the system can resolve an approximate step function through the center of the sphere along the lateral (x) dimension, it is not a step function function along the axial (z) direction. This is caused by the limited extent of the measurable volume of 3D k-space (i.e., the limited bandpass).
The "missing cone" of information, primarily surrounding the k z axis, creates a noticeably wide point-spread function along z, which leads to its distinct sinc shape. Various methods are available to computationally fill in the missing cone using prior sample information [49,50].
For our third observation, 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 attempt to holographically refocus its complex 2D solution. We obtain this solution using the same number of images (q = 675) and with the FP procedure in [32], 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, quite 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). FPT thus appears more accurate here than these two alternatives.
For our third and final quantitative test, we verify the axial resolution of FPT along z. Here, we prepare a sample containing two closely separated layers. Each layer contains 2 µm microspheres (n s = 1.59) distributed across the surface of a glass slide, which we sandwich together with oil in between (n o = 1.515). The 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, as diagrammed in Fig. 5(a)). The 3.9 µm center-to-center distance is close to the expected axial resolution limit of 3.7 µm for the FPT microscope.
Conventional microscope images of the sample, using the center LED for illumination, are in Not only is each layer clearly distinguishable (as predicted theoretically), but we now also have quantitative information about the sample's complex refractive index.

B. Biological experiments
Next, we use FPT to reconstruct the 3D complex refractive index of two different thick biological specimens. Since the exact composition of these specimens is unknown, it is challenging to verify the quantitative accuracy of our reconstructions here, especially given the accuracy of first Born approximation is only guaranteed up to a total phase shift of approximately one wavelength.
However, we will point out the multiple qualitative benefits of FPT in these examples.
Our first tomographic reconstruction is of a trichinella spiralis parasite (Fig. 6). Here, since the worm extends along a larger distance than the width of our detector, we performed FPT twice, shifting the FOV between to capture the left and right side of the worm with 10% overlap between. We then merged each tomographic reconstruction together with a simple averaging operation (matching that from FP [32]). The total captured volume here is 0.8 mm×0.4 mm×110 µm. If we were to replace our current digital detector with one that occupied the entire microscope FOV, we would increase our fixed imaging volume to 2.2 mm×2.2 mm×110 µm, and obtain tomograms that each contain approximately 10 9 voxels per acquisition.
We display a thresholded 3D scattering potential reconstruction of the parasite at the top of  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. 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 one worm section (black dash) in three insets. The intensity contrast drops by over a factor of 2 at infocus locations, which will prevent the success of depth segmentation techniques (e.g., focal stack deconvolution [4]). Since FPT effectively offers phase contrast, points along the parasite within its reconstruction voxels instead show maximum contrast, which enables direct segmentation via thresholding, as achieved in the top plot.
For our second 3D biological example, we tomographically reconstruct a starfish embryo at its larval stage (Fig. 7). Here, we again show three different closely spaced z-slices of the reconstructed scattering potential (Re[∆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 ∆z = 0 plane, which is a developing stomach, nearly completely disappears in the ∆z = −3.3 µm plane. Now at this z-slice, however, small structures, which we expect to be the developing vibratile celia as well as epithelia cells at various locations [51], clearly appear in the lower right.
Both the particular plane of the developing stomach and even the presence of the vibratile celia 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 (x, y) resolution with respect to FPT.

V. CONCLUSION
The FPT method performs diffraction tomography using intensity measurements, captured with a standard microscope and an LED illuminator. Its reconstruction algorithm extends previous work with FP to now operate in 3D. 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, which leads to approximately one giga-voxel of complex sampling points per acquisition if imaging over the total microscope FOV (2.2 × 2.2 mm). We demonstrated quantitative measurement of the complex index of refraction within thick biological specimens.
We believe that FPT can be significantly improved with additional experimental development.
First, a better LED array geometry will enable a higher angle of illumination to improve resolution.
Second, we set the number of captured images here to match previously determined data redundancy requirements [52]. However, we have observed that reconstructions are successful with much fewer images than otherwise expected. Along with using a multiplexed illumination strategy [53], this may help significantly speed up tomogram capture time. Third, 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 future setups.
Alternative computational approaches may also improve FPT. Here, we list a number of possible directions. First, we adopted the well-known alternating projections method (i.e., the ePIE algorithm [54]) for ptychographic update. Other methods, such as convex-based approaches [55], can perform better in the presence of noise. Second, alternative approximations are also available to solve the first Born approximation [56]. Third, a big detriment to resolution is currently the missing cone in 3D k-space, and various methods are available to fill this cone in, e.g., by assuming the sample is positive-only, sparse, or of a finite spatial support [49,50]. Finally, there are already suggested methods to solve for the full Born series, which take into account the effects of multiple scattering [44]. Connections between this type of multiple scattering solver, recent methods applying the multi-slice approximation [34,35,57], and FPT may lead to successful reconstruction of increasingly turbid samples.