## Abstract

We develop a novel algorithm for large-scale holographic reconstruction of 3D particle fields. Our method is based on a multiple-scattering beam propagation method (BPM) combined with sparse regularization that enables recovering dense 3D particles of high refractive index contrast from a single hologram. We show that the BPM-computed hologram generates intensity statistics closely matching with the experimental measurements and provides up to 9× higher accuracy than the single-scattering model. To solve the inverse problem, we devise a computationally efficient algorithm, which reduces the computation time by two orders of magnitude as compared to the state-of-the-art multiple-scattering based technique. We demonstrate the superior reconstruction accuracy in both simulations and experiments under different scattering strengths. We show that the BPM reconstruction significantly outperforms the single-scattering method in particular for deep imaging depths and high particle densities.

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

## 1. Introduction

Single-shot holographic particle 3D imaging is fundamental to many important applications, such as flow cytometry [1], biological sample characterization [2,3], and flow measurement [4]. The technique works by first recording a 2D intensity measurement from a particle volume under coherent illumination and then estimating the 3D refractive index distribution. The problem is challenging because of the dimensional mismatch between the single 2D measurement and the unknown 3D object and the “phaseless” measurement, both of which make the inverse problem severely ill-posed, especially for large-scale problems. Furthermore, multiple scattering effects become significant as the particle density and the refractive index contrast increase, which necessitates a nonlinear forward model to accurately describe the image formation process.

Traditional holographic particle 3D reconstruction is based on the linear single-scattering model derived from the first Born approximation method (FBM) [5]. The most widely used algorithm is known as the “back-propagation method” [4], which is equivalent to apply the pseudo-inverse of the (3D-to-2D) single-scattering forward operator to the captured hologram. To improve the reconstruction accuracy, compressive holography [6] has been developed that supplements the FBM forward model with a sparsity constraint and solves the 3D reconstruction problem by an iterative optimization procedure. This approach is particularly effective in alleviating the twin-image artifacts arising [6] and improving the robustness to unaccounted multiple-scattering effects at high scattering densities [7]. However, these methods are fundamentally limited by the underlying single-scattering assumption, which is valid only for weakly scattering objects. The model gradually breaks down as the particle density and/or the refractive index contrast increase. The multiple-scattering effects result in a model mismatch, which in turn limits the reconstruction accuracy by these techniques.

Several multiple-scattering models have been developed, such as iterative Born series [8,9], contrast source inversion [10], discrete-dipole approximation [11–14], and series expansion with accelerated gradient descent on Lippmann–Schwinger equation [15]. However, computational challenges typically restrict these methods to be used only for small-scale problems. This is because they involve iteratively estimating the internally scattered fields within the object volume for modeling the multiple-scattering process, which incurs a computation cost that grows rapidly as the object size, and thus is not ideal for our intended applications containing on the order of 100 million voxels in the object volume. Recently, the multi-slice beam propagation method (BPM) [16–19] has emerged as an attractive, computationally efficient multiple-scattering model. The utility of the BPM has been demonstrated for reconstructing 3D refractive index distributions using tomographic measurements from multiple interferometric complex-field measurements [16] or intensity-only measurements [17–19]. Our goal here is to critically examine the utility of the BPM to handle the inversion from a single in-line hologram, which inherently suffers from more severe missing-cone artifacts [20,21] and greater dimensional mismatch as compared to the previous tomography studies.

We demonstrate a novel 3D reconstruction algorithm that combines a BPM multiple-scattering forward model and a sparsity prior. The accuracy of our forward model is quantified by comparing the simulation with the experiments. We show that the BPM-computed intensity statistics closely match with the experimental data at different scattering densities, and provide up to 9$\times$ higher accuracy in predicting the hologram contrast than that from the FBM. Benefited from the high computational efficiency, the BPM reduces the computational complexity by $N_z$ times (where $N_z$ is the number of axial slices of the object volume), as compared to the state-of-the-art multiple-scattering model-based holographic particle reconstruction method [9]. We demonstrate our proposed method’s superior 3D imaging performance in both simulation and experiments under different scattering densities and refractive index contrast. To quantify the improvement by the multiple-scattering model over the single-scattering model, we compare our method with the compressive holography method. We show that the BPM reconstruction significantly outperforms the FBM in particular at deep imaging depths and for high particle densities.

