Single-shot holographic 3D particle-localization under multiple scattering

We introduce a novel framework that incorporates multiple scattering for large-scale 3D particle-localization using single-shot in-line holography. Traditional holographic techniques rely on single-scattering models which become inaccurate under high particle-density. We demonstrate that by exploiting multiple-scattering, localization is significantly improved. Both forward and back-scattering are computed by our method under a tractable recursive framework, in which each recursion estimates the next higher-order field within the volume. The inverse scattering is presented as a nonlinear optimization that promotes sparsity, and can be implemented efficiently. We experimentally reconstruct 100 million object voxels from a single 1-megapixel hologram. Our work promises utilization of multiple scattering for versatile large-scale applications.


Introduction
Three dimensional (3D) particle-localization using in-line holography is fundamental to many applications, such as biological sample characterization, 1, 2 flow cytometry, 3,4 fluid mechanics, 5,6 and optical measurement. [7][8][9] Reconstructing dense samples, however, still remains challenging. 10 Standard back-propagation method (BPM) can only handle low particle density. 10 Compressive holography based on the first Born approximation significantly improves upon BPM by imposing sparsity constraints. 11,12 However, it suffers from an underlying single scattering assumption, which becomes invalid at high particle densities where multiple scattering effects become significant. In this work, we propose, for the first time to our knowledge, a framework that accounts for multiple scattering in in-line digital holography, and enables accurate 3D particle localization at high density in a computationally efficient fashion.
Multiple scattering induces a nonlinear relation between the permittivity contrast and the scattered field, making it difficult to invert. 13 Many algorithms have been proposed to solve the inverse multiple scattering problem and demonstrated improved performance over single-scattering methods, such as iterative Born series, [14][15][16][17][18] contrast source inversion, [19][20][21][22] modified gradient, 23,24 Series Expansion with Accelerated Gradient Descent on Lippmann-Schwinger Equation (SEAGLE), 25,26 and hybrid methods. [27][28][29][30] However, computational challenges restrict them to be demonstrated only for small-scale problems. This is because modeling multiple scattering necessitates computing the internal scattered field within the object volume. Furthermore, the effectiveness of existing multiple scattering methods has been demonstrated only under multi-shot tomography. While multiple measurements do alleviate the ill-posedness of the inverse problem, they also increase acquisition time and system complexity. Therefore, it is of interest to investigate if one can exploit multiple scattering using only a single-shot measurement. In this work, we demonstrate successful inverse multiple scattering for large-scale problems and reconstruct 100 million voxels from a single 1megapixel in-line hologram. We show that even under such highly ill-posed conditions, inversion of multiple scattering is still possible and can be used to improve results compared to single scattering techniques.
To calculate multiple scattering, we build our model based on the Born series expansion. 13 To make it computationally efficient, we take a multislice approximation by discretizing the 3D object volume into a series of 2D thin axial slices. At each slice, each object voxel takes a uniform refractive index value. Between neighboring slices, the uniform background medium is assumed.
By adjusting the voxel size and inter-slice distance, our model allows to flexibly trade off computational complexity for model accuracy. At the limit when the voxel size equals the inter-slice distance, our discretization reduces to the existing approaches in. [14][15][16] Our computational structure closely resembles the multislice model (i.e. beam propagation method). [31][32][33] However, our model has the benefit of computing both forward and backward scattering, whereas the latter only accounts for forward multiple scattering.
To compute multiple scattering, we introduce a novel 3D-to-3D operator to efficiently evaluate the internal scattered fields within the volume. The computational framework discretizes the 3D object as a set of 2D slices, and multiple scattering is modeled as recursive propagation among them. Starting from the initial field, each subsequent recursion estimates the next higher-order scattering term within the object volume. This process can be carried out up-to an arbitrary order until the field converges to a steady-state. To evaluate the convergence, we adapt a metric derived from the residual error of the internal field. 19,25 Next, we devise a 3D-to-2D operator that computes the external scattered fields by propagating the multiply scattered internal 3D field to the 2D sensor plane. Finally, the intensity measured by the hologram is the interference between the scattered and the unscattered fields. This further complicates the model by introducing the "twinimage" problem. 34 If only single scattering is considered, our model reduces to linear compressive holography. 11 As a result of multiple scattering computation, the hologram encodes information about the high-angle scattering within the volume, which is otherwise ignored in single scattering based methods. We show that this extra information leads to better recovery of the scatterers, in particular at larger depths.
To solve the inverse scattering problem, we derive an optimization procedure that iteratively minimizes the data-fidelity term measuring the difference between the estimated and measured holograms, and imposes a sparsity-promoting regularization on the object. The overall structure of the algorithm follows the proximal-gradient method. 35 The key ingredient is the gradient com-putation of the data-fidelity term. Conveniently, our recursive forward model leads to a similarly structured recursive gradient computation. Further exploiting the convolution structure in the scattering operators, the algorithm is implemented using efficient FFT-based computations.
Distinct from prior Born-series based models, [14][15][16][17][18] we do not directly measure the full complex field. The effect of sparsity regularization to the twin-image artifact has been studied using single scattering, 11,36 and 2D multiple scattering models. 27 Here, we show that the sparsity is also effective in suppressing twin-image artifacts under 3D multiple scattering models.
An important feature of our multislice-based framework is that the 3D object can be flexibly estimated with any desired number of axial slices, as set by the targeted resolution. In particular, we show that it is possible to use much fewer axial slices to achieve high localization accuracy while still exploiting the extra information contained in the multiple scattering. This allows us to handle much larger scale problems with reduced computational cost as compared to existing techniques that are often limited by fine sampling requirements.
Single scattering based methods tend to under-estimate the refractive index contrast. This under-estimation can be mitigated by incorporating multiple scattering. 17,25,37,38 We show this effect using our multislice-based approach in single-shot in-line holography, and demonstrate improved particle localization and axial resolution under multiple scattering.
Next, we demonstrate the localization accuracy of our method by imaging 3D distributed particles in water at various densities in both simulation and experiment. To facilitate quantitative comparison of different methods, we use a classification framework and use the receiver operating characteristic (ROC) curve to determine each method's best performance. At low particle density, our multiple scattering model converges to the single scattering solution as expected, since the information is dominated by the first order scattering. At high particle density, our model largely improves the accuracy since multiple-scattering becomes more significant. We observe that the localization accuracy is highly depth-dependent. Following the classification framework, we use the dice coefficient 39 to quantify the localization result slice-by-slice. We show that our multiple scattering model provides greater improvement at larger depths.

