1. Introduction

Digital inline holography (DIH) is a diffractive imaging method in which volumetric information is encoded and subsequently extracted from a 2D image [1]. The ability to resolve the position of objects in 3D naturally leads to temporal tracking [2, 3] with applications in particle dynamics [4], microorganism swimming [5], insect flight [6], measurements of 3D flow fields [7], and flow-structure interaction [8], among many others. The primary advantage of DIH particle tracking velocimetry (DIH-PTV) is that it is able to capture time-resolved 3D velocities using only a single camera. DIH-PTV is substantially less expensive than methods such as tomographic PTV or traditional PIV due to the need for only one camera and compatibility with low cost laser sources. The low cost and hardware simplicity of DIH has enabled multiple in situ applications [9–12].

Despite these advantages, there are several drawbacks of DIH-PTV that have limited broad application of the method. Since the inception of DIH-PTV, poor longitudinal (i.e., out of the plane of the image) resolution has consistently been the greatest challenge [2, 13]. Known as the depth-of-focus (DOF) problem, the apparent elongation of reconstructed particles is caused by the finite sensor extent and the intensity integration effect of physical pixels [1, 2]. Additional imaging noise further reduces the ability to resolve the critical high frequency fringes. Because of this limitation, some applications of DIH-PTV (namely [13, 14]) have primarily focused on the analysis of the much more accurate in-plane velocities. When high particle concentrations are used, cross-interference and other noise sources (twin image, out-of-focus particles) result in a low signal-to-noise ratio (SNR) [15]. Malek et al. showed that the reconstruction quality depends on both the shadow density and the depth of the sample [15]. The shadow density is

Many approaches to overcoming these drawbacks focus on hardware design to improve the recorded image quality or encode more information in the recording. The most common is the use of multiple viewing angles (tomographic DIH), using the lateral accuracy and some depth information from each view to more accurately localize the particles [16–18]. This method requires only two cameras (or one by using mirrors [16]) compared to four required for conventional tomography. Other approaches specific to DIH-PTV seek to reduce the effective shadow density by illuminating only a limited volume of interest [19] or using localized particle seeding [20]. Due to their mechanical and optical complexity, these methods are non-trivial to implement and care must be taken to avoid flow disturbance. Off-axis holography is also commonly used as it does not have the twin image problem and separates the reference beam from the object (reducing contamination) [11, 21]. However, this requires precision alignment and higher laser coherence. Collectively, these methods requiring multiple optical paths, viewing angles, and calibration thereof negate the principle advantages of DIH-PTV, namely ease of use and hardware simplicity.

Many authors have focused on improving the numerical processing of DIH-PTV images. Much of this work has been focused on automatic detection of the object focal plane [9, 22]. However, these do not address the problem of low SNR and usually assume that accurate 2D segmentation is trivial which is only true when the particle concentration is very low. Iterative phase retrieval methods have been shown to solve the twin image problem and improve the reconstructed SNR [23, 24], but have not been applied for PTV. Holographic deconvolution [13, 25–27] borrows a method from optical microscopy to treat the apparent blurring of point objects in the 3D reconstruction as convolution of the true object with a blurring point spread function (PSF). However, the dependence of deconvolution on a 3D Fourier transform makes this method memory intensive and windowing may be needed to process large holograms (more than ${10}^8$ voxels). The point object assumption also limits the range of applications suitable for deconvolution.

A more recent approach to holographic reconstruction is the inverse method [28, 29]. Inverse methods, rather than reconstructing the object from the image, instead find the optimal object that would produce the observed image while satisfying some physical constraints. One of the first inverse problem formulations was proposed by Soulez et al. [28] who performed a 4D parametric optimization to find the 3D location and radius of spherical particles. Their stated linear time dependence on the number of particles makes this method unsuitable for fluid flow applications where thousands of particles must be tracked for hundreds of frames. Furthermore, the assumption of spherical particles restricts the scope of potential applications. The use of the term "compressive holography" (CH) to refer to the inverse problem was introduced by Brady et al. in 2009 who borrowed concepts from the field of compressive sensing for holographic processing [29]. They used a total variation regularized approach to produce in-focus images of two dandelion seed parachutes recorded concurrently at two different focal planes. Denis et al. [30] used a similar approach with a simpler sparsity-based regularization. Both approaches show a significant reduction in the out-of-focus noise, twin images, and other noise. Recent work has used a set of physically meaningful constraints (including sparsity, smoothness, and non-amplification) in a unified parametric and CH framework to achieve excellent reconstructions of absorption and phase of individual evaporating particles and their evaporation tails [31, 32]. However, these methods operate on a limited number of planes and cannot reconstruct the 3D shape of complex objects.

