## Abstract

We propose an accurate and computationally efficient 3D scattering model, multi-layer Born (MLB), and use it to recover the 3D refractive index (RI) of thick biological samples. For inverse problems recovering the complex field of thick samples, weak scattering models (e.g., first Born) may fail or underestimate the RI, especially with a large index contrast. Multi-slice (MS) beam propagation methods model multiple scattering to provide more realistic reconstructions; however, MS does not properly account for highly oblique scattering, nor does it model backward scattering. Our proposed MLB model uses a first Born model at each of many slices, accurately capturing the oblique scattering effects and estimating the backward scattering process. When used in conjunction with an inverse solver, the model provides more accurate RI reconstructions for high-resolution phase tomography. Importantly, MLB retains a reasonable computation time that is critical for practical implementation with iterative inverse algorithms.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## Corrections

Michael Chen, David Ren, Hsiou-Yuan Liu, Shwetadwip Chowdhury, and Laura Waller, "Multi-layer Born multiple-scattering model for 3D phase microscopy: publisher’s note," Optica**7**, 923-923 (2020)

https://www.osapublishing.org/optica/abstract.cfm?uri=optica-7-8-923

27 May 2020: A typographical correction was made to Eq. (6).

## 1. INTRODUCTION

Phase imaging methods that reconstruct the 3D refractive index (RI) of semi-transparent objects generally rely on two key components: a light scattering model and a phase retrieval algorithm. Optical diffraction tomography (ODT), for example, scans the angle of incident light and captures many 2D complex fields using interferometry [1–4] or pupil plane modulation [5–7]. Intensity-based 3D phase methods use through-focus or color images from off-the-shelf microscopes [8–12]. 3D RI can be recovered via 2D phase retrieval followed by a filtered back propagation or a 3D deconvolution [10,13,14]. All of these methods use a weakly scattering approximation (first Born or Rytov), which works well for thin index-matched samples (e.g., single cells), but causes errors for thick samples that incur multiple scattering.

For 3D phase imaging with optically dense samples, a forward model that describes light propagation beyond the single-scattering regime is needed. The multi-slice (MS) beam propagation method accounts for multiple-scattering effects in a computationally efficient manner, by approximating the sample as a sequence of infinitesimally thin planar slices along the optical axis, each modeled as a 2D transmission layer separated by a uniform medium. MS has been used with both interferometric and intensity-only measurements to reconstruct multiple-scattering samples [15–17]. However, MS has two key drawbacks: 1) quantitative accuracy of MS reconstructions may decrease as the angle of the illumination beam increases. Heuristic approaches have been proposed to mitigate this issue [18,19], but they are either unstable or greatly increase computation cost. 2) MS does not model back-scattering, which contains valuable high-resolution information from a sample’s missing cone region [20].

Scattering models have been developed to generate forward and backward scattering fields close to the analytic solution of the Helmholtz equation via a series of convolutions with a 3D Green’s function. These include the recursive Born [21–23], contrast source inversion [24], coupled dipole [25], hybrid method [26], and series expansion with accelerated gradient descent on the Lippmann–Schwinger equation (SEAGLE) [27]. Even after significant advances in improving the rate of convergence and memory requirement of SEAGLE [28,29], it still requires orders-of-magnitude more computation and memory than weakly scattering or MS models, due to operations on 3D arrays, which limits application to 2D slice reconstructions or small-volume 3D phase tomography with iterative inverse algorithms [30,31].

Here, we introduce a new multiple-scattering model that calculates both the forward and backward scattering fields accurately, even for high-angle illumination, in a computationally efficient manner. Our multi-layer Born (MLB) model decomposes the 3D object into a series of 3D slabs with finite thickness, then sequentially applies the first Born scattering process to each slab. This removes the paraxial assumption of the MS method, achieves accurate prediction for high-angle illumination, and also enables computation of backward scattered light from the sample. When used with an inverse algorithm, MLB recovers quantitative RI with lower error than and computation time similar to MS. These properties make MLB suitable for imaging large 3D volumes with iterative inverse algorithms. Here, we implement MLB in a phase tomography framework using intensity-only images acquired on a LED array microscope [10,15,32,33].

## 2. MULTI-LAYER BORN SCATTERING MODEL FOR 3D PHASE TOMOGRAPHY

Generally, 3D phase tomography methods solve an inverse problem by minimizing the difference between the measurements and the expected measurements, given an estimate of the sample’s 3D complex RI distribution. Iterative algorithms compute the expected measurements at each iteration via a forward model that describes how light propagates through the sample. The forward model is critical, as its accuracy affects the reconstruction quality, and its computational efficiency determines compute time, since it must be invoked dozens or hundreds of times before the solution converges. Below we describe our MLB forward model and the inverse algorithm that we use to demonstrate it.

