Translator Disclaimer
18 June 2019 Holographic particle localization under multiple scattering
Author Affiliations +
We introduce a computational framework that incorporates multiple scattering for large-scale three-dimensional (3-D) particle localization using single-shot in-line holography. Traditional holographic techniques rely on single-scattering models that become inaccurate under high particle densities and large refractive index contrasts. Existing multiple scattering solvers become computationally prohibitive for large-scale problems, which comprise millions of voxels within the scattering volume. Our approach overcomes the computational bottleneck by slicewise computation of multiple scattering under an efficient recursive framework. In the forward model, each recursion estimates the next higher-order multiple scattered field among the object slices. In the inverse model, each order of scattering is recursively estimated by a nonlinear optimization procedure. This nonlinear inverse model is further supplemented by a sparsity promoting procedure that is particularly effective in localizing 3-D distributed particles. We show that our multiple-scattering model leads to significant improvement in the quality of 3-D localization compared to traditional methods based on single scattering approximation. Our experiments demonstrate robust inverse multiple scattering, allowing reconstruction of 100 million voxels from a single 1-megapixel hologram with a sparsity prior. The performance bound of our approach is quantified in simulation and validated experimentally. Our work promises utilization of multiple scattering for versatile large-scale applications.



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.79 Reconstructing dense samples, however, remains challenging.10 Standard backpropagation 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 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.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,1418 contrast source inversion,1922 modified gradient,23,24 series expansion with accelerated gradient descent on Lippmann–Schwinger equation (SEAGLE),25,26 and hybrid methods.2730 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.13 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. 1415.16. Our computational structure closely resembles the multislice model (i.e., beam propagation method).3133 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.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 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,1418 we do not directly measure the full complex field. The effect of sparsity regularization on the twin-image artifact has been studied using single-scattering11,36 and 2-D multiple-scattering models.27 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 coefficient39 to quantify the localization result slice-by-slice. We show that our multiple-scattering model provides greater improvement at larger depths.


Theory and Method


Forward Model

Consider the imaging geometry in Fig. 1(a). An in-line hologram Im at the measurement plane can be written as