Inverse methods have only been used to reconstruct objects for which the axial location is known either a priori or from a conventional reconstruction method. There have been no applications of compressive holography for DIH-PTV of flows containing thousands of tracer particles. The primary barrier preventing such application is the high computational cost. For illustration, we consider a case with 1000 particles per image and 100 images each sized 1024 × 1024 pixels with 1024 reconstruction planes (${10}^9$ voxels). The best reported speeds for parametric methods is approximately 4 seconds per particle (4.6 days for our example) [33]. A previous GPU-accelerated compressive holography implementation can reconstruct a volume up to $1024\times 1024\times 10$ voxels in 7.6 seconds (22 hours for our example) [34]. Other methods have even longer extrapolated times including 1000 days [31] and 400 days (on modern hardware) [29]. In addition to the time required for processing, memory requirements for CH place severe restrictions on the size of hologram that can be processed. Storage of the hypothetical test case hologram would require approximately 8 GB to store in memory (complex floating-point values, 8 bytes each). Several additional variables of this size are needed for CH algorithms. However, contemporary GPU memory is limited to approximately 12 GB for consumer hardware and memory transfers from the GPU to RAM are slow. Therefore, current application of CH must either limit the volume size to take advantage of the speed increase of GPU-acceleration or rely on the much larger RAM available on most desktop computers and rely on much slower CPU processing.

In the present study, we first summarize the fundamentals of CH. We then introduce our proposed method using fused lasso regularization and a sparse storage structure to enable processing of very large images in a realistic time (55 hours for the 100 image sequence of large holograms described above). We then provide several synthetic and experimental evaluation cases to demonstrate the quality and performance of the proposed method.

2. Methodology

We formulate the 3D reconstruction of the object volume as an inverse optimization problem, following the method of Brady, Endo, and others [29, 34]. The optimization problem formulation, Eq. (2), seeks to find the object field (x) that minimizes the difference between the observed hologram (b) and the estimated hologram produced by propagating the object to the imaging plane ($\widehat{b}=Hx$). A regularization term $\lambda R\left(x\right)=g\left(x\right)$ is included in Eq. (2) to enforce physically realistic properties (principally sparsity and smoothness) and avoid trivial solutions. The size of the regularization parameter λ determines how strongly these properties are enforced. To ensure solution convergence, we use a linearized form of the forward model for hologram formation, Eq. (3), which implicitly treats any nonlinear terms (i.e., twin image and cross interference) as noise.

Here k₀ is the wavenumber, h_z is the Rayleigh-Sommerfeld diffraction kernel, and convolution is performed in the Fourier domain. The discrete Fourier domain representation of h_z is given by Eq. (4) [2] where m andn are the discrete coordinates, M and N are the total number of elements, and s is the pixel size.

The only difference between forward propagation (hologram formation) and back-propagation (reconstruction) is the sign of the depth, z. Note that Eq. (3) models the object field as a series of discrete planes in z. The true object must fall on one or more of these planes in order to be correctly modeled.

We solve the inverse problem using FISTA (Fast Iterative Shrinkage-Thresholding Algorithm) [35] as implemented for FASTA [36]. This method is selected due to its high convergence rate, relative simplicity, and similarity to prior work [29, 34]. FISTA has two steps: a shrinkage step using the proximal operator (Eq. (5)) and an accelerated update using the previous estimate (a simple element-wise operation not included here for brevity).

In the context of holography, the gradient term $\nabla f\left(x_{k-1}\right)$ is the reconstructed volume from the residual hologram ($\widehat{b}-b$). The scalar factor L is determined automatically using backtracking [35]. The proximal operator (Eq. (6)) can be interpreted as a gradient descent step with step size L [37]. The reader is referred to the references for further details on FISTA.

The form of the regularization function $R\left(x\right)$ determines which properties of the solution will be enforced. The ${\mathcal{l}}^1$ norm (Eq. (7)) enforces a sparse solution (i.e., one with few non-zero elements) which for DIH-PTV means that the object volume fraction is low. This sparsity-based regularization has been demonstrated for holography by Denis et al. [30] and Endo et al. [34].

The Total Variation (TV) norm (Eq. (8)) is the sum of the first-order gradients over the image (size $N_x\times N_y$). These gradients are computed using a circular boundary condition (i.e., $x_{0,0}=x_{N_x,N_y}$) which is valid when the image is zero-padded to ensure no objects are near the borders. TV regularization enforces a smooth solution (small gradients) and is naturally extensible to higher dimensions.

The TV approach has been used by Brady et al. [29] and Endo et al. [34] who demonstrate that it is superior to the ${\mathcal{l}}^1$ regularization for sufficiently large objects. The reason for this benefit is that the smooth solution enforced by the TV regularization results in reconstructed objects that are not over-segmented. However, we will see that TV regularization is substantially more computationally demanding that the ${\mathcal{l}}^1$ method.

