High-throughput intensity diffraction tomography with a computational microscope

We demonstrate a motion-free intensity diffraction tomography technique that enables direct inversion of 3D phase and absorption from intensity-only measurements for weakly scattering samples. We derive a novel linear forward model, featuring slice-wise phase and absorption transfer functions using angled illumination. This new framework facilitates flexible and efficient data acquisition, enabling arbitrary sampling of the illumination angles. The reconstruction algorithm performs 3D synthetic aperture using a robust, computation and memory efficient slice-wise deconvolution to achieve resolution up to the incoherent limit. We demonstrate our technique with thick biological samples having both sparse 3D structures and dense cell clusters. We further investigate the limitation of our technique when imaging strongly scattering samples. Imaging performance and the influence of multiple scattering is evaluated using a 3D sample consisting of stacked phase and absorption resolution targets. This computational microscopy system is directly built on a standard commercial microscope with a simple LED array source add-on, and promises broad applications by leveraging the ubiquitous microscopy platforms with minimal hardware modifications.


Introduction
Quantitative characterization of thick biological samples is a challenging task. Unstained biological cells appear transparent when imaged under a standard brightfield microscope, since only the sample's absorption information is directly visible. Though techniques based on exogenous labels (e.g. dyes and fluorophores) have been developed [1], they suffer from the need for staining or labeling using external contrast agents which may alter cellular behavior [2,3]. Here, we develop a new label-free phase tomography technique, which provides 3D cellular information with intrinsic structural sensitivity. Our technique is fast, motion-free, and easy to implement with a computational microscope platform, in which a standard commercial microscope is modified with an LED array source [4][5][6][7][8][9][10][11].
Intensity diffraction tomography (IDT) refers to a class of 3D imaging techniques that employ tomographic phase reconstruction from intensity-only measurements. One IDT approach [29][30][31][32] combines a defocus-based phase contrast technique [38] and diffraction tomography model [39] to recover 3D phase. The measurement involves taking multiple defocused images while rotating the sample or the illumination. Recent works further incorporate partially coherent illumination using symmetric [33][34][35][36] or asymmetric [37] light source to achieve up to 2× resolution improvement in the recovered phase. However, changing focus not only requires mechanical scanning, but also increases the acquisition time and data size, both of which are undesirable for high-throughput applications.
An alternative IDT approach extracts 3D phase information using angled illumi-nation without mechanical scanning [7][8][9]. In [7], 3D phase contrast is computed with an algorithm combining lightfield refocusing [40] and differential phase contrast [41]; however, the results suffer from low-resolution as diffraction effects are neglected. In [8], a multislice model is proposed to incorporate diffraction and multiple forward scattering, and achieves high-resolution 3D recovery. However, since the multislice model is nonlinear, it necessitates an iterative reconstruction algorithm, which is non-ideal for time-constrained applications. In [9], an iterative algorithm combining Fourier ptychography [5] and the first Born approximation [39] was proposed. However, their model ignored the interference term between the scattered and unscattered fields. Here, we show that this term is the primary source of 3D phase contrast in IDT measurements.
In this work, we develop a novel linear IDT model that relates the sample's 3D permittivity contrast to intensity measurements using angled illumination (Fig. 1). Previous efforts [33][34][35][36][37] formulate the phase-intensity mapping in the 3D Fourier space; this approach unfortunately suffers from stringent sampling requirement for the measurement, resulting in hundreds of defocused images needed in practice. In addition, the reconstruction requires computation and memory intensive 3D deconvolution. Our approach overcomes all these limitations by employing a slice-based framework. The 3D sample is first modeled as a series of 2D slices along the axial direction. We then derive the slice-wise (2D) phase and absorption transfer functions (TF) at different depths for each illumination angle. We show that this framework enables flexible and efficient data acquisition, allowing using arbitrary patterning of the illumination angles and much fewer images required compared to other techniques [9,[33][34][35][36][37]. Our model fully accounts for the interference between the scattered and unscattered fields. The linearization is achieved via the first Born approximation that considers single scattering and neglects higher order nonlinear effects. Our linear model enables non-iterative 3D reconstruction directly from intensity-only measurements. We demonstrate volumetric reconstructions with closed-form Tikhonovregularized solutions, implemented using a computation and memory efficient 2D FFT-based algorithm. We demonstrate our technique on both stained and unstained thick biological samples. The imaging performance and limitation is further investigated in the presence of strong multiple scattering. Experiments demonstrate that our technique is robust even for samples with large permittivity contrast.

