Holographic particle localization under multiple scattering

Waleed Tahir; Ulugbek S. Kamilov; Lei Tian

doi:10.1117/1.AP.1.3.036003

Abstract

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.

Keywords

digital holography multiple scattering particle localization

1 Introduction

Three-dimensional (3-D) particle-localization using in-line holography is fundamental to many applications, such as biological sample characterization,¹^,² flow cytometry,³^,⁴ fluid mechanics,⁵^,⁶ and optical measurement.⁷^–⁹ 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.¹¹^,¹² 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,¹⁴^–¹⁸ contrast source inversion,¹⁹^–²² modified gradient,²³^,²⁴ series expansion with accelerated gradient descent on Lippmann–Schwinger equation (SEAGLE),²⁵^,²⁶ and hybrid methods.²⁷^–³⁰ 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).³¹^–³³ However, our model has the benefit of computing both forward and backward scattering, whereas the latter only accounts for forward multiple scattering.

Sign up for Advanced Photonics TOC. Get the latest issue of Advanced Photonics delivered right to you！Sign up now

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.¹⁹^,²⁵ 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,¹⁴^–¹⁸ 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¹¹^,³⁶ 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.¹⁷^,²⁵^,³⁷^,³⁸ 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.

2 Theory and Method

2.1 Forward Model

Consider the imaging geometry in Fig. 1(a). An in-line hologram $I_{m}$ at the measurement plane can be written as $I_{m} (x_{0}) = {| u_{in} (x_{0}, 0) + E (x_{0}, 0) |}^{2} = {| u_{in} |}^{2} + 2 u_{in} Re {E} + {| E |}^{2},$ (1)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, and $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. The scattered field and its “twin-image” contribution is related to the measured hologram after background removal by $Re {E (x_{0}, 0)} = \frac{[E (x_{0}, 0) + E^{*} (x_{0}, 0)]}{2} = \frac{I_{m} (x_{0}) - {| u_{in} |}^{2}}{2 u_{in}} .$ (2)