We propose using the Fused Lasso (FL) regularization method (Eq. (9)) which is a combination of the TV and ${\mathcal{l}}^1$ norms ("fusion" and "lasso" being alternative terms for TV and ${\mathcal{l}}^1$ respectively) [38]. Solutions to the FL problem are both smooth and sparse while having some characteristics that make it less computationally demanding than TV.

The shrinkage step of FISTA (Eq. (5)) is by far the most computationally complex and depends on the choice of $R\left(x\right)$. As such, the computational cost of FISTA is closely linked to that of evaluating the proximal operator. For the ${\mathcal{l}}^1$ regularizer, the proximal operator has a simple closed-form solution as soft thresholding:

However, the proximal for the TV function does not have a closed-form solution and requires an iterative solution, for which we use the gradient projection method of Beck & Teboulle [39]. This method requires storage of each directional derivative for the duration of the iterations which may require a substantial amount of memory. The FL regularization function has the useful property that it is separable and can be computed by soft thresholding the solution to the TV problem (i.e., with ${\lambda }_{{\mathcal{l}}^1}=0$

). Because the non-sparse TV solution must first be computed before soft thresholding to produce the sparse FL solution, high memory requirements of the TV proximal still apply within each FISTA step even though the result is sparse. We limit our TV regularization to 2D planes which can be computed independently, reducing the memory requirements to those of a single plane. It is worth noting that prior compressive holography methods using TV regularization have reported only the 2D variant.

Because FISTA is an iterative solution method, the computational time required to process a single image may be relatively high. PTV requires processing thousands of large, well-resolved volumes. Therefore, it is crucial to reduce the processing time to amanageable level to enable application to real flow studies. We utilize a CUDA/C++ GPU implementation of our algorithm to accelerate the processing. Key components such as the fast Fourier transform (FFT) already have efficient GPU library implementations while the reconstruction kernel and proximal operator largely use highly parallelizable pixel-wise operations. One of the principle challenges in implementing this approach is memory usage. A reasonably sized reconstruction volume ($1024\times 1024s1024$ voxels) would require approximately 8 GB to store in memory (complex floating-point values, 8 bytes each). Several additional variables of this size are needed for FISTA. However, contemporary GPU memory is limited to approximately 12 GB which would place limits on the type of holograms which could be processed. In order to circumvent this challenge, we exploit the sparsity of the ${\mathcal{l}}^1$ and FL regularized solutions to dramatically reduce the memory requirements. We use the coordinate (commonly, COO) sparse matrix format to store all volume data during iterations. The COO format stores the indices (row and column) and value for each non-zero element in a plane. Because data is accessed per plane for both the forward and reverse propagation as well as the 2D TV proximal operator, each plane is independently indexed. The total storage for each non-zero element is thus 24 bytes (8 bytes per index, 8 bytes for complex floating-point value) compared to 8 bytes per element for a non-sparse structure. Thus, memory usage should be reduced as long as the data sparsity (number of zero elements divided by total) exceeds 67%. Experience suggests that most PTV holograms have sparsity exceeding 99% [15, 27].

The primary advantage of the compressive holography approach is that it produces very high SNR reconstructions that are more easily segmented for particle localization. In one sense, the sparsity regularization inherently separates objects (non-zero voxels) from the background (zero voxels), negating the need for complex volume normalization and SNR enhancement such as that used by Toloui et al. [13]. While these directly thresholded results are reasonable, we have found that two additional filters greatly reduce the instances of over-segmentation. The first is a very low intensity threshold on the order of $1/{256}^{th}$ the maximum intensity of the image. This value is selected as any values below it are indistinguishable from zero when using a min-max scaling and 8 bit discretization for visualization. The second filter is a minimum object volume. This must be adjusted slightly depending on the size of the particles, noise level, and apparent elongation length. At this time, it is not directly linked to the true particle volume. Usually, objects of 5 voxels or fewer are treated as noise. While crucial for counting the number of particles in a single hologram, these parameters have minimal effect when applied to a sequence of images for which the particles are tracked because over-segmentation noise rarely persists for multiple frames.

Because inverse approaches depend on matching the recorded image using a model, any image features that are not modeled (e.g., background intensity variations, imperfections in the optical path) will reduce the quality of the fit. To rectify this, we preprocess the recorded images by removing the constant background using the method of [4] (subtraction followed by division by the square root). The background is computed as the mean of a time sequence of holograms in which the objects of interest have significant motion.