#### A. Forward Scattering Model

A thick sample can be described by its 3D scattering potential, $ V(x,y,z) $, which is related to its RI $ n(x,y,z) $, via the following expression:

where $ (x,y,z) $ are the 3D coordinates, $ {k_o} = \frac{{2\pi }}{\lambda } $, $ \lambda $ is the wavelength of incident light in vacuum, and $ {n_{\text{b}}} $ denotes background RI of the surrounding media. As light propagates through the sample, the incident field gets scattered; this scattered field interferes with the incident field to form the total field. Modeling such dynamics, including both forward and backward scattering interactions, is complicated and highly nonlinear; hence, we often make simplifying approximations.The first Born approximation assumes weak scattering, resulting in a model where the scattered field is linearly related to the scattering potential by a Green’s function, $ G(r) = \frac{{ - \exp ( {\text{i}{k_o}{n_{\text{b}}}\text{r}} )}}{{4\pi \text{r}}} $, where $ r = \sqrt {{x^2} + {y^2} + {z^2}} $ [21]. For an incident field, $ {U_{\text{in}}} $, the total field, $ {U_{\text{tot}}} $, can be written as

Our MLB scattering model includes multiple-scattering effects by dividing the 3D object into multiple slabs and sequentially applying the first Born scattering process on each. Each slab has a thin but finite thickness, $ \Delta z $ (see Fig. 1). Therefore, Eq. (2) describes the field after the wave propagates through the $ {n^{\text{th}}} $ layer. This serves as the incident field of the $ {(n + 1)^{\text{th}}} $ layer. By applying the first Born scattering process on each layer recursively, the total field computed at the last layer is the multiply scattered exit field. We let $ {U^n} $ and $ U_{\text{s}}^n $ denote the incident field and the scattered field, respectively, within the $ {n^{\text{th}}} $ layer. For mathematical simplicity, we define the $ {n^{\text{th}}} $ layer as the one that occupies the space $ z \in [ {(n - \frac{1}{2})\Delta z,(n + \frac{1}{2})\Delta z} ] $. The recursive formula can be written as

The first term on the right-hand side of Eq. (3) is the incident field at the center of the $ {n^{\text{th}}} $ layer propagated by one slab thickness $ \Delta z $. Given the 2D Fourier transform operator $ {\cal F}\{ \cdot \} $, the spatial frequency spectrum of the incident field $ {\tilde U^n}(\textbf{u}) = {\cal F}\{ {U^n}({\boldsymbol \rho },n\Delta z)\} $, and the angular spectrum propagation kernel $ P(\textbf{u},z) = {e^{\text{i}2\pi z\sqrt {{{({n_{\text{b}}}/\lambda )}^2} - \parallel \textbf{u}{\parallel ^2}} }} $ with the 2D spatial frequency space coordinates vector $ \textbf{u} $, we get

To simplify Eq. (4), we find it useful to use the 2D Fourier spectrum of the Green’s function [34]:

Equations (3), (5), and (9) make up the forward scattering part of our MLB model. Equation (9) is computationally expensive; however, when $ \Delta z $ is chosen to be sufficiently small, the sinc function is approximately one, and the integration becomes a convolution, which can be efficiently evaluated via fast Fourier transform (FFT).

Usually, the total field after the light passes through all the layers is imaged by a low-pass imaging system, where the spatial frequency bandwidth is set by the numerical aperture (NA) of the objective lens. Hence, we digitally refocus the output field to the image plane and apply a low-pass filter. With the ideal circular low-pass filter $ C(\textbf{u}) = \text{circ}(\textbf{u}\lambda /\text{NA}) $, the measured forward scattered field can be efficiently computed using FFTs, as shown in Algorithm 1.

Both the MS and MLB scattering models break the object into layers and model the multiple-scattering effects as light propagates through the sample. One major difference between the two models is that MLB considers non-paraxial 3D scattering effects within each layer, while MS simplifies it with a 2D transmission function. Hence, MLB is more accurate for high-angle illumination, or for highly scattering samples where more of the propagating light is oblique. MLB is also capable of modeling the backward scattering field (see derivation in Supplement 1).

#### B. Inverse Problem for 3D Intensity-Based Phase Tomography