Theory and method 2.1 Forward model
The scattering of a sample can be characterized by its scattering potential V ( r) = 1 4π k 2 0 ∆ ( r) [42], where ∆ = ∆ Re +i∆ Im = − 0 is the permittivity contrast between the sample and the surrounding medium 0 . The real part ∆ Re characterizes the phase effect, and the imaginary part ∆ Im describes absorption. k 0 = 2π/λ is the wave number in free space, where λ the wavelength of the illumination. r = ( x, z) denotes the 3D spatial coordinates, with transverse coordinates x and axial position z.
We employ the first Born approximation [42] to model the light-sample interaction. Given the incident field f i , the total field f after propagating through the sample is given by where G( r) = exp(ik| r|)/| r| is the outgoing Green's function and k = √ 0 k 0 . The integral term represents the field f s scattered off the sample. The plane-wave illumination is modeled as f i ( r| u i ) = S( u i )e −i( u i · x+η i z) with transverse frequency u i and axial frequency η i = k 2 − | u i | 2 . S( u i ) represents the intensity of the i th LED, which we model as a point source [4][5][6][7][8][9][10][11]. The notation a(·| u i ) indicates the illumination direction for the field a. As the field propagates through the microscope, it is filtered by the pupil P ( u), corresponding to a convolution with the coherent point spread function h( x). The resulting intensity at the back focal plane of the microscope where * denotes 2D convolution. For notational simplicity, we have neglected the microscope's magnification. The total field at the front focal plane (z = 0) is in which the scattered field is simplified using the Weyl expansion of the Green's function [42] as where · represents the 2D Fourier transform (FT); ∆ ( u, z) is the 2D FT of the slice ∆ ( x, z) at depth z; F −1 {·} denotes the 2D inverse FT (IFT); u is the transverse frequency variable, and η( u) = k 2 − | u| 2 the axial frequency. (??) calculates the total scattered field as the coherent superposition of the scattered field from each slice; inter-slice scattering is ignored as a consequence of the first Born approximation. Importantly, (??) contains the interference information between the unscattered and scattered field, and can be expanded into four terms: represents squared modulus of the scattered field, which is negligible if the permittivity contrast is small, and can thus be dropped [43]. Hence, the phase and absorption information is mainly contained only in the two cross-terms, is . This is also highlighted by the Fourier spectrum of the intensity measurement in Fig. 1. The relation between the sample's permittivity contrast and the intensity spectrum is thus linear; and the phase and absorption terms are decoupled: where H Re and H Im are the angle-dependent phase and absorption TFs for each sample slice at depth z, respectively, and are Equations (??-??) define the forward model of our technique. From these equations, the computational complexity of our model is set by the N z evaluations of where N x , N y , N z correspond to the dimension of the data in 3D. The TFs have the following essential properties for scan-free 3D reconstruction. (a) The transverse frequency of the incident field u i determines the off-axis shift of the pupil function. Due to intensity-only measurement, a pair of shifted pupils, shifting to opposite directions, are super-imposed in the TFs (akin to the "twin-image" in holography). This is illustrated in the computed phase and absorption TFs in Fig. 1(c), and validated experimentally by visualizing the intensity spectrum of the raw images in Fig. 1(b). (b) The depth information is encoded in the phase term, which enables scan-free 3D reconstruction. Specifically, the linear phase term e iη i z corresponds to a geometrical shift in the real space, which was accounted for previously using the lightfield model [7]; the propagation term e iη( u+ u i )z models the diffraction effects. (c) At the focal plane, assuming a real and symmetric pupil function (i.e. an ideal unaberrated microscope), the phase TF is imaginary and anti-symmetric, whereas the absorption TF is real and symmetric. In general, because the phase and absorption information are encoded asymmetrically in the Fourier space [ Fig. 1(c)], it leads to measurable intensity contrast difference for phase and absorption features, and allows inversion of both quantities simultaneously. (d) The support of the combined TFs sets the resolution limit of our technique, which can be analyzed following the same framework as [33][34][35][36][37]. Since the forward model relies on the existence of a strongly unscattered field, it is only valid for brightfield measurements. As a result, we use illumination angles up to the objective NA. This allows us to achieve resolution equivalent to the incoherent limit. In the transverse direction, the support is uniformly 4NA/λ. In the axial direction, the system suffers from the missing cone problem [43,44], resulting in the axial elongation present in the reconstruction. The axial Fourier coverage varies with sample's feature size, and is up to 2 − 2 1 − NA 2 /λ.