## 2. Method

Our reconstruction algorithm solves an $\ell _1$-regularized least-squares problem, in which the data fidelity term measures the difference between the BPM-estimated and captured holograms. Below, we describe our forward model and the reconstruction algorithm.

#### 2.1 Forward model

Our imaging geometry is shown in Fig. 1(a). A plane wave passes through a particle volume and the interference pattern resulting from the scattered fields are captured as the hologram. This multiple-scattering image-formation process is modeled as,

where $\hat {\textbf {I}}\in \mathbb {R}^M$ denotes the vectorized, BPM computed intensity hologram containing $M$ pixels. $\textbf {w}\in \mathbb {R}^N = \textbf {n}-n_0$ represents the vectorized 3D refractive index contrast distribution that contains $N$ voxels, and is defined by the difference between the particle’s index distribution $\textbf {n}$ and the constant background index $n_0$. The particle’s refractive index is assumed to be real-valued since the absorption is negligible. $\textbf {S}: \mathbb {R}^N\rightarrow \mathbb {C}^M$ is the BPM operator that maps the 3D refractive index to the 2D complex field at the sensor plane.The BPM is computed by a recursive sequence, as illustrated in Fig. 1(b). It approximates the 3D object as $N_z$ infinitesimally thin axial planar slices along the optical axis, each modeled as a 2D phase mask separated equally by a uniform medium. The field is then computed by a sequence of diffraction and refraction operations. Specifically, the field exiting the $j\textrm {th}$ slice $\textbf {S}_j(\textbf {w})$ is related to the field from the $(j-1)\textrm {th}$ slice $\textbf {S}_{j-1}(\textbf {w})$ by

The implementation of the BPM requires an initial condition, for which we set the initial field as a constant-valued on-axis plane wave $\textbf {S}_0(\textbf {w}) = \mathbf {1}$.

The diffraction operator for a propagation distance $\Delta z$ is denoted by $\textbf {H}$ and computed by $\textbf {H} = \textbf {F}^H\textrm {diag}(\textbf {h})\textbf {F}$, where $\textbf {F}$ and $\textbf {F}^H$ are the 2D discrete Fourier and inverse Fourier transform matrices respectively, $\large {\cdot }^H$ denotes matrix Hermitian transpose, $\textrm {diag}(\textbf {h})$ is a diagonal matrix formed by the vector $\textbf {h}$ as the main diagonal and is implemented by element-wise multiplication, and $\textbf {h}$ represents the vectorized 2D angular spectrum transfer function $h(p,q)$:

where $p,q$ index the 2D wave vector in the $k$-space, $k=k_0n_0$ is the wave-number in the background medium, $k_0={2\pi }/{\lambda _0}$, $\lambda _0$ is the wavelength in vacuum, and $\Delta k$ is the sampling interval in the $k$-space. $\mathcal {P}$ represents the low-pass filter due to the evanescent cutoff.The refraction operator is implemented by element-wise multiplications between the diffracted field and the accumulated phase $\textbf {p}_j(\textbf {w}_j)=e^{jk_0\Delta z\textbf {w}_j}$, where $\textbf {w}_j$ denotes the discretized 2D refractive index map in the $j\textrm {th}$ slice.

The hologram is the intensity of the exit field from the last slice at $N_z$ low-pass filtered by the pupil function of the imaging system.

The computation of the BPM thus requires implementing Eq. (2) $N_z$ times. The computational complexity of the BPM is $O(N\textrm {log}(N_x N_y))$, as set approximately by $2N_z$ times evaluations of 2D FFT, where the total number of object voxels is $N=N_xN_yN_z$, and $N_x$ and $N_y$ are the number of pixels in the $x$ and $y$ directions, respectively. For comparison, the computational complexity of the Born series model is $O(N_zN\textrm {log}(N_x N_y))$ [9]. Thus, the BPM reduces the computational complexity by $N_z$ times, which becomes significant when the object depth is large.