While compressive holography and inverse methods have existed for over a decade, this is the first application to 3D PTV. Previous uses of a parametric inverse method for particle tracking have tracked fewer than 10 particles concurrently [28, 33, 40, 41]. Furthermore, CH is usually used with a small number ($\sim 10$) of reconstruction planes with a large spacing ($\sim 1$ mm) between planes. Here we demonstrate the ability to reconstruct volumes with over 1000 planes with both cubic and elongated reconstructed voxels. The largest volume reconstructed by Endo et al. [34] contained ${10}^7$ voxels while our sparse representation enables reconstruction of volumes containing more than ${10}^9$ voxels on a desktop computer. The use of the fused lasso regularization to enforce both smoothness and sparsity has not been previously demonstrated for compressive holography. To emphasize these distinctions, we refer to our method as a Regularized Inverse Holographic Volume Reconstruction (RIHVR, pronounced "river"). RIHVR dramatically increases the SNR of the reconstructed volume. This enables processing of high noise and high particle concentration holograms (both traits are common in DIH-PTV applications) that could not be reliably processed using existing methods. Because RIHVR does not assume a size or shape of the object, it can be used when the imaged particles are polydisperse or non-spherical. We next present several practical examples to demonstrate these capabilities.

3. Demonstration cases

To demonstrate that the proposed method is applicable to a variety of DIH-PTV cases, we present the results for processing four classes of holograms: an isolated nanowire, simulated tracer particles in isotropic turbulence, swimming microorganisms, and an experimental T-junction flow seeded with microfibers. The first case, the isolated nanowire, demonstrates improved 3D reconstruction of a continuous object with a significant 3D shape. Simulated holograms then provide a realistic flow case for which ground truth exists for the particle locations. The RIHVR method is evaluated against deconvolution (with inverse iterative particle extraction where applicable) which has been previously validated against conventional PIV and shown to provide substantial improvement over other DIH-PTV approaches [13]. A simple reconstruction method following the approach of Pan & Meng [42] (global thresholding followed by peak intensity depth localization) is also shown for comparison. Experimental holograms of swimming microorganisms and microfibers in a T-junction flow represent real measurement domains for which some flow behaviors are known from prior studies. These later cases demonstrate that RIHVR can be applied to broad measurement domains where other DIH-PTV methods fail.

 

Fig. 1 (a) Recorded holographic image of a 90 nm Ag nanowire. (b) Hologram after image enhancement. (c) Volumetric reconstruction of the sample using the deconvolution method. (d) Reconstruction using RIHVR with sparsity (${\mathcal{l}}^1$) regularization. (e)Reconstruction using RIHVR with fused lasso (FL) regularization. For visualization, (c) uses the intensity as the transparency alpha value while (d) and (e) show all non-zero values at equal intensity.

Download Full Size | PDF

3.1. Isolated nanowire

A qualitative evaluation of the proposed inverse reconstruction method uses a silver nanowire in suspension. This is an example of a continuous object with significant extent in all three reconstruction dimensions. The length of the wire is not known a priori. As such, a parametric inverse model such as the one used by Soulez et al. [43] is unsuitable. The sample is a suspension of 90 nm diameter Ag nanowires in isopropyl alcohol. The illumination source is a 450 nm fiber-coupled laser diode (QPhotonics QFLD-450-10S), collimated using a Nikon CFI Plan Fluor 10X objective lens. A Nikon CFI Apo TIRF 100X oil immersion microscopic objective and video camera (Andor Zyla 5.5 sCMOS) are used to image the sample. The recorded pixel size is 70 nm. The recorded image (2560 × 2160 pixels) is cropped to a 1024 × 1024 pixel region around a selected nanowire to ensure that only a single object is in the image and to reduce unnecessary computational cost. Reconstruction is performed at 70 nm intervals (equal to the lateral pixel pitch) for a depth of 42 μm (600 planes). Measurement of similar samples using DIH has been undertaken by Dixon et al. [25] who measured the diffusion of nanowires and Kempkes et al. [44] who demonstrated a 2° accuracy for the orientation of microfibers. Unlike the prior methods, our approach does not assume a linear fiber and is suitable for measuring non-rigid wires.

The raw and enhanced holograms are shown in Figs. 1(a) and 1(b) respectively alongside renderings of the reconstructed volumes produced using deconvolution, RIHVR with sparsity regularization (${\lambda }_{{\mathcal{l}}^1}=0.27$), and RIHVR with fused lasso regularization (${\lambda }_{{\mathcal{l}}^1}=0.15$, ${\lambda }_{TV}=0.12$). Figures 1(c)–1(e) show that both regularization methods substantially reduce the DOF of the reconstruction. Measured as the width at half the measured intensity averaged along the wire length, the DOF decreases from 1.97 μm using deconvolution to 0.63 μm and 0.89 μm using the sparsity and fused lasso regularization methods respectively. Similarly, a 99% decrease in the segmented volume and 90% decrease in the segmented cross-sectional area is observed between deconvolution and RIHVR (with similar reduction for both RIHVR regularizers).