Figure 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.⁴⁰ (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 $I_{br} (x_{0}) = [I_{m} (x_{0}) - {| u_{in} |}^{2}] / 2 u_{in}$ . To model the hologram resulting from multiple scattering up-to the $K$ ’th order, we apply the framework of Born series expansion,¹³^,¹⁶ which gives two coupled equations: $E (x_{0}, 0) = ∭_{Ω} f (x^{'}, z^{'}) u_{K} (x^{'}, z^{'}) h (x_{0} - x^{'}, - z^{'}) d^{2} x^{'} d z^{'},$ (3) $u_{k} (x, z) = u_{in} (x, z) + ∭_{Ω} f (x^{'}, z^{'}) u_{k - 1} (x^{'}, z^{'}) h (x - x^{'}, z - z^{'}) d^{2} x^{'} d z^{'},$ (4)where $h$ is the 3-D Green’s function and $u_{K} (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^{'}) = \frac{k^{2}}{4 π} [n^{2} (x^{'}, z^{'}) - n_{med}^{2}]$ ,¹³ 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 3-D support, $(x, z) \in Ω$ , $(x^{'}, z^{'}) \in Ω$ , whereas the hologram is measured outside the support, $(x_{0}, z_{0}) \notin Ω$ . To compute higher order scattering, the initial condition is $u_{0} (x, z) = 0$ , and $k = 1,2, \dots, 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: $E = H diag (u_{K}) f,$ (5) $u_{k} = u_{in} + G diag (u_{k - 1}) f,$ (6)for $k = 1,2, \dots, K$ . Here, $f$ and $u_{k}$ have dimension $(N_{x} N_{y} N_{z}) \times 1$ , and $E$ is $(N_{x} N_{y}) \times 1$ . $N_{x}$ , $N_{y}$ , and $N_{z}$ 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 = K^{H} Q_{0} B$ , where $B = bldiag (K, .., K)$ is a block-diagonal matrix. $K$ and $K^{H}$ are the 2-D DFT and the inverse 2-D DFT matrices, respectively, each with dimension $(N_{x} N_{y}) \times (N_{x} N_{y})$ . $Q_{0} = [\begin{matrix} L_{0 - 1} & L_{0 - 2} & \dots & L_{0 - N_{z}} \end{matrix}]$ , 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 Ref. 11. $G$ is the multiple scattering operator that performs propagation from the object volume to within itself (Fig. 2). $G = B^{H} RB$ , where $R = [Q_{1}, Q_{2}, \dots, Q_{N_{z}}]$ and $Q_{m} = [\begin{matrix} L_{m - 1} & L_{m - 2} & \dots & L_{m - N_{z}} \end{matrix}]$ . $R$ has the dimension of $(N_{x} N_{y} N_{z}) \times (N_{x} N_{y} N_{z})$ 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 dimensions $(N_{x} N_{y}) \times (N_{x} N_{y})$ , the direct and the angular spectrum methods.⁴¹ We use the direct method, in which the Green’s function $h (r) = \exp (\dot{ι} k | r |) / | r |$ is sampled in the spatial domain, followed by slicewise 2-D FFT, where $r = (x, y, z)$ .

Figure 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 $u_{k - 1}$ ; it is then propagated to every other slice within the volume. (b) This computed scattered-field $u_{k}^{s}$ 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.

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.

2.2 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 algorithm⁴⁴^,⁴⁵ that minimizes a cost function $C (f)$ to compute the estimated object $\hat{f}$ as follows: $\hat{f} = \underset{f \in F}{argmin} {C (f)} = \underset{f \in F}{argmin} {D (f) + τ {‖ f ‖}_{TV}},$ (7)where $D (f) = \frac{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)], $‖ \cdot ‖$ 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} \overset{Δ}{=} \sum_{n = 1}^{N} {‖ {[Df]}_{n} ‖}_{l_{2}} = \sum_{n = 1}^{N} \sqrt{{| {[D_{x} f]}_{n} |}^{2} + {| {[D_{y} f]}_{n} |}^{2} + {| {[D_{z} f]}_{n} |}^{2}},$ (8)where $D : R^{N} \to R^{N \times 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,⁴⁶ in which the $t$ ’th iteration is written as $f^{t} \leftarrow {prox}_{τ TV} [f^{t - 1} - α \frac{\partial D (f)}{\partial f}],$ (9)where ${prox}_{τ TV} (g) \overset{Δ}{=} \underset{f \in F}{argmin} {\frac{1}{2} {‖ 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 $f^{0} = 0$ .

Similar to the forward model, the gradient computation is also a $K$ -order recursion: $\frac{\partial D (f)}{\partial f} = Re {diag ({\bar{u}}_{K}) H^{H} r + {[\frac{\partial u_{K}}{\partial f}]}^{H} diag (f) H^{H} r},$ (10) ${[\frac{\partial u_{k}}{\partial f}]}^{H} a = diag ({\bar{u}}_{k - 1}) G^{H} a + {[\frac{\partial u_{k - 1}}{\partial f}]}^{H} diag (f) G^{H} a,$ (11)for $k = 1,2, \dots, K$ . Here, $r = Re {E_{est}} - I_{br}$ is the residual; $A^{H}$ and $\bar{A}$ represent the Hermitian and complex conjugate of the matrix $A$ , respectively. The recursion is initialized with $\partial u_{0} / \partial 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.

3 Results

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 \times$ 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 \pm 0.021 μ m$ (Thermofisher Scientific 4009A) suspended in deionized water. The suspension is held in a quartz-cuvette with inner dimensions $40 mm \times 40 mm \times 0.5 mm$ . 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 / n_{water} = 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 $λ / (1 - \sqrt{1 - {NA}^{2}}) = 5.7 μ m$ . During the computation, $2 \times$ 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 \times 176 μ m \times 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 \geq 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 \times 512 \times 50$ grid were 58 and 257 s, respectively. All reconstructions were run for 100 iterations. In what follows, we present our findings.

3.1 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.¹⁷^,²⁵^,³⁷ 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 \times 44 μ m \times 6.4 μ m$ , discretized as a $256 \times 256 \times 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 \times 6.2 μ m \times 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.

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

3.2 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 \times 88 μ m \times 250 μ m$ , discretized on a $512 \times 512 \times 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 $R_{g}$ , which corresponds to the fraction of the hologram area directly occluded by the scatterers, defined as $R_{g} = \frac{total cross-sections of all scatterers}{area of the hologram} \approx \frac{N_{p} π r^{2}}{N_{x} N_{y} δ x δ y},$ (12)where $N_{p}$ 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 $R_{g}$ , we thus expect greater contributions of multiple scattering. From the signal processing perspective, $R_{g}$ 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 $R_{g}$ tested are 0.01, 0.02, 0.05, 0.1, 0.2, and 0.4, corresponding to $N_{p} = 100$ , 200, 500, 1000, 2000, 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 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} \leq 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.

Figure 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 $R_{g}$ tested in the simulation. For $R_{g} \leq 0.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,¹⁹^,²⁵^,³⁰ as $e = ‖ {Au}_{K} - u_{0} ‖,$ (13)where $A = I - G diag (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.

$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 (Rg≤0.02), and better than single-scattering for 0.02<Rg≤0.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.$

Figure 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 ( $R_{g} \leq 0.02$ ), and better than single-scattering for $0.02 < R_{g} \leq 0.1$ . Reconstruction fails for very high concentration ( $R_{g} > 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 $R_{g} = 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) $SNR = 20 \log_{10} (\frac{‖ f_{true} ‖}{‖ f_{true} - \hat{f} ‖}),$ (14)where $f_{true}$ and $\hat{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 ( $R_{g} \leq 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} \leq 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 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 = \frac{2 \sum_{i = 1}^{N} p_{i} g_{i}}{\sum_{i = 1}^{N} p_{i}^{2} + \sum_{i = 1}^{N} g_{i}^{2}},$ (15)where $p_{i}$ and $g_{i}$ 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)].

Figure 6.Reconstruction performance as a function of depth. (a) Segmentation maps of reconstructed slices (zoomed-in $51 μ m \times 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.

3.3 Large-Scale Experimental Validation

Finally, we demonstrate our method on a set of large-scale experiments. We reconstruct over 100 million object voxels ( $1024 \times 1024 \times 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 $R_{g}$ 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 $R_{g}$ , depicting degradation of reconstruction with increase in particle density.

Figure 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 $R_{g} \leq 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 CR.

Next, we assess the estimated particle density. Our multiple-scattering model consistently performs better than the single-scattering model for $R_{g} \leq 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.⁵¹ At the critical $R_{g} = 0.1$ concentration, the CR is around 0.335.

Finally, we closely examine the 3-D 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.

Figure 8.A 3-D visualization of the localized particles under different concentrations from our experiment and their $200 \times 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.

4 Conclusion

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⁵³^–⁵⁵ may be explored in the future.

References

[1] I. Moon et al. Automated three-dimensional identification and tracking of micro/nanobiological organisms by computational holographic microscopy. Proc. IEEE, 97, 990-1010(2009).

[2] T.-W. Su, L. Xue, A. Ozcan. High-throughput lensfree 3D tracking of human sperms reveals rare statistics of helical trajectories. Proc. Natl. Acad. Sci. U. S. A., 109, 16018-16022(2012).

[3] F. C. Cheong et al. Flow visualization and flow cytometry with holographic video microscopy. Opt. Express, 17, 13071-13079(2009).

[4] J. Garcia-Sucerquia et al. Digital in-line holographic microscopy. Appl. Opt., 45, 836-850(2006).

[5] J. Sheng, E. Malkiel, J. Katz. Digital holographic microscope for measuring three-dimensional particle distributions and motions. Appl. Opt., 45, 3893-3901(2006).

[6] 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, 1549-1554(2010).

[7] P. Picart et al. Tracking high amplitude auto-oscillations with digital Fresnel holograms. Opt. Express, 15, 8263-8274(2007).

[8] U. Schnars, W. P. Jüptner. Digital recording and numerical reconstruction of holograms. Meas. Sci. Technol., 13, R85(2002).

[9] F. Soulez et al. Inverse-problem approach for particle digital holography: accurate location based on local optimization. J. Opt. Soc. Am. A, 24, 1164-1171(2007).

[10] W. Chen et al. Empirical concentration bounds for compressive holographic bubble imaging based on a Mie scattering model. Opt. Express, 23, 4715-4725(2015).

[11] D. J. Brady et al. Compressive holography. Opt. Express, 17, 13040-13049(2009).

[12] Y. Rivenson, A. Stern, B. Javidi. Compressive Fresnel holography. J. Disp. Technol., 6, 506-509(2010).

[13] M. Born, E. Wolf. Scattering from inhomogeneous media. Principles of Optics, 695-734(2003).

[14] W. C. Chew, Y. M. Wang. Reconstruction of two-dimensional permittivity distribution using the distorted Born iterative method. IEEE Trans. Med. Imaging, 9, 218-225(1990).

[15] F.-C. Chen, W. C. Chew. Experimental verification of super resolution in nonlinear inverse scattering. Appl. Phys. Lett., 72, 3080-3082(1998).

[16] U. S. Kamilov et al. A recursive Born approach to nonlinear inverse scattering. IEEE Signal Process. Lett., 23, 1052-1056(2016).

[17] G. A. Tsihrintzis, A. J. Devaney. Higher-order (nonlinear) diffraction tomography: reconstruction algorithms and computer simulation. IEEE Trans. Image Process., 9, 1560-1572(2000).

[18] R. Pierri, G. Leone. Inverse scattering of dielectric cylinders by a second-order born approximation. IEEE Trans. Geosci. Remote Sens., 37, 374-382(1999).

[19] P. M. van den Berg, R. E. Kleinman. A contrast source inversion method. Inverse Prob., 13, 1607-1620(1997).

[20] M. T. Bevacquad et al. Non-linear inverse scattering via sparsity regularized contrast source inversion. IEEE Trans. Comput. Imaging, 3, 296-304(2017).

[21] P. M. van den Berg, A. Van Broekhoven, A. Abubakar. Extended contrast source inversion. Inverse Prob., 15, 1325-1344(1999).

[22] P. Van den Berg, A. Abubakar. Contrast source inversion method: state of art. Prog. Electromagn. Res., 34, 189-218(2001).

[23] R. E. Kleinman, P. M. van den Berg. A modified gradient method for two-dimensional problems in tomography. J. Comput. Appl. Math., 42, 17-35(1992).

[24] R. E. Kleinman, P. M. van den Berg. An extended range-modified gradient technique for profile inversion. Radio Sci., 28, 877-884(1993).

[25] H.-Y. Liu et al. SEAGLE: sparsity-driven image reconstruction under multiple scattering. IEEE Trans. Comput. Imaging, 4, 73-86(2018).

[26] Y. Ma et al. Accelerated image reconstruction for nonlinear diffractive imaging. IEEE Int. Conf. on Acoustics, Speech, and Signal Processing (ICASSP), 6473-6477(2018).

[27] T.-A. Pham et al. Versatile reconstruction framework for diffraction tomography with intensity measurements and multiple scattering. Opt. Express, 26, 2749-2763(2018).

[28] J. Lim et al. Beyond Born–Rytov limit for super-resolution optical diffraction tomography. Opt. Express, 25, 30445-30458(2017).

[29] 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).

[30] E. Mudry et al. Electromagnetic wave imaging of three-dimensional targets using a hybrid iterative inversion method. Inverse Prob., 28, 065007(2012).

[31] L. Tian, L. Waller. 3D intensity and phase imaging from light field measurements in an LED array microscope. Optica, 2, 104-111(2015).

[32] U. S. Kamilov et al. Optical tomographic image reconstruction based on beam propagation and sparse regularization. IEEE Trans. Comput. Imaging, 2, 59-70(2016).

[33] U. S. Kamilov et al. Learning approach to optical tomography. Optica, 2, 517-522(2015).

[34] M. Guizar-Sicairos, J. R. Fienup. Understanding the twin-image problem in phase retrieval. J. Opt. Soc. Am. A, 29, 2367-2375(2012).

[35] A. Beck, D. P. Palomar, Y. C. Eldar, M. Teboulle. Gradient-based algorithms with applications to signal-recovery problems. Convex Optimization in Signal Processing and Communications, 42-88(2009).

[36] W. Zhang et al. Twin-image-free holography: a compressive sensing approach. Phys. Rev. Lett., 121, 093902(2018).

[37] R. Ling et al. High-throughput intensity diffraction tomography with a computational microscope. Biomed. Opt. Express, 9, 2130-2141(2018).

[38] J. Lim et al. Learning tomography assessed using Mie theory. Phys. Rev. Appl., 9, 034027(2018).

[39] L. Rokach, O. Maimon. Clustering Methods, 321-352(2005).

[40]

[41] D. Mas et al. Fast algorithms for free-space diffraction patterns calculation. Opt. Commun., 164, 233-245(1999).

[42] A. D. Yaghjian. Electric dyadic Green’s functions in the source region. Proc. IEEE, 68, 248-263(1980).

[43] J. Van Bladel. Some remarks on Green’s dyadic for infinite space. IRE Trans. Antennas Propag., 9, 563-566(1961).

[44] A. Fannjiang. Tv-min and greedy pursuit for constrained joint sparsity and application to inverse scattering. Math. Mech. Complex Syst., 1, 81-104(2013).

[45] A. Fannjiang. Compressive sensing theory for optical systems described by a continuous model(2015).

[46] A. Beck, M. Teboulle. Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems. IEEE Trans. Image Process., 18, 2419-2434(2009).

[47] U. S. Kamilov. A parallel proximal algorithm for anisotropic total variation minimization. IEEE Trans. Image Process., 26, 539-548(2017).