Strategies for threedimensional particle tracking with holographic video microscopy
Abstract.
The video stream captured by an inline holographic microscope can be analyzed on a framebyframe basis to track individual colloidal particles' threedimensional motions with nanometer resolution. In this work, we compare the performance of two complementary analysis techniques, one based on fitting to the exact LorenzMie theory and the other based on phenomenological interpretation of the scattered light field reconstructed with RayleighSommerfeld backpropagation. Although LorenzMie tracking provides more information and is inherently more precise, RayleighSommerfeld reconstruction is faster and more general. The two techniques agree quantitatively on colloidal spheres' inplane positions. Their systematic differences in axial tracking can be explained in terms of the illuminated objects' light scattering properties.
§ I. Introduction
Particle tracking based on quantitative analysis of video microscopy images (1) has become a standard technique in several branches of science and engineering. Timeresolved particle trajectories obtained from sequences of video images have been used to study the interactions between colloidal particles (2), to address fundamental questions in statistical physics (3), to probe the viscoelastic properties of soft media (4) and to measure the dynamical properties of single polymers (5), among many other applications. Most of these studies have focused on motion parallel to the microscope's focal plane because threedimensional tracking with conventional techniques requires extensive calibrations and accesses only a limited axial range (6). Recent advances in inline holographic video microscopy (7); (8); (9); (10) promise to overcome these limitations, providing access to threedimensional trajectory data with nanometer resolution over very large volumes, with only the overall magnification of the optical train requiring calibration.
Here, we report on a headtohead comparison between two complementary approaches to analyzing holographic video snapshots. The first involves fitting each image to predictions of the LorenzMie theory of light scattering (9); (10). It provides a wealth of information with outstanding resolution (9); (10); (11), but is computationally very demanding and works only for particles whose lightscattering properties are described well by the theory. The second is more general and less computationally expensive, but relies on phenomenological interpretation to yield particletracking information. This approach uses Fourier diffraction theory (12) to approximately reconstruct the threedimensional light field responsible for the recorded hologram (7); (8); (13) and associates features in this reconstruction with the structure (14) and orientation (7); (8); (14) of the sample. Using measurements on micrometerscale colloidal spheres as a test case, we demonstrate that these two approaches agree quantitatively, except for a systematic offset in the measured axial positions, which can be accounted for quantitatively with LorenzMie optics.
§ II. Inline Holographic Video Microscopy
Our inline holographic microscope, depicted schematically in Fig. 1(a), is based on a commercial inverted light microscope (Nikon TE 2000U). The collimated beam from a 10 HeNe laser (Uniphase, ) replaces the conventional incandescent illuminator. The resulting irradiance of roughly is comparable to that of conventional microscope illumination and is too weak to exert measurable forces on the illuminated objects or to heat the sample appreciably. An object at position scatters a small portion of the incident plane wave, as depicted in Fig. 1. The scattered light then propagates to the focal plane of the microscope, where it interferes with the unscattered portion of the laser beam. The resulting interference pattern is magnified by the microscope's objective lens (Nikon Plan Apo, , numerical aperture 1.4, oil immersion) and projected by a video eyepiece () onto a CCD camera (NEC TI324AII), which records its intensity at 30 frames/s with an exposure time of 1 . The resulting video signal is recorded as uncompressed digital video with a digital video recorder (Pioneer DVR520H) before being analyzed.
This system provides a calibrated magnification of 0.135 pixel. When used with conventional illumination and analyzed with standard methods of digital video microscopy (1), it typically yields a micrometerscale sphere's inplane position to within 1/10 pixel, or roughly 10 . With appropriate calibration, a sphere's axial position also can be measured to within 100 (1). The axial range, however, is restricted to roughly by the limited depth of focus of the highnumericalaperture objective lens. Illuminating the sample instead with coherent light provides access to more information at much higher resolution and over a much larger axial range, without requiring detailed calibrations. Extracting this information, however, requires a substantially more sophisticated analysis strategy.
The unprocessed holographic image of a colloidal sphere in Fig. 1(b) illustrates some of the challenges of inline holographic microscopy. In addition to the sphere's lightscattering pattern, this image is marred by speckle and is obscured by fringes due to all of the other surfaces and particles along the optical train. Many of these distracting features could be mitigated by reducing the correlation length of the illumination (15), for example with a rotating diffuser (16). This, however, would reduce the benefits of holographic imaging relative to conventional incoherent imaging. Most of the unwanted interference fringes could be scaled out of the spatial bandwidth of the imaging system by illuminating the sample with a diverging laser beam. Collimated illumination, however, presents a much simpler scattering geometry that lends itself to precision analysis.
§ III. Holographic Particle Tracking
We model the incident illumination as a plane wave,
(1) 
propagating along with wavenumber in a medium of refractive index . Both the transverse amplitude profile and the polarization are assumed to be independent of . The wave scattered by the sample,
(2) 
propagates in three dimensions with complex amplitude and spatially varying polarization . The measured intensity at point in the focal plane is due to the superposition of the incident and scattered waves,
(3)  
(4) 
If we assume that most of the static defects in the recorded image arise from variation in the illumination, then they can be removed by normalizing with a measured background image (8); (10); (11); (14). The success of this procedure can be judged by the normalized hologram in Fig. 1(c), whose visual quality approaches that of a conventional microscope image. Assuming, furthermore, that the object's height above the focal plane is greater than its size, the geometric rotation of the scattered polarization will be comparatively small, so that . This approximation can be improved by deliberately defocusing the microscope. An optically isotropic sample's normalized hologram then is described by
(5) 
where the reduced scattered field is .
Equation (5) is the point of departure for the two approaches to holographic particle tracking that we will compare. The first involves fitting the recorded images to predictions of the exact LorenzMie theory for light scattering. The second involves tracking features in the threedimensional light field reconstructed from a recorded hologram using the RayleighSommerfeld diffraction integral.
§ III.1. LorenzMie Fitting
The scattered field at position in the focal plane due to an object at relative to the center of the focal plane may be written as
(6) 
where is the normalized scattering function. In this case, Eq. (5) reduces to
(7) 
Given an analytic form for , a normalized hologram may be fit to Eq. (7) for the position of the particle as well as any free parameters in the scattering function (9).
LorenzMie theory describes the scattering function as an expansion (17)
(8) 
in the vector spherical harmonics, and , whose coefficients, and , depend on the size, shape, composition, and orientation of the object and on the structure of the illuminating field. For a homogeneous isotropic sphere of radius and refractive index illuminated by a linearly polarized plane wave, the expansion coefficients are expressed in terms of spherical Bessel functions and spherical Hankel functions as (17)
(9)  
(10) 
where and where primes denote derivatives with respect to the argument. This form assumes , which is appropriate for our microscope. The scattering coefficients fall off rapidly with order , and the series is found to converge after a number of terms (17); (18), which typically is less than 30 for micrometerscale latex spheres in water. To compute in practice, we use the accurate but numerically intensive continued fraction algorithm due to Lentz for the expansion coefficients (19) and the more efficient recurrence algorithm due to Wiscombe for the vector spherical harmonics (20).
The images in Fig. 2 are the normalized image of a 1.51 diameter polystyrene sphere in water (Polysciences Lot 526826) and a pixelbypixel nonlinear leastsquares fit to Eqs. (7) through (10) for , and . Fits are performed with the MPFIT implementation (21) of the LevenbergMarquardt algorithm (22) whose rigorous estimates for the adjustable parameters' uncertainties suggest that the sphere's center has been identified with a precision of 1 in all three dimensions. This estimate of the measurement resolution has been confirmed by independent measurements of similar colloidal particles' dynamics (9); (10). Similar partperthousand resolution is achieved for the particle's radius and refractive index. The quality of the fit can be seen in the azimuthally averaged results that also are plotted in Fig. 2. In this particular example, the particle is found to have a radius of , a purely real refractive index of , and is centered at above the focal plane.
Fits to LorenzMie theory can be this successful because the interference fringes in each holographic image can be defocused enough to cover thousands of pixels. In conventional brightfield imaging, by contrast, a particle at such a large height above the focal plane would be effectively invisible.
If the illumination were not sufficiently collimated, the incident wavefronts' curvature would introduce positiondependent distortions into the holograms that would not be accounted for by Eq. (7). We find the fit values for a particle's size and refractive index to be independent of its position in the microscope's field of view. Wavefront curvature therefore should not affect estimates for particles' positions.
§ III.2. RayleighSommerfeld BackPropagation
The hologram of a colloidal particle also can be used to reconstruct the threedimensional light field that created the hologram. Although such reconstructions may be treated as volumetric renderings of the scattered object (7); (8); (23), they are more properly interpreted as renderings of the scattered field, or its intensity, . Such volumetric reconstructions can be computed from any hologram without knowledge of the scattering object's shape, size or composition. They also are much less computationally expensive than LorenzMie fits. Features in the reconstructed light field may be associated with the positions and orientations of objects in the scattering volume (7); (8); (14). Volumetric reconstructions therefore promise highspeed modelfree measurements of multiple objects' positions and motions in three dimensions. How objects' actual positions are related to features in their backprojected holograms has not been assessed quantitatively in previous studies, however. To address this, we compare the precision measurements of the previous section with threedimensional particle tracking results obtained with RayleighSommerfeld back propagation.
The wave scattered by a particle spreads as it propagates. Its amplitude, consequently, should be substantially smaller than the illumination's once it reaches the focal plane. In that case, the third term in Eq. (5) is negligible compared with the other two and
(11) 
If the complex scattered field were completely specified in the focal plane, it could be reconstructed at height above the focal plane as the convolution,
(12) 
of the scattered amplitude in the focal plane with the RayleighSommerfeld propagator (12)
(13) 
where . The sign convention for accounts for the object's position upstream of the focal plane. Equation (12) may be rewritten with the Fourier convolution theorem as
(14) 
where
(15) 
(16) 
To use this formalism to reconstruct the scattered field, we note that the Fourier transform of is
(17) 
where is the Fourier transform of . From this,
(18) 
may be recognized to be the superposition of the scattered field at height above the focal plane and a spurious field due to the object's mirror image in the focal plane, which is known as the twin image. The twin image's influence on the reconstructed field is reduced by defocusing, which also improves the approximations underlying Eq. (11). In the further approximation that the illumination is uniform, the reconstructed scattered field is
(19) 
The associated intensity, , is an estimate for the scattered light's intensity at height above the focal plane. The example in Fig. 3(a) shows a typical volumetric reconstruction in 0.135 slices for a 1.5 diameter polystyrene sphere in water, the reconstructed rays converging to the diffractionlimited bright red spot roughly 10 above the focal plane.
Because the projection operation in Eq. (12) assumes the medium to be homogeneous and isotropic, it cannot account for light's distortion by the object itself. Consequently, the reconstructed field upstream of an object is degraded by artifacts, as is the case for in Fig. 3(a). Reconstruction artifacts overlap the images of objects that are located upstream of a scattering particle. Nevertheless, the volumetric image in Fig. 3(b) (Media 1) of the intensity scattered by a collection 1.58 diameter silica spheres in water (Duke Scientific Lot 24169) shows clearly resolved maxima, even when spheres are occluded by their neighbors (8). These 22 spheres were arranged with holographic optical traps (26); (27); (28) into a bodycentered cubic (bcc) lattice with a 10.8 lattice constant. The colored regions in Fig. 3(b) (Media 1) are isosurfaces of the brightest 1 percent of the reconstructed voxels and clearly suggest the positions of all of the spheres in the lattice. Interpreting such a complicated hologram with fits to LorenzMie theory would be challenging.
Whereas LorenzMie fitting requires a specific model for each object's light scattering properties, RayleighSommerfeld backpropagation also can be applied to objects whose structures are not known, or for which no numerically stable scattering function is available. To illustrate this, the normalized holograms in Figures 4(a) and (c) were obtained with two 1.58 diameter silica spheres, the second of which had a 40 thick cap of permalloy (80% Ni and 20% Fe) applied to roughly on octant. Asymmetric structure in the latter particle's hologram arises from the the superposition of two contributions, the symmetric scattering pattern due to the sphere and another less well characterized contribution from its cap. Asymmetries in the hologram change form and direction as the capped sphere rotates. This structure also is apparent in the volumetric reconstructions in Figs. 4(b) and (c). Whereas the uncapped sphere gives rise to a rotationally symmetric intensity pattern, the capped sphere's volumetric image has obvious structure. Both the location of the capped sphere and an indication of its orientation are captured by the volumetric reconstruction in Figure 4, whose intensity maximum contrasts with those of ordinary spheres in that it is asymmetric.
It is tempting to associate the positions (1) of local maxima in the reconstructed intensity distribution with the actual positions of objects in the sample (7); (8). The approximations used in deriving Eq. (19) coupled with the nonlinear phenomenological centroid identification (1) make it difficult to assess this method's accuracy and resolution a priori. Comparison with results of LorenzMie analysis addresses this need.
§ IV. Comparison of Strategies
Figure 5(a) shows 10 in the 10 duration trajectory of a 1.5 diameter polystyrene sphere diffusing in water sampled at 1/30 and analyzed with both LorenzMie fitting and RayleighSommerfeld backpropagation. The two approaches agree quantitatively on the particle's inplane position, as is shown by the histogram of differences in Fig. 5(b). This distribution has peaks at 0.35 pixel offsets because brightnessweighted centroid detection is biased toward pixels' centers unless particular care is taken to match the detection window to the size of the object (1). Because the resolution of the LorenzMie fits is consistently better than 1 over the entire trajectory, the difference can be ascribed entirely to errors in interpreting the RayleighSommerfeld reconstructions.
Figure 5(c) shows the distribution of differences of axial positions measured with RayleighSommerfeld (RS) back propagation and LorenzMie (LM) fitting as a function of . These differences are normally distributed about and vary only slightly with axial position.
This systematic offset can be ascribed to the fundamental difference in what the two methods track. The LorenzMie approach uses the exact theory for light scattering by a sphere. It therefore yields an estimate for the actual threedimensional position of the sphere's center relative to the center of the focal plane. The RayleighSommerfeld approach, by contrast, seeks the center of brightness in the reconstructed scattering pattern, which generally falls downstream of the sample's center. This distinction has been noted in passing by previous studies (7); (8). How it varies with the sample's position and physical properties has not previously been addressed.
Although the axial offset plotted in Fig. 5(c) varies only weakly with axial position, the data in Fig. 5(d) for colloidal silica and polystyrene spheres of varying radii show a much stronger dependence on both composition and size. This means that backpropagation cannot reliably be used to measure threedimensional separations in polydisperse and heterogeneous samples. LorenzMie analysis accounts explicitly for individual spheres' size and refractive index and so should provide accurate estimates for axial positions regardless of size and composition.
Trends in the observed axial offsets constitute a selfconsistency test both for the analyses and also for our instrumentation. The solid curves in Fig. 5 represent the distance from the center of a sphere to the point of maximum intensity in its forward scattering pattern computed with Eqs. (8) through (10). This distance should correspond to the difference in axial position obtained with RayleighSommerfeld and LorenzMie analyses. Quantitative agreement between these predictions and our measurements on two different materials over an order of magnitude in particle size supports both the correctness of our analysis and also the efficacy of our experimental implementation.
The axial offset is comparable to the measurement error in the RayleighSommerfeld method for the smallest spheres in Fig. 5(d). This observation suggests that objects with features smaller than the wavelength of light may be tracked in three dimensions at high resolution using RayleighSommerfeld back propagation, even if the shape and composition of the object is not known. This, in turn, explains why RayleighSommerfeld backpropagation has proved successful at tracking the threedimensional positions and orientations of diffusing metal oxide nanorods, two of whose dimensions are smaller than the wavelength of light (14).
Although RayleighSommerfeld backpropagation provides less information than LorenzMie fitting, it is more flexible and substantially faster. Computing each reconstructed axial slice requires roughly a dozen floatpoint array operations and a Fourier transform. Depending on the desired axial range and resolution, a few hundred slices may be computed for a single volume. By contrast, each evaluation of the LorenzMie scattering function with Eq. (8) requires well over a hundred floatpoint array operations, and each fit requires as many as a hundred evaluations. Ten full RayleighSommerfeld reconstructions consequently can be computed in the time required for a single LorenzMie fit. Identifying and locating maxima in the reconstructions does not add substantially to the processing time (1).
We reduce the computational overhead of both analysis methods by performing array calculations on the graphical processing unit (GPU) of a highend graphics card using the GPUlib library (10); (29). Hardwareaccelerated analysis typically yields ten volumetric reconstructions per second or one LorenzMie fit, which is fast enough for nearrealtime characterization of individual colloidal particles (9); (10); (11); (30), visualization of threedimensional flows (10), measurement of rheological properties (31) and performance of medical diagnostic tests (10).
§ V. Conclusion
Our measurements on colloidal spheres demonstrate that RayleighSommerfeld backpropagation can be a fast, flexible and effective method for tracking objects in three dimensions. Axial positions obtained by tracking intensity maxima in the reconstructed light field differ systematically from the true position for spheres larger than the wavelength of light. These deviations depend on the spheres' size and refractive index, and thus can lead to errors in estimates for relative separations among dissimilar objects in three dimensions. Fortunately, particles smaller than the wavelength of light have negligibly small offsets, suggesting that RayleighSommerfeld back propagation may be most useful for analyzing the motions of small objects whose shapes sizes and compositions are unknown.
Our comparison of these techniques over an order of magnitude in particle size suggests that the measurement error in RayleighSommerfeld tracking can be better than 100 in all directions over the entire accessible axial range. Although this compares poorly with the nanometer resolution of LorenzMie tracking, RayleighSommerfeld tracking appears to be more effective for manybody systems and can be applied to objects of any shape (14).
Alternatives to LorenzMie theory such as the discrete dipole approximation (DDA) or finite difference timedomain (FDTD) calculations may be able to provide the basis for highresolution tracking and characterization for more highly structured samples than isotropic homogeneous spheres. In all such cases, RayleighSommerfeld backpropagation can provide a useful starting point for the fitting process, provided that account is taken of the size and compositiondependent axial offset.
§ VI. Acknowledgements
This work was supported in part by the National Science Foundation through Grant Number DMR0922680 and in part through the MRSEC program of the NSF through Grant Number DMR0820341.
References

