3D intensity and phase imaging from light field measurements in an LED array microscope

Realizing high resolution across large volumes is challenging for 3D imaging techniques with high-speed acquisition. Here, we describe a new method for 3D intensity and phase recovery from 4D light field measurements, achieving enhanced resolution via Fourier ptychography. Starting from geometric optics light field refocusing, we incorporate phase retrieval and correct diffraction artifacts. Further, we incorporate dark-field images to achieve lateral resolution beyond the diffraction limit of the objective ( 5 × larger NA) and axial resolution better than the depth of field, using a low-magnification objective with a large field of view. Our iterative reconstruction algorithm uses a multislice coherent model to estimate the 3D complex transmittance function of the sample at multiple depths, without any weak or single-scattering approximations. Data are captured by an LED array microscope with computational illumination, which enables rapid scanning of angles for fast acquisition. We demonstrate the method with thick biological samples in a modified commercial microscope, indicating the technique ’ s versatility for a wide range of applications.


INTRODUCTION
3D imaging techniques are crucial for thick samples and in situ studies; however, high-speed acquisition with high resolution remains a challenging task. Confocal microscopy [1] and multiphoton microscopy [2] are extremely popular for their high resolution, but require slow point scanning through the 3D volume. Wide-field 3D imaging generally involves tomography, in which projections of the sample are taken at many angles. Here, we consider phase imaging, which provides stain-free and labelfree intrinsic contrast for transparent biological samples.
Here, we describe an alternative method in which only a single intensity image is captured for each angle. This is possible because data with angular diversity provide both 3D information and phase contrast. 2D phase of thin samples can be computed from images taken at multiple illumination angles [20][21][22][23][24] because of asymmetry introduced in the pupil plane [25]. All of these approaches assume a thin sample; in thick samples, angle-dependent data usually represent tomographic information. Here, instead of choosing either 2D phase imaging of thin samples or 3D recovery of thick samples, we achieve both.
Further, we use angles of illumination greater than that allowed by the numerical aperture (NA) of the objective, so the resulting dark-field images contain subresolution feature information. Using Fourier ptychography [23,24], we use these to build up a larger effective NA, limited by the sum of the illumination and objective NAs. Thus, a low-magnification objective having a large field of view (FoV) can recover highresolution gigavoxel-sized 3D intensity and phase images.
We take a holistic approach to the inverse problem, where an optimization procedure is used to recover 3D intensity and phase from the captured data, and remove aberrations. The algorithm is inspired by 3D ptychography [26,27], in which a multislice approach models the sample as a series of thin slices [28]. The wave field propagates through the sample from slice to slice, with each slice modulating the field, and the objective NA limiting the resolution of the captured images. To solve the inverse problem, we iteratively update the 3D complex transmittance function for each illumination angle, effectively implementing a nonlinear 3D deconvolution [29] that removes out-of-plane blur. Since the multislice model makes no weak or single-scattering approximations, it can correct for multiple scattering.
Our work is best understood by its relation to light fields, which use space and angle to parameterize rays in 3D. Light field microscopes [30] capture projections across a 2D range of angles, whereas tomography generally only scans angles in 1D. Thus, light fields describe 4D data that can be thought of as limited-angle tomography with multiple directions of rotation [31,32]. The standard algorithm for light field digital refocusing fully incorporates 3D geometric effects [33,34]. However, when wave-optical effects become more prominent (e.g., at smaller feature sizes), the lateral resolution in the digitally refocused images deteriorates due to unaccounted for diffraction [35]. Using the light field refocused result as an initial guess, our algorithm can be thought of as a diffraction correction routine for light fields, with embedded phase retrieval.
To collect 4D space-angle data, we use an LED array microscope, which scans through illumination angles quickly and with no moving parts. Similar LED array illumination was demonstrated previously for on-chip lensless imaging techniques [36,37]. Our system is built on a commercial microscope in which the illumination unit has been replaced by a programmable LED array. This simple, inexpensive hardware modification enables not only 4D light field capture [35,38], but also dark field [38,39], phase contrast [35,39], Fourier ptychography [23,24], and digital aberration removal [40].