To illustrate the multiple scattering effects modeled by the BPM, we show an example $yz$ cross-section of the intensity distribution computed from a particle volume in Fig. 1(b). In the zoom-in region, complex interference patterns are visible due to multiple particles in the close proximity, which introduces strong multiple scattering effects [22].

#### 2.2 Inverse problem

### 2.2.1 Problem formulation

We formulate the 3D reconstruction as a minimization problem:

The regularization term is the $\ell _1$ norm of the refractive index distribution

which promotes the sparsity of the reconstructed object.The minimization Eq. (4) is not a trivial task. The primary difficulty stems from that the data fidelity term $\mathcal {D}$ involves a nonlinear forward operator $\textbf {S}$ and the regularization term $\mathcal {R}$ is non-smooth. We next present a novel algorithm based on the proximal-gradient descend technique, which extends [16] to intensity-only measurement.

### 2.2.2 Computation of the gradient

The crucial step is the gradient computation for $\mathcal {D}$ with respect to $\textbf {w}$, as summarized in Algorithm 1 and briefly explained below. Additional details are provided in Supplement 1.

First, we take the derivative of $\mathcal {D}(\textbf {w})$ with respect to $\textbf {w}$ and rearrange it into a column vector

To derive a tractable algorithm for computing Eq. (8), we apply the recursive BPM formula in Eq. (2) and compute the local gradient at the $j$th slice:

Since the input plane wave does not depend on $\textbf {w}$, so for the initial condition $j=0$, we have

Based on Eq. (9) and the initial condition in Eq. (10), we obtain a practical implementation of Eq. (8) to calculate the gradient of the data fidelity term, as summarized in Algorithm 1.

Intuitively, the gradient computation is similar to the “error backpropagation” algorithm used in deep neural networks [16]. The algorithm iterates between two major steps, including update the intermediate field gradient slice-by-slice (first term in Eq. (9)) and backpropagate the slice-wise field residual (second term in Eq. (9)). Since this gradient computation takes the same recursive procedure as the forward BPM model, its computation complexity is also $O(N\textrm {log}(N_x N_y))$.

### 2.2.3 Reconstruction algorithm

Our algorithm is summarized in Algorithm 2, that reconstructs the refractive index $\textbf {w}$ from a single hologram based on the proximal gradient algorithm. Conceptually, this algorithm is similar to the fast iterative shrinkage/thresholding algorithm (FISTA) [23], which is widely used to minimize objective function that consist of the sums between a smooth and a non-smooth term.

A major component of this algorithm is the proximal operator for the $\ell _1$-regularizer

The parameters in Algorithm 2 are set as follows. We set the initial guess $\hat {\textbf {w}}^0 = \mathbf {0}$, and the step size $\gamma =5\times 10^{-6}$. The stopping criterion is the maximum iteration number to be 200-300. The regularization parameter $\tau$ is tuned under different scattering conditions. When scattering is stronger, $\tau$ is set larger. The reconstruction is found to be insensitive to the fine-tuning of $\tau$.

## 3. Evaluation of BPM forward model accuracy in large scale

#### 3.1 Intensity statistics analysis

We first validate the forward model accuracy by comparing the BPM-computed holograms with experimental measurements. In practice, we assess the accuracy based on analyzing the intensity statistics under different scattering conditions, since the ground-truth particle positions in the experiments are not known.