Forward model
Consider the imaging geometry in Fig. 1(a). An in-line hologram I m at the measurement plane can be written as where E is the scattered field on the measurement plane, u in is the incident plane wave and is assumed to be real on the hologram plane (z 0 = 0) without loss of generality, x 0 = (x 0 , y 0 ) represents the transverse spatial coordinates on the hologram plane. The self-interference term of the scattered field |E| 2 is ignored; the validity of this assumption is discussed in Sec. 3

.2. From
Eq. (2), the scattered field and its "twin-image" contribution is related to the measured hologram after background removal by The background-removed hologram thus represents the real component of the scattered field at the measurement plane, and is given as To model the hologram resulting from multiple scattering up-to the K th order, we apply the framework of Born series This computed scattered-field u s k is added to the incident-field u in to obtain the next higher-order Born-field u k . This process is recursively applied to compute the multiply scattered field within the volume. expansion, 13,16 which gives two coupled equations: where h is the 3D Green's function and u K (x , z ) is the K th order multiply scattered field within the volume Ω. This mathematical model is based the scalar Helmholtz equation, 13 and polarization effects are neglected. This internal field is computed recursively within the support Ω using the 13 where n med is the index of the homogeneous background medium, and k is the wave number in free space. For simplicity, f is assumed to be real valued and absorption effects are ignored. We note that the spatial coordinates associated with the object are within the 3D support: (x, z) ∈ Ω, (x , z ) ∈ Ω; whereas the hologram is measured outside the support: To compute higher order scattering, the initial condition is u 0 (x, z) = 0, and k = 1, 2, ..., K indexes the scattering order. When K = 1, u 1 (x, z) = u in (x, z) is the incident field [from Eq. (4)], and Eq. (3) reduces to the first-Born approximation that linearly relates the object to the singly scattered field. When K = 2, this relation becomes nonlinear and the second order multiple scattering is taken into account via modeling of the additional interaction between the object and the field within the volume. For larger K, the approximation becomes more accurate by accounting for K-order multiple scattering.
Equations (3) and (4) can be discretized to get the following recursive forward model and N z are the number of pixels along the x, y, and z within Ω, respectively. diag(u) represents a diagonal matrix with the vector u on the main diagonal.
We consider the 3D volume to be a set of discrete 2D slices along the longitudinal axis. H and G are the scattering operators which represent propagation among the object slices and to the hologram plane. H propagates the field from the object support to the hologram plane. H = where L z −z is a diagonal matrix representing the discrete transfer function that performs propagation between two slices, from the slice z to z . This treatment of the H operator is similar to that in. 11 G is the multiple scattering operator that performs propagation from the object volume to within itself (Fig. 2).
, and contains transfer functions that propagate the field from each slice to every slice within the support. There are two methods of computing the elements in the transfer function L z −z having dimension (N x N y ) × (N x N y ), the direct and the angular spectrum methods. 41 We use the direct method, in which the Green's function h(r) = exp(ιk|r|) |r| is sampled in the spatial domain, followed by slice-wise 2D FFT, where r = (x, y, z).
An important numerical treatment to h(r) is around the singularity at |r| = 0. We adapt the technique in 42,43 and consider a spherical exclusion zone around |r| = 0 of radius a, inside which the Green's function is assumed to take a constant value. Effectively, this assigns an "averaged" Green's function value around the singular region. Empirically, we found that at low refractive index, the choice of a does not significantly affect the result as long as the center voxel (at |r| = 0) is excluded. For high refractive index, a highly affects the convergence of the forward model. We set a to match the largest expected radius of the particles. This means that the strong multiple scattering inside each individual particle cannot be reliably modeled at high contrast; hence it is ignored. Only particle-particle interactions are modeled. Correspondingly, during the inversion, a sets the largest particle size that can be recovered by our model for high index particles.