When comparing the results generated using the sparsity and fused lasso regularization methods in Figs. 1(d) and 1(e), the smoothing effect of the fused lasso is apparent. The fused lasso regularized results show fewer gaps in the wire profile and an overall more contiguous object. However, this comes at the cost of some expansion of the object and a slightly larger DOF. Interpolated cross-sections normal to the wire axis (insets in Fig. 1) illustrate that both RIHVR approaches approximate the true circular shape of the wire. Conversely, Fig. 1(c) shows that deconvolution produces an X-shaped cross-section characteristic of simple holographic reconstructions. RIHVR also demonstrates robustness to image noise. The raw image in Fig. 1(a) has a substantial amount of background noise and even enhancement via background removal does not produce a noise free image. Additional fringe patterns – caused by vibrations, fluctuations in illumination intensity, and out-of-view objects – are visible in Fig. 1(b) (the enhanced image) but do not result in artifacts in the reconstructed volume when using RIHVR.

3.2. Synthetic turbulent flow

Turbulent flows represent the most challenging case for 3D flow measurements as they are highly three-dimensional and involve velocity fluctuations across a broad dynamic range of scales. Here we assess the accuracy and limitations of our method using simulated holograms of a homogeneous isotropic turbulent flow. The simulated tracer particle trajectories are determined by querying the forced isotopic turbulence data from the Johns Hopkins Turbulence Database with Lagrangian particle tracking [45–47]. The simulation domain is scaled to $5\times 5\times 5$ mm${ }^3$ and sampled with a nondimensional time step of 0.012 (60 DNS time steps) to capture 100 instants (image frames). The Reynolds number based on the domain size is 23,000 and the Kolmogorov length and time scales (smallest scales of turbulent fluctuations) are 67 μm and 3.7 frames respectively. The RMS velocity is 6.7 μm/frame. Maintaining the Reynolds number and using a low viscosity fluid ($\nu ={10}^{-7}$ m${ }^2$/s), this corresponds to a frame rate of 75 kHz which is achievable with modern cameras. One thousand particles with a diameter of 7.5 μm are initially randomly spatially distributed throughout the 3D domain and their positions at subsequent time steps are determined using a Lagrangian tracking method [47]. A periodic boundary condition is applied to the particles to ensure that the number of objects in the field of view is constant (this is ignored during processing). The simulated holograms are 512 × 512 pixel images with a 10 μm pixel size and 632 nm illumination wavelength. A sample hologram is shown if Fig. 2.

 

Fig. 2 Sample simulated hologram containing 1000 particles

Download Full Size | PDF

The reconstruction plane spacing is equal to the lateral pixel spacing (10 μm). The total size of the reconstructed volume is $512\times 512\times 700$ voxels ($1.8\times {10}^8$). 100 FISTA iterations are used with regularization parameters ${\lambda }_{{\mathcal{l}}^1}=0.5$ and ${\lambda }_{TV}=0.2$. For segmentation, the minimum intensity threshold is set to $2/{256}^{th}$ the maximum and the minimum volume is 5 voxels. The particle positions are estimated using the weighted centroid of each connected component and the particles are tracked using the method of Crocker & Grier [48] with a maximum per frame displacement of 70 μm and a minimum tracked duration of 10 frames. Two alternative hologram reconstruction methods are presented for comparison. The first is a simple reconstruction method following the approach of Pan & Meng [42]. The second is the deconvolution method of Toloui & Hong [27] with two passes of the inverse iterative particle extraction step. The particle trajectories for all methods are smoothed with a total variation filter. The tracked results using RIHVR are shown in Fig. 3(a).

 

Fig. 3 (a) Smoothed 3D particle trajectories extracted from a synthetic hologram using RIHVR. (b) Localization error of tracked particles relative to their true locations. (c) Error in the RMS velocity components of the three test methods compared to ground truth.

Download Full Size | PDF

To evaluate the localization error, extracted particles are matched to their true location using a nearest neighbor method [48]. The resulting error distributions in each dimension are summarized in Fig. 3(b). For all three methods, the error in x and y is very small (smaller than the pixel size). However, the error in z is substantially greater, demonstrating the DOF problem. Comparing the three methods, the ${75}^{th}$ percentile decreases from 11.5 voxels using reconstruction to 6 voxels using deconvolution and 3.5 voxels with RIHVR. The same trends are seen at the other percentiles as well. Thus, RIHVR produces a 40% improvement in longitudinal localization over the prior best method and a 70% improvement over simple reconstruction. These trends also hold when evaluating the mean error of the unfiltered instantaneous velocity measurements. The error in the axial velocity, w, decreases from 13.4 voxels/frame with reconstruction to 3.8 voxels/frame with deconvolution and 2.4 voxels/frame with RIHVR.

 