Eq. (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

Eq. (2)


Fig. 1

In-line holography with multiple scattering. (a) A plane-wave is incident on a 3-D object containing distributed scatterers. The field undergoes multiple scattering within the volume and then propagates to the image plane. A hologram is recorded, which is then used to estimate the unknown scatterers’ distribution. (b) An inline holography setup is used that consists of a collimated laser for illumination and a 4F system for magnification. (c) The raw data are a single hologram. (d) The reconstruction implements a nonlinear inverse multiple scattering algorithm.40 (e) The output estimates the 3-D distribution of the scatterers.


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:

Eq. (3)


Eq. (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,13 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],13 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:

Eq. (5)


Eq. (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=[L01L02L0Nz], where Lzz 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=[Lm1Lm2LmNz]. 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 Lzz having dimensions (NxNy)×(NxNy), 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 slicewise 2-D FFT, where r=(x,y,z).

Fig. 2

Illustration of the 3-D internal scattered field operator G in Eq. (6). (a) Each object slice f is first voxelwise multiplied by the lower order scattered field uk1; it is then propagated to every other slice within the volume. (b) This computed scattered-field uks is added to the incident-field uin to obtain the next higher-order Born-field uk. This process is recursively applied to compute the multiply scattered field within the volume.


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.


Inverse Problem

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 algorithm44,45 that minimizes a cost function C(f) to compute the estimated object f^ as follows:

Eq. (7)

where D(f)=12IbrRe{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, fTV imposes a penalty on the total variation (TV) of the object and is defined as

Eq. (8)

where D:RNRN×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,46 in which the t’th iteration is written as

Eq. (9)

where proxτTV(g)=ΔargminfF{12fg2+τfTV} is the proximal operator for TV minimization,47 and α is the step size set via backtracking line search.48 The initialization is f0=0.

Similar to the forward model, the gradient computation is also a K-order recursion:

Eq. (10)


Eq. (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.5  nm. This satisfies the sampling requirement in the medium, where the wavelength is λ=630  nm/nwater=473.7  nm. We set the voxel size along axial direction δz=172.5  nm, 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 λ/(11NA2)=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,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., K10, 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.


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 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.49 We therefore expect that incorporating multiple scattering will improve object estimation from the hologram.

Fig. 3

Small-scale multiple-scattering inversion. (a) An accurate 3-D forward model is used to simulate the hologram. (b) Multislice 3-D reconstruction is performed from a single simulated measurement using our method. The number of slices in the inverse reconstruction can be flexibly chosen. (c) Full 3-D inversion is performed by reconstructing all axial slices in the original object using our method. The multiple-scattering method outperforms the single-scattering method by providing both more accurate permittivity contrast estimation and improved optical sectioning. (d) Our multislice approach enables 3-D reconstruction using a much reduced number of slices while still maintaining the benefit of incorporating multiple scattering. Reconstruction using only three slices is compared to demonstrate the improved localization capability by our method.


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.


Large-Scale Inversion of Multiple-Scattering: Simulation

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.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 Rg, which corresponds to the fraction of the hologram area directly occluded by the scatterers, defined as

Eq. (12)

Rg=total cross-sections of all scatterersarea of the hologramNpπr2NxNyδxδy,
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,51 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 Rg0.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.

Fig. 4

Effect of particle density on the scattered intensity term |E|2 contribution in the hologram. (a) Contribution is negligible compared to the hologram for low particle densities and becomes gradually important as the particle density increases. (b) The ratio between the total intensity of the hologram and the |E|2 terms for all values of Rg tested in the simulation. For Rg0.1, the total intensity of the hologram is at least an order of magnitude larger than the |E|2 term.


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

Eq. (13)

where A=IGdiag(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.

Fig. 5

Validation of our multiple-scattering method on large-scale simulation. (a) Convergence properties of the forward model are studied under varying particle densities. Higher-order scattering is generally required for convergence when the object is strongly scattering. In most cases studied, two orders of scattered field sufficiently capture the majority of the contribution. (b) For higher refractive index contrast (δn=0.19), multiple-scattering performs similarly to single-scattering for low concentration (Rg0.02), and better than single-scattering for 0.02<Rg0.1. Reconstruction fails for very high concentration (Rg>0.1), i.e., when the SNR drops below an empirically chosen value of 1 dB. The error in the predicted versus the ground truth particle concentrations also shows a similar trend. (c) For lower contrast (δn=0.01), multiple scattering contributions are negligible and both methods give similar performance. (d) A 3-D rendering depicting localized particles is shown for δn=0.19 and Rg=0.1. Both methods have similar performance for slices close to the image plane, but our multiple-scattering model performs better at increased depths.


Next, we evaluate the reconstruction accuracy by measuring the signal-to-noise ratio (SNR)

Eq. (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 (Rg0.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<Rg0.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.10 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

Eq. (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)].

Fig. 6

Reconstruction performance as a function of depth. (a) Segmentation maps of reconstructed slices (zoomed-in 51μm×51  μm regions) at different depths (true positive, white; true negative, black; false positive, green; false negative, pink). For object slices close to the hologram, both multiple and single scattering methods provide high accuracy. 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 slicewise 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.



Large-Scale Experimental Validation

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.

Fig. 7

Experimental validation of our method in large-scale. (a) The multiple-scattering model converges to a lower cost than the single-scattering model for all concentrations indicating better fit to the cost function. (b) The reconstructed particle density follows a trend similar to the simulation where multiple-scattering performs better than the single-scattering method for Rg0.1; both methods fail for Rg>0.1. (c) As Rg increases, the hologram gradually resembles speckle patterns, as quantified by the CR.


Next, we assess the estimated particle density. Our multiple-scattering model consistently performs better than the single-scattering model for Rg0.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.51 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.

Fig. 8

A 3-D visualization of the localized particles under different concentrations from our experiment and their 200×200 lateral cross sections at different depths. 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.




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 expansion52 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 approaches5355 may be explored in the future.


The authors thank Anne Sentenac, Alex Matlock, and Benedict Diederich for helpful discussions regarding multiple scattering, Yujia Xue for discussions on the simulations, and Quentin Lin for help on the experiment. Disclosures: there are no conflicts of interest related to this article. Funding: National Science Foundation (NSF) (Grant Nos. 1813848 and 1813910).



I. Moon et al., “Automated three-dimensional identification and tracking of micro/nanobiological organisms by computational holographic microscopy,” Proc. IEEE, 97 (6), 990 –1010 (2009). IEEPAD 0018-9219 Google Scholar


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. U. S. A., 109 (40), 16018 –16022 (2012). Google Scholar


F. C. Cheong et al., “Flow visualization and flow cytometry with holographic video microscopy,” Opt. Express, 17 (15), 13071 –13079 (2009). OPEXFF 1094-4087 Google Scholar


J. Garcia-Sucerquia et al., “Digital in-line holographic microscopy,” Appl. Opt., 45 (5), 836 –850 (2006). APOPAI 0003-6935 Google Scholar


J. Sheng, E. Malkiel and J. Katz, “Digital holographic microscope for measuring three-dimensional particle distributions and motions,” Appl. Opt., 45 (16), 3893 –3901 (2006). APOPAI 0003-6935 Google Scholar


L. Tian et al., “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). APOPAI 0003-6935 Google Scholar


P. Picart et al., “Tracking high amplitude auto-oscillations with digital Fresnel holograms,” Opt. Express, 15 (13), 8263 –8274 (2007). OPEXFF 1094-4087 Google Scholar


U. Schnars and W. P. Jüptner, “Digital recording and numerical reconstruction of holograms,” Meas. Sci. Technol., 13 (9), R85 (2002). MSTCEP 0957-0233 Google Scholar


F. Soulez et al., “Inverse-problem approach for particle digital holography: accurate location based on local optimization,” J. Opt. Soc. Am. A, 24 (4), 1164 –1171 (2007). JOAOD6 0740-3232 Google Scholar


W. Chen et al., “Empirical concentration bounds for compressive holographic bubble imaging based on a Mie scattering model,” Opt. Express, 23 (4), 4715 –4725 (2015). OPEXFF 1094-4087 Google Scholar


D. J. Brady et al., “Compressive holography,” Opt. Express, 17 (15), 13040 –13049 (2009). OPEXFF 1094-4087 Google Scholar


Y. Rivenson, A. Stern and B. Javidi, “Compressive Fresnel holography,” J. Disp. Technol., 6 (10), 506 –509 (2010). IJDTAL 1551-319X Google Scholar


M. Born and E. Wolf, “Scattering from inhomogeneous media,” Principles of Optics, 695 –734 7 ed.Cambridge University Press, Cambridge (2003). Google Scholar


W. C. Chew and Y. M. Wang, “Reconstruction of two-dimensional permittivity distribution using the distorted Born iterative method,” IEEE Trans. Med. Imaging, 9 218 –225 (1990). ITMID4 0278-0062 Google Scholar


F.-C. Chen and W. C. Chew, “Experimental verification of super resolution in nonlinear inverse scattering,” Appl. Phys. Lett., 72 (23), 3080 –3082 (1998). APPLAB 0003-6951 Google Scholar


U. S. Kamilov et al., “A recursive Born approach to nonlinear inverse scattering,” IEEE Signal Process. Lett., 23 1052 –1056 (2016). IESPEJ 1070-9908 Google Scholar


G. A. Tsihrintzis and A. J. Devaney, “Higher-order (nonlinear) diffraction tomography: reconstruction algorithms and computer simulation,” IEEE Trans. Image Process., 9 (9), 1560 –1572 (2000). IIPRE4 1057-7149 Google Scholar


R. Pierri and G. Leone, “Inverse scattering of dielectric cylinders by a second-order born approximation,” IEEE Trans. Geosci. Remote Sens., 37 (1), 374 –382 (1999). IGRSD2 0196-2892 Google Scholar


P. M. van den Berg and R. E. Kleinman, “A contrast source inversion method,” Inverse Prob., 13 1607 –1620 (1997). INPEEY 0266-5611 Google Scholar


M. T. Bevacquad et al., “Non-linear inverse scattering via sparsity regularized contrast source inversion,” IEEE Trans. Comput. Imaging, 3 (2), 296 –304 (2017). Google Scholar


P. M. van den Berg, A. Van Broekhoven and A. Abubakar, “Extended contrast source inversion,” Inverse Prob., 15 (5), 1325 –1344 (1999). INPEEY 0266-5611 Google Scholar


P. Van den Berg and A. Abubakar, “Contrast source inversion method: state of art,” Prog. Electromagn. Res., 34 (11), 189 –218 (2001). PELREX 1043-626X Google Scholar


R. E. Kleinman and P. M. van den Berg, “A modified gradient method for two-dimensional problems in tomography,” J. Comput. Appl. Math., 42 (1), 17 –35 (1992). JCAMDI 0377-0427 Google Scholar


R. E. Kleinman and P. M. van den Berg, “An extended range-modified gradient technique for profile inversion,” Radio Sci., 28 877 –884 (1993). RASCAD 0048-6604 Google Scholar


H.-Y. Liu et al., “SEAGLE: sparsity-driven image reconstruction under multiple scattering,” IEEE Trans. Comput. Imaging, 4 (1), 73 –86 (2018). Google Scholar


Y. Ma et al., “Accelerated image reconstruction for nonlinear diffractive imaging,” 6473 –6477 (2018). Google Scholar


T.-A. Pham et al., “Versatile reconstruction framework for diffraction tomography with intensity measurements and multiple scattering,” Opt. Express, 26 (3), 2749 –2763 (2018). OPEXFF 1094-4087 Google Scholar


J. Lim et al., “Beyond Born–Rytov limit for super-resolution optical diffraction tomography,” Opt. Express, 25 (24), 30445 –30458 (2017). OPEXFF 1094-4087 Google Scholar


F. Simonetti, “Multiple scattering: the key to unravel the subwavelength world from the far-field pattern of scattered wave,” Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 73 036619 (2006). Google Scholar


E. Mudry et al., “Electromagnetic wave imaging of three-dimensional targets using a hybrid iterative inversion method,” Inverse Prob., 28 065007 (2012). INPEEY 0266-5611 Google Scholar


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). Google Scholar


U. S. Kamilov et al., “Optical tomographic image reconstruction based on beam propagation and sparse regularization,” IEEE Trans. Comput. Imaging, 2 (1), 59 –70 (2016). Google Scholar


U. S. Kamilov et al., “Learning approach to optical tomography,” Optica, 2 (6), 517 –522 (2015). Google Scholar


M. Guizar-Sicairos and J. R. Fienup, “Understanding the twin-image problem in phase retrieval,” J. Opt. Soc. Am. A, 29 (11), 2367 –2375 (2012). JOAOD6 0740-3232 Google Scholar


A. Beck, M. Teboulle, “Gradient-based algorithms with applications to signal-recovery problems,” Convex Optimization in Signal Processing and Communications, 42 –88 Cambridge University Press, Cambridge (2009). Google Scholar


W. Zhang et al., “Twin-image-free holography: a compressive sensing approach,” Phys. Rev. Lett., 121 (9), 093902 (2018). PRLTAO 0031-9007 Google Scholar


R. Ling et al., “High-throughput intensity diffraction tomography with a computational microscope,” Biomed. Opt. Express, 9 (5), 2130 –2141 (2018). BOEICL 2156-7085 Google Scholar


J. Lim et al., “Learning tomography assessed using Mie theory,” Phys. Rev. Appl., 9 (3), 034027 (2018). PRAHB2 2331-7019 Google Scholar


L. Rokach and O. Maimon, Clustering Methods, 321 –352 Springer US, Boston, Massachusetts (2005). Google Scholar


D. Mas et al., “Fast algorithms for free-space diffraction patterns calculation,” Opt. Commun., 164 (4–6), 233 –245 (1999). OPCOB8 0030-4018 Google Scholar


A. D. Yaghjian, “Electric dyadic Green’s functions in the source region,” Proc. IEEE, 68 (2), 248 –263 (1980). IEEPAD 0018-9219 Google Scholar


J. Van Bladel, “Some remarks on Green’s dyadic for infinite space,” IRE Trans. Antennas Propag., 9 (6), 563 –566 (1961). IRAPAG 0096-1973 Google Scholar


A. Fannjiang, “Tv-min and greedy pursuit for constrained joint sparsity and application to inverse scattering,” Math. Mech. Complex Syst., 1 (1), 81 –104 (2013). Google Scholar


A. Fannjiang, “Compressive sensing theory for optical systems described by a continuous model,” (2015) Google Scholar


A. Beck and M. Teboulle, “Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems,” IEEE Trans. Image Process., 18 (11), 2419 –2434 (2009). IIPRE4 1057-7149 Google Scholar


U. S. Kamilov, “A parallel proximal algorithm for anisotropic total variation minimization,” IEEE Trans. Image Process., 26 (2), 539 –548 (2017). IIPRE4 1057-7149 Google Scholar


N. Andrei, “An acceleration of gradient descent algorithm with backtracking for unconstrained optimization,” Numer. Algorithms, 42 (1), 63 –73 (2006). NUALEG 1017-1398 Google Scholar


M. Azimi and A. Kak, “Distortion in diffraction tomography caused by multiple scattering,” IEEE Trans. Med. Imaging, 2 (4), 176 –195 (1983). ITMID4 0278-0062 Google Scholar


S. Bolte and F. Cordelieres, “A guided tour into subcellular colocalization analysis in light microscopy,” J. Microsc., 224 (3), 213 –232 (2006). JMICAR 0022-2720 Google Scholar


A. Ishimaru, Wave Propagation and Scattering in Random Media, 2 Academic Press, New York (1978). Google Scholar


G. Osnabrugge, S. Leedumrongwatthanakun and I. M. Vellekoop, “A convergent born series for solving the inhomogeneous Helmholtz equation in arbitrarily large media,” J. Comput. Phys., 322 113 –124 (2016). JCTPAH 0021-9991 Google Scholar


Y. Sun, Z. Xia and U. S. Kamilov, “Efficient and accurate inversion of multiple scattering with deep learning,” Opt. Express, 26 (11), 14678 –14688 (2018). OPEXFF 1094-4087 Google Scholar


S. Li et al., “Imaging through glass diffusers using densely connected convolutional networks,” Optica, 5 (7), 803 –813 (2018). Google Scholar


Y. Li, Y. Xue and L. Tian, “Deep speckle correlation: a deep learning approach toward scalable imaging through scattering media,” Optica, 5 (10), 1181 –1190 (2018). Google Scholar


Waleed Tahir is a PhD student working under supervision of Lei Tian at the Computational Imaging Systems Lab in the Electrical and Computer Engineering Department at Boston University. He received his BE degree in 2013 from the National University of Sciences and Technology, Pakistan, after which he worked as a design engineer at Burqstream Technologies, Pakistan from 2014 to 2016. His research is geared toward computational imaging and machine learning for biomedical imaging.

Ulugbek S. Kamilov is a director of the Computational Imaging Group at Washington University in St. Louis. He received his BSc and MSc degrees in communication systems, and his PhD in electrical engineering from EPFL, Switzerland, in 2008, 2011, and 2015, respectively. From 2015 to 2017, he was a research scientist at Mitsubishi Electric Research Laboratories, Cambridge, Massachusetts, USA. His main area of research is computational imaging, with an emphasis on the mathematical and computational aspects of image reconstruction and analysis.

Lei Tian directs the Computational Imaging Systems Lab in the Electrical and Computer Engineering Department at Boston University. He received his PhD in 2013 and MS in 2010 from Massachusetts Institute of Technology. He was a postdoctoral associate in the EECS Department at University of California, Berkeley, from 2013 to 2016. His research focuses on computational microscopy, neurophotonics, imaging in complex media, and machine learning for biomedical microscopy.

© The Authors. Published by SPIE and CLP under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Waleed Tahir, Ulugbek S. Kamilov, and Lei Tian "Holographic particle localization under multiple scattering," Advanced Photonics 1(3), 036003 (18 June 2019).
Received: 8 March 2019; Accepted: 20 May 2019; Published: 18 June 2019

Back to Top