Inverse problem
To estimate the object f from the holographic measurement, we need to solve Eqs. (5,6). Unlike traditional digital holography, this problem is nonlinear when K > 1. We devise an inverse scattering algorithm 44, 45 that minimizes a cost function C(f ) to compute the estimated objectf as followŝ where D(f ) = 1 2 I br − Re{E est } 2 is the data fidelity term in which I br is the real valued measured hologram after background removal, and E est is the complex-valued scattered field estimate from our model [Eq. (5)]. Re{E est } provides an estimate for I br [Eq.
(2)], · represents the L2-norm, F is the convex set that constrains the object to be nonnegative, and τ is the regularization parameter that is empirically tuned. f TV imposes a penalty on the total variation (TV) of the object, and is defined as where D : R N → R N ×3 is the discrete gradient operator with matrices D x , D y , and D z denoting the finite difference operations along x-, y-, and z-directions, respectively.
The minimization in Eq. (7) is implemented via the proximal-gradient method, 46 in which the t th iteration is written as where prox τ TV g) 47 and α is the step size set via backtracking line search. 48 The initialization is f 0 = 0.
Similar to the forward model, the gradient computation is also a K-order recursion.
for k = 1, 2, ..., K. Here, r = Re{E est } − I br is the residual, A H and A represent the Hermitian and complex conjugate of the matrix A, respectively. The recursion is initialized with ∂u 0 /∂f = 0.
Brute-force evaluation of the gradient is highly computationally intensive for large scale problems, with each vector having more than a few million elements. We devise a computationally efficient implementation by making use of the FFT-based structures in G and H operators. This algorithm extends the framework in 16 on small-scale 2D to large-scale 3D problems, and further demonstrates reconstruction from intensity-only as opposed to full-field measurements.