We perform simulations using parameters that match with the experiments. More details about the experiments are provided in Supplement 1. Holograms are computed at four particle densities $\rho$, including $1.60, 3.20, 6.41, 12.82\times 10^{4} /\mathrm{\mu} \textrm {L}$, corresponding to $250, 500, 1000, 2000$ particles in a $176.64 \times 176.64 \times 500~\mathrm{\mu} \textrm {m}^3$ volume. $1~\mathrm{\mu} \textrm {m}$ particles are randomly placed in 3D positions. The refractive index of the particle $n$ and the background medium (water) $n_0$ is 1.59 and 1.33, respectively, and the contrast is $\Delta n=0.26$. To simulate 3D refractive index contrast distribution $\textbf {w}$, the voxels inside the particles are assigned with the constant $\Delta n$, and the rest of the background voxels are assigned with zero. The lateral sampling size $\Delta x, \Delta y$ are both $0.1725\mathrm{\mu} \textrm {m}$, and the axial sampling size $\Delta z$ is $\lambda /16=0.0297\mathrm{\mu} \textrm {m}$, where $\lambda =\lambda _0/n_0$ and $\lambda _0=0.632\mathrm{\mu} \textrm {m}$. The resulting object size in our forward model is $1024\times 1024\times 16840$ voxels. Example computed holograms at four particle densities are shown in Fig. 2(a). As expected, characteristic fringe patterns from individual particles are still visible at the lowest density case. As the density increases, the holograms gradually become partially developed speckle patterns.

First, we analyze the BPM’s accuracy by comparing the histograms of the computed hologram with the captured hologram at each particle density. As shown in Fig. 2(b), the histograms match well at all four densities.

Next, we assess the BPM’s accuracy in the spatial frequency domain. We calculate the normalized spectra of the computed and captured holograms. As shown in Fig. 2(c), the frequency components of the computed holograms closely match with the experimental measurements within the NA of the imaging system.

Finally, we compare the accuracy of the BPM and FBM based on the hologram contrast $\mathcal {C} = \frac {\sigma }{\mathrm{\mu} }$, where $\sigma$ and $\mathrm{\mu}$ are the standard deviation and mean of the hologram, respectively. Briefly, the FBM is a linear model that assumes a weakly scattering object, and the resulting scattered field is linearly related to the scattering potential $V(x,y,z) = \frac {1}{4\pi }k_0^2[n(x,y,z)^2-n_0^2]$. Given the plane-wave incident field $U_{in}(x,y,z) = e^{in_0k_0z}$, the total field $U_{\textrm {FBM}}$ can be approximated as

As shown in Fig. 2(d), the contrast from the BPM-computed holograms agree well with the experiments. The BPM slightly under-estimates the contrast. A possible reason is that the BPM approximates the multiple forward scattering by computing the complex field slice-by-slice but ignores backscattering. As a result, the BPM may slightly under-estimates the high frequency information, which reduces the contrast in the calculated holograms. As a comparison, the contrast from the FBM-computed holograms are consistently higher than the experimental data. The contrast discrepancy increases as the particle density. At the highest density, the discrepancy for the BPM is 0.0048, whereas the FBM is 0.0422, representing a $9\times$ improvement by the BPM.

Overall, these studies show that the BPM can accurately model multiple scattering in a dense 3D particle field and significantly outperforms the single-scattering model.

#### 3.2 Sampling distance $\Delta z$ and scattering strength effect in a BPM forward model

The BPM accuracy is primarily influenced by the axial sampling distance $\Delta z$ and the scattering strength of the 3D object. In the following, we quantitatively evaluate the model accuracy under different axial sampling and scattering conditions (by changing the particle density $\rho$ and refractive index contrast $\Delta n$), while fixing the lateral sampling $\Delta x = \Delta y = 0.1725\mathrm{\mu} \textrm {m}$.

In Fig. 3(a), the accuracy of BPM is plotted for different $\Delta z$ under the same scattering condition. We compare the holograms computed with different $\Delta z$ with the reference $\hat {\textbf {I}}_0$ using $\Delta z=\lambda /16$. The difference is quantified by the signal-to-noise ratio (SNR), $\textrm {SNR}=10\log _{10}\frac {\|\hat {\textbf {I}}_0\|_2}{\|\hat {\textbf {I}}_0-\mathbf {\hat {I}}\|_2}$, where $\hat {\textbf {I}}$ is the computed hologram with axial down-sampling. The accuracy of the BPM drops rapidly when $\Delta z$ is smaller than the particle diameter ($1\mathrm{\mu} \textrm {m}$) and reduces slowly as $\Delta z$ further increases. This indicates that to accurately compute the inner-particle multiple scattering using the BPM, dense axial sampling is generally needed. Nevertheless, even under coarse axial sampling up to $\Delta z=16\lambda$, the BPM is still more accurate than the FBM computed with the reference dense axial sampling $\Delta z=\lambda /16$.

