Multi-layer Born multiple-scattering model for 3D phase microscopy

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 theOSAOpen Access Publishing Agreement


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][2][3][4] or pupil plane modulation [5][6][7]. Intensitybased 3D phase methods use through-focus or color images from off-the-shelf microscopes [8][9][10][11][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][16][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 backscattering, 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][22][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].

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 = 2π λ , λ is the wavelength of incident light in vacuum, and n 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 ) = − exp(ik o n b r) 4π r , where r = x 2 + y 2 + z 2 [21]. For an incident field, U in , the total field, U tot , can be written as This linear assumption holds when the magnitude of V is small (i.e., a weakly scattering object). For the MS method, the first Born model is equivalent to assuming that each slice of the sample interacts only with the unscattered incident beam, without accounting for any scattering that occurred before light reached each slice. Hence, Eq. (2) does not account for multiple scattering. 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, z (see Fig. 1). Therefore, Eq. (2) describes the field after the wave propagates through the n th layer. This serves as the incident field of the (n + 1) 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 n s denote the incident field and the scattered field, respectively, within the n th layer. For mathematical simplicity, we define the n th layer as the one that occupies the space z ∈ [(n − 1 2 ) z, (n + 1 2 ) z]. The recursive formula can be written as where ρ = (x , y ) is the 2D vector in real space coordinates; the out-of-domain functions, U n and U n s with their z > (n + 1 2 ) z, are convenient expressions and should be understood as those at z = (n + 1 2 ) z propagated in uniform background medium, and Intensity measurements with spatially coherent illumination from different angles are captured on an optical microscope and fed in as inputs to our 3D phase tomography algorithm. By solving a nonlinear optimization problem with our MLB scattering forward model, the 3D refractive index of a multiple-scattering object is recovered.
The first term on the right-hand side of Eq. (3) is the incident field at the center of the n th layer propagated by one slab thickness z. Given the 2D Fourier transform operator F{·}, the spatial frequency spectrum of the incident field U n (u) = F{U n (ρ, n z)}, and the angular spectrum propagation kernel P (u, z) = e i2π z √ (n b /λ) 2 − u 2 with the 2D spatial frequency space coordinates vector u, we get To simplify Eq. (4), we find it useful to use the 2D Fourier spectrum of the Green's function [34]: where The assumption made below is that the scattering potential does not vary axially within each layer, which is generally valid when the thickness z is small (i.e., V (ρ, n z + ζ ) = V n (ρ) for ζ ∈ [− z/2, z/2]). Based on Eqs. (4)-(6) and the convolution theorem, the spatial frequency spectrum of the scattered field, whereṼ n (u) = F{V n (ρ)} denotes the 2D spatial frequency spectrum of the scattering potential at the n th layer. The integration along the ζ dimension involves only two propagation kernels: where sinc(x ) = sin(π x )/(π x ). Note that | z − ζ | = z − ζ in the integration interval. Plugging the integration result into Eq. (7), we obtainŨ Equations (3), (5), and (9) make up the forward scattering part of our MLB model. Equation (9) is computationally expensive; however, when 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 (u) = circ(uλ/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 intensitybased 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 L(V ): Here, N LED is the total number of LEDs used, which determines the number of illuminations. The first term of Eq. (10) is the data fidelity term, which gives a metric for how well the current object estimate fits the measurements (through the Euclidean distance between the estimated amplitudes and the square root of measured intensities I measure ). The second term is a regularization term that enforces prior knowledge on the 3D scattering potential. τ is a tunable parameter that balances the strength of data fidelity and regularization. Because the minimum of Eq. (10) cannot be computed analytically, an iterative algorithm is required. Due to the large-scale nature of 3D problems, we choose the proximal gradient method for minimizing Eq. (10), which has relatively lower memory requirements and computation complexity, as compared to the alternative direction method of multipliers (ADMM) or second-order Newton's methods. To implement the inverse algorithm, the gradient of the data fidelity term with respect to V , {V n grad (ρ)} N n=1 , is required at each iteration (see Algorithm 2), where {·} * stands for the complex conjugate operation, and R, Q, U † image , and U n, † are the auxiliary variables for gradient computation. Typically, the algorithm converges faster if V is sequentially updated over different angles of incidence [15,16,39], as opposed to summing the gradients computed from all measurements and refining V . The penalty function helps regularize the updated V to avoid over-fitting due to the nonlinearity of the model. The following proximal operator is applied after the gradient steps: 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

Algorithm 2. Gradient Computation
Input: Measured intensity I measure , predicted total field on the image plane U image , incident fields at each layer {U n (ρ, n z)} N n=1 , current estimated 3D scattering potential {V n (ρ)} N n=1 , layer thickness z, and the refocusing distance f . for j ← 1 to N LED do sequential gradient descent 3: : t k ← 1 10: Start k + 1 iteration 11: end if 12: : 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 α is set within a range from 0.5-5, and τ 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, α is manually tuned to achieve the best rate of convergence, and τ is adjusted to mitigate most of the noise artifacts without degrading physical structures.

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 × 15 × 7.5 µ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 × 450 × 150 voxels with 0.05 µm resolution. Assuming an illumination wavelength of 532 nm and objective NA obj = 0.8, we calculate the amplitude images using forward and backward scattered light under on-axis (NA illu = 0.00) and off-axis (NA 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 ∼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 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π 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 (NA illu ≈ 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 NA illu close to NA 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.

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 × 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 × 0.65NA objective, and 46 mm for the 60× 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 µm pixel size) at the front port of the microscope, where an additional 2× 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 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 µm polystyrene bead (Sigma-Aldrich, n = 1.6010 at λ = 514 nm ) immersed in two different index-matching oils (Cargille), resulting in RI contrasts of n = 0.0113 and n = 0.0452. We used a 40 × 0.65 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 × 0.80 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 highfrequency 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× in each dimension. Figure 5(a) shows orthonormal cross sections of the recovered 3D RI (after halo-correction [44]). The reconstruction volume contains 612 × 612 × 120 voxels with voxel size 0.108 × 0.108 × 0.540 µ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 ∼17 µ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 (∼0.64 µ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 fibercoupled LED (λ = 530 nm center wavelength) is combined with a 2D scanning mirror to achieve programmable angle-scanning illumination. Two high-NA oil-immersion objective lenses (effective 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 × 1200 × 100 voxels, with voxel size 0.12 × 0.12 × 0.35 µm 3 . The full length of the C. elegans (∼1 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.

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.