Regularized Inverse Holographic Volume Reconstruction for 3D Particle Tracking

TL;DR

This paper introduces a regularized inverse holographic reconstruction method that enhances 3D particle tracking by improving resolution, handling higher concentrations, and reducing computational costs through GPU acceleration.

Contribution

The paper presents a novel inverse problem approach with fused lasso regularization for full volumetric particle reconstruction, enabling better accuracy and efficiency in DIH.

Findings

Accurate 3D particle reconstruction demonstrated on synthetic and experimental data.

Significant reduction in computational cost using GPU processing.

Extended DIH capabilities to high concentration microorganism and microfiber applications.

Abstract

The key limitations of digital inline holography (DIH) for particle tracking applications are poor longitudinal resolution, particle concentration limits, and case-specific processing. We utilize an inverse problem method with fused lasso regularization to perform full volumetric reconstructions of particle fields. By exploiting data sparsity in the solution and utilizing GPU processing, we dramatically reduce the computational cost usually associated with inverse reconstruction approaches. We demonstrate the accuracy of the proposed method using synthetic and experimental holograms. Finally, we present two practical applications (high concentration microorganism swimming and microfiber rotation) to extend the capabilities of DIH beyond what was possible using prior methods.

Figures

Equations

$$s_d=n_{s}L{d}^2$$

$$\hat{x}=\underset{x}{\mathrm{argmin}}\left\{{\Vert Hx-b\Vert }_2^2+\lambda R(x)\equiv f(x)+g(x)\right\}$$

$$R(x)={\Vert x\Vert }_1=\sum _{i=1}^N\left|x_i\right|$$

$$R(x)={\Vert x\Vert }_{TV}=\sum _{i=1}^{N_x}\sum _{j=1}^{N_y}\sqrt{(x_{i,j}-x_{i-1,j}{)}^2+(x_{i,j}-x_{i,j-1}{)}^2}$$

$$g(x)=\lambda R(x)={\lambda }_{{\ell }^1}{\Vert x\Vert }_1+{\lambda }_{TV}{\Vert x\Vert }_{TV}$$

$$pro{x}_L(v)=\underset{x}{\mathrm{argmin}}\left(g(x)+\frac{1}{2L}{\Vert x-v\Vert }_2^2\right)\approx v-L\nabla g(x)$$

$$x_k=pro{x}_{\lambda L}\left(x_{k-1}-L\nabla f(x_{k-1})\right)$$

$$pro{x}_{\lambda L}(v)={\left(1-\frac{\lambda L}{\left|v\right|}\right)}_{+}sign(v)$$

$${\dot{p}}_i={\mathbf{\Omega }}_{ij}{p}_j+({\mathbf{S}}_{ij}{p}_j-p_i{p}_j{\mathbf{S}}_{jk}{p}_k)$$

Full text

††journal: oe

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], measurements of 3D flow fields [6], and flow-structure interaction [7, 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, 10, 11].

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 [12, 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 [2, 1]. Additional imaging noise further reduces the ability to resolve the critical high frequency fringes. Because of this limitation, some applications of DIH-PTV (namely [14, 13]) 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, 16]. 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

[TABLE]

where $n_{s}$ is particle number concentration, $L$ is sample depth, and $d$ is the particle diameter. Methods for improving DIH-PTV (see next paragraph) often require increasing the optical complexity, extensive process tuning by an expert, and expensive computations.

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 [17, 18, 19, 20]. This method requires only two cameras (or one by using mirrors [18]) 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 [21, 22] or using localized particle seeding [23]. 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) [24, 25]. 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 [26, 27]. 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 [28, 29], but have not been applied for PTV. Holographic deconvolution [30, 31, 32, 33, 13] 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 [34, 35]. 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. [34] 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 [35]. 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. [36] 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. More recently, a set of physically meaningful constraints (including sparsity, smoothness, and non-amplification) have been used to achieve excellent reconstructions of absorption and phase of individual evaporating particles and their evaporation tails [37].

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\times 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) [38]. 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) [39]. Other methods have even longer extrapolated times including 1000 days [37] and 400 days (on modern hardware) [35]. 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 [35, 39]. The optimization problem formulation (2) seeks to find the object field ($x$) that minimizes the difference between the observed hologram ($b$) and the estimated hologram produced from the object ($\hat{b}=Hx$). To ensure that the solution converges, we use a linearized form of the forward model for hologram formation which implicitly treats any nonlinear terms (i.e. twin image and cross interference) as noise.