Fig. 4 (a) Extraction rate for each method for increasing particle concentrations. (b) The number of particles which can be accurately extracted is higher for RIHVR than the other methods. Solid line is 100% EP, dashed line is 50%.

Download Full Size | PDF

For turbulence measurements, it is common to measure Reynolds stresses, which are velocity fluctuation statistics [7]. Here we present measurements of the root-mean-square (RMS) velocity in Fig. 3(c). This is comparable to Reynolds stress when the mean is zero (as it is for this flow) while maintaining intuitive meaning for applications other than flow measurement. For this flow in the period during which particles are simulated, the true RMS velocities averaged over the whole volume are (6.3, 7.6, 6.2) voxels/frame in the u, v, and w directions respectively. The trajectory smoothing produces a 3% error in the $u_{RMS}$ and $v_{RMS}$ measurements but significantly reduces spurious fluctuations in z. Using reconstruction, the measured $w_{RMS}$ differs from the true value by over 800%. This is reduced to 30% using deconvolution. However, this is still unacceptably high for real measurements. The error using RIHVR is only 7% which is substantially better and only 2× greater than the error in the u and v measurements.

As the velocity vector spacing in PTV is directly related to the particle spacing, the maximum particle concentration is a critical concern for many PTV measurements. The quality of recorded holograms depends on several factors including the particle concentration, size, volume depth, wavelength, magnification, and numerical aperture. Here we focus on only the particle concentration to compare the tested methods. Note that the results may vary with changes to the other parameters but the relative comparisons between the methods should remain consistent. In general, the extraction rate (Ep, number of correctly extracted particles divided by true particle count) decreases with increasing concentration. Here we use the number of particles per pixel to scale the concentration because it allows for the most direct comparison with the literature. It has previously been shown that the commonly used shadow density (Eq. (1)) does not completely explain the extraction rate in all situations [15, 27]. For comparison, Toloui et al. [13] performed measurements with a concentration of 0.0035 particles/pixel while other sources used significantly lower concentrations. Figure 4(a) shows that concentrations of up to 0.035 particles/pixel can be processed using RIHVR while maintaining $Ep>60\%$. The increased number of extractable particles enabled by RIHVR, as shown in Fig. 4(b), enables increased resolution in velocimetry applications and higher particle concentrations in other applications including studies of biological flows and fluid-particle interaction where high concentration may be crucial for the sample being studied. An example of such acase is given in section 3.3.

3.3. Swimming algae

One practical application of DIH-PTV is the study of microorganism swimming behaviors. Previous studies have used small sample volumes ($\sim 0.05\mu$L) in order to measure the large cell concentration present in cultures ($\sim {10}^6$ cells/mL) [5, 21]. Here we demonstrate that RIHVR is superior to prior DIH algorithms for these experiments. We also demonstrate the ability to record and process much larger sample volumes ($\sim 1\mu$L) which could enable new scientific studies.

The alga Dunaliella primolecta is a unicellular species which can be used for biofuel production [49]. Cells have a length of 10 μm and swim using two flagella [50]. In this study, D. primolecta is grown at 37°C in a growth medium. Manual concentration measurements using a microscope indicate that the sample has a concentration of $1.8\times {10}^6$ cells/mL. The sample container is a $10\times 30\times 1$ mm${ }^3$ glass cuvette. Holograms are recorded at 100 Hz using a 2048 × 1088 pixel sensor (Flare 2M360-CL). The sensor pixel size is 5 μm and the exposure time is 50 μs. A 5x microscopic objective is used, resulting in a recorded sample volume of $2.05\times 1.09\times 1$ mm³). For simplicity and speed, the recorded image is cropped to a size of 1024×1024 pixels (1×1 mm²). The light source is a 532 nm diode laser (Thorlabs CPS532) which is expanded and filtered with a spatial filter, as shown in Fig. 5(a). While the number of particles per pixel is relatively low for this sample (0.002), the particles are large enough that the shadow density (Eq. (1)) becomes significant, $s_d=18\%$. The maximum shadow density used by Toloui et al. [13] was 10.5% using deconvolution, while Malek et al. [15] achieved an extraction rate of only 20% for $s_d=10\%$. Reducing the measurement depth can enable holograms to be processed using conventional methods [50], but risks introducing wall effects that influence the behavior. Similarly, we have found that dilution of the sample changes the cell swimming behavior. Therefore, high concentration holograms – which can only be processed using RIHVR – are important to these microbiological studies. Studies of microorganism behavior using non-holographic methods are challenging because their 3D motion leads to low residence time in a microscopic depth of field and size constraints make multi-camera imaging difficult.

 

Fig. 5 (a) DIH imaging system, recorded hologram, XY projection of processed reconstruction. Scale bars are 100 μm. (b) Subset (25%) of tracked particles. (c) PDF of velocity fluctuations in each direction. The w distribution has longer tails and a sharper peak but is not substantially wider than the other two components.