(1)
J. C. Crocker and D. G. Grier, “Methods of digital video microscopy for colloidal studies,” J. Colloid Interface Sci. 179, 298–310 (1996).

(2)
J. C. Crocker and D. G. Grier, “Microscopic measurement of the pair interaction potential of chargestabilized colloid,” Phys. Rev. Lett. 73, 352–355 (1994).

(3)
G. M. Wang, E. M. Sevick, E. Mittag, D. J. Searles, and D. J. Evans, “Experimental demonstration of violations of the second law of thermodynmaics for small systems and short time scales,” Phys. Rev. Lett. 89, 050601 (2002).

(4)
T. G. Mason, K. Ganesan, J. H. van Zanten, D. Wirtz, and S. C. Kuo, “Particle tracking microrheology of complex fluids,” Phys. Rev. Lett. 79, 3282–3285 (1997).

(5)
T. T. Perkins, D. E. Smith, R. G. Larson, and S. Chu, “Stretching of a single tethered polymer in a uniform flow,” Science 268, 83–87 (1995).

(6)
C. Gosse and V. Croquette, “Magnetic tweezers: Micromanipulation and force measurement at the molecular level,” Biophys. J. 82, 3314–3329 (2002).

(7)
J. Sheng, E. Malkiel, and J. Katz, “Digital holographic microscope for measuring threedimensional particle distributions and motions,” Appl. Opt. 45, 3893–3901 (2006).