A. Light Field Refocusing from LED Array Measurements
By placing the LED array sufficiently far above the sample that the illumination is considered spatially coherent, we can treat every LED's illumination as a plane wave from a unique angle. Sequentially turning on each LED in the 2D array, while capturing images, therefore builds up a 4D data set of two spatial and two angular variables, similar to a light field measurement.
For 3D samples, light field refocusing is intuitively understood as a compensation for the geometric shift that occurs upon propagation. Figure 1 illustrates how off-axis illumination causes the intensity to shift from its original position in the plane of focus by a distance Δx x i Δz∕z i . With higher angles or larger Δz, the rays shift further across the plane, as shown in Fig. 2(b). The slope of the line created by each feature is determined by its depth, while the x location is defined by its θ x 0 crossing. The light field refocusing routine undoes the shift and sums over all angles to synthesize the intensity image.
However, diffraction and phase effects can cause the light to deviate from the straight lines predicted by geometrical optics. This is evidenced by the diffraction rings surrounding each dark line in Fig. 2(b). While light field refocusing corrects for geometric shifts, additional wave-optical effects degrade the resolution with defocus. Our algorithm starts from the light field refocused result, which captures most of the energy redistribution. We then iteratively estimate the phase and diffraction effects.
To achieve resolution beyond the diffraction limit of the objective, we need dark-field illumination from LEDs at high angles. For thin samples, each illumination angle shifts the sample spectrum around in Fourier space, with the objective aperture selecting out different sections. Thus, by scanning through different angles, many sections of Fourier space are captured. These can be stitched together with synthetic aperture approaches [41,42] to create a high-resolution image in real space. The caveat is that phase is required, which the Fourier ptychography algorithm [23,24] provides by performing translational diversity phase retrieval [43][44][45][46] in Fourier space.
When the sample is thick, each angle of illumination takes a different path through the sample. Thus, the Fourier spectrum of each illumination angle's exit field is different, but all of these data are interrelated by the multislice model that we use here. Combining many angle-dependent low-resolution images can therefore still achieve enhanced resolution at all slices, limited by the sum of the illumination and objective NAs.

B. Multislice Forward Model
Our forward model assumes that the illumination from the nth LED is a tilted plane wave f n 1 r expi2πu n · r, where spatial frequency is related to illumination angle by u n sin θ x;n ∕λ; sin θ y;n ∕λ and λ is wavelength.
The field propagating through the thick sample is modeled by a multislice approximation [26,28] that splits the 3D sample into a series of N thin slices, each having a complex transmittance function o m r (m 1; 2; …; N ), where r x; y denotes the lateral coordinates and m indexes the slices. As Research Article light passes through each slice, the field is first multiplied by the 2D transmittance function of that slice, and then propagated to the next slice. The spacing between neighboring slices is modeled as a uniform medium (e.g., air) of thickness Δz m . Thus, the field exiting the sample can be calculated using a series of multiply-and-propagate operations: where the superscripts and subscripts denote the indices of the LED and the slice, respectively. f and g are the fields before and after each slice, respectively, and P Δz f·g F −1 exp−i2πΔz ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1∕λ 2 − juj 2 p F f·g denotes propagation by a distance Δz, with F f·g and F −1 f·g being the 2D Fourier transform and its inverse, respectively. After passing through the thick sample, the exiting complex field is then imaged to the camera plane of the microscope, a process that involves a low-pass filtering by the objective's pupil function. The Fourier spectrum of the field at the camera plane, C n u, is thus where G n u denotes the spectrum of the exit field from the last slice, Pu the pupil function including aberrations, and is an additional defocus term assuming the last slice is Δz N distance away from the actual focal plane. The final intensity image from the nth LED illumination is I n r jF fC n grj 2 :