[TABLE]

To avoid trivial solutions, a constraint must be applied to ensure a physically realistic solution. This constraint is implemented as the additional regularization term in (2), $\lambda R(x)=g(x)$. The form of this regularization function determines which properties of the solution will be enforced. The $\ell^{1}$ norm (3) enforces a sparse solution (i.e. one with few non-zero elements). This sparsity-based regularization has been demonstrated for holography by Denis et al. [36] and Endo et al. [39] who showcase its utility when the fraction of the sample volume occupied by objects is very small.

[TABLE]

The Total Variation (TV) norm (4) is the sum of the 1st order gradients over the image (size $N_{x}\times N_{y}$). It is naturally extensible to higher dimensions. TV regularization enforces a smooth solution (small gradients).

[TABLE]

The TV approach has been used by Brady et al. [35] and Endo et al. [39] who demonstrate that it is superior to the $\ell^{1}$ regularization for sufficiently large objects. However, we will see that TV regularization is substantially more computationally demanding that the $\ell^{1}$ method. We propose using the Fused Lasso (FL) regularization method (5) which is a combination of the TV and $\ell^{1}$ norms ("fusion" and "lasso" being alternative terms for TV and $\ell^{1}$ respectively) [40].

[TABLE]

Solutions to the FL problem are both smooth and sparse while having some characteristics that make it less computationally demanding than TV.

We solve the inverse problem using FISTA (Fast Iterative Shrinkage-Thresholding Algorithm) [41] as implemented for FASTA [42]. This method is selected due to its high convergence rate, relative simplicity, and similarity to the approaches of Brady et al. [35] and Endo et al. [39]. FISTA is a proximal gradient method which makes use of the proximal operator [43] which can be interpreted as a gradient step with step size $L$ (6).

[TABLE]

FISTA has two steps: a shrinkage step using the proximal operator and an accelerated update using the previous estimate. The reader is referred to the references for further details on FISTA. The shrinkage step of FISTA is by far the most computationally complex (the accelerated update uses only basic arithmetic), taking the form:

[TABLE]

As such, the computational cost of FISTA is closely linked to that of evaluating the proximal operator. For the $\ell^{1}$ regularizer, the proximal operator has a simple closed-form solution as soft thresholding