Download Full Size | PDF

The holographic volume is reconstructed with 270 planes, separated by 3 μm, with the volume confirmed to include both walls of the cuvette. The regularization parameters are ${\lambda }_{{\mathcal{l}}^1}=0.1$ and ${\lambda }_{TV}=0.1$ with 100 FISTA iterations and 5 gradient projection steps to evaluate the TV proximal operator. A sequence of 2000 frames is analyzed. A sliding window of 151 images (1.5 sec) is used to compute the mean background in order to reduce the effect of cells starting or stopping their motion. Because RIHVR uses a complex object model, it is equally suited to measuring amplitude objects (e.g., nanowire and synthetic objects) as well as phase objects such as the D. primolecta cells. More complex regularization methods (with associated higher computation cost) could be used to further enhance phase objects [31].

RIHVR detects and tracks an average of 294 objects per frame. This is dramatically lower than the expected count of 2000 cells/frame from the concentration measurement. However, a substantial number of particles are seen to remain stationary on the two walls. These are treated as background noise and are removed during the image enhancement. A selection of 3D tracks is shown in Fig. 5(b). For clarity, only a subset of 25% of the tracked data is shown in Fig. 5 while the full density is shown in Fig. 6. The cell trajectories have been smoothed using a Savitzky-Golay filter of 20 frames. The frame rate is sufficiently high that this filter does not suppress any real motions. Under the assumption that cell swimming motions are isotropic, the probability distribution functions (PDF) of velocity fluctuations (normlized by the RMS velocity) are expected to coincide for each component. Figure 5(c) shows that while the u and v velocities are in good agreement, the same is not true for w, even after smoothing. This indicates that the DOF problem is not entirely eliminated for this extremely noisy case. However, gross motions in the longitudinal direction are visible and fine scale complex behaviors such has helical swimming can be seen from an enlarged view of the sample in Fig. 6.

 

Fig. 6 2D view of the reconstructed cell trajectories showing complex behaviors Also illustrates the true cell concentration of processed volume. Scale bars are 100 μm.

Download Full Size | PDF

3.4. Rotating rods in flow

In addition to improvements to the measurement accuracy and seeding density limits, RIHVR enables the measurement of complex shapes as previously illustrated in Fig. 1. Here we present a flow case where the seeding particles are rods rather than the usual spherical tracers. Using RIHVR, we are able to extract both the location and orientation of each rod and track their evolution in the flow. This type of multimodal measurement using a single camera has not been previously reported.

To demonstrate this measurement, we use the T-junction flow of the type studied by [51] which occurs frequently in industrial and biological flows. The rotation and alignment of fibers in flow have been extensively studied for target applications including paper manufacturing and microorganism alignment (see for example [52–54]). The fibers used for the present study are marketed as an additive to strengthen composite materials, where the alignment of the fibers may have an impact on the material properties. Prior experimental work has either been restricted to 2D measurements [52, 53] or multi-camera 3D measurements of individual fibers [54]. Holography is a valuable alternative when the motion is three-dimensional, seeding density is high, or optical access is restricted.

The experimental channel has a square cross-section with a side length of 1 mm. The junction is at a right angle and all three branches (inlet and two outlets) have the same geometry. The inlet flow rate is 1000 mL/hr which corresponds to a Reynolds number $Re=290$. The seeding particles are 7 μm diameter SiC (ρ = 3 g/cm${ }^3$) microfibers (Haydale Technologies) with an aspect ratio of 10. The response time of the particle, computed using the equivalent diameter, is ${\tau }_p=14\mu$s. The characteristic time of the flow is ${\tau }_f=1.9$ ms. The resultant Stokes number is $S{t}_k=0.007$ indicating that the particles will trace the flow. A high speed video camera (NAC Memrecam HX-5) is used to record the holograms at 6000 Hz with an exposure time of 50 μs. A microscopic objective (Edmund Optics, 10x, NA=0.45) is used to image the sample, resulting in a 1024 × 1024 image with a pixel size of 0.91 μm. The light source is a spatially filtered HeNe laser (λ = 632 nm). The FL regularization method is used with ${\lambda }_{{\mathcal{l}}^1}=0.1$ and ${\lambda }_{TV}=0.12$. 110 reconstruction planes are used with a spacing of 9.1 μm ($10\times$ the pixel pitch). An intensity-weighted principle component analysis is used to determine the orientation of the fibers (similar to the method of [44]).

For validation of the flow field, the flow (absent any particles) is simulated using ANSYS Fluent (ANSYS, Inc.), with the results found to be in agreement with the simulations of [51] and the experimental particle pathlines. The fiber rotation rate is modeled using the Jeffery equation in the limit where the particle aspect ratio is $\gg$1 [55]:

Here p is a unit vector aligned with the particle axis, Ω is the rotation tensor, and S is the strain rate tensor. Because the particle rotation rate is coupled to the orientation, the rotation rates for the simulation in Fig. 7(c) assume that the particles are initially aligned with the inlet flow direction.

 

Fig. 7 (a) Visualization of the measured rod trajectories. Solid lines show the measured 3D orientation of the rods. Colors indicate individual particles. Vorticity isosurfaces (ω = 3000s⁻¹) are from the CFD simulation. (b) View of the experimental fiber orientations in the yz plane. (c) Contour map of the particle rotation rate (s⁻¹) expected from the simulation. (d) Measured 3D particle rotation rate.

Download Full Size | PDF

The experimental fiber orientations are shown in Fig. 7(a) along with vorticity isosurfaces from the simulation which illustrate that the principle flow features predicted by the simulation (two vortices aligned with x) are present in the experiment. A 2D projection of the fibers is shown in Fig. 7(b). The optical reconstruction direction corresponds to the crossflow (z) axis in Figs. 7(a)–7(d). The clear appearance of the two counter-rotating vortices demonstrates that RIHVR sufficiently reduces the DOF to enable the recovery of this 3D flow feature. Additionally, changes in the orientation of the fibers can be seen. The measured rotation rate ($\left|\dot{p}\right|$) is higher near the centers of the vortices in Fig. 7(d) which matches the expected behavior from the simulation shown in Fig. 7(c). The peak rotations rates are accurate to 30% while the location of the peaks are accurate within 0.1 mm. Some discrepancies between the measured and predicted rotation rates can be attributed to misalignment between the two domains. Because the rotation rate is dependent on the velocity gradients (which have substantial variation) and on the particle location and orientation (which have some measurement uncertainty), even small misalignment of the two domains would cause deviations between the two results. Non-ideal flow conditions, including unsteadiness and geometrical imperfections, could cause differences in the locations of the vortices and the peak rotation magnitude. The accuracy of the machining process used to make the channel is $\pm 0.05$ mm which is comparable to the peak location error. Since the flow rate is constant, uncertainty in the channel geometry also produces uncertainty in the inlet flow velocity. Finally, the Jeffery equation (Eq. (11)) used to predict the particle rotation rate assumes non-inertial ellipsoidal particles in Stokes flow. The true particle motion is expected to deviate slightly from this idealization. Given these uncertainties, we conclude that the agreement between the simulation and experimental results is adequate for demonstrating that RIHVR enables direct 3D particle rotation rate measurements.

4. Discussion and conclusions

The processing speed of RIHVR is suitable for many applications. Further adjustments by the user can enable even faster processing for specific cases. If axial uncertainty is less important, the number of reconstruction planes can be reduced, provided the planes are still sufficient to model the object field. The use of pure sparsity (${\mathcal{l}}^1$) regularization further improves speed at the cost of poorer SNR. Finally, the number of solution iterations can be reduced at the cost of a less converged solution.

While RIHVR is broadly applicable, there remain some limitations. The enhancement approach using a time-averaged background is only suitable for moving objects. The solution to the minimization problem (Eq. (2)) depends on the choice of the regularization parameters ${\lambda }_{{\mathcal{l}}^1}$ and ${\lambda }_{TV}$ which are user-specified. For most experimental cases, parameter values near (${\lambda }_{{\mathcal{l}}^1}=0.1$ and ${\lambda }_{TV}=0.1$ are appropriate. These can each be increased (or decreased) to reduce noise (${\lambda }_{{\mathcal{l}}^1}$) or increase the smoothness (${\lambda }_{TV}$) of the result. Automatic selection of these parameters is an area of ongoing research. The assumption of a sparse and smooth object field may not be true in all cases, for example if the medium has fluctuations in the index of refraction. Finally, FISTA assumes that the regularized minimization problem is convex. This assumption breaks down for large objects when a large number of reconstruction planes are used. These restrictions generally do not apply to PTV applications.

We have demonstrated the application of CH techniques to volumetric reconstruction, using the presented RIHVR approach for the reconstruction and tracking of 3D particle fields. The reconstructed volumes are both sparse and smooth, assumptions that apply equally for most particle tracking applications. The use of GPU processing and sparse storage enable the reconstruction of volumes containing over ${10}^9$ voxels which is orders of magnitude larger than previously reported for any CH method. RIHVR provides a substantial improvement in the longitudinal position and displacement measurement accuracy in addition to an increase in the particle concentration limit. These improved capabilities have allowed the extension of DIH-PTV to the tracking of a dense culture of microorganisms and measuring the orientation of microfibers in 3D flow. RIHVR is a broadly applicable approach capable of enabling low-cost 3D measurements for wide-ranging applications.