To demonstrate the use of our MLB model, we incorporate it as the forward model of an iterative inverse problem that uses intensity-only images to reconstruct 3D phase (RI) information. An intuitive way to probe the 3D structure of an object is by rotating the object [35,36]; however, this generally requires additional hardware and limits the types of objects that can be imaged [14,37]. Alternatively, we can illuminate the object at various angles [20], as in many ODT systems that use either scanning galvanometers or spatial light modulators (SLMs) [1,6,17,38]. For intensity-based 3D phase tomography techniques, the requirement of spatial and temporal coherence of the light source is less stringent, and the illumination scanning process can be realized using an inexpensive LED array [9,12,15,32,33]. Combining the illumination scanning scheme and the proposed MLB scattering model, we develop a 3D intensity-based phase tomography algorithm that applies to multiple-scattering objects. Similar to previously proposed multiple-scattering methods [15,16,27], an iterative algorithm is needed to recover the 3D RI of the sample because the measurements are nonlinearly related to the scattering potential. First, we formulate 3D intensity-based phase tomography as an optimization problem with an objective function $ {\cal L}(V) $:

In this paper, we choose total variation (TV) as the regularization function because it has an efficient proximal operator [40]. Nesterov’s acceleration method is used at the end of every iteration to further speed up convergence [41], and a momentum restarting mechanism [42] promotes stability. Detailed implementation of our proposed 3D intensity-based phase tomography with MLB is summarized in Algorithm 3. Once the 3D complex scattering potential is recovered, the real-valued RI can also be retrieved from Eq. (1), as shown in our simulations and experiments (Section 3 and Section 4). Based on Algorithms 1–3, our 3D phase imaging framework using MLB is summarized in Fig. 1.

Note that our inverse problem and forward scattering model can be interchanged with others. For complex-field measurements from holographic phase tomography (instead of intensity-only datasets), we can incorporate our MLB model by simply changing the data fidelity term in Eq. (10). Our inverse Algorithm 3 is also compatible with other scattering models by using a different forward model and deriving its gradient (Algorithms 1,2). In both simulations and experiments, step size $ \alpha $ is set within a range from 0.5–5, and $ \tau $ is chosen to be around $ {10^{ - 2}} $ for the first Born, Rytov, and MLB models. In contrast, these two parameters are two to three orders of magnitude smaller when MS is applied, due to the large local Lipschitz constant of the objective function when using MS, in which RI is solved instead of the scattering potential. For each model, $ \alpha $ is manually tuned to achieve the best rate of convergence, and $ \tau $ is adjusted to mitigate most of the noise artifacts without degrading physical structures.

## 3. SIMULATIONS

#### A. Comparison of Forward Models

To analyze the accuracy of our MLB forward model and compare with other scattering models used in 3D phase imaging (first Born, Rytov, and MS) independent of the inverse solver, we first simulate the amplitude measurements of a 3D cell phantom with different methods. Since we do not have ground truth for the measurements, we compare against those generated by SEAGLE [27], which should be the most accurate method, but requires very long computation times. Here, since we are using only SEAGLE to compare forward model accuracy, we need not run it for the entire large-volume 3D reconstruction.

The RI contrast of the 150-layer cell phantom ensures multiple scattering: it has a $ 15 \times 15 \times 7.5\,\unicode{x00B5} {{\rm m}^3} $ ellipsoid body, which contains cytoplasm (RI $ n = 1.35 $), a nucleus ($ n = 1.33 $), two nucleoli ($ n = 1.36 $), several small organelles ($ n = 1.39 $), and a thin plasma membrane ($ n = 1.37 $) enclosing the body. The cell is surrounded by medium of $ n = 1.33 $, within a volume of $ 450 \times 450 \times 150 $ voxels with $ 0.05\,\,\unicode{x00B5}\text{m} $ resolution. Assuming an illumination wavelength of 532 nm and objective $ {\text{NA}_{\text{obj}}} = 0.8 $, we calculate the amplitude images using forward and backward scattered light under on-axis ($ {\text{NA}_{\text{illu}}} = 0.00 $) and off-axis ($ {\text{NA}_{\text{illu}}} = 0.75 $) illumination.

Processing was done on a NVIDIA TITAN X GPU installed on a desktop computer (Intel i7-5960X CPU). The forward model computation time for one illumination angle with first Born, Rytov, MS, and MLB were 0.22, 0.22, 0.25, and 0.32 s, respectively. The memory requirement for SEAGLE was too large to run on the GPU, and the computation took ${\sim}15 $ min for a single illumination angle on our CPU.

Ground truth (SEAGLE) in-focus amplitudes at two different illumination angles are shown in the top row of Fig. 2(a), and the error maps (difference from ground truth) for each forward model are in the bottom four rows. Since the weakly scattering approximation does not hold, both first Born and Rytov methods produce large error. In contrast, MS and MLB scattering models account for multiple scattering within the object, greatly increasing accuracy. For off-axis illumination, however, the MS model incurs error due to the paraxial approximation. We show both forward and backward scattering predictions of the amplitude (the MS model does not predict back scatter). Overall, the MLB is the most accurate among the computationally efficient models.