Next, we study the BPM’s accuracy for different $\Delta z$ for different particle densities $\rho$ and refractive index contrasts $\Delta n$ (Fig. 3(b) and (c)). Our study shows that the shape of the curve remain the same for different scattering conditions, which suggests that it is only determined by the sampling distance $\Delta z$.

Next, we study how the scattering strength affects the BPM’s accuracy for a fixed $\Delta z=\lambda /4$. We compute holograms at different particle densities $\rho$, refractive index contrasts $\Delta n$, and volume thicknesses $Z$. We then compare these holograms with the corresponding reference ($\Delta z=\lambda /16$). The results are summarized in Fig. 3(d)-(f). In general, the accuracy decreases as the scattering becomes stronger. The SNR is found to satisfy the following scaling law:

## 4. 3D imaging results

In this section, we report 3D particle imaging results in both simulation and experiments.

#### 4.1 Simulation results

Given the high accuracy of the BPM model validated in Sec. 3., we next quantify the accuracy of our reconstruction algorithm in simulation. To implement the reconstruction algorithm, we first select a suitable $\Delta z$. Due to memory limitations, it is not feasible to perform reconstructions using the same dense $\Delta z=\lambda /16$ as the forward simulation. Based on our quantitative study in Fig. 3, we select $\Delta z=6.2\lambda$ for all the reconstructions, which balances the model accuracy and computational cost. The corresponding object volume contains around 177 million voxels.

In Fig. 4(a), we simulated a hologram from a particle field ($\rho =6.41\times 10^4/\mathrm{\mu} \textrm {L}, \Delta n =0.26$, particle diameter:$1\mathrm{\mu} \textrm {m}$) using the densely sampled BPM model ($\Delta z = \lambda /16$). The 3D reconstruction result is visualized in Fig. 4(b). For better visualization, we show the $x$-projection of the central $69\times 69\times 500\mathrm{\mu} \textrm {m}^3$ region in Fig. 4(c). As expected, the reconstructed particles are elongated along the $z$-axis due to the missing cone problem. To further evaluate the reconstruction, we extracted the centroids of the reconstructed particles and then overlay them onto the ground-truth. The $z$- and $x$-projections are shown in Fig. 4(d) and (e). As shown in the $z$-projection, the reconstruction errors mostly occurred in the peripheral FOV region. Further inspecting the the $x$-projection, the localized centroids agree well with the ground truth.

To further quantify the reconstruction accuracy, we repeat the simulation at seven different refractive index contrasts and four different particle densities. We then use Jaccard index (JI), lateral root mean square error (RMSE), and axial RMSE for quantification [24]:

As shown in Fig. 5(a), our method achieves JI > 0.9 for imaging particles at densities $\rho \le 6.41\times 10^{4} /\mathrm{\mu} \textrm {L}$ (1000 particles in the $176.74\times 176.74\times 500 \mathrm{\mu} \textrm {m}^3$ volume) with $\Delta n = 0.26$. Our method also achieves localization accuracy better than the diffraction limit. The lateral RMSE is less than 0.25$\mathrm{\mu} \textrm {m}$ and the axial RMSE is less than 3.5$\mathrm{\mu} \textrm {m}$ in all cases. Figure 5(b) shows representative $z$-projections of the centorids overlaid onto the ground-truth for the central $69\times 69\times 500\mathrm{\mu} \textrm {m}^3$ region. These results show that our method performs well for particle densities as high as $6.41\times 10^4 /\mathrm{\mu} \textrm {L}$ and refractive index contrast as high as 0.32.

