Three-dimensional (3-D) 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–9 Reconstructing dense samples, however, remains challenging.¹⁰ Standard backpropagation method (BPM) can only handle low particle density.¹⁰ 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 a framework that accounts for multiple scattering in in-line digital holography and enables accurate 3-D 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.¹³ 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–18 contrast source inversion,^19–22 modified gradient,^23,24 series expansion with accelerated gradient descent on Lippmann–Schwinger equation (SEAGLE),^25,26 and hybrid methods.^27–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 multishot 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 whether 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 1-megapixel 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.¹³ To make it computationally efficient, we take a multislice approximation by discretizing the 3-D object volume into a series of two-dimensional (2-D) 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 interslice distance, our model allows us to flexibly trade computational complexity for model accuracy. At the limit when the voxel size equals the interslice distance, our discretization reduces to the existing approaches in Refs. 14–16. Our computational structure closely resembles the multislice model (i.e., beam propagation method).^31–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 3-D-to-3-D operator to efficiently evaluate the internal scattered fields within the volume. The computational framework discretizes the 3-D object as a set of 2-D 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 3-D-to-2-D operator that computes the external scattered fields by propagating the multiply scattered internal 3-D field to the 2-D 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 “twin-image” problem.³⁴ If only single scattering is considered, our model reduces to linear compressive holography.¹¹ 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.³⁵ The key ingredient is the gradient computation 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–18 we do not directly measure the full complex field. The effect of sparsity regularization on the twin-image artifact has been studied using single-scattering^11,36 and 2-D multiple-scattering models.²⁷ Here, we show that the sparsity is also effective in suppressing twin-image artifacts under 3-D multiple-scattering models.

An important feature of our multislice-based framework is that the 3-D 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 underestimate the refractive index contrast. This underestimation 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 3-D 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³⁹ to quantify the localization result slice-by-slice. We show that our multiple-scattering model provides greater improvement at larger depths.

Consider the imaging geometry in Fig. 1(a). An in-line hologram Im at the measurement plane can be written as Im(x0)=|uin(x0,0)+E(x0,0)|2=|uin|2+2uinRe{E}+|E|2,(1)where E is the scattered field on the measurement plane, uin is the incident plane wave and is assumed to be real on the hologram plane (z0=0) without loss of generality, and x0=(x0,y0) 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. The scattered field and its “twin-image” contribution is related to the measured hologram after background removal by Re{E(x0,0)}=[E(x0,0)+E*(x0,0)]2=Im(x0)−|uin|22uin.(2)