[TABLE]

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 [44]. 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_{\ell^{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 a manageable 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. However, a reasonably sized reconstruction volume ($1024\times 1024\times 1024$ 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 $\ell^{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% [33, 15].

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), thus 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}$ of the maximum intensity of the image. This value is selected as any values below it would be 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.

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 ([45, 46, 47, 48, 38]) have tracked fewer than 10 particles concurrently. 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 cubic reconstructed voxels. The largest volume reconstructed by Endo et al. [39] 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 [49] (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.

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. [50] 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\times 2160$ pixels) is cropped to a $1024\times 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 $\mu$m (600 planes). Measurement of similar samples using DIH has been undertaken by Dixon et al. [30] who measured the diffusion of nanowires and Kempkes et al. [51] who demonstrated a $2^{\circ}$ 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 Figures 1a and 1b respectively alongside renderings of the reconstructed volumes produced using deconvolution, RIHVR with sparsity regularization, and RIHVR with fused lasso regularization. The figure shows 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 $\mu$m using deconvolution to 0.63 $\mu$m and 0.89 $\mu$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, the smoothing effect of the fused lasso is apparent (Figure 1d and e). 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 Figure 1) illustrate that both RIHVR approaches approximate the true circular shape of the wire. Conversely, deconvolution (Figure 1c) produces an X-shaped cross-section characteristic of simple holographic reconstructions. RIHVR also demonstrates robustness to image noise. The raw image (Figure 1a) has a substantial amount of background noise and even enhancement via background removal does not produce a noise free image (Figure 1b). Additional fringe patterns – caused by vibrations, fluctuations in illumination intensity, and out-of-view objects – are visible in 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 [52, 53, 54]. The simulation domain is scaled to $5\times 5\times 5$ mm3 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 $\mu$m and 3.7 frames respectively. The RMS velocity is 6.7 $\mu$m/frame. Maintaining the Reynolds number and using a low viscosity fluid ($\nu=10^{-7}$ m2/s), this corresponds to a frame rate of 75 kHz which is achievable with modern cameras. The particles are initially randomly spatially distributed throughout the 3D domain and their positions at subsequent time steps are determined using a Lagrangian tracking method [54]. 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\times 512$ pixel images with a 10 $\mu$m pixel size and 632 nm illumination wavelength.

The reconstruction plane spacing is equal to the lateral pixel spacing (10 $\mu$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_{\ell^{1}}=0.5$ and $\lambda_{FL}=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 [55] with a maximum per frame displacement of 70 $\mu$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 [49]. The second is the deconvolution method of Toloui & Hong [33] 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 Figure 2a.

To evaluate the localization error, extracted particles are matched to their true location using a nearest neighbor method [55]. The resulting error distributions in each dimension are summarized in Figure 2b. 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.

For turbulence measurements, it is common to measure Reynolds stresses which are velocity fluctuation statistics [6]. Here, we present measurements of the root-mean-square (RMS) velocity (Figure 2c). 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$\times$ 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, and image resolution. 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 1 does not completely explain the extraction rate in all situations [15, 33]. For comparison, Toloui et al. [13] performed measurements with a concentration of 0.0035 particles/pixel while other sources used significantly lower concentrations. Using RIHVR, concentrations of up to 0.035 particles/pixel can be processed while maintaining $Ep>60\%$ (Figure 3a). The increased number of extractable particles enabled by RIHVR (Figure 3b) 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 a case 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) [56, 5]. 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 [57]. Cells have a length of 10 $\mu$m and swim using two flagella [58]. In this study, D. primolecta is grown at $37^{\circ}$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$ mm3 glass cuvette. Holograms are recorded at 100 Hz using a $2048\times 1088$ pixel sensor (Flare 2M360-CL). The sensor pixel size is 5 $\mu$m and a 5x microscopic objective is used. The recorded sample volume is $2.05\times 1.09\times 1$ mm3. For simplicity and speed, the recorded image is cropped to a size of $1024\times 1024$ pixels ($1\times 1$ mm2). The light source is a 532 nm diode laser (Thorlabs CPS532) which is expanded and filtered with a spatial filter (see Figure 4a). While the number of particles per pixel is relatively low for this sample (0.002), the particles are large enough that the shadow density (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 [58], 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.

The holographic volume is reconstructed with 270 planes, separated by 3 $\mu$m, with the volume confirmed to include both walls of the cuvette. The regularization parameters are $\lambda_{\ell^{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. Images are preprocessed by removing the mean of the sequence using the method of [4] (subtraction followed by division by the square root). 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.

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 Figure 4b. For clarity, only a subset of 25% of the tracked data is shown in Figure 4 while the full density is shown in Figure 5. 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 4c 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 (Figure 5).

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 Figure 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 [59] 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 [60, 61, 62]). 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 [63, 61] or multi-camera 3D measurements of individual fibers [62]. 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 2 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 $\mu$m diameter SiC ($\rho=3$ g/cm3) 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 $St_{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. A microscopic objective (Edmund Optics, 10x, NA=0.45) is used to image the sample, resulting in a $1024\times 1024$ image with a pixel size of 0.91 $\mu$m. The light source is a spatially filtered HeNe laser ($\lambda=632$ nm). The FL regularization method is used with $\lambda_{\ell^{1}}=0.1$ and $\lambda_{TV}=0.12$. 110 reconstruction planes are used with a spacing of 9.1 $\mu$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 [51]).

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 [59] 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 [64, 65]:

[TABLE]

Where $\boldsymbol{p}$ is a unit vector aligned with the particle axis, $\boldsymbol{\Omega}$ is the rotation tensor, and $\boldsymbol{S}$ is the strain rate tensor. Because the particle rotation rate is coupled to the orientation, the rotation rates for the simulation (Figure 6c) assume that the particles are initially aligned with the inlet flow direction.

The experimental fiber orientations are shown in Figure 6a 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 Figure 6b. The optical reconstruction direction corresponds to the crossflow ($z$) axis in these figures. 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\lvert\dot{p}\right\rvert$) is higher near the centers of the vortices (Figure 6d) which matches the expected behavior from the simulation (Figure 6c). 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 (9) 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 Conclusions

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.

Funding

National Science Foundation (NSF) (1427014); University of Minnesota Informatics Institute.

Acknowledgments

We acknowledge the assistance of Santosh Kumar and Jiaqi You in performing measurements presented here.