Holographic time-resolved particle tracking by means of three-dimensional volumetric deconvolution

Holographic particle image velocimetry allows tracking particle trajectories in time and space by means of holography. However, the drawback of the technique is that in the three-dimensional particle distribution reconstructed from a hologram, the individual particles can hardly be resolved due to the superimposed out-of-focus signal from neighboring particles. We demonstrate here a three-dimensional volumetric deconvolution applied to the reconstructed wavefront which results in resolving all particles simultaneously in three-dimensions. Moreover, we apply the three-dimensional volumetric deconvolution to reconstructions of a time-dependent sequence of holograms of an ensemble of polystyrene spheres moving in water. From each hologram we simultaneously resolve all particles in the ensemble in three dimensions and from the sequence of holograms we obtain the time-resolved trajectories of individual polystyrene spheres.


Introduction
Holographic particle image velocimetry (HPIV) is an established technique for tracking particle trajectories in time and space by means of optical holography [1][2][3][4]. Off-axis [5][6][7][8][9], inline  and hybrid schemes [12,32] of experimental arrangements have been proposed allowing a three-dimensional reconstruction of the wave scattered by the particle field. One of the critical issues in HPIV is assigning locations of individual particles from the reconstructed scattered wave. While the lateral position of a particle in the plane orthogonal to the optical axis can relatively easy be determined, the precise position of the particle along the optical axis (z-axis) is difficult to allocate due to the extended depth-of-focus signal along z-direction. In an optical system, the image of an ideal point source appears as an Airy disk with the central maximum remaining pronounced over an extended defocusing distance along the zaxis. In fact, 80% of the maximal intensity of the central peak is still apparent at the defocus distance [2,[32][33]: where  is the effective angular aperture (half-angle) of the hologram and  being the wavelength. The depth of focus of an ideal point source thus amounts to 2. For realistic sizelimited particles, the depth of focus is even larger and can reach 40 times of the particle size [15,20]. The axial position of the reconstructed particle can roughly be assigned to the maximum of the intensity of the reconstructed wavefront [13-14, 19, 28] or by a threshold of the reconstructed intensity [17][18]34]. Several techniques have been proposed for trying to find the exact position of a particle in axial direction. It has been found that the amplitude of the reconstructed complex-valued wave exhibits oscillations in axial directions with the minimum at the exact particle position [35]. Pan et al demonstrated that the particle position can also be found from the minimum of the variance of the reconstructed imaginary part of the scattered wave [15]. Fournier et al. proposed to implement a window-like function into the reconstruction algorithm which reduces the number of oscillations in the reconstructed amplitude profile and thus allows finding the particle position more precisely [16]. Also, the entropy method [23] and Hough transform [24] were proposed for finding the axial positions of particles. Alternatively, the position of the particle and its shape can be recovered by fitting the holographic image with the interference pattern predicted by Mie theory [34,36]. In practice, these methods are time-consuming and require the analysis of the wave front distribution of each particle. Recently, we demonstrated that all the positions of particles distributed in three dimensions can be retrieved at once from a holographic reconstruction by three-dimensional volumetric deconvolution with the point-spread function of the optical system [25]. Here we show how this three-dimensional deconvolution can be applied to obtain the positions of an ensemble of moving particles.

Experimental
The experimental setup for a HPIV measurement is shown in Fig. 1. A drop of aqueous solution of 10 m diameter polystyrene spheres is immersed into a cuvette filled with distilled water and the moving spheres are imaged by inline holography. The quartz glass cuvette, placed right adjacent to the microscope objective, exhibits an interior width of 10 mm and a wall thickness of 1.25 mm. The hologram size projected by the microscope objective onto the screen amounts to roughly 25 mm in diameter and the field of view is about 625 m × 625 m. The camera allows capturing digital 10 bit images of 1000 × 1000 pixels corresponding to about 16 pixels for each 10 m diameter polystyrene sphere. The concentration of the spheres in the drop was empirically selected such that the presence of several flowing spheres could be observed at any time in the scene. For obtaining time-resolved trajectories of the moving spheres, a sequence of 500 frames was recorded with an exposure time of 20 ms per frame and a time interval between two successive images of 110 ms. For the particle-tracking analysis presented here 30 frames were selected.