Next, we quantify the improvement of our multiple-scattering BPM model as compared to the single-scattering FBM model. To do so, we performed reconstructions using the compressive holography algorithm that uses the FBM forward model and a total variation (TV) regularization to enforce sparsity [6,7]. In Fig. 6, we compare results across four different particle densities with $\Delta n = 0.26$. To quantify the depth-dependent reconstruction accuracy, we calculated JI by dividing the entire depth into 17 segments and then calculating JI for each segment. The depth-wise JIs under different scattering densities are shown in Fig. 6(a). The results show that the BPM-based algorithm consistently detects the particles across all the depths for particle densities $\rho \le 6.41\times 10^{4} /\mathrm{\mu} \textrm {L}$. In contrast, the FBM-based algorithm suffers from rapid degradation as the depth increases. In particular, when the depth is around $500 \mathrm{\mu} \textrm {m}$, the JI of BPM is $>34$ times higher than the FBM when $\rho \le 6.41\times 10^{4} /\mathrm{\mu} \textrm {L}$. Representative $z$-projections of the reconstructed centroids overlaid onto the ground truth across different depths are shown in Fig. 6(b) for the central $69\times 69\times 500 \mathrm{\mu} \textrm {m}^3$ region. These results highlight that our multiple-scattering algorithm significantly outperforms the traditional single-scattering method, in particular for reconstructing particles at deep depths.

#### 4.2 Experimental results

We demonstrate our reconstruction algorithm on experimentally captured holograms at 4 different densities. In Fig. 7, 3D visualizations of example reconstruction results are shown in depth-color coded 3D renderings and $x$- and $z$-projections. The axial elongations are visible in all cases, which are resulted from the missing spatial frequency information limited by the single-view holographic measurement. We further observe that more particles are reconstructed at the depths closer to the hologram plane. This is expected since, for particles closer to the hologram plane, a greater amount of scattered information are captured by the finite-sized image sensor.

Finally, we compare our method with the FBM algorithm in Fig. 8. Similar to our observations in the simulation, the BPM algorithm can better reconstruct particles at deeper depths than the FBM algorithm, as highlighted by the $x$- and $z$-projections for different depth regions.

## 5. Conclusion

We have developed a new reconstruction algorithm for large-scale holographic particle 3D imaging based on a multiple-scattering beam propagation model. Our forward model demonstrates superior accuracy for multiple-scattering particle fields as compared to the traditional first Born-based single scattering model. The computational efficiency of our iterative algorithm allows reducing the computational complexity by more than 2 orders of magnitude for reconstructing volumes containing morn than 100 million voxels, as compared to the iterative Born series based method. The accuracy of our proposed algorithm is demonstrated in both simulation and experiment. In particular, we show that our method is particularly effective to improve the imaging performance for particles at deep depths. Together, these advances may open up new exciting opportunities for large-scale holograph particle 3D imaging in various applications.

Though our algorithm achieved start-of-the-art performance, it is still limited by the severe missing cone problem, which elongates the reconstructed particles. A promising future direction is to combine our multiple-scattering model and advanced deep-learning priors to further improve the imaging performance [25,26].

## Funding

National Science Foundation (1813848, 1846784).

## Acknowledgments

We thank Boston University Shared Computing Cluster for the computing resources.

## Disclosures

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

## Data availability

Data underlying the results presented in this paper are available in Ref. [27].

## Supplemental document

See Supplement 1 for supporting content.

## References

**1. **F. C. Cheong, B. S. R. Dreyfus, J. Amato-Grill, K. Xiao, L. Dixon, and D. G. Grier, “Flow visualization and flow cytometry with holographic video microscopy,” Opt. Express **17**(15), 13071–13079 (2009). [CrossRef]

**2. **I. Moon, M. Daneshpanah, B. Javidi, and A. Stern, “Automated three-dimensional identification and tracking of micro/nanobiological organisms by computational holographic microscopy,” Proc. IEEE **97**(6), 990–1010 (2009). [CrossRef]

**3. **T.-W. Su, L. Xue, and A. Ozcan, “High-throughput lensfree 3d tracking of human sperms reveals rare statistics of helical trajectories,” Proc. Natl. Acad. Sci. **109**(40), 16018–16022 (2012). [CrossRef]