C. Reconstruction Algorithm
Using the multislice forward model, we develop an iterative reconstruction routine that makes explicit use of the light field result as an initial guess. Light field refocusing predicts the intensity image I Δz at Δz from the actual focal plane to be [35,38] I Δz x; y X where the coordinates of the nth LED are defined by x n ; y n . Our initial guess of the 2D transmittance function at the corresponding sample slice is then o Δz x; y ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi I Δz x; y p . We then improve the estimate for the sample's intensity and phase at each slice, as well as an estimate of the pupil function aberrations [40,[44][45][46], by an iterative Fourier ptychography (c) Reconstructions using our multislice method are compared to light field refocusing and physically changing the microscope focus. Diffraction effects severely blur the light field results, whereas our multislice method is able to recover the full diffraction-limited resolution with improved image contrast, while also removing out-of-focus blur.
Research Article [23,24] reconstruction process, combined with the multislice inversion procedure in [26]. The reconstruction procedure aims to minimize the difference between the actual and estimated intensity measurements in a least-square sense: min fo m rg;Pu X n X r jI n r − jF fC n grj 2 j 2 : At each iteration, the sample's transmittance function is updated for each illumination angle as follows: (1) Starting from the current guess of the multislice transmittance function, we use our forward model to generate the current estimate of the Fourier spectrum of the field at the camera plane, C n u, when illuminating with the nth LED.
(2) We update C n u by replacing the estimated intensity with the actual measurement [47,48]: F fC n g jF fC n gj ; where· denotes an updated variable.
(3) We update the exit field's spectrum and pupil function: G n u UG n ; P; C n ∕H ;Ĉ n ∕H ; Pu UP; G n ; C n ∕H ;Ĉ n ∕H : Since the two functions come as a product, we use a gradient descent procedure, described by U, to separate the two updates [24,44]. The procedure is general for updating anyψ from the product of its previous estimate ψ with another function ϕ, for example, β ψϕ: ψ Uψ; ϕ; β;β ψ jϕjϕ β − β jϕj max · jϕj 2 δ ; where δ is a regularization constant to ensure numerical stability. The updated exit fieldĝ n N is then g n N r F −1 fĜ n gr: (4) The field is back-propagated through the 3D sample and the following steps are repeated until the first slice. At the mth slice, the transmittance function o n m and the incident field f n m of this slice are updated using the same procedure as Eq. (11): o n m r Uo n m r; f n m r; g n m r;ĝ n m r; f n m r Uf n m r; o n m r; g n m r;ĝ n m r: The updated exit field of the previous sliceĝ n m−1 is related tô f n m by back-propagation: g n m−1 r P −Δz m−1 ff n m rg: (5) At the first slice, the incident field is kept unchanged as the original illumination:f n 1 r f n 1 r.
After looping through the images from each of the LEDs, we check the convergence of the current sample estimate by computing the mean squared difference between the measured and estimated intensity images from each angle. The algorithm converges reliably within only a few iterations in all cases we tested.
Although iterative methods like the one used here often get stuck in local minima [48], we find that the close initial guess provided by the light field result helps to avoid this problem. Our observation agrees well with recent phase retrieval theory for similar iterative methods [49,50], which guarantees convergence to a global solution, given proper initialization. Further, the data contain significant redundancy (4D data for 3D reconstruction) and diversity (from angular variation), providing a highly constrained solution space. Similar robustness has also been shown with translational diversity data in real-space ptychography [44][45][46]. Thus, when our estimate correctly predicts the captured data and returns a good convergence criterion, we can be confident that the result is correct.

EXPERIMENTAL RESULTS
Our experimental setup (Fig. 1) consists of a 32 × 32 custommade LED array (4 mm spacing, central wavelength λ 643 nm, with a 20 nm bandwidth) placed z LED 74 mm above the sample, replacing the microscope's standard illumination unit (Nikon TE300 inverted). The LED array is controlled by an ARM microcontroller and is synchronized with the camera (PCO.edge) to scan through the LEDs at cameralimited speeds. Our camera is capable of 100 frames per second at full frame (2560 × 2160 pixels) and with 16-bit data, although we use longer exposure times for dark-field images. Thus, the acquisition time can be easily traded off for image quality or resolution. Each LED has a d 120 μm × 120 μm square emitting area, resulting in an illumination coherence area of A c 400 μm × 400 μm (according to the van Cittert-Zernike theorem, A c λz LED ∕d ). Thus, our coherent plane wave illumination assumption holds as long as we reconstruct the image in patches that have an area smaller than the coherence area. The final full FoV reconstruction is obtained by stitching together all the patches.