The background-removed hologram thus represents the real component of the scattered field at the measurement plane and is given as Ibr(x0)=[Im(x0)−|uin|2]/2uin. To model the hologram resulting from multiple scattering up-to the K’th order, we apply the framework of Born series expansion,^13,16 which gives two coupled equations: E(x0,0)=∭Ωf(x′,z′)uK(x′,z′)h(x0−x′,−z′)d2x′dz′,(3)uk(x,z)=uin(x,z)+∭Ωf(x′,z′)uk−1(x′,z′)h(x−x′,z−z′)d2x′dz′,(4)where h is the 3-D Green’s function and uK(x′,z′) is the K’th-order multiply scattered field within the volume Ω. This mathematical model is based on the scalar Helmholtz equation,¹³ and polarization effects are neglected. This internal field is computed recursively within the support Ω using the Born series [Eq. (4)]. The permittivity contrast f is related to the refractive index by f(x′,z′)=k24π[n2(x′,z′)−nmed2],¹³ where nmed 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 3-D support, (x,z)∈Ω, (x′,z′)∈Ω, whereas the hologram is measured outside the support, (x0,z0)∉Ω. To compute higher order scattering, the initial condition is u0(x,z)=0, and k=1,2,…,K indexes the scattering order. When K=1, u1(x,z)=uin(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: E=Hdiag(uK)f,(5)uk=uin+Gdiag(uk−1)f,(6)for k=1,2,…,K. Here, f and uk have dimension (NxNyNz)×1, and E is (NxNy)×1. Nx, Ny, and Nz are the number of pixels along the x, y, and z within Ω, respectively. Here, diag(u) represents a diagonal matrix with the vector u on the main diagonal.

We consider the 3-D volume to be a set of discrete 2-D slices along the longitudinal axis. H and G are the scattering operators that represent propagation among the object slices and to the hologram plane. H propagates the field from the object support to the hologram plane. H=KHQ0B, where B=bldiag(K,..,K) is a block-diagonal matrix. K and KH are the 2-D DFT and the inverse 2-D DFT matrices, respectively, each with dimension (NxNy)×(NxNy). Q0=[L0−1L0−2…L0−Nz], where Lz′−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 Ref. 11. G is the multiple scattering operator that performs propagation from the object volume to within itself (Fig. 2). G=BHRB, where R=[Q1,Q2,…,QNz] and Qm=[Lm−1Lm−2…Lm−Nz]. R has the dimension of (NxNyNz)×(NxNyNz) 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 Lz′−z having dimensions (NxNy)×(NxNy), the direct and the angular spectrum methods.⁴¹ 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 slicewise 2-D FFT, where r=(x,y,z).

An important numerical treatment to h(r) is around the singularity at |r|=0. We adapt the technique from Refs. 42 and 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 a 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.

To estimate the object f from the holographic measurement, we need to solve Eqs. (5) and (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 object f^ as follows: f^=argminf∈F{C(f)}=argminf∈F{D(f)+τ‖f‖TV},(7)where D(f)=12‖Ibr−Re{Eest}‖2 is the data fidelity term in which Ibr is the real valued measured hologram after background removal, and Eest is the complex-valued scattered field estimate from our model [Eq. (5)]. Re{Eest} provides an estimate for Ibr [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. Here, ‖f‖TV imposes a penalty on the total variation (TV) of the object and is defined as ‖f‖TV=Δ∑n=1N‖[Df]n‖l2=∑n=1N|[Dxf]n|2+|[Dyf]n|2+|[Dzf]n|2,(8)where D:RN→RN×3 is the discrete gradient operator with matrices Dx, Dy, and Dz denoting the finite difference operations along x-, y-, and z-directions, respectively.

The minimization in Eq. (7) is implemented via the proximal-gradient method,⁴⁶ in which the t’th iteration is written as ft←proxτTV[ft−1−α∂D(f)∂f],(9)where proxτTV(g)=Δargminf∈F{12‖f−g‖2+τ‖f‖TV} is the proximal operator for TV minimization,⁴⁷ and α is the step size set via backtracking line search.⁴⁸ The initialization is f0=0.

Similar to the forward model, the gradient computation is also a K-order recursion: ∂D(f)∂f=Re{diag(u¯K)HHr+[∂uK∂f]Hdiag(f)HHr},(10)[∂uk∂f]Ha=diag(u¯k−1)GHa+[∂uk−1∂f]Hdiag(f)GHa,(11)for k=1,2,…,K. Here, r=Re{Eest}−Ibr is the residual; AH and A¯ represent the Hermitian and complex conjugate of the matrix A, respectively. The recursion is initialized with ∂u0/∂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 Ref. 16 from small-scale 2-D to large-scale 3-D problems and further demonstrates reconstruction from intensity-only as opposed to full-field measurements.

We test our model on both simulations and experiments. In our experiment, the inline holography setup uses a linearly polarized HeNe laser (632.8 nm, 500:1 polarization ratio, Thorlabs HNL210L) that is collimated for illumination ^[Fig. 1(b)]. A 4F system with a 20× objective lens (0.4 NA, CFI Plan Achro) and a 200-mm tube lens is used to collect the scattered field with the Nyquist sampling requirement satisfied. A CMOS sensor (FLIR GS3-U3-123S6M-C) is used to capture the holograms. The object consists of polystyrene microspheres with nominal diameter 0.994±0.021μm (Thermofisher Scientific 4009A) suspended in deionized water. The suspension is held in a quartz-cuvette with inner dimensions 40mm×40mm×0.5mm. We are interested in localizing the individually suspended scatterers. A shutter speed of 5 ms was used and found to be sufficiently fast to capture the holograms without any motion artifacts from suspended particles. The illumination beam diameter is less than the width of the cuvette and larger than the CMOS sensor area to avoid edge artifacts. The front focal plane of the objective lens was set just outside the inner wall of the cuvette for hologram recording.

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/nwater=473.7nm. We set the voxel size along axial direction δz=172.5nm, such that 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 λ/(1−1−NA2)=5.7μm. During the computation, 2× zero-padding is used in all FFTs to avoid boundary artifacts. We demonstrate large-scale inverse scattering that reconstructs 100 million voxels in a 176μm×176μm×500μm volume.

For large-scale simulation, we model the system parameters to approximately match the physical setup. On such a scale, rigorous solutions such as FDTD are computationally prohibitive, and sample complexity makes analytical solutions such as Mie theory nontrivial. We first study the effect of multiple scattering on simulated holograms using 3-D SEAGLE,²⁵ 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 s, respectively. All reconstructions were run for 100 iterations. In what follows, we present our findings.

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 3-D 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μm×44μm×6.4μm, discretized as a 256×256×37 object, containing eight spheres in water, each with refractive index n=1.43 and diameter 1μm. In Fig. 3(a), we depict the central 6.2μm×6.2μm×6.4μm 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.⁴⁹ 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 first-, second-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 first-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 first-order model. The underestimation is mitigated when second-order scattering is included in the inverse model.

In the second case, we estimate the object using only three slices to perform the inverse scattering reconstruction ^[Fig. 3(d)]. The reconstruction in this case approximates the 3-D object comprising of three discrete slices. We observe that our method is able to detect the eight spheres as disks at the correct axial locations corresponding to the centers of particles. When using only three slices in the reconstruction, the model has a smaller 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 first-order result, we observe smaller contrast and worse axial sectioning. The second-order multiple scattering improves the contrast as well as axial sectioning, resulting in better localization capability.

Next, we demonstrate the inversion of multiple scattering from single-shot measurement in large-scale. For this purpose, we design a simulation that 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μm×88μm×250μm, 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 first- and second-order models. The particle density is estimated for each reconstructed volume using the ImageJ 3-D objects counter toolbox.⁵⁰ 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 Rg, which corresponds to the fraction of the hologram area directly occluded by the scatterers, defined as Rg=total cross-sections of all scatterersarea of the hologram≈Npπr2NxNyδxδy,(12)where Np represents the total number of scatterers in the volume. This metric is valid for scattering domains that are not very thick, as in our case. For a collection of identical particles suspended in a homogeneous medium, the geometric and scattering cross-sections are directly related,⁵¹ in which the latter is a direct measure of the fraction of the incident light scattered by an object. For higher values of Rg, we thus expect greater contributions of multiple scattering. From the signal processing perspective, Rg also measures the sparsity of the problem as it approximates the ratio between the number of nonzero unknowns and the number of measurements. The values of Rg tested are 0.01, 0.02, 0.05, 0.1, 0.2, and 0.4, corresponding to Np=100, 200, 500, 1000, 2000, and 4000, respectively. We simulate five random object volumes for each value of Rg and refractive index contrast and report the mean statistics of the reconstructions.

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 Rg≤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 second-order scattering is sufficient for most of the cases studied. Our convergence metric e is defined by the residual error of the field within the 3-D volume,^19,25,30 as e=‖AuK−u0‖,(13)where A=I−Gdiag(f). This convergence metric essentially measures the self-consistency of the total internal field.²⁵ 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) SNR=20log10(‖ftrue‖‖ftrue−f^‖),(14)where ftrue and f^ 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 ^{[Figs. 5(b) and 5(c)]}. In general, the reconstruction performance from both methods drops 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 worsening the ill-posedness of the problem.

For higher contrast (δn=0.19), at low particle density (Rg≤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<Rg≤ 0.1, our method outperforms the single scattering method, providing a better estimate of the actual particle density. For Rg>0.1, the SNR drops below 1 dB, and we empirically consider that 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)]. 3-D renderings and cross-sectional reconstructions at different depths are shown in Fig. 5(d).

The depth-dependent performance is highlighted in Figs. 6(a) and 6(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.¹⁰ This allows us to evaluate the localization accuracy slice-by-slice, whose statistics are accumulated by five different object volumes for each particle density. The statistic we use for comparison is the Dice coefficient that is used to gauge the similarity of two samples and is defined as D=2∑i=1Npigi∑i=1Npi2+∑i=1Ngi2,(15)where pi and gi each represent a voxel from the predicted and ground truth binary segmentation volumes, respectively, and i indexes the voxels of each 3-D 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(b)]. In addition, the area under each ROC provides a direct measure of the algorithm’s overall classification performance. Our results indicate that the multiple-scattering model consistently outperforms the single scattering method ^[Fig. 6(c)].

Finally, we demonstrate our method on a set of large-scale experiments. We reconstruct over 100 million object voxels (1024×1024×100) from each 1-megapixel hologram. Our multiple-scattering model significantly improves the 3-D localization accuracy as compared to the BPM and single-scattering methods. Notably, our experimental results closely match the simulation.

We prepare polystyrene microsphere suspensions, ranging from dense to sparse concentrations via successive dilution, with corresponding Rg values of 0.4, 0.2, 0.1, 0.05, 0.0125, and 0.0063. Five 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 methods, 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 Rg, 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 Rg≤0.1 ^[Fig. 7(b)]. For Rg>0.1, reconstruction fails for both methods as also found in our simulation. Hence, we use Rg=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.⁵¹ At the critical Rg=0.1 concentration, the CR is around 0.335.

Finally, we closely examine the 3-D reconstructions for Rg=0.0063 and Rg=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.

We have presented a new computational framework for utilizing multiple scattering in in-line holography for large-scale 3-D 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 depth-dependent artifacts and improving the 3-D localization accuracy compared to traditional methods. Our method may open 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 its application to particle density higher than 0.1 geometric cross-section. Recent work on convergent Born series expansion⁵² provides a promising avenue to extend our model to higher scattering scenarios.

The 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.