**4. **L. Tian, N. Loomis, J. A. Domínguez-Caballero, and G. Barbastathis, “Quantitative measurement of size and three-dimensional position of fast-moving bubbles in air-water mixture flows using digital holography,” Appl. Opt. **49**(9), 1549–1554 (2010). [CrossRef]

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

**6. **D. J. Brady, K. Choi, D. L. Marks, R. Horisaki, and S. Lim, “Compressive holography,” Opt. Express **17**(15), 13040–13049 (2009). [CrossRef]

**7. **W. Chen, L. Tian, S. Rehman, Z. Zhang, H. P. Lee, and G. Barbastathis, “Empirical concentration bounds for compressive holographic bubble imaging based on a mie scattering model,” Opt. Express **23**(4), 4715–4725 (2015). [CrossRef]

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

**9. **W. Tahir, U. S. Kamilov, and L. Tian, “Holographic particle localization under multiple scattering,” Adv. Photonics **1**(03), 1 (2019). [CrossRef]

**10. **P. M. Van Den Berg and R. E. Kleinman, “A contrast source inversion method,” Inverse Probl. **13**(6), 1607–1620 (1997). [CrossRef]

**11. **B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for scattering calculations,” J. Opt. Soc. Am. A **11**(4), 1491–1499 (1994). [CrossRef]

**12. **K. D. Unger, P. C. Chaumet, G. Maire, A. Sentenac, and K. Belkebir, “Versatile inversion tool for phaseless optical diffraction tomography,” J. Opt. Soc. Am. A **36**(11), C1–C8 (2019). [CrossRef]

**13. **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**(6), 065007 (2012). [CrossRef]

**14. **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**(6), 609–612 (2016). [CrossRef]

**15. **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. Imaging **4**(1), 73–86 (2018). [CrossRef]

**16. **U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “Optical tomographic image reconstruction based on beam propagation and sparse regularization,” IEEE Trans. Comput. Imaging **2**(1), 59–70 (2016). [CrossRef]

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

**18. **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**(9), 1211–1219 (2019). [CrossRef]

**19. **M. Chen, D. Ren, H.-Y. Liu, S. Chowdhury, and L. Waller, “Multi-layer born multiple-scattering model for 3d phase microscopy,” Optica **7**(5), 394–403 (2020). [CrossRef]

**20. **K. Tam and V. Perez-Mendez, “Tomographical imaging with limited-angle input,” J. Opt. Soc. Am. **71**(5), 582–592 (1981). [CrossRef]

**21. **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**(13), 16933–16948 (2015). [CrossRef]

**22. **J. Lim, A. Goy, M. H. Shoreh, M. Unser, and D. Psaltis, “Learning tomography assessed using mie theory,” Phys. Rev. Appl. **9**(3), 034027 (2018). [CrossRef]

**23. **A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci. **2**(1), 183–202 (2009). [CrossRef]

**24. **D. Sage, T.-A. Pham, H. Babcock, T. Lukes, T. Pengo, J. Chao, R. Velmurugan, A. Herbert, A. Agrawal, S. Colabrese, A. Wheeler, A. Archetti, B. Rieger, R. Ober, G. M. Hagen, J.-B. Sibarita, J. Ries, R. Henriques, M. Unser, and S. Holden, “Super-resolution fight club: assessment of 2d and 3d single-molecule localization microscopy software,” Nat. Methods **16**(5), 387–395 (2019). [CrossRef]

**25. **Z. Wu, Y. Sun, A. Matlock, J. Liu, L. Tian, and U. S. Kamilov, “Simba: Scalable inversion in optical tomography using deep denoising priors,” IEEE J. Sel. Top. Signal Process. **14**(6), 1163–1175 (2020). [CrossRef]

**26. **A. Matlock and L. Tian, “Physical model simulator-trained neural network for computational 3d phase imaging of multiple-scattering samples,” arXiv preprint arXiv:2103.15795 (2021).

**27. **H. Wang, W. Tahir, J. Zhu, and L. Tian, “Large-Scale-3D-Holographic-Imaging-with-Beam-Propagation,” GitHub (2021), https://github.com/bu-cisl/Large-Scale-3D-Holographic-Imaging-with-Beam-Propagation.