A. Improved Light Field Refocusing
We first demonstrate improvement over light field refocusing with a 10× objective (0.25 NA), using only bright-field LEDs. This corresponds to a 69 LED circle at the center of our array. The two-slice test sample consists of two resolution targets, one placed above the focal plane and the other placed below the focal plane and rotated relative to the first one. Some raw images are shown in Fig. 2(a) for LEDs illuminating with a varying angle θ x . The image shifts with illumination angle as shown in the space-angle plots in Fig. 2(b). In Fig. 2(c), we compare reconstructions from light field refocusing and our multislice method with images captured from physical focus (with all bright-field LEDs on). The resolution in the physically refocused images is 0.78 μm (Group 9, Element 3); however, the light field refocused image cannot resolve such small features, due to diffraction blurring. Multislice reconstructions

Research Article
with a two-slice model, however, do recover diffraction-limited resolution at both depths, providing a significant improvement over the light field refocused result.
Since our result represents the sample transmittance function, out-of-plane blur from the other resolution target is largely removed, unlike in physical focusing. In addition, our results have better high-frequency contrast. This is because the physically focused images are taken with all bright-field LEDs on, so the incoherent optical transfer function implies a 2× larger frequency cutoff than the coherent case, but with decreased response at higher spatial frequencies [51]. In our result, we synthesize a coherent transfer function (CTF) that has a uniform frequency response within the passband.

B. Multislice Fourier Ptychography
Next, we demonstrate multislice Fourier ptychography for obtaining resolution beyond the objective's diffraction limit by including dark-field LEDs up to 0.41 illumination NA. To do this, we switch to a 4× objective (0.1 NA), then use our method to recover lateral resolution with an effective NA of 0.51 (Fig. 3). We can resolve Group 9, Element 4, giving a resolution of 0.69 μm (five times better NA than the objective), as expected. The added benefit is that the FoV (1.8 mm × 2.1 mm) of the 4× objective is bigger, resulting in a large volume reconstruction. Note that the physically focused images now display significant out-of-plane blur, since the small NA provides a large depth of field. Our multislice reconstruction successfully mitigates most of this blur, resulting in a clean image at each depth.
Finally, we demonstrate our method on a continuous thick Spirogyra algae sample (Carolina Biological) having both absorption and phase effects (Fig. 4). Here, we use the 10× objective (0.25 NA) with LEDs that provide a best possible lateral resolution of 0.59 μm (Rayleigh criteria with an effective NA of 0.66). The sample has a total thickness of ∼100 μm, which we split into 11 slices spaced by 10 μm, representing a step size midway between the axial resolution of the objective and our predicted axial resolution [Eq. (17)]. Although the sample is continuous through the entire depth range, our multislice method will recover slices that only contain parts of the sample within the axial resolution range around each corresponding depth.