Inverse problem
The reconstruction algorithm combines all the intensity images to estimate the sample's complex permittivity contrast (i.e. phase and absorption) in 3D. Since our forward model enables inferring the full field information, the reconstruction essentially stitches all the Fourier components coherently, akin to the 3D synthetic aperture [45,46]. To facilitate efficient reconstruction, we propose a slice-wise deconvolution algorithm, in which the phase and absorption are reconstructed slice-by-slice.
To implement our algorithm, we first discretize the 3D sample into a stack of 2D slices, equivalently replacing the integral in (??) by a discrete sum over the slice index. Next, each intensity image is pre-processed to perform background normalization and removal: g = (I − I i )/I i . The corresponding discretized and normalized TFs are: H Re = H Re /I i and H Im = H Im /I i . This normalization also removes the unknown scaling in the source intensity. Direct deconvolution is known to suffer from noise amplifications [47]; therefore we employ Tikhonov regularization that imposes an energy constraint to suppress these artifacts. The closed-form solutions for phase and absorption are

Results
Our system consists of a Nikon TE 2000-U microscope with a programmable red (central wavelength 630nm) LED array source (specifications same as [11]) placed 79mm above the sample [ Fig. 1(a)]. An sCMOS camera (Pco.Edge 5.5) is used for image acquisition, which is synchronized in real-time with the LED source via a microcontroller. Since our technique completely removes any mechanical scanning, the acquisition speed is only limited by the camera's frame rate (up to 50Hz). To achieve the incoherent resolution limit, we acquire data using brightfield LEDs fully covering the NA of the objective used in each experiment. Angle-varying intensity images are captured by sequentially turning on one LED at a time. We conduct experiments using 10× (0.25 NA, CFI Plan Achro) and 40× (0.65 NA, CFI Plan Achro) microscope objectives (MO). For the 10× MO, 89 LEDs are within the brightfield region; whereas 697 LEDs for the 40× MO. We only use a small subset of the LEDs, as detailed in each experiment section. Different LED sampling patterns are explored and their reconstruction results are compared experimentally, demonstrating the flexibility in data acquisition enabled by our IDT framework. In all the images, the reconstructed permittivity contrast has been converted to the real (phase) and imaginary (absorption) part of the refractive index. The background medium index for all the biological samples are assumed to be 1.33 (i.e. close to water).

Imaging of stained 3D sample
We first demonstrate our technique to image a stained spirogyra sample (Fisher Scientific S68786) using the 10× MO. The sample contains both highly absorbing features (e.g. chloroplasts) and "phase" features (e.g. filaments) orientated in 3D, as seen in the full field-of-view (FOV) brightfield image in Fig. 2(a). We used all of the 89 brightfield images to perform the IDT reconstruction. To demonstrate the best possible axial resolution, the slice spacing is set to be 5µm in the reconstruction, corresponding to approximately 2× oversampling axially. We have reconstructed 25 phase and absorption slices equally spaced between −20µm and 200µm. The number of images used is about twice the number of unknowns, the underlying linear problem is thus over-determined, albeit ill-conditioned due to the missing axial frequency information that in turn sets the axial resolution. The reconstruction performance is illustrated on features with different length scales. In Fig. 2(b), we zoom in on a region containing four clustered spirogyras. Our reconstruction demonstrates that axial sectioning can be successfully obtained using our IDT algorithm -features at different depths are clearly distinguished in the reconstruction. In Fig. 2(c), we turn to a spiral structure on a single spirogyra. The helical structure that presents much finer axial features, is successfully reconstructed. Finally, we look at the filaments in Fig. 2(d), which is a "phase sample" that does not provide high contrast in the brightfield image. Using our technique, the filaments are clearly resolved with high contrast in the phase reconstruction. As expected, the absorption reconstruction does not provide much contrast.

Imaging of unstained dense cell clusters
Next, we image unstained MCF-7 cancer cells fixed in formaldehyde solution using the 40× MO (Fig. 3). The sample contains both monolayer cells and dense cell clusters. A qualitative visualization of the sample is performed via phase contrast (PhC) mode (MO: 0.65 NA, CFI Plan Achro), shown in Fig. 3(a). Features imparting longer optical path delay generally produce darker contrast in the PhC image. Halos are present at the boundaries of thick cell regions.
We first use 153 images to reconstruct 22 phase and absorption slices with 1µm spacing from −6µm to 15µm. To demonstrate the versatility of our IDT technique, Fig. 3(b) shows the phase reconstructions of a few cell regions covering both thin and thick features, along with their PhC images for comparison. Subcellular features are correctly reconstructed, matching the PhC images. Features appear darker contrast in PhC are reconstructed with higher refractive index values.
Next, we zoom-in on a region containing a dense cluster of 5 cells with partial axial overlapping [ Fig. 3(c)]. Phase reconstructions are shown at 6 slices along with the corresponding PhC images, captured by mechanically adjusting the focus. Since both techniques use high angle illumination that results in similar Fourier coverage, we expect that they should provide similar lateral resolution and axial sectioning capability. This can be qualitatively verified in Fig. 3(c).
Finally, we demonstrate the flexibility and robustness of our technique when different illumination patterns and number of LEDs are used for data acquisition. In Fig. 3(d), we investigate illumination strategies using symmetric and pseudorandom patterns. Each symmetric pattern uses LEDs that are equally spaced along both azimuthal and radial directions. Each pseudorandom pattern is designed such that each quadrant contains the same number of randomly located LEDs. In all cases, we reconstruct the same 22 phase and absorption slices (i.e. in-total 44 unknown slices). With 37 images, the number of measurement is slightly smaller than the number of unknowns. Although the problem is under-determined, the reconstruction only degrades slightly as compared to the ones from using more images (e.g. 105 and 153). When we further reduce the number of images, the reconstruction further degrades, but is still able to recover cellular features. Using the pseudorandom pattern does not make significant difference when using a large number of LEDs (e.g. 105 and 153). As the number of LEDs is reduced, the pseudorandom pattern is observed to achieve slightly better results, likely because it can provide more uniform Fourier coverage.

Imaging of strongly scattering sample
Finally, we evaluate our technique when imaging a strongly scattering sample, which contains both highly absorbing and strong phase features with high permittivity contrast. The goal is to investigate the performance of our technique when multiple scattering effects become strong. The influence of multiple scattering to ODT has been investigated extensively [48,49]. A general observation is that the unaccountedfor multiple scattering leads to under-estimation of the permittivity contrast. The stronger the multiple scattering, the more severe the under-estimation [48]. In IDT, the data is further limited by the intensity-only information. Our forward model only considers single scattering; intensity variations due to multiple scattering essentially act as non-random noise in the inversion, so they can only be partially suppressed by the Tikhonov regularization.
The sample consists of an absorption target (58-198, Edmund Optics) placed underneath a phase target (QPT, Benchmark Technologies), separated axially by ∼ 790µm [Fig. 4(a)]. The absorption target contains chromium patterns on a glass substrate. The phase target consists of structures of different sizes and heights with refractive index 1.52, which are raised above a glass substrate. All the experiments are performed in air, so the permittivity contrast is large (¿1), violating the requirement of the first Born approximation. To investigate how the reconstruction degrades as the strength of scattering increases, we perform controlled experiments on phase patterns with increasing height (50nm, 100nm, and 200nm), while keeping the same region of the absorption target underneath. To further investigate the influence of sample depths, we repeat the experiments at three different focal planes. The first dataset is taken with the phase target in focus [(A) in Fig. 4]; the second with focus 300nm below (A) [(B) in Fig. 4]; the third with the absorption target in focus [(C) in Fig. 4]. In each experiment, we take 89 brightfield images with the 10× MO. A sample brightfield image from the on-axis LED is shown in Fig.4(b).
Since the sample consists of only two layers of interest, the reconstruction can be efficiently performed just on the two slices, without the need to compute any slices in-between. To implement the reconstruction, we use an axial sampling ∆z = 10µm, matching the theoretical resolution. Since ∆z is much larger than the pattern height, the reconstructed permittivity contrast ∆ rec should be interpreted as the average over the slice thickness. Given the refractive index of the phase target n ph , we can convert the averaged permittivity contrast to an estimate of the pattern height by h rec = ∆ rec ∆z/(n 2 ph − 1). The reconstructed pattern heights for the three different cases are shown in Fig. 4(c). Due to the intensity-only measurement, the phase's DC component is always lost, and thus all the height estimates are centered around zero. Under-estimation of the heights is observed in all cases due to the strong multiple scattering that is unaccounted for in our model. As we image a taller pattern, the multiple scattering becomes stronger; correspondingly, the reconstruction shows a larger error in the estimation. The strong multiple scattering also leads to the hazy background in each reconstructed image.
To quantitatively validate these results, we perform simulations on a similar twolayer sample [ Fig. 4(d)]. The multislice model [8] is used to simulate the intensity measurements since it can accurately model multiple scattering at this length scale [50].
We then use our IDT model to perform the reconstruction. As seen in the recovered images, the amount of height under-estimation in each case matches well with the corresponding experimental results. A similar level of background noise is also observed.
Using our IDT model, we expect to achieve the same lateral resolution at all slices, set by the bandwidth 4NA/λ. As shown in cases (A), when the phase target is near focus, we successfully achieve the theoretically predicted resolution (1.54µm based on the Rayleigh Criterion). However, as the focal position is defocused from the target, the resolution degrades, [(B) and (C) in Fig. 4(c)]. The same trend is observed for the absorption target [ Fig. 4(e)]. A similar resolution degradation was also observed in [8] using a similar setup, though the experiments reported here span ∼ 8× wider defocus range (e.g. 800µm here vs 110µm in [8]). We expect that this degradation is likely due to LED position mis-calibration. The amount of geometrical shift scales linearly with the defocus distance due to the linear phase term in the TFs. Thus, the same amount of angular mis-calibration would result in a larger error in data with farther defocus. Notably, despite this resolution degradation, the recovered phase values remain consistent regardless of the focal positions, as shown in the cutlines of the estimated heights in Fig. 4(c).

Conclusion
We have demonstrated a new computational microscopy technique that enables scanfree, 3D phase and absorption reconstruction using intensity-only measurements. Our slice-based intensity diffraction tomography framework allows flexible illumination patterning for data acquisition and a direct inversion algorithm that has both low computational complexity and is memory efficient. We have experimentally validated this technique on dense 3D biological samples. Although our model only accounts for single scattering, the reconstruction is robust for samples with high permittivity contrast. The effects of multiple scattering are examined via both simulation and experiment, and their results have been found to be in good agreement. The system is simple, and built directly on a commercial microscope with an LED array to enable rapid illumination-angle scanning. This is particularly attractive for the wide adoption of our system to existing microscopy facilities and may open up many biomedical imaging applications.