Single-shot holographic particle 3D imaging is fundamental to many important applications, such as flow cytometry [1], biological sample characterization [2,3], and flow measurement [4]. The technique works by first recording a 2D intensity measurement from a particle volume under coherent illumination and then estimating the 3D refractive index distribution. The problem is challenging because of the dimensional mismatch between the single 2D measurement and the unknown 3D object and the "phaseless" measurement, both of which make the inverse problem severely ill-posed, especially for large-scale problems. Furthermore, multiple scattering effects become significant as the particle density and the refractive index contrast increase, which necessitates a nonlinear forward model to accurately describe the image formation process.
Traditional holographic particle 3D reconstruction is based on the linear single-scattering model derived from the first Born approximation method (FBM) [5]. The most widely used algorithm is known as the "back-propagation method" [4], which is equivalent to apply the pseudo-inverse of the (3D-to-2D) single-scattering forward operator to the captured hologram. To improve the reconstruction accuracy, compressive holography [6] has been developed that supplements the FBM forward model with a sparsity constraint and solves the 3D reconstruction problem by an iterative optimization procedure. This approach is particularly effective in alleviating the twin-image artifacts arising [6] and improving the robustness to unaccounted multiple-scattering effects at high scattering densities [7]. However, these methods are fundamentally limited by the underlying single-scattering assumption, which is valid only for weakly scattering objects. The model gradually breaks down as the particle density and/or the refractive index contrast increase. The multiple-scattering effects result in a model mismatch, which in turn limits the reconstruction accuracy by these techniques.
Several multiple-scattering models have been developed, such as iterative Born series [8,9], contrast source inversion [10], discrete-dipole approximation [11][12][13][14], and series expansion with accelerated gradient descent on Lippmann-Schwinger equation [15]. However, computational challenges typically restrict these methods to be used only for small-scale problems. This is because they involve iteratively estimating the internally scattered fields within the object volume for modeling the multiple-scattering process, which incurs a computation cost that grows rapidly as the object size, and thus is not ideal for our intended applications containing on the order of 100 million voxels in the object volume. Recently, the multi-slice beam propagation method (BPM) [16][17][18][19] has emerged as an attractive, computationally efficient multiple-scattering model. The utility of the BPM has been demonstrated for reconstructing 3D refractive index distributions using tomographic measurements from multiple interferometric complex-field measurements [16] or intensity-only measurements [17][18][19]. Our goal here is to critically examine the utility of the BPM to handle the inversion from a single in-line hologram, which inherently suffers from more severe missing-cone artifacts [20,21] and greater dimensional mismatch as compared to the previous tomography studies.
We demonstrate a novel 3D reconstruction algorithm that combines a BPM multiple-scattering forward model and a sparsity prior. The accuracy of our forward model is quantified by comparing the simulation with the experiments. We show that the BPM-computed intensity statistics closely match with the experimental data at different scattering densities, and provide up to 9× higher accuracy in predicting the hologram contrast than that from the FBM. Benefited from the high computational efficiency, the BPM reduces the computational complexity by N z times (where N z is the number of axial slices of the object volume), as compared to the state-of-the-art multiple-scattering model-based holographic particle reconstruction method [9]. We demonstrate our proposed method's superior 3D imaging performance in both simulation and experiments under different scattering densities and refractive index contrast. To quantify the improvement by the multiple-scattering model over the single-scattering model, we compare our method with the compressive holography method. We show that the BPM reconstruction significantly outperforms the FBM in particular at deep imaging depths and for high particle densities.
Our reconstruction algorithm solves an ℓ 1 -regularized least-squares problem, in which the data fidelity term measures the difference between the BPM-estimated and captured holograms. Below, we describe our forward model and the reconstruction algorithm.
Our imaging geometry is shown in Fig. 1(a). A plane wave passes through a particle volume and the interference pattern resulting from the scattered fields are captured as the hologram. This multiple-scattering image-formation process is modeled as,
where Î ∈ R M denotes the vectorized, BPM computed intensity hologram containing M pixels. w ∈ R N = nn 0 represents the vectorized 3D refractive index contrast distribution that contains N voxels, and is defined by the difference between the particle's index distribution n and the constant background index n 0 . The particle's refractive index is assumed to be real-valued since the absorption is negligible. S : R N → C M is the BPM operator that maps the 3D refractive index to the 2D complex field at the sensor plane. The BPM is computed by a recursive sequence, as illustrated in Fig. 1(b). It approximates the 3D object as N z infinitesimally thin axial planar slices along the optical axis, each modeled as a 2D phase mask separated equally by a uniform medium. The field is then computed by a sequence of diffraction and refraction operations. Specifically, the field exiting the jth slice S j (w) is related to the field from the (j -1)th slice S j-1 (w) by S j (w) = diag(p j (w j ))HS j-1 (w).
(
where the S j : R M → C M is the local BPM operator that maps the jth phase screen w j to the 2D complex field after the jth slice. The implementation of the BPM requires an initial condition, for which we set the initial field as a constant-valued on-axis plane wave S 0 (w) = 1.
The diffraction operator for a propagation distance ∆z is denoted by H and computed by H = F H diag(h)F, where F and F H are the 2D discrete Fourier and inverse Fourier transform matrices respectively, • H denotes matrix Hermitian transpose, diag(h) is a diagonal matrix formed by the vector h as the main diagonal and is implemented by element-wise multiplication, and h represents the vectorized 2D angular spectrum transfer function h(p, q):
where p, q index the 2D wave vector in the k-space, k = k 0 n 0 is the wave-number in the background medium, k 0 = 2π/λ 0 , λ 0 is the wavelength in vacuum, and ∆k is the sampling interval in the k-space. P represents the low-pass filter due to the evanescent cutoff.
The refraction operator is implemented by element-wise multiplications between the diffracted field and the accumulated phase p j (w j ) = e jk 0 ∆zw j , where w j denotes the discretized 2D refractive index map in the jth slice.
The hologram is the intensity of the exit field from the last slice at N z low-pass filtered by the pupil function of the imaging system.
The computation of the BPM thus requires implementing Eq. ( 2) N z times. The computational complexity of the BPM is O(Nlog(N x N y )), as set approximately by 2N z times evaluations of 2D FFT, where the total number of object voxels is N = N x N y N z , and N x and N y are the number of pixels in the x and y directions, respectively. For comparison, the computational complexity of the Born series model is O(N z Nlog(N x N y )) [9]. Thus, the BPM reduces the computational complexity by N z times, which becomes significant when the object depth is large.
To illustrate the multiple scattering effects modeled by the BPM, we show an example yz cross-section of the intensity distribution computed from a particle volume in Fig. 1(b). In the zoom-in region, complex interference patterns are visible due to multiple particles in the close proximity, which introduces strong multiple scattering effects [22].
We formulate the 3D reconstruction as a minimization problem:
where D is the data fidelity term and R is the regularization term. The parameter τ>0 controls the amount of regularization. The data fidelity term is given by
where I is the captured hologram and ∥ • ∥ 2 is the ℓ 2 norm. The regularization term is the ℓ 1 norm of the refractive index distribution
which promotes the sparsity of the reconstructed object. The minimization Eq. ( 4) is not a trivial task. The primary difficulty stems from that the data fidelity term D involves a nonlinear forward operator S and the regularization term R is non-smooth. We next present a novel algorithm based on the proximal-gradient descend technique, which extends [16] to intensity-only measurement.
The crucial step is the gradient computation for D with respect to w, as summarized in Algorithm 1 and briefly explained below. Additional details are provided in Supplement 1.
First, we take the derivative of D(w) with respect to w and rearrange it into a column vector
where r ≜ Î -I is the residual between the estimated and captured holograms. The estimated measurement can be written as
, where • denotes taking the complex conjugate. Thus, the gradient of the data fidelity term is
To derive a tractable algorithm for computing Eq. ( 8), we apply the recursive BPM formula in Eq. ( 2) and compute the local gradient at the jth slice:
Since the input plane wave does not depend on w, so for the initial condition j = 0, we have
Based on Eq. ( 9) and the initial condition in Eq. ( 10), we obtain a practical implementation of Eq. ( 8) to calculate the gradient of the data fidelity term, as summarized in Algorithm 1.
]︁ H r N z using the following iterative procedure for slices
Intuitively, the gradient computation is similar to the "error backpropagation" algorithm used in deep neural networks [16]. The algorithm iterates between two major steps, including update the intermediate field gradient slice-by-slice (first term in Eq. ( 9)) and backpropagate the slice-wise field residual (second term in Eq. ( 9)). Since this gradient computation takes the same recursive procedure as the forward BPM model, its computation complexity is also O(Nlog(N x N y )).
Our algorithm is summarized in Algorithm 2, that reconstructs the refractive index w from a single hologram based on the proximal gradient algorithm. Conceptually, this algorithm is similar to the fast iterative shrinkage/thresholding algorithm (FISTA) [23], which is widely used to minimize objective function that consist of the sums between a smooth and a non-smooth term.
A major component of this algorithm is the proximal operator for the ℓ 1 -regularizer
where a and γ are explained in Algorithm 2. This proximal operator has a closed-form solution, known as soft-thresholding, see Supplement 1.
The parameters in Algorithm 2 are set as follows. We set the initial guess ŵ0 = 0, and the step size γ = 5 × 10 -6 . The stopping criterion is the maximum iteration number to be 200-300. The regularization parameter τ is tuned under different scattering conditions. When scattering is stronger, τ is set larger. The reconstruction is found to be insensitive to the fine-tuning of τ.
)︁ s t ← ŵt + ((q t-1 -1)/q t )( ŵtŵt-1 ) t ← t + 1 Until stopping criterion Return estimate of the refractive index ŵt
We first validate the forward model accuracy by comparing the BPM-computed holograms with experimental measurements. In practice, we assess the accuracy based on analyzing the intensity statistics under different scattering conditions, since the ground-truth particle positions in the experiments are not known.
We perform simulations using parameters that match with the experiments. More details about the experiments are provided in Supplement 1. Holograms are computed at four particle densities ρ, including 1.60, 3.20, 6.41, 12.82 × 10 4 /µL, corresponding to 250, 500, 1000, 2000 particles in a 176.64 × 176.64 × 500 µm 3 volume. 1 µm particles are randomly placed in 3D positions. The refractive index of the particle n and the background medium (water) n 0 is 1.59 and 1.33, respectively, and the contrast is ∆n = 0.26. To simulate 3D refractive index contrast distribution w, the voxels inside the particles are assigned with the constant ∆n, and the rest of the background voxels are assigned with zero. The lateral sampling size ∆x, ∆y are both 0.1725µm, and the axial sampling size ∆z is λ/16 = 0.0297µm, where λ = λ 0 /n 0 and λ 0 = 0.632µm. The resulting object size in our forward model is 1024 × 1024 × 16840 voxels. Example computed holograms at four particle densities are shown in Fig. 2(a). As expected, characteristic fringe patterns from individual particles are still visible at the lowest density case. As the density increases, the holograms gradually become partially developed speckle patterns.
First, we analyze the BPM's accuracy by comparing the histograms of the computed hologram with the captured hologram at each particle density. As shown in Fig. 2(b), the histograms match well at all four densities.
Next, we assess the BPM's accuracy in the spatial frequency domain. We calculate the normalized spectra of the computed and captured holograms. As shown in Fig. 2(c), the frequency components of the computed holograms closely match with the experimental measurements within the NA of the imaging system.
Finally, we compare the accuracy of the BPM and FBM based on the hologram contrast C = σ µ , where σ and µ are the standard deviation and mean of the hologram, respectively. Briefly, the FBM is a linear model that assumes a weakly scattering object, and the resulting scattered field is linearly related to the scattering potential V(x, y, z)
. Given the plane-wave incident field U in (x, y, z) = e in 0 k 0 z , the total field U FBM can be approximated as where the Green's function G(r) = exp(ik 0 n 0 r)/r, n(x, y, z) is the refractive index at location (x, y, z) and r = √︁
To perform the computation, the 3D volume is first discretized with voxels in size ∆x × ∆y × ∆z. In our simulation, ∆x = ∆y = 0.1725µm and ∆z = λ/16. The scattering potential is calculated as 1 4π k 2 0 ∆x∆y∆z(n 2n 2 0 ) for voxels within the particle, and is zero for the rest of the voxels. To compute the scattered field, we treat the discretized volume as a series of 2D slabs. The scattered field from a given slab can be efficiently calculated by first multiplying the incident field at the given depth with the scattering potential of the slab and then convolve with the Green's function implemented with the fast Fourier transform (FFT)-based algorithm. The total field U FBM is obtained by summing the incident field at the hologram plane and the total scattered field from all the slabs.
As shown in Fig. 2(d), the contrast from the BPM-computed holograms agree well with the experiments. The BPM slightly under-estimates the contrast. A possible reason is that the BPM approximates the multiple forward scattering by computing the complex field slice-by-slice but ignores backscattering. As a result, the BPM may slightly under-estimates the high frequency information, which reduces the contrast in the calculated holograms. As a comparison, the contrast from the FBM-computed holograms are consistently higher than the experimental data. The contrast discrepancy increases as the particle density. At the highest density, the discrepancy for the BPM is 0.0048, whereas the FBM is 0.0422, representing a 9× improvement by the BPM.
Overall, these studies show that the BPM can accurately model multiple scattering in a dense 3D particle field and significantly outperforms the single-scattering model.
The BPM accuracy is primarily influenced by the axial sampling distance ∆z and the scattering strength of the 3D object. In the following, we quantitatively evaluate the model accuracy under different axial sampling and scattering conditions (by changing the particle density ρ and refractive index contrast ∆n), while fixing the lateral sampling ∆x = ∆y = 0.1725µm.
In Fig. 3(a), the accuracy of BPM is plotted for different ∆z under the same scattering condition. We compare the holograms computed with different ∆z with the reference Î0 using ∆z = λ/16. The difference is quantified by the signal-to-noise ratio (SNR), SNR = 10 log 10
, where Î is the computed hologram with axial down-sampling. The accuracy of the BPM drops rapidly when ∆z is smaller than the particle diameter (1µm) and reduces slowly as ∆z further increases. This indicates that to accurately compute the inner-particle multiple scattering using the BPM, dense axial sampling is generally needed. Nevertheless, even under coarse axial sampling up to ∆z = 16λ, the BPM is still more accurate than the FBM computed with the reference dense axial sampling ∆z = λ/16.
Next, we study the BPM's accuracy for different ∆z for different particle densities ρ and refractive index contrasts ∆n (Fig. 3(b) and (c)). Our study shows that the shape of the curve remain the same for different scattering conditions, which suggests that it is only determined by the sampling distance ∆z.
Next, we study how the scattering strength affects the BPM's accuracy for a fixed ∆z = λ/4. We compute holograms at different particle densities ρ, refractive index contrasts ∆n, and volume thicknesses Z. We then compare these holograms with the corresponding reference (∆z = λ/16). The results are summarized in Fig. 3(d)-(f). In general, the accuracy decreases as the scattering becomes stronger. The SNR is found to satisfy the following scaling law:
where we define a scaling parameter A(∆z) to describe the effects of ∆z whose values are proportional to that shown in Fig. 3(a). X, Y, Z are the lateral and axial sizes of the object volume and P is the total number of particles. Intuitively, Eq. ( 13) shows that the SNR is approximately inversely proportional to the total scattering potential V (where V =
)P, n is the refractive index of the particles, R = 0.5µm is the radius of the particles) of the particle volume.
In this section, we report 3D particle imaging results in both simulation and experiments. For visualization, the extracted centroid of each particle is enlarged to a 1µm disc. Note the definition of the z-direction is the same as that in Fig. 1, in which the left sides of (c) and (e) lie closest to the hologram plane.
Given the high accuracy of the BPM model validated in Sec. 3., we next quantify the accuracy of our reconstruction algorithm in simulation. To implement the reconstruction algorithm, we first select a suitable ∆z. Due to memory limitations, it is not feasible to perform reconstructions using the same dense ∆z = λ/16 as the forward simulation. Based on our quantitative study in Fig. 3, we select ∆z = 6.2λ for all the reconstructions, which balances the model accuracy and computational cost. The corresponding object volume contains around 177 million voxels.
In Fig. 4(a), we simulated a hologram from a particle field (ρ = 6.41 × 10 4 /µL, ∆n = 0.26, particle diameter:1µm) using the densely sampled BPM model (∆z = λ/16). The 3D reconstruction result is visualized in Fig. 4(b). For better visualization, we show the x-projection of the central 69 × 69 × 500µm 3 region in Fig. 4(c). As expected, the reconstructed particles are elongated along the z-axis due to the missing cone problem. To further evaluate the reconstruction, we extracted the centroids of the reconstructed particles and then overlay them onto the groundtruth. The z-and x-projections are shown in Fig. 4(d) and (e). As shown in the z-projection, the reconstruction errors mostly occurred in the peripheral FOV region. Further inspecting the the x-projection, the localized centroids agree well with the ground truth. To further quantify the reconstruction accuracy, we repeat the simulation at seven different refractive index contrasts and four different particle densities. We then use Jaccard index (JI), lateral root mean square error (RMSE), and axial RMSE for quantification [24]:
where TP, FP, and FN denote the number of true positive, false positive, and false negative particles, respectively, δx, δy, and δz measure the distance between the centroids of the reconstructed and the matching ground-truth particle, and B is the set of all TP particles, see more details in Supplement 1.
As shown in Fig. 5(a), our method achieves JI > 0.9 for imaging particles at densities ρ ≤ 6.41 × 10 4 /µL (1000 particles in the 176.74 × 176.74 × 500µm 3 volume) with ∆n = 0.26. Our method also achieves localization accuracy better than the diffraction limit. The lateral RMSE is less than 0.25µm and the axial RMSE is less than 3.5µm in all cases. Figure 5(b) shows representative z-projections of the centorids overlaid onto the ground-truth for the central 69 × 69 × 500µm 3 region. These results show that our method performs well for particle densities as high as 6.41 × 10 4 /µL and refractive index contrast as high as 0.32. Next, we quantify the improvement of our multiple-scattering BPM model as compared to the single-scattering FBM model. To do so, we performed reconstructions using the compressive holography algorithm that uses the FBM forward model and a total variation (TV) regularization to enforce sparsity [6,7]. In Fig. 6, we compare results across four different particle densities with ∆n = 0.26. To quantify the depth-dependent reconstruction accuracy, we calculated JI by dividing the entire depth into 17 segments and then calculating JI for each segment. The depth-wise JIs under different scattering densities are shown in Fig. 6(a). The results show that the BPM-based algorithm consistently detects the particles across all the depths for particle densities ρ ≤ 6.41 × 10 4 /µL. In contrast, the FBM-based algorithm suffers from rapid degradation as the depth increases. In particular, when the depth is around 500µm, the JI of BPM is >34 times higher than the FBM when ρ ≤ 6.41 × 10 4 /µL. Representative z-projections of the reconstructed centroids overlaid onto the ground truth across different depths are shown in Fig. 6(b) for the central 69 × 69 × 500µm 3 region. These results highlight that our multiple-scattering algorithm significantly outperforms the traditional single-scattering method, in particular for reconstructing particles at deep depths.
We demonstrate our reconstruction algorithm on experimentally captured holograms at 4 different densities. In Fig. 7, 3D visualizations of example reconstruction results are shown in depth-color coded 3D renderings and x-and z-projections. The axial elongations are visible in all cases, which are resulted from the missing spatial frequency information limited by the single-view holographic measurement. We further observe that more particles are reconstructed at the depths closer to the hologram plane. This is expected since, for particles closer to the hologram plane, a greater amount of scattered information are captured by the finite-sized image sensor. Finally, we compare our method with the FBM algorithm in Fig. 8. Similar to our observations in the simulation, the BPM algorithm can better reconstruct particles at deeper depths than the FBM algorithm, as highlighted by the x-and z-projections for different depth regions.
We have developed a new reconstruction algorithm for large-scale holographic particle 3D imaging based on a multiple-scattering beam propagation model. Our forward model demonstrates superior accuracy for multiple-scattering particle fields as compared to the traditional first Born-based single scattering model. The computational efficiency of our iterative algorithm allows reducing the computational complexity by more than 2 orders of magnitude for reconstructing volumes containing morn than 100 million voxels, as compared to the iterative Born series based method. The accuracy of our proposed algorithm is demonstrated in both simulation and experiment. In particular, we show that our method is particularly effective to improve the imaging performance for particles at deep depths. Together, these advances may open up new exciting opportunities for large-scale holograph particle 3D imaging in various applications.
Though our algorithm achieved start-of-the-art performance, it is still limited by the severe missing cone problem, which elongates the reconstructed particles. A promising future direction is to combine our multiple-scattering model and advanced deep-learning priors to further improve the imaging performance [25,26].
Funding. National Science Foundation (1813848, 1846784).