(8)
S.H. Lee and D. G. Grier, “Holographic microscopy of holographically trapped threedimensional structures,” Opt. Express 15, 1505–1512 (2007).

(9)
S.H. Lee, Y. Roichman, G.R. Yi, S.H. Kim, S.M. Yang, A. van Blaaderen, P. van Oostrum, and D. G. Grier, “Characterizing and tracking single colloidal particles with video holographic microscopy,” Opt. Express 15, 18,275–18,282 (2007).

(10)
F. C. Cheong, B. Sun, R. Dreyfus, AmatoGrill, K. Xiao, L. Dixon, and D. G. Grier, “Flow visualization and flow cytometry with holographic video microscopy,” Opt. Express 17, 13,071–13,079 (2009).

(11)
F. C. Cheong, K. Xiao, and D. G. Grier, “Characterization of individual milk fat globules with holographic video microscopy,” J. Dairy Sci. 92, 95–99 (2009).

(12)
J. W. Goodman, Introduction to Fourier Optics, 3rd ed. (McGrawHill, New York, 2005).

(13)
J. GarciaSucerquia, W. Xu, S. K. Jericho, P. Klages, M. H. Jericho, and H. J. Kreuzer, “Digital inline holographic microscopy,” Appl. Opt. 45, 836–850 (2006).