Hologram normalization
For hologram normalization [37][38] a background image is needed. It is provided by summing up all the recorded holographic frames which smears out the signal from the moving spheres but maintains pronounced signals from static object such as scratches on the cuvette surface, see Fig. 2. Normalized holograms were subsequently obtained by dividing the experimental holograms with the background image.

Reconstruction of holograms
The reconstruction of the normalized holograms is based on applying the Kirchhoff scalar diffraction theory with the Fresnel-Sommerfeld solution. The Huygens-Fresnel principle, as predicted by the Fresnel-Sommerfeld solution, can be expressed as follows [39]: where  is the wavelength, r denotes the vector from point P 1 to point P 0 , and n denotes the vector normal to the hologram surface,  the angle between vectors n and r , U(P 0 ) the reconstructed complex-valued field, U(P 1 ) the complex-valued field at the hologram, and the integration is performed over the hologram surface. In the experiments presented here, the Fresnel approximation is not fulfilled; here z is the distance between sample and detector, and (x,y) and (X,Y) are the coordinates in the object and detector plane respectively. Since the Fresnel approximation is not fulfilled, we apply the angular spectrum theory [39] for a plane wave propagating short distances. The complex-valued wave scattered by the object U Object (x,y,z) is reconstructed from the hologram by back-propagation of the optical field from the hologram plane (X,Y) to the object plane (x,y) [39]: where FT and 1 FT  are the Fourier transform and inverse Fourier transform, respectively, and (f x , f y ) are the spatial frequencies.  There is significant background noise that we associate with diffuse scattering in the solvent and possible interference between internally reflected light from the cuvette surfaces. It manifests itself also in the very noisy appearance of the normalized holograms, as shown in Fig. 2 and Fig. 3, which does not allow fitting the hologram of each sphere using Mie theory.

Three-dimensional volumetric deconvolution
The reconstruction of a hologram results in the three-dimensional distribution of the scattered wave. The positions of individual scatterers can be obtained by applying a three-dimensional deconvolution of the reconstructed wave front with the point-spread function (PSF) of the optical system The latter is a complex-valued three-dimensional wave scattered by a point-like object and imaged by the optical system [25]. In [25] we proposed two deconvolution methods: instant and iterative deconvolution and demonstrated their application to experimental holograms of microspheres distributed three-dimensionally. As a result, we were able to obtain the three-dimensional distribution of particle positions at increased resolution in both axial and lateral directions. The algorithm of instant three-dimensional volumetric deconvolution was applied by Dixon et al [27] for tracking of a few objects in the scene. In the work presented here we apply the method of instant three-dimensional deconvolution of the intensity of the reconstructed wavefront with the intensity of the PSF of the optical system. Each hologram in the experimental sequence was resampled with 200 × 200 pixels and reconstructed at a series of distances from the hologram ranging from 1 to 7 mm with a step width of 0.3 mm to obtain a total of 200 reconstructions. Thus, from each hologram we obtain the reconstructed three-dimensional volumetric complex-valued field sampled with 200 × 200 × 200 pixels. The number of pixel is limited by the performance of the MATLAB fftn routine employed for three-dimensional deconvolution. The deconvolution was performed on 64-bit operating system with 16 Gb RAM.

a. PSF
The three-dimensional volumetric deconvolution method does not require any calibration experiments and the PSF of the system is obtained by reconstructing either a simulated hologram of a point scatterer or an experimental hologram of a single scatterer (a cutout of an experimental hologram where just a single scatterer is observed) [25]. In the experiments presented here, the PSF was obtained as the three-dimensional complex-valued wavefront reconstructed from a simulated hologram of a point scatterer, see Fig. 4. The hologram for obtaining the PSF is simulated as follows: a point scatterer is placed at z = 4 mm from the detector plane and its hologram is calculated with the following formula: where (x,y) is the point-scatterer described by a delta-function and assumed to be of the size of one pixel for the numerical simulation. The parameters of the simulated hologram are identical to those of the experiment: a wavelength = 532 nm and an imaged area of 625 m × 625 m. The simulated hologram is reconstructed by back propagating the wavefront calculated by Eq. (3): and the reconstructed complex-valued wavefront U PSF (x,y,z) constitutes the PSF.