Results
We test our model on both simulations and experiments. In our experiment, the inline holography setup uses a linearly polarized HeNe laser (632.8nm, 500:1 polarization ratio, Thorlabs HNL210L) that is collimated for illumination [ Fig. 1 Importantly, Eq. (6) requires computation of high-angle multiply scattered field propagating within the volume; thus the internal field needs to be sampled at the Nyquist rate λ/2. In our system, the camera's pixel-size is 3.45µm, and the effective lateral sampling size after magnification is δx = δy = 172.5nm. This satisfies the sampling requirement in the medium, where the wavelength is λ = 630nm/n water = 473.7nm. We set the voxel size along axial direction δz = 172.5nm, such that the the voxels are cubic. The spacing between slices is assumed to contain uniform background medium, and is set to be 5µm, approximately matching the system's axial resolution of  25 which is an accurate forward model that incorporates multiple scattering, including scattering within each particle. It is based on a rigorous optimization procedure that solves the Lippmann Schwinger equation. We further simulate the hologram at high particle densities using our model with a sufficiently high scattering order, e.g. K ≥ 10, such that the model converges and the simulated field closely estimates the actual. In order to validate the convergence, we present an evaluation metric and show that the model converged within the first few scattering orders for all tested scenarios. Improvement of our method compared to single scattering is presented quantitatively.
The Boston University Shared Computing Cluster (SCC) was used for all computations. The average times of computing one iteration for the single and multiple scattering models on a 512 × 512 × 50 grid were 58 and 257 seconds, respectively. All reconstructions were run for 100 itera- tions. In what follows, we present our findings.

Effect of multiple scattering in small-scale inversion: A multislice-based approach
It has been shown that in the presence of strong multiple scattering, the single scattering models underestimate the permittivity contrast. 17,25,37 Here we validate our model on a small-scale simulation and make similar observation by showing that the underestimation is mitigated as multiple scattering is incorporated in the inversion.
The utility of our multislice-based computational approach is also demonstrated, in which the number of axial slices can be arbitrarily chosen in the inverse reconstruction. Effectively, we approximate the 3D object with a fixed number of slices, such that the computation is tractable when expanding to large-scale problems.
We simulate a volume of 44 × 44 × 6.4µm 3 , discretized as a 256 × 256 × 37 object, containing 8 spheres in water, each with refractive index n =1.43 and diameter 1 µm. In Fig. 3(a), we depict the central 6.2 × 6.2 × 6.4µm 3 region of this object. The multiple scattering is significant in the presence of occluding geometry along the optical axis, and a refractive index contrast of δn = 0.1. It is known that occlusion causes strong axial field coupling via multiple scattering between scatterers, which is ignored by the first-order model. 49 We therefore expect that incorporating multiple scattering will improve object estimation from the hologram.
An inline hologram is simulated at 5µm from the front slice using the SEAGLE. The hologram is then inverted using our multislice-based method incorporating 1 st , 2 nd order scattering. The scattered intensity |E| 2 is included when simulating the hologram. During the inversion, this term is ignored following the procedure in Sec. 2. Our results indicate that even under this approximation, our model suppresses the underestimation artifacts by incorporating multiple scattering.
In order to test the utility of our multislice-based approach, we perform reconstruction for two cases. In the first case, we reconstruct all 37 slices for the object [ Fig. 3(c)]. The reconstruction based on the 1 st order scattering underestimates the refractive indices. We attribute this artifact to the strong axial coupling via multiple scattering between the occluding particles, which is ignored by the 1 st order model. The underestimation is mitigated when 2 nd order scattering is included in the inverse model.
In the second case, we estimate the object using only 3 slices to perform the inverse scattering reconstruction [ Fig. 3(d)]. The reconstruction in this case approximates the 3D object comprising of 3 discrete slices. We observe that our method is able to detect the 8 spheres as disks at the cor-rect axial locations corresponding to the centers of each particles. When using only 3 slices in the reconstruction, the model has less number of slices to create the same effect at the measurement plane as the 37-slice ground truth object. In order to compensate for this, the reconstructed scattering density can be approximated as the integrated permittivity contrast along the optical axis, while still correctly localizing the particles. In the 1 st order result, we observe smaller contrast and worse axial sectioning. The 2 nd order multiple scattering improves the contrast as well as axial sectioning, resulting in better localization capability.

Large-scale inversion of multiple-scattering: simulation
Next, we demonstrate the inversion of multiple scattering from single-shot measurement in largescale. For this purpose, we design a simulation which involves estimating the concentration of particles in a suspension from its inline hologram. We show that our multiple scattering model improves the accuracy in estimating the particle density, particularly at larger depths.
The simulated volume is 88 × 88 × 250µm 3 , discretized on a 512 × 512 × 50 grid, in which disk-shaped scatterers of 1 µm diameter and constant refractive index are suspended randomly in water, in varying densities. The disk-shape may not represent actual spheres, as would be in the real application, however it is taken as an approximation due to stringent sampling requirements for such a large volume. We consider two values of the refractive index contrast δn = 0.01 and 0.19. For each volume, we first simulate holograms using 20 th order scattering, followed by the reconstruction using 1 st and 2 nd order models. The particle density is estimated for each reconstructed volume using the ImageJ 3D objects counter toolbox. 50 The optimal threshold parameter used for calculating the density is determined using the ROC.
As a measure of particle density, we consider the geometric cross-section R g , which corre- is negligible compared to the hologram for low particle densities, and becomes gradually important as the particle density increases. (b) The ratio of between the total intensity of the hologram and the |E| 2 terms for all values of R g tested in the simulation. For R g ≤ 0.1, the total intensity of the hologram is at least an order of magnitude larger than the |E| 2 term.
sponds to the fraction of the hologram area directly occluded by the scatterers, defined as R g = total cross-sections of all scatterers area of the hologram where N p represents the total number of scatterers in the volume. This metric is valid for scattering and 4000, respectively. We simulate five random object volumes for each value of R g and refractive index contrast, and report the mean statistics of the reconstructions in Fig. 5. In Sec. 2.1, we assumed that the intensity of the scattered field |E| 2 is negligible in the forward model. In this study, this assumption holds true when R g ≤ 0.1, where the contribution of |E| 2 is at least an order of magnitude smaller than the total intensity of the hologram [ Fig. 4]. For higher particle density, |E| 2 becomes increasingly significant, which leads to greater model error.
For the series expansion approach used in our model, it is important to evaluate its convergence.
In Fig. 5(a), we present the convergence properties of the forward model under our experimental conditions. While in general higher-order terms are required for convergence under stronger scattering, the 2 nd order scattering is sufficient for most of the cases studied. Our convergence metric At larger depths, the accuracy deteriorates for both methods. Our multiple scattering method performs notably better at larger depths for higher particle densities. (b) The slice-wise dice coefficient plotted as a function of slice depth also indicates that the multiple scattering model provides improved segmentation accuracy, especially at greater depth. (c) The particle localization accuracy is quantified using the ROC curve. The curves corresponding to the multiple scattering solutions consistently have larger areas underneath, indicating better localization accuracy as compared to the single scattering method in all cases studied.
e is defined by the residual error of the field within the 3D volume, 19,25,30 as where A = I − Gdiag(f ). This convergence metric essentially measures the self-consistency of the total internal field. 25 For K-order scattering, it computes the norm of the residual contribution from (K + 1)-order scattering, which must approach zero in the case of convergence.
Next, we evaluate the reconstruction accuracy by measuring the signal to noise ratio (SNR) where f true andf are the true and estimated objects, respectively. For higher index contrast (δn = 0.19), our multiple-scattering model performs consistently better than the single-scattering model for all densities tested [ Fig. 5(b,c)]. Generally, the reconstruction performance from both methods drop as the density increases. We attribute this to stronger higher-order scattering and decrease in the object sparsity. The stronger scattering introduces higher order nonlinearity through Eq. (6), making the problem harder to invert. The decrease in object sparsity leads to an effective smaller measurement-to-unknown ratio, further worsen the ill-posedness of the problem.
For higher contrast (δn = 0.19), at low particle density (R g ≤ 0.02), single and multiple scattering methods perform similarly. This is expected since multiple scatterings are weak due to the small scattering cross section. For 0.02 < R g ≤ 0.1, our method outperforms the single scattering method, providing a better estimate of the actual particle density. For R g > 0.1, the SNR drops below 1dB, and we empirically consider the reconstruction has failed [ Fig. 5(b)].
For lower index contrast (δn = 0.01), both multiple and single scattering methods perform almost identically for all densities tested, which indicates that the contribution from multiple scattering is negligible for low refractive index contrast [ Fig. 5(c)]. 3D renderings and cross-sectional reconstructions at different depths are depicted in Fig. 5(d).
The depth-dependent performance is highlighted in Fig. 6(a,b). Close to the hologram plane (z = 5µm), single and multiple scattering reconstructions are similar, and match the ground truth.
At larger depth (z = 190µm), the single-scattering reconstruction degrades and results in a large number of missing particles [ Fig. 6(a)]. Our multiple scattering model improves the localization at larger depths. By treating particle localization as a binary classification problem, we use the ROC curve to determine the optimal segmentation threshold when quantifying the voxels reconstructed by each method. 10 This allows us to evaluate the localization accuracy slice-by-slice, whose statistics are accumulated by 5 different object volumes for each particle density. The statistic we use for comparison is the dice coefficient which is used to gauge the similarity of two samples and is defined as where p i and g i each represents a voxel from the predicted and ground truth binary segmentation volumes respectively and i indexes the voxels of each 3D volume. Consistent with the visual inspection in Fig. 6(a), the dice coefficient clearly indicates improvement at larger depth using our multiple scattering model [ Fig. 6  The reconstructed particle density follows a trend similar to the simulation where multiple-scattering performs better than the single scattering method for R g ≤ 0.1; both methods fail for R g > 0.1. (c) As R g increases, the hologram gradually resembles speckle patterns, as quantified by the contrast ratio (CR). For low density, both multiple and single scattering methods perform similarly. For high density, the underestimation of particles from the single scattering method is clearly visible, especially at increased depth. Our multiple scattering model mitigates the underestimation as it accounts for the intercoupling between particles whose strength increases as the depth. The traditional BPM is effective for low density, but completely fails for high density and the reconstruction resembles speckles throughout the volume.
holograms are recorded at each concentration, and then used for reconstructions and density estimation. Background subtraction is performed on each hologram as a preprocessing step to remove static artifacts. The inversion is performed using our method with second-order multiple scattering (K = 2), the single scattering method, and BPM.
First, we evaluate the results based on the optimization convergence cost [ Fig. 7(a)]. The multiple scattering method converges to a lower value than single scattering for all densities, indicating better fit to the cost function C(f ). The cost increases for both methods with R g , depicting degradation of reconstruction with increase in particle density.
Next, we assess the estimated particle density. Our multiple scattering model consistently performs better than the single scattering model for R g ≤ 0.1 [ Fig. 7(b)]. For R g > 0.1, reconstruction fails for both methods as also found in our simulation. Hence, we use R g = 0.1 as an empirical performance bound of our method for this application.
Evidently, the recorded holograms gradually resemble speckle patterns as the particle density increases [ Fig. 7(c)]. We quantify the hologram's contrast ratio (CR) at each density, which can be used as an alternative metric. The CR is calculated as the ratio of the standard deviation to the mean. 51 At the critical R g = 0.1 concentration, the CR is around 0.335.
Finally, we closely examine the 3D reconstructions for R g = 0.0063 and R g = 0.1 [ Fig. 8].
For the low density case, single and multiple scattering methods perform similarly due to weak multiple scattering. For the high density case where multiple scattering becomes significant, our method outperforms the single scattering model. Particle localization degrades with increased depth; our multiple-scattering method provides a more uniform estimation and better localization at increased depth, matching our observations in the simulation. While BPM is able to reconstruct individual particles at the low density, it completely fails for high density, and resembles speckles extended across the object volume.

Conclusion
We have presented a new computational framework for utilizing multiple scattering in in-line holography for large-scale 3D particle localization. Our model recursively computes both forward and backward multiple scattering in a computationally efficient manner. Both simulations and experiments demonstrate the significance of modeling multiple scattering in alleviating depthdependent artifacts and improving the 3D localization accuracy compared to traditional methods.
Our method may opens up new opportunities for large-scale imaging applications utilizing multiple scattering.
Our model is currently limited by the convergence regime of the classical Born series expansion, preventing it to be applied to particle density higher than 0.1 geometric cross-section. Recent work on convergent Born series expansion 52 provides a promising avenue to extend our model to higher scattering scenarios.
The novel multislice structure proposed in our model provides a flexible framework for trading computational cost for model accuracy. Still, higher order scattering calculation necessitates longer computational times, which is less appealing for applications requiring real-time reconstructions.
To facilitate rapid volumetric estimation without sacrificing accuracy, recent machine learning based inverse scattering approaches 53-55 may be explored in the future.

Disclosures
There are no conflicts of interest related to this article.