(14)
F. C. Cheong and D. G. Grier, “Rotational and translational diffusion of copper oxide nanorods measured with holographic video microscopy,” Opt. Express 18, 6555–6562 (2010).

(15)
F. Dubois, L. Joannes, and J. C. Legros, “Improved threedimensional imaging with a digital holography microscope with a source of partial spatial coherence,” Appl. Opt. 38, 7085–7094 (1999).

(16)
J. GarciaSucerquia, J. H. Ramírez, and R. Castaneda, “Incoherent recovering of the spatial resolution in digital holography,” Opt. Commun. 260, 62–67 (2006).

(17)
C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley Interscience, New York, 1983).

(18)
M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption and Emission of Light by Small Particles (Cambridge University Press, Cambridge, 2001).

(19)
W. J. Lentz, “Generating Bessel functions in Mie scattering calculations using continued fractions,” Appl. Opt. 15, 668–671 (1976).

(20)
W. J. Wiscombe, “Improved Mie scattering algorithms,” Appl. Opt. 19, 1505–1509 (1980).

(21)
C. B. Markwardt, “Nonlinear least squares fitting in IDL with MPFIT,” in Astronomical Data Analysis Software and Systems XVIII, D. Bohlender, P. Dowler, and D. Durand, eds. (Astronomical Society of the Pacific, San Francisco, 2009).

(22)
J. Moré, “The LevenbergMarquardt algorithm: Implementation and theory,” in Numerical Analysis, G. A. Watson, ed., vol. 630, p. 105 (SpringerVerlag, Berlin, 1977).

(23)
B. Rappaz, P. Marquet, E. Cuche, Y. Emery, C. Depeursinge, and P. J. Magistretti, “Measurement of the integral refractive index and dynamic cell morphometry of living cells with digital holographic microscopy,” Opt. Express 13, 9361–9373 (2005).

(24)
G. C. Sherman, “Application of the convolution theorem to Rayleigh's integral formulas,” J. Opt. Soc. Am. 57, 546–547 (1967).

(25)
U. Schnars and W. P. O. Jüptner, “Digital recording and reconstruction of holograms,” Meas. Sci. Tech. 13, R85–R101 (2002).

(26)
E. R. Dufresne and D. G. Grier, “Optical tweezer arrays and optical substrates created with diffractive optical elements,” Rev. Sci. Instrum. 69, 1974–1977 (1998).

(27)
D. G. Grier, “A revolution in optical manipulation,” Nature 424, 810–816 (2003).

(28)
M. Polin, K. Ladavac, S.H. Lee, Y. Roichman, and D. G. Grier, “Optimized holographic optical traps,” Opt. Express 13(15), 5831–5845 (2005).

(29)
P. Messmer, P. J. Mullowney, and B. E. Granger, “GPULib: GPU computing in highlevel languages,” Comp. Sci. Engin. 10, 70–73 (2008).

(30)
K. Xiao and D. G. Grier, “Multidimensional optical fractionation with holographic verification,” Phys. Rev. Lett. 104, 028302 (2010).

(31)
F. C. Cheong, S. Duarte, S.H. Lee, and D. G. Grier, “Holographic microrheology of polysaccharides from Streptococcus mutans biofilms,” Rheol. Acta 48, 109–115 (2009).