b. Results of three-dimensional volumetric deconvolution
Three-dimensional volumetric deconvolution is calculated by applying the following formula to the three-dimensional distributions [ x y z U x y z  , and  is a small constant added to avoid division by zero. In this three-dimensional deconvolution, a numerical apodization window function A F (X,Y,Z) is applied in the Fourier domain: where N F is the cut-off frequency in pixels; we used N F = 70. A F (X,Y,Z) cuts off higher order frequencies that are mostly due to noise. The second numerical filter A O (x,y,z) is applied in the object domain: whereby N O denotes the cut-off limit in pixels and we used N O = 95. This filter allows avoiding accumulation of artefacts at the edges of the three-dimensional reconstruction volume. This three-dimensional volumetric deconvolution results in the three-dimensional distribution of the positions of the individual scatterers. In Fig. 5, the reconstructions of three consecutive holograms before and after deconvolution are shown. Fig. 5. Three successive holograms out of the sequence and their threedimensional reconstruction before and after applying the three-dimensional volumetric deconvolution are shown. The area of the reconstructed volume amounts to 625 m × 626 m × 6000 m.
In Fig. 6, an intensity profile of a single scatterer before and after three-dimensional deconvolution is displayed. Obviously, the resolution in both, axial in lateral directions is greatly improved. As evident from Fig. 6(b), the lateral positions of scatterers before deconvolution are represented as blurred spots. After deconvolution they appear sharp and limited in size to just 2 or 3 pixels, as shown in Fig. 6(c). Another advantage of the threedimensional deconvolution with the three-dimensional PSF is that it allows distinguishing between the intensity maximum caused by the particle and its caustics. Caustics are maxima of intensity formed by the superposition of the waves scattered by neighboring particles, they do not possess the out-of-focus distribution of a particle, see Fig. 4. The deconvoluted distribution exhibits a maximum only when the signal scattered by a particle perfectly matches the PSF being the signal from a perfect point scatterer including its out-of focus distribution, see Fig. 4. The caustics will thus not result in maxima in the three-dimensional deconvolution. As evident from Fig 6, some of the intensity maxima are missing after deconvolution; they are attributed to caustics formed by the superposition of waves scattered by the neighboring particles, as for example indicated by the green arrow in Fig. 6(b) and 6(c). The three-dimensional volumetric deconvolution is thus a much more accurate technique than other methods for determining particle positions. The best condition for a particle to be resolved by three-dimensional deconvolution is when the particle is not too close to an edge of the reconstruction volume since part of its three-dimensional signal is then cut off. In this case the match between the three-dimensional PSF and three-dimensional particle reconstruction is only partial and the deconvolution might miss it. Reconstructed intensity distributions summed up along axial direction before a three-dimensional deconvolution and (c) after three-dimensional deconvolution.

Time-resolved particle tracking
A time-dependent sequence of holograms allows observing the evolution in the threedimensional distribution of particles. Such time-resolved holographic particle tracking has various applications. When only a few particles are present in the scene, their movement can easily be observed by the naked eye and matching of each particle position in different time frames can be done manually. In this manner, for example, evaporation process of droplets in jets [31] and velocity profiles in a tube showing parabolic distribution [19] have been studied. When a larger number of scatterers is flowing in a solvent, as for example when studying flow dynamics, the velocity distribution is typically obtained by the cross-correlation between particle positions [8]. However, the method of cross-correlation requires the position of each particle as an input, and those positions are difficult to find when the number of particles is large. Conventionally, positions of individual particles are extracted and stored as coordinates (x,y,z). Then, a three-dimensional volume is created where pixels corresponding to particle positions are marked. This creates a three-dimensional representation of the particle distribution. In our method, after a three-dimensional deconvolution, single particles are clearly resolved and their trajectories can be followed individually, see Fig. 5. As an example, we have selected four spheres which are close to each other but exhibit different sedimentation velocities, as illustrated in Fig. 7. The sedimentation velocity of each particle is estimated from the linear fitting. The pairs of closely spaced particles, red-green (34 m distance in xz-plane) and blue-cyan (300 m distance in xz-plane) exhibit similar velocities.

Conclusions
We have demonstrated that the three-dimensional volumetric deconvolution can be successfully applied to simultaneously resolve individual positions within an ensemble of particles moving in a solvent. Good reconstructions can be achieved despite the interference between numerous particles and the background due to the solvent. The three-dimensional volumetric deconvolution allows separating between the signal originating from a particle and the caustics. After deconvolution, particle positions are clearly resolved and their tracking allows following their time-resolved trajectories.