C. Analysis of Resolution
The stacked resolution targets provide a convenient way to experimentally characterize lateral resolution at multiple depths. Figure 5 plots simulated (theoretical) versus experimentally measured resolution for the two-slice situation using several methods and varying defocus depths, where the defocus distance refers to the relative distance of the test target from the physical focus plane of the microscope. We define resolution according to the closest set of bars that can be discriminated.
First, we examine the light field refocusing result. As the defocus distance z increases, the lateral resolution degrades due to diffraction effects, as predicted by theory [35]. Note that this is a different type of diffraction effect than that pointed out for the original light field microscope [30,52]. Our multislice method recovers resolution back to the full diffraction limit of the system. This provides considerable improvement in resolution over the light field refocusing case as the defocus distances get larger.
Using our multislice Fourier ptychography method, we expect to achieve lateral resolution at all slices that is limited by the sum of the NAs of the illumination and objective. When the target is at or near focus, we successfully achieve the maximum resolution expected of this system (0.59 μm). However, as the sample plane moves away from focus the resolution degrades, although this is not predicted by simulations (see Fig. 5). We believe this error to be due to LED position Fig. 3. Experimental results using multislice Fourier ptychography to achieve enhanced resolution at two depths simultaneously. Using a 4× objective, we achieve a resolution of 0.69 μm (five times better NA than the objective). (Top) Low-resolution raw image. (Bottom) Zoom-in at two depths comparing our multi-slice recovery to physical refocusing. miscalibration, since a higher defocus is more sensitive to angle error. Thus, accurate calibration of LED positions will be crucial for extending our work to higher magnifications.
To analyze both lateral resolution and depth (z) sectioning theoretically, we use the 3D CTF of the imaging system [53][54][55][56][57], with the caveat that this theory assumes a single-or weak scattering approximation (e.g., Born or Rytov model) [58]. Analytical theory for multiple scattering is not available, but many procedures start from single scattering and apply it recursively [58], so this should provide a good starting point for resolution analysis. The 3D CTF of our LED array microscope is sketched in Fig. 6, where the thick arcs describe the frequency coverage in the u x − u z space. As expected, the lateral bandwidth Δu x is determined by the sum of the objective NA (NA obj ) and illumination NA (NA illum ): Thus, the lateral resolution when using Fourier ptychography is equivalent to that which would be achieved by an objective having the sum of the two NAs. The axial resolution, however, does not follow this trend. Its bandwidth, Δu z , is neither that of the objective, nor is it that would result from using an objective having the sum of the two NAs. Instead, it is somewhere in between (detailed analysis in Supplement 1), and can be calculated from Fig. 6 as This equation shows that the depth resolution improvement achieved in our experiments will not be as large as the lateral resolution improvement factor. For the 0.25 NA objective and  17) is δz 1∕Δu z 5.4 μm. This estimation sets a lower bound for an achievable axial resolution, since multiple scattering may reduce it further, as observed in simulations (see Supplement 1). In practice, we may suffer further loss at large defocus distances, due to LED miscalibration.

DISCUSSION
The multislice Fourier ptychography method we present here is similar in concept to 2D Fourier ptychography [23] in the sense that it recovers resolution beyond the band limit of the objective. However, the extension to 3D becomes more similar to light fields or tomography. Although our inverse algorithm is quite similar in procedure to real-space 3D ptychography [26,27], the interpretation of the data is different. In our case, we collect real-space images for different angles of illumination, and thus different projection paths through the sample. This leads to the light field refocusing being a close initial guess. The two situations can be related in phase space by a 90 degree rotation [59]. Ptychography collects the same data as a spectrogram [60,61]-Fourier space intensity with real-space aperture scanning [62], whereas Fourier ptychography collects real-space intensity with Fourier space scanning. Thus, the light field refocusing initial guess that we use here could in principle be applied to real-space 3D ptychography. Interestingly, the connection to phase space also predicts the connection between ray and wave optics [63], where the shearing of a light field is analogous to the shearing of a Wigner phase-space function.
One of the key factors of success of our method is the large amount of data collected, because it provides combined data redundancy and diversity, which will improve the convergence of the algorithm. While it is difficult to theoretically quantify the number of images necessary, we find empirically that we always need to collect significantly more pixels of data than we reconstruct. For example, in the experiment in Fig. 4, we reconstruct 11 slices of 2D intensity and phase, plus a single complex pupil function for digital aberration removal. Our data set has 225 captured images, giving a factor of 10× more data collected than recovered. These numbers are comparable to the ratios typically used in 2D ptychography to achieve reliable convergence [64], although multiplexing has been shown to significantly reduce the captured data [24]. Future work will explore the limits of data requirements for the 3D case.

CONCLUSION
We have presented a new method for multislice 3D Fourier ptychography that recovers 3D sample intensity and phase with resolution beyond the diffraction limit of the microscope objective used. Our data are captured by an LED array microscope, which is particularly attractive for commercial microscopy, since it can achieve rapid scanning of angles by LED array illumination. The method is label-free and stain-free, and so has wide application in the biological imaging of live samples.