There is a trade-off between the model accuracy and the maximum phase within each layer when using MLB, since the first Born approximation for each slab becomes less accurate as the phase contrast in each slice becomes larger. Figure 2(b) illustrates this trade-off by plotting the evolution of mean squared error (MSE) versus the maximum phase change in each layer. By binning the neighboring layers and adapting a larger $ \Delta z $ along the $ z $ axis, we effectively increase the maximum phase in each layer. The MSE of the first Born results is not displayed because its value is an order of magnitude larger than the others. In the forward scattering case, the MSE of the MLB results grows linearly with the phase in each layer, while the MSEs of Rytov and MS remain fairly flat. The MLB model provides the most accurate predictions for off-axis illumination when the phase within each layer is smaller than $ 0.04\pi $ for the cell phantom. For backward scattering, MSEs are smaller due to weaker backward scattered light, but the error quickly increases as the phase change in each layer grows. Similarly, the MLB provides significant accuracy improvement when thinner slices are used, since the weakly scattering assumption is valid within each layer.

#### B. Comparison of 3D RI Reconstructions

To understand how the accuracy of the forward scattering model used in Algorithm 3 affects 3D RI reconstructions, we compare the results from the four different models described in Section 3.A. SEAGLE was again used to simulate the measurements, where 104 amplitude images were generated with illuminations ($ {\text{NA}_{\text{illu}}} \approx 0.75 $) from an annular region on the source plane. The rest of the imaging parameters (NA, pixel size, wavelength) remain the same as listed in Section 3.A. Measurements contain high spatial frequency content and rich phase contrast when the object is illuminated by oblique light with $ {\text{NA}_{\text{illu}}} $ close to $ {\text{NA}_{\text{obj}}} $ [43], as seen in the first row of Fig. 3. Hence, the high-resolution 3D RI information of the sample is well encoded in the measurements.

We compare the ground truth RI of the cell phantom with the results computed using different forward models and our inverse solver. Accuracy is quantified by the MSE of each reconstruction, which are: first Born 0.0159, Rytov 0.0076, MS 0.0088, and MLB 0.0054. Since the first Born method overestimates the amplitude values on the image plane [(see Fig. 2(a)], the RI retrieved using it is underestimated. The Rytov model results in a more accurate reconstruction due to a less strict approximation, but suffers from incorrect background RI and less uniform RI inside the object. In contrast, MS and MLB recover RIs with more uniform background, which are similar to the ground truth. However, the 3D phantom recovered using MS introduces undesired contrast on the $ x - y $ cross section (elongation indicated by yellow arrows in the MS panel of Fig. 3), and demonstrates consistently higher RI values in the internal cytoplasm, compared to that of the ground truth. This is why the MSE with MS is slightly higher than with the Rytov model, though MS provides better visual appearance. Despite suffering from the missing cone problem and sacrificing the axial resolution for low spatial frequencies, our algorithm with MLB recovers the most quantitative 3D RI because of the high accuracy of the scattering model.

Total reconstruction times using our algorithm with 100 iterations were 0.67 (first Born), 1.19 (Rytov), 1.50 (MS), and 1.67 (MLB) hours. The differences are caused mainly by the various computation efficiencies of each scattering model. The computation time grows linearly with the product of number of layers, number of iterations, and number of measurements, but barely changes as the number of pixels in each layer increases. This is because the most computationally expensive operation in each layer, the FFT, is efficiently parallelized on a GPU.

## 4. EXPERIMENTS

Experiments in Sections 4.A and 4.B use an LED array microscope to realize our proposed 3D intensity-based phase tomography. A planar $ 32 \times 32 $ programmable LED array is mounted on a Nikon TE300 inverted microscope, which generates 514 nm light from different angles of incidence [15]. Since the LEDs are far away from the sample relative to the field of view (FoV) (68 mm for the $ 40 \times 0.65\text{NA} $ objective, and 46 mm for the $ 60 \times $ objective), the incident light from each LED can be treated as a plane wave propagating at an angle determined by the LED position. All images from the LED array microscope are captured by a monochrome sCMOS camera (PCO.edge 5.5, $ 6.5\,\,\unicode{x00B5}\text{m} $ pixel size) at the front port of the microscope, where an additional $ 2 \times $ magnification is provided. To demonstrate high-NA imaging capability of the proposed method, an oil-immersion microscope is used (see system design details in Section 4.C). For each illumination with both the LED array microscope and high-NA oil-immersion microscope, we record two intensity images, with and without the sample. By dividing the former by the latter, the normalized intensity, $ {I_{\text{measure}}} $, serves as the input to the inverse problem in Algorithm 3.

#### A. 3D Imaging of Polystyrene Beads

In order to quantify accuracy for different scattering models in experimental conditions, we start by imaging a known sample—a $ 5\,\,\unicode{x00B5}\text{m} $ polystyrene bead (Sigma-Aldrich, $ n = 1.6010 $ at $ \lambda = 514\,\,\text{nm }$) immersed in two different index-matching oils (Cargille), resulting in RI contrasts of $ \Delta n = 0.0113 $ and $ \Delta n = 0.0452 $. We used a $ 40 \times 0.65\,\,\text{NA} $ objective lens (Nikon CFI Plan Achromat), with its front focal plane aligned at the center of the bead. One hundred images with angles of incidence ranging from approximately 23° to 44° were captured at 6.66 frames per second (fps). Reconstructions using each of the four forward models (all with a positivity constraint and TV regularization in the proximal step) are shown in Fig. 4.

For the low index contrast case [Fig. 4(b)], the sample can be considered weakly scattering, and so all of the models are valid; hence, the single-scattering and multiple-scattering models yield similar results. We note that low SNR and partial coherence effects occur when using an LED array for intensity-based imaging of a weakly scattering sample. These effects are model mismatches that contribute to quantitative inaccuracies, as shown when comparing the quantitative RI profiles of the tomographic reconstruction models to the ideal RI profile.

In the case of high index contrast, the resulting multiple scattering causes RI accuracy to vary significantly depending on the tomographic reconstruction framework [Fig. 4(c)]. The first Born approximation fails by underestimating the RI value and missing content in the center of the bead. The Rytov model mitigates the RI underestimation, but results in a corrupted shape axially. Reconstructions with MS and MLB both account for multiple scattering, but the reconstructed bead using MS has an elongated shape and slightly less quantitative contrast, similar to our simulation. The recovered polystyrene bead using MLB has higher quantitative accuracy and more isotropic resolution. This demonstrates that MLB is the best forward model for imaging objects of high RI contrast and multiple scattering.

#### B. 3D Imaging of Weakly Scattering Sample

Next, we look at a weakly scattering biological sample, a fixed 3T3 cell. We capture 100 intensity measurements at 2 fps using a $ 60 \times 0.80\,\text{NA} $ objective lens (Nikon CFI Achromat), with illumination angles ranging from 49° to 51°. Since the illumination NA is close to the objective NA, we simultaneously encode high-frequency and low-frequency phase contrast of the 3D sample in the measurements. This is critical for high-resolution phase tomography. The images at the camera plane are over-sampled due to the additional magnification provided by the tube lens, so we down-sample the raw data by a factor of $ 2 \times $ in each dimension.

Figure 5(a) shows orthonormal cross sections of the recovered 3D RI (after halo-correction [44]). The reconstruction volume contains $ 612 \times 612 \times 120 $ voxels with voxel size $ 0.108 \times 0.108 \times 0.540\,\unicode{x00B5}{{\rm m}^3} $. In this case, neither the positivity constraint nor TV regularization were used. We can observe the components within the cell, such as the nucleus, indicated by arrows. The physical thickness of the cell is ${\sim}17\,\,\unicode{x00B5}\text{m} $. We also render the 3T3 cell with false colors labeling different RIs for visualization in a 3D volume [see Fig. 5(b) and Visualization 1].

Besides providing quantitative density of the sample, having 3D phase information also gives better contrast when comparing a slice of the 3D reconstruction to a conventional 2D phase contrast image. Figure 5(c) shows several regions of interest (ROIs) and depths for three different imaging modalities: 2D phase contrast from assymetric half-circle illumination, quantitative phase using differential phase contrast (DPC) with four measurements [43,45], and a slice of our recovered 3D RI. While all theoretically have similar depths of field (${\sim}0.64\,\,\unicode{x00B5}\text{m} $), which helps reveal different structures at each plane (e.g., vesicles indicated by the arrows), their contrasts vary dramatically. Out-of-focus information is integrated with in-focus cell content in the 2D phase image, making it hard to distinguish one from the other. Although the phase contrast images highlight the in-focus components of the cell, it is not possible to interpret the actual structure from only phase contrast images. In contrast, 3D RI reconstruction naturally displays quantitative optical density at each depth and builds up high contrast contours that separate distinct organelles.

#### C. 3D Imaging of Multiple-Scattering Sample

MLB describes the multiple-scattering process based on scalar wave theory in the object space without any assumptions on the imaging system. Hence, it naturally applies to many different microscope setups in addition to the LED array microscope. Moreover, unlike the MS or beam propagation method, MLB considers non-paraxial wave interaction.

To demonstrate, we conduct intensity-based phase tomography on a
custom-built high-NA optical microscope [17]. A fiber-coupled LED ($ \lambda =
530\,\,\text{nm} $ center wavelength) is combined with a
2D scanning mirror to achieve programmable angle-scanning
illumination. Two high-NA oil-immersion objective lenses (effective $ \text{NA} = 1.45
$) serve as the condenser and the
imaging lenses, which enables sub-wavelength resolution. A fixed
*C. elegans* sample is sandwiched between
two coverslips and intensity images are acquired for each illumination
angle. We collect 120 images with the illumination angles scanned on a
spiral trajectory. A computational self-calibration algorithm is used
to accurately measure the illumination angles [46]. Algorithm 3 is applied to reconstruct quantitative RI within a volume of $ 1200 \times 1200 \times
100 $ voxels, with voxel size $ 0.12 \times 0.12 \times
0.35\,\unicode{x00B5}{{\rm m}^3} $. The full length of the *C. elegans* (${\sim}1\,\,\text{mm
}$) does not fit in a single FoV, so 12
phase tomography datasets at different regions are captured and
digitally combined to visualize the entire worm.

Large differences among RI results using the first Born, MS, and MLB
models can be observed in Fig. 6. As expected, the RI contrast of the body of the *C. elegans* is large and the sample is thick, so
the weakly scattering approximation (first Born) fails to capture most
of the low-frequency content. The RI values recovered using MS are
slightly higher than RIs in the MLB case. This agrees well with our
simulation results. In addition, the halo artifacts pointed out by the
arrows in Fig. 6 are mitigated
by changing the scattering model from MS to MLB. All of the above
suggest that MLB works well with high-NA imaging systems and
multiple-scattering samples. Additionally, the computation and memory
costs of MLB are similar to MS, which is an efficient method to
evaluate 3D scattering. Hence, we are able to achieve 2 giga-voxels 3D
RI retrieval of the whole *C. elegans*,
shown at the bottom of Fig. 6,
on a desktop computer with GPUs.

## 5. CONCLUSION

We have introduced the MLB scattering model, an efficient and accurate forward model for 3D phase (RI) microscopy. MLB demonstrates superior accuracy for multiple-scattering objects, as compared to widely used single-scattering methods (first Born, Rytov) and the multiple-scattering MS method. Its moderate computation complexity and highly parallelizable operations, as compared to other multiple-scattering models (e.g., SEAGLE, FDTD), are important for making large-scale inverse problems feasible, since they require evaluation of the forward model at each of many iterations, which becomes computationally expensive for large datasets. Our MLB model does not rely on the paraxial approximation, so is suitable for high-NA imaging, and also predicts backscattered light, improving accuracy. Together, these advances open up 3D phase microscopy to an entirely new range of biological samples that are thicker and more highly scattering than what was previously possible.

## Funding

STROBE: A National Science Foundation Science Technology Center (DMR 1548924); Gordon and Betty Moore Foundation (GBMF4562); Office of Naval Research (N00014-17-1-2401); David and Lucile Packard Foundation; Chan Zuckerberg Initiative; Kirschstein National Research Service Award (F32GM129966).

## Acknowledgment

The authors would like to thank Nicole Repina, Fan Wu, and Angus Lee for the preparation of biological samples.

## Disclosures

The authors declare that there are no conflicts of interest related to this paper.

See Supplement 1 and Visualization 3 for supporting content.

## REFERENCES

**1. **Y. Sung, W. Choi, C. Fang-Yen, K. Badizadegan, R. R. Dasari, and M. S. Feld, “Optical diffraction tomography
for high resolution live cell imaging,” Opt.
Express **17**,
266–277 (2009). [CrossRef]

**2. **Y. Cotte, F. Toy, P. Jourdain, N. Pavillon, D. Boss, P. Magistretti, P. Marquet, and C. Depeursinge, “Marker-free phase
nanoscopy,” Nat. Photonics. **7**, 113–117
(2013). [CrossRef]

**3. **K. Kim, Z. Yaqoob, K. Lee, J. W. Kang, Y. Choi, P. Hosseini, P. T. C. So, and Y. Park, “Diffraction optical tomography
using a quantitative phase imaging unit,” Opt.
Lett. **39**,
6935–6938 (2014). [CrossRef]

**4. **S. Shin, K. Kim, J. Yoon, and Y. Park, “Active illumination using a
digital micromirror device for quantitative phase
imaging,” Opt. Lett. **40**, 5407–5410
(2015). [CrossRef]

**5. **T. Kim, R. Zhou, M. Mir, S. D. Babacan, P. S. Carney, L. L. Goddard, and G. Popescu, “White-light diffraction
tomography of unlabelled live cells,” Nat.
Photonics. **8**,
256–263 (2014). [CrossRef]

**6. **S. Chowdhury, W. J. Eldridge, A. Wax, and J. Izatt, “Refractive index tomography
with structured illumination,” Optica **4**, 537–545
(2017). [CrossRef]

**7. **D. Jin, R. Zhou, Z. Yaqoob, and P. T. C. So, “Dynamic spatial filtering
using a digital micromirror device for high-speed optical diffraction
tomography,” Opt. Express **26**, 428–437
(2018). [CrossRef]

**8. **G. Gbur and E. Wolf, “Diffraction tomography without
phase information,” Opt. Lett. **27**, 1890–1892
(2002). [CrossRef]

**9. **C. Zuo, J. Sun, J. Zhang, Y. Hu, and Q. Chen, “Lensless phase
microscopy and diffraction tomography with multi-angle
and multi-wavelength illuminations using a LED
matrix,” Opt. Express **23**, 14314–14328
(2015). [CrossRef]

**10. **M. Chen, L. Tian, and L. Waller, “3D differential phase contrast
microscopy,” Biomed. Opt.
Express **7**,
3940–3950 (2016). [CrossRef]

**11. **J. M. Soto, J. A. Rodrigo, and T. Alieva, “Label-free quantitative 3D
tomographic imaging for partially coherent light
microscopy,” Opt. Express **25**, 15699–15712
(2017). [CrossRef]

**12. **R. Horstmeyer, J. Chung, X. Ou, G. Zheng, and C. Yang, “Diffraction tomography with
Fourier ptychography,” Optica **3**, 827–835
(2016). [CrossRef]

**13. **J. Lim, K. Lee, K. H. Jin, S. Shin, S. Lee, Y. Park, and J. C. Ye, “Comparative study of iterative
reconstruction algorithms for missing cone problems in optical
diffraction tomography,” Opt. Express **23**, 16933–16948
(2015). [CrossRef]

**14. **M. H. Jenkins and T. K. Gaylord, “Three-dimensional quantitative
phase imaging via tomographic deconvolution phase
microscopy,” Appl. Opt. **54**, 9213–9227
(2015). [CrossRef]

**15. **L. Tian and L. Waller, “3D intensity and phase imaging
from light field measurements in an LED array
microscope,” Optica **2**, 104–111
(2015). [CrossRef]

**16. **U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “Learning approach to optical
tomography,” Optica **2**, 517–522
(2015). [CrossRef]

**17. **S. Chowdhury, M. Chen, R. Eckert, D. Ren, F. Wu, N. Repina, and L. Waller, “High-resolution 3D refractive
index microscopy of multiple-scattering samples from intensity
images,” Optica **6**,
1211–1219 (2019). [CrossRef]

**18. **P. Bon, B. Wattellier, and S. Monneret, “Modeling quantitative phase
image formation under tilted illuminations,”
Opt. Lett. **37**,
1718–1720 (2012). [CrossRef]

**19. **X. Ma, W. Xiao, and F. Pan, “Optical tomographic
reconstruction based on multi-slice wave propagation
method,” Opt. Express **25**, 22595–22607
(2017). [CrossRef]

**20. **S. S. Kou and C. J. Sheppard, “Image formation in holographic
tomography: high-aperture imaging conditions,”
Appl. Opt. **48**,
H168–H175 (2009). [CrossRef]

**21. **M. Born and E. Wolf, *Principles of Optics: Electromagnetic
Theory of Propagation, Interference and Diffraction of Light*,
7th ed. (Cambridge
University, 1999).

**22. **U. S. Kamilov, D. Liu, H. Mansour, and P. T. Boufounos, “A recursive born approach to
nonlinear inverse scattering,” IEEE Signal
Process. Lett. **23**,
1052–1056 (2016). [CrossRef]

**23. **W. Tahir, U. S. Kamilov, and L. Tian, “Holographic particle
localization under multiple scattering,” Adv.
Photon. **1**,
1–12 (2019). [CrossRef]

**24. **A. Dubois, K. Belkebir, and M. Saillard, “Retrieval of inhomogeneous
targets from experimental frequency diversity data,”
Inverse Probl. **21**,
S65–S79 (2005). [CrossRef]

**25. **P. C. Chaumet and K. Belkebir, “Three-dimensional
reconstruction from real data using a conjugate gradient-coupled
dipole method,” Inverse Probl. **25**, 024003 (2009). [CrossRef]

**26. **E. Mudry, P. C. Chaumet, K. Belkebir, and A. Sentenac, “Electromagnetic wave imaging
of three-dimensional targets using a hybrid iterative inversion
method,” Inverse Probl. **28**, 065007 (2012). [CrossRef]

**27. **H.-Y. Liu, D. Liu, H. Mansour, P. T. Boufounos, L. Waller, and U. S. Kamilov, “SEAGLE: sparsity-driven image
reconstruction under multiple scattering,”
IEEE Trans. Comput. Imag. **4**,
73–86 (2018). [CrossRef]

**28. **E. Soubies, T.-A. Pham, and M. Unser, “Efficient inversion of
multiple-scattering model for optical diffraction
tomography,” Opt. Express **25**, 21786–21800
(2017). [CrossRef]

**29. **Y. Ma, H. Mansour, D. Liu, P. T. Boufounos, and U. S. Kamilov, “Accelerated image
reconstruction for nonlinear diffractive imaging,” in
*IEEE International Conference on Acoustics, Speech and Signal
Processing (ICASSP)* (2018),
p. 6473–6477.

**30. **T.-A. Pham, E. Soubies, A. Goy, J. Lim, F. Soulez, D. Psaltis, and M. Unser, “Versatile reconstruction
framework for diffraction tomography with intensity measurements and
multiple scattering,” Opt. Express **26**, 2749–2763
(2018). [CrossRef]

**31. **T. Zhang, C. Godavarthi, P. C. Chaumet, G. Maire, H. Giovannini, A. Talneau, M. Allain, K. Belkebir, and A. Sentenac, “Far-field diffraction
microscopy at λ/10 resolution,”
Optica **3**,
609–612 (2016). [CrossRef]

**32. **G. Zheng, R. Horstmeyer, and C. Yang, “Wide-field, high-resolution
Fourier ptychographic microscopy,” Nat.
Photonics **7**,
739–745 (2013). [CrossRef]

**33. **Z. Liu, L. Tian, S. Liu, and L. Waller, “Real-time brightfield,
darkfield, and phase contrast imaging in a light-emitting diode array
microscope,” J. Biomed. Opt. **19**, 106002 (2014). [CrossRef]

**34. **E. Wolf, “Three-dimensional structure
determination of semi-transparent objects from holographic
data,” Opt. Commun. **1**, 153–156
(1969). [CrossRef]

**35. **H. N. Chapman, A. Barty, S. Marchesini, A. Noy, S. P. Hau-Riege, C. Cui, M. R. Howells, R. Rosen, H. He, J. C. Spence, U. Weierstall, T. Beetz, C. Jacobsen, and D. Shapiro, “High-resolution ab initio
three-dimensional x-ray diffraction microscopy,”
J. Opt. Soc. Am. A **23**,
1179–1200 (2006). [CrossRef]

**36. **M. Dierolf, A. Menzel, P. Thibault, P. Schneider, C. M. Kewish, R. Wepf, O. Bunk, and F. Pfeiffer, “Ptychographic x-ray computed
tomography at the nanoscale,” Nature **467**, 436–439
(2010). [CrossRef]

**37. **M. Habaza, B. Giboa, Y. Roichman, and N. T. Shaked, “Tomographic phase microscopy
with 180° rotation of live cells in suspension by holographic optical
tweezers,” Opt. Lett. **40**, 1881–1884
(2015). [CrossRef]

**38. **K. Lee, K. Kim, G. Kim, S. Shin, and Y. Park, “Time-multiplexed structured
illumination using a DMD for optical diffraction
tomography,” Opt. Lett. **42**, 999–1002
(2017). [CrossRef]

**39. **H. Erdogan and J. A. Fessler, “Ordered subsets algorithms for
transmission tomography,” Phys. Med.
Biol. **44**,
2835–2851 (1999). [CrossRef]

**40. **A. Beck and M. Teboulle, “Fast gradient-based algorithms
for constrained total variation image denoising and deblurring
problems,” IEEE Trans. Image Process. **18**, 2419–2434
(2009). [CrossRef]

**41. **Y. Nesterov, “A method of solving a convex
programming problem with convergence rate O
(1/k^{2}),” Sov. Math.
Doklady **27**,
372–376
(1983).

**42. **B. O’Donoghue and E. Candes, “Adaptive restart for
accelerated gradient schemes,” Found. Comput.
Math. **15**,
715–732
(2015). [CrossRef]

**43. **L. Tian and L. Waller, “Quantitative differential
phase contrast imaging in an LED array microscope,”
Opt. Express **23**,
11394–11403 (2015). [CrossRef]

**44. **M. E. Kandel, M. Fanous, C. Best-Popescu, and G. Popescu, “Real-time halo correction in
phase contrast imaging,” Biomed. Opt.
Express **9**,
623–635 (2018). [CrossRef]

**45. **S. B. Mehta and C. J. Sheppard, “Quantit. phase-gradient
imaging at high resolution with asymmetric illumination-based
differential phase contrast,” Opt.
Lett. **34**,
1924–1926 (2009). [CrossRef]

**46. **R. Eckert, Z. F. Phillips, and L. Waller, “Efficient illumination angle
self-calibration in Fourier ptychography,”
Appl. Opt. **57**,
5434–5442 (2018). [CrossRef]