Side scattered light from micrometer-sized particles is recorded using an off-axis digital holographic setup. From holograms, a volume is reconstructed with information about both intensity and phase. Finding particle positions is non-trivial, since poor axial resolution elongates particles in the reconstruction. To overcome this problem, the reconstructed wavefront around a particle is used to find the axial position. The method is based on the change in the sign of the curvature around the true particle position plane. The wavefront curvature is directly linked to the phase response in the reconstruction. In this paper we propose a new method of estimating the curvature based on a parametric model. The model is based on Chebyshev polynomials and is fit to the phase anomaly and compared to a plane wave in the reconstructed volume. From the model coefficients, it is possible to find particle locations. Simulated results show increased performance in the presence of noise, compared to the use of finite difference methods. The standard deviation is decreased from 3–39 μm to 6–10 μm for varying noise levels. Experimental results show a corresponding improvement where the standard deviation is decreased from 18 μm to 13 μm.
© 2017 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
Digital holography is a suitable technique for tracking micro-particles since it records 3D information. Holography was originally proposed by Gabor in 1948 , and its digital version is usually attributed to Goodman . Extracting information about particle fields is a natural application of digital holography because 3D information about position, size, and morphology is recorded in the hologram. Applications include studies of flow phenomena using holographic particle velocimetry (HPIV) , tracking microorganisms , and characterizing aerosols [5–7] and the dynamic behavior of combustion fuel particles [8,9] and air bubbles , to mention a few. There are two main configurations of holographic setups used for particle analysis. In in-line holography, the illumination wave is also used as a reference wave. Here the scattered light from particles interferes with the un-scattered light in the forward direction. This is a very simple and compact setup, but it has disadvantages in that only the forward direction can be used, and a twin image and zero-order terms disturb the reconstruction. The other main configuration is off-axis holography, which uses a separate reference wave with an off-axis angle. This setup is more complex, but with it, it is possible to remove the twin image and zero-order content from the image of interest. Any imaging direction can be used with this setup, not only the forward direction. In this investigation, the off-axis configuration is utilized because it allows for a general registration direction. In digital holography, the hologram is recorded on a CCD or CMOS detector and numerically reconstructed to obtain 3D information [11,12]. For low numerical apertures, the resolution in the reconstruction along the optical axis is poor compared to that of the lateral one ; therefore, particles appear elongated. To find the correct axial particle position, a metric of some kind is needed. Previously reported metrics are based on intensity [14–16], which are most common, or they are based on the whole complex field [17,18]. It is also possible to solve the inverse problem directly from the recorded hologram [19,20]. This is, however, more time consuming and, therefore, unsuitable for large particle sets.
We have, in a previous paper, presented a metric that is based on the wavefront curvature of the scattered light . This metric utilizes the fact that scattered light from a particle converges toward its true position and thereafter diverges. The change from a converging to a diverging wave is found by tracking the wavefront curvature. One significant advantage with a phase-based metric is that it does not have an offset that is present in intensity-based metrics . The wavefront curvature in our previous paper was estimated from phase gradients that were computed using finite differences. This was a simple way to estimate the curvature; unfortunately, in the presence of noise, the precision decreases rapidly due to the noise enhancement. In this paper, we propose a method that parameterizes the phase anomaly compared to a plane wave, without using finite differences. The method consists of fitting the phase anomaly to a three-dimensional polynomial model, and from the polynomial coefficients the particle position is determined. By using a model, we reduce the effect of noise. The rest of the paper is structured as follows: in Section 2, theory regarding particle scattering, digital hologram formation, and reconstruction are presented. The proposed polynomial model is also outlined. In Section 3, the setup is described along with the simulations and measurements performed. In Section 4, the results from both simulations and measurements are presented and discussed. Finally, the paper is summarized in Section 5, and the conclusions are presented.
Portions of this work were presented at the Digital Holography & 3-D Imaging conference in 2017, Paper M4A.3.
A. Particle Holograms
Light incident on a particle generates an angularly dependent outwards propagating scattered field. Consider the geometry sketched in Fig. 1. The particle is located at the origin of a right handed Cartesian coordinate system and is illuminated by a plane wave propagating along the -axis. Detection is done along the vicinity of the -axis, which is denoted as the optical axis. We denote the Jones vectors of the incident light and the scattered light . For spherical micrometer-sized particles, the Lorenz-Mie theory is a suitable scattering model . The theory is most commonly expressed with respect to the scattering plane; this is the plane spanned by the illumination and observation directions, denoted and , respectively, see Fig. 1. The incident and scattered fields are then expressed in polarization components parallel and perpendicular to the scattering plane. These components are denoted , for the incident light and , for the scattered light, respectively. The relationship between these components is as follows:1) and the full expression for the scattering coefficients, we refer the reader to .
In digital particle holography, the object light is given by Eq. (1), and a reference wave is given by11]. In our setup, the light is recorded by the two-lens telecentric imaging system, depicted in Fig. 2. The magnification between the object and image space is then given by , where and are the focal lengths of lens and lens , respectively. We choose to denote the coordinate system in the object space by and the coordinates in the image plane on the detector by . The imaging system then maps from the object plane to the image plane. The coordinates are related by the magnification of the system as , , and , where and are the defocus distances in object and image space, respectively. The reference wave is an off-axis point source in the aperture plane that is reshaped by the lens to become a plane wave. In these notations, the intensity on the detector for an arbitrary polarization becomes 3(a) shows the recorded intensity on the detector for a simulated hologram. The particle is defocused at distance .
The hologram captured by the detector has to be reconstructed before further analysis is performed. It is possible to extract individual interference terms when using off-axis holography, since they are separated into different lobes in the Fourier domain. The first four terms in Eq. (3), that is, the intensity components, are located in a central lobe. The four interference terms each have a separate side lobe in a well configured system. The interference term of interest is extracted by applying a rectangular mask around the corresponding lobe, removing all content outside of it. If the reference wave is a tilted plane wave, it will scale the object wave and add a phase gradient proportional to the off-axis angle. It is possible to remove the phase gradient by re-centering the lobe after the filtering process. By doing so, reduces to a scalar depending on the reference wave intensity. Because of this, we drop the from now on; hence, the complex amplitude is known in Fourier space for the component of interest, where is the axial position of the detector in image space. It is now possible to compute the field at an axial position from the detector with the use of the angular spectrum method , and the field becomes
The basic principle of our metric is that the reconstructed wave from a particle converges toward its position and thereafter diverges. Converging and diverging waves have different signs of curvature, so the plane where the curvature changes sign is then the plane where the particle is located. The wavefront curvature can be obtained from the phase in the reconstructed volume. The phase in the volume oscillates quickly along the optical axis, much faster than the step length in the reconstruction. This is problematic, since it then contains discontinues in the z-direction, making it unsuitable to fit it to a polynomial model. Figure 3(b) shows the reconstructed phase along the z-direction along the center of the particle. Aliasing effects are present, causing the oscillations to appear slower than they really are. By instead comparing the phase to a plane wave situated in the same plane, it becomes possible to compare the data between two planes. This phase anomaly is calculated as follows:
C. Chebyshev Polynomial Model
It is assumed that the phase anomaly around a single particle only depends on the response from that particle. Only a small volume around a potential particle is therefore used in the parameterization. This assumption holds in general if the particle distribution is sparse. Potential particles are found from a coarse intensity-based search where voxels with an intensity over a defined threshold are chosen.
For a single particle the phase anomaly, , is fit to a model based on Chebyshev polynomials of the first kind . Chebyshev polynomials are chosen as base functions for the parameterization for two reasons. The first is that they are orthonormal in a Cartesian coordinate system and, therefore, suitable for computations in a rectangular volume. Second, terms of different orders are associated with properties in the imaging process when using small to moderate numerical apertures. For instance, first-order terms are related to a spatial shift in direction, and second-order terms are related to defocus. The modeled phase anomaly in three dimensions is defined from one-dimensional Chebyshev polynomials. The polynomials are defined in the interval , and the coordinates are rescaled to normalized coordinates in this domain. The 3D polynomial is then as follows:6) can be rewritten in matrix form as
The objective is then to find the set of estimated coefficients that best represents the measured phase anomaly . We minimize the squared error between the measured data, and the estimated data, , from the polynomial model. The coefficients are then the solution to the following minimization problem:6), represent the phase anomaly correctly, a sufficiently high polynomial degree is needed. To find this value, scattered light from a single particle is simulated using Eqs. (1) and (3), respectively. For the full description of the simulation implementation, we refer the reader to . The simulated hologram is reconstructed according to Section 3. In the reconstructed volume, the parameterization is performed and the mean squared error (MSE) per voxel is evaluated for different values of . The result is shown in Fig. 4. A good fit is obtained for , where the mean error is approximately . Using higher orders does not significantly improve the fit. Therefore, a fit with is assumed to be sufficient and will be used further on in this paper. The phase anomaly corresponding to Fig. 3(b) is shown in Fig. 3(c). The markers are the simulated data, and the solid line is the Chebyshev polynomial fit. Figure 3(d) shows phase anomaly cross-sections in the y-direction for different z-positions. Here the change in wavefront curvature becomes apparent.
The three-dimensional data is now represented by the coefficient array . The relationship between the elements in the array and information about the particle is now investigated. A single particle is simulated, and coefficients are estimated for different axial positions by translating the computational box along the optical axis. The dimension of the box used is voxels. All coefficients are investigated, and as an example, the results for the coefficients and are shown in Fig. 5. The coefficient has a linear behavior around the particle and becomes zero at the particle position. is therefore used as a metric to find the axial position of the particle. This coefficient corresponds to the quadratic term in the direction perpendicular to the scattering plane, so it describes the curvature of the wavefront in that direction. When becomes zero, the wavefront changes curvature. This behavior is consistent with the findings in . The corresponding coefficient , which represents the quadratic term in the direction parallel to the scattering plane, is also shown in Fig. 5. It shows a similar behavior, but it is slightly shifted compared to . Therefore, only will be used as a metric.
To use the metric on a set of particles recorded in a single hologram, some practical considerations have to be made. First, we need to segment particles from each other, which is done by applying a threshold on the intensity in the volume. From the binary volume, connected voxels are clustered and labeled. The point with maximum intensity has a slight offset from the particle position, but it is still in the vicinity . This point is therefore a good starting point when searching for the plane where becomes zero. A small 3D domain is formed around the maximum intensity point, and the Chebyshev model is evaluated in this domain. The position where becomes zero is then the position of the particle. This position is found by computing the zero-crossing using linear interpolation. It is important that the particles are sparse, so that only scattered light from an individual particle appear in the 3D domain. The whole procedure for obtaining the axial position of particles from a hologram can therefore be summarized in the following steps:
- 1. Reconstruct a three-dimensional volume using the angular spectrum method.
- 2. Apply a threshold on the intensity in the volume to obtain a binary volume.
- 3. Automatic segmentation of each particle from the binary volume.
- 4. Find the point of maximum intensity for each particle.
- 5. Evaluate the coefficient in a small domain around the point of maximum intensity for each particle.
- 6. The axial location where the coefficient becomes zero is stored as the position of each particle.
Consider the setup shown in Fig. 6. A 532 nm Nd:YAG continuous laser with a power of 85 mW is used for illumination. The light is polarized parallel to the scattering plane. This polarization is chosen to reduce Rayleigh scattering from inhomogeneous parts in the sample. The light is split into two parts using a 50/50 beam splitter (BS). The first part of the light is expanded by the lenses L1 and L2, with focal lengths and , respectively. The expanded light illuminates the sample, which consists of a silicon cube with spherical particles molded into it. The particles have an average diameter of 10 μm, and the refractive indices of the silicon and particles are 1.4 and 1.5, respectively. The concentration of particles is low so that the imaged scene is sparse, and only a few particles are present. The cube is placed so that only scattered light from the particles is imaged. There could be internal variations of the refractive index in the silicone cube that affect the experimental results compared to if the particles were imaged in air. These will then show up as an increase in random error. Scattered light from the particles is recorded using a telecentric imaging system, formed by lenses L3 and L4 at a 90° angle. Lenses L3 and L4 have focal lengths and , respectively. This corresponds to a system magnification of . A aperture is used to band limit the image. It is placed in the joint focal plane between lenses L3 and L4, making the system telecentric. The other light part from the BS is used as the reference wave. It is focused into a polarization maintaining fiber using a mirror (M) and a fiber port (FP). The fiber end is placed in the aperture plane. By doing so, we ensure that the reference light becomes a point source in the Fourier domain. In the back focal plane of lens L4, a detector (Sony XCL 5005) is placed to record the hologram. It has pixels and a pixel pitch of 3.45 μm, giving the detector a total size of . The images are stored in an 8-bit format. The fiber has a numerical aperture of 0.12, which is sufficient to cover the whole detector. The numerical aperture of this system is in the object plane and in the image plane. Since the sample is a silicon cube, the effective numerical aperture is even smaller, since the silicon-air boundary will refract the scattered light. The effective numerical becomes
B. Simulations and Measurements
To compare how the parametric Chebyshev method compares to the previously reported finite difference method, a large set of simulations was performed. Once again we use the simulation routine described in , that is, an implementation of Eqs. (1) and (3), respectively. The simulation properties are set as described in Section 3.A. Some assumptions are made: the illumination is a perfect plane wave, and particles are identical perfect spheres with a diameter of 10 μm. Then, 10 holograms each containing 20 particles located at random lateral and axial positions are simulated. Since the axial positions are known, it is easy to compute the error for each particle. The standard deviation, false positive ratio, and false negative ratio are computed. This is repeated for different noise levels ranging from 0 to 4 percent of the maximum intensity.
For experimental measurements, the backside of the silicon cube is placed at approximately in front of the imaging system, so that the focus plane is inside the sample cube. Unlike in simulations, the true particle position is not known, and there are several holograms of stationary particles recorded in the series. Particle positions are evaluated with the method outlined in Section 4 from the holograms. A set of positions for five different particles are then obtained. For each particle, the axial mean position is estimated and subtracted from the recorded positions. The associated errors are therefore computed.
4. RESULTS AND DISCUSSION
The resulting standard deviations from the simulations are shown in Fig. 7 as a function of the noise ratio. Data for both the Chebyshev model and the finite difference method are presented. It is clear that by using the Chebyshev model, the standard deviation drastically decreases for larger noise ratios. At the highest simulated noise ratio of 4%, the standard deviation is reduced from 39 μm to 10 μm, an improvement of about 4 times. This can be explained by the fact that the finite difference is influenced by the noise in two voxels, while the Chebyshev polynomial method uses far more data points. Since the noise is assumed to be zero mean and independent, it will have a smaller effect with the use of more data points. For no or low noise ratios, the proposed polynomial model performed worse than the finite difference method. This is likely due to the extension from 2D to 3D computational domain. In a volume, it is more likely that light from other particles influences the estimation. The false positive ratio is around 1%, so it is very unlikely that noise is mistaken for a particle. This is due to the fact that noise rarely originates from a geometrical point, and hence, it does not have the behavior of an outwards propagating wavefront. Therefore, it does not have the same change in wavefront curvature as a particle.
The false negative ratio increased to around 5% when using the Chebyshev model, compared to 3% when using finite difference. A closer study of the undetected particles shows two main reasons for the increase. The first reason is inaccuracy in the course estimation. The model searches for a zero-crossing of close to the point obtained from the course estimation. So if that estimation is very inaccurate, the model never finds the zero-crossing. The second reason is that particles located near the edges of the reconstructed volume do not use the same number of data points, since parts of the computational domain are shifted outside the reconstruction. This can cause the model to behave inaccurately.
The performance is also compared to the total sum of gradients method, which is commonly used to determine particles axial positions [25,26]. Applying the method to the simulated data set shows that the sum of the gradients method results in a small offset from the true axial position, around 3 μm, which is still inside the particle. The standard deviation for the different noise levels showed a slight increase ranging from 7–12 μm compared to 5–10 μm for the polynomial fit. The two methods, hence, show comparable noise susceptibility.
The results from the experimental measurements are presented in Fig. 8. This figure shows the error histogram after outliers have been removed. The error seems to have a Gaussian profile, where the standard deviation is , which is an improvement compared to the finite difference method, where the standard deviation was 18 μm . Even though the standard deviation decreased, it is still higher than for the simulations. From Fig. 7, the expected standard deviation would be around 7 μm. This can depend on several factors. In the simulations, we assume only single scattering and ignore any effects from the finite size of the silicon cube. For a single particle, scattered light from other particles can be viewed as noise. It is possible that these effects influence the results. Also, aspects like mechanical and thermal stability in the experimental setup influence the results. Still, the standard deviation is in the same order of magnitude as the diameter of the particle. It is, therefore, likely that the measured position is inside or in the close vicinity of the particle.
Finally, a brief discussion of the parametrization process is in order. Consider the errors presented in Fig. 4. For , the error is negligible. However, since the coefficient used as a metric is a second-order term, higher terms are redundant due to the orthonormal properties of the base functions. This means adding higher order terms should not change lower order coefficients. The biggest error reduction is from to . This is not a surprise, since a spherical wave over a small aperture can be approximated by a second-order term. This is also why a second-order term is used as a metric. From Fig. 5, it is clear that only can be used as a metric, since is shifted along the optical axis. These coefficients represent quadratic terms parallel and perpendicular to the scattering plane, respectively. The wavefront parallel to the scattering plane is modulated by Lorenz-Mie scattering coefficients, while the perpendicular direction is not. The modulation distorts the wavefront in the parallel direction, causing to become zero at a shifted position, while is unaffected. This effect has previously been reported by Pu and Meng . They related the modulation effects to aberrations for different observation directions. These effects are important to consider when interpreting results from a certain observation direction. A future extension of the proposed polynomial method to include non-spherical particles could be possible. Non-spherical particles do, in general, depend on the perpendicular direction as well , so the modulation effects need to be handled. This is something that could be addressed in future work. Modulation effects can also be introduced by the reference wave. If the reference is a curved wave instead of a plane wave, the curvature will modulate the phase response. It is therefore important to ensure the quality of the reference wave.
In this paper, a parameterized model based on Chebyshev polynomials is introduced to estimate the wavefront curvature. It is shown that the phase anomaly around a scattering particle can be modeled by a three-dimensional Chebyshev polynomial of degree . To find this degree, the mean error per voxel is evaluated. The polynomial coefficients hold information about the particle position, and by using simulations, the behavior of all obtained coefficients was studied. It is found that the coefficient becomes zero in the axial plane where the particle is located. This coefficient corresponds to the quadratic term perpendicular to the scattering plane. By evaluating this coefficient at different positions and searching for the axial plane where it becomes zero, it is possible to find the true particle position. From simulations, this method has standard deviations between 6 and 10 μm. This should be compared to the standard deviation using the finite difference method previously reported, which were 3–39 μm for the same conditions. The Chebyshev polynomial model performs better in higher noise ratios, since the noise in individual voxels does not affect the result to the same degree. For no or very low noise, the worse performance is due to the extension of the computational domain from 2D to 3D. In a volume, light from other particles is more easily introduced in the computation. Experimental results also show improvement when using this method, from 18 μm to 13 μm. This improvement is not as large as the improvement expected from simulations; however, ideal conditions in simulations are rarely obtainable.
Vetenskapsrådet (VR) (621-2014-4906).
1. D. Gabor, “A new microscopic principle,” Nature 161, 777–778 (1948). [CrossRef]
2. J. W. Goodman and R. W. Lawrence, “Digital image formation from electronically detected holograms,” Appl. Phys. Lett. 11, 77–79 (1967). [CrossRef]
3. H. Meng, G. Pan, Y. Pu, and S. H. Woodward, “Holographic particle image velocimetry: from film to digital recording,” Meas. Sci. Technol. 15, 673–685 (2004). [CrossRef]
4. J. Sheng, E. Malkiel, and J. Katz, “Digital holographic microscope for measuring three-dimensional particle distributions and motions,” Appl. Opt. 45, 3893–3901 (2006). [CrossRef]
5. M. J. Berg and S. Holler, “Simultaneous holographic imaging and light-scattering pattern measurement of individual microparticles,” Opt. Lett. 41, 3363–3366 (2016). [CrossRef]
6. M. J. Berg, N. R. Subedi, and P. A. Anderson, “Measuring extinction with digital holography: nonspherical particles and experimental validation,” Opt. Lett. 42, 1011–1014 (2017). [CrossRef]
7. O. Kemppinen, Y. Heinson, and M. Berg, “Quasi-three-dimensional particle imaging with digital holography,” Appl. Opt. 56, F53–F60 (2017). [CrossRef]
8. Y. Wu, M. Brunel, R. Li, L. Lan, W. Ao, J. Chen, X. Wu, and G. Gréhan, “Simultaneous amplitude and phase contrast imaging of burning fuel particle and flame with digital inline holography: model and verification,” J. Quant. Spectrosc. Radiat. Transfer 199, 26–35 (2017). [CrossRef]
9. Y. Wu, L. Yao, X. Wu, J. Chen, G. Gréhan, and K. Cen, “3D imaging of individual burning char and volatile plume in a pulverized coal flame with digital inline holography,” Fuel 206, 429–436 (2017). [CrossRef]
10. S. I. Satake, Y. Yonemoto, T. Kikuchi, and T. Kunugi, “Detection of microbubble position by a digital hologram,” Appl. Opt. 50, 5999–6005 (2011). [CrossRef]
11. J. W. Goodman, Introduction to Fourier Optics, 2nd ed. (McGraw-Hill, 1996).
12. N. Verrier and M. Atlan, “Off-axis digital hologram reconstruction: some practical considerations,” Appl. Opt. 50, H136–H146 (2011). [CrossRef]
13. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light, 7th ed. (Cambridge University, 1999).
14. Y. S. Bae, J. I. Song, and D. Y. Kim, “Volumetric reconstruction of Brownian motion of a micrometer-size bead in water,” Opt. Commun. 309, 291–297 (2013). [CrossRef]
15. J. K. Abrantes, M. Stanislas, S. Coudert, and L. F. A. Azevedo, “Digital microscopic holography for micrometer particles in air,” Appl. Opt. 52, A397–A409 (2013). [CrossRef]
16. L. Wilson and R. Zhang, “3D localization of weak scatterers in digital holographic microscopy using Rayleigh-Sommerfeld back-propagation,” Opt. Express 20, 16735–16744 (2012). [CrossRef]
17. G. Pan and H. Meng, “Digital holography of particle fields: reconstruction by use of complex amplitude,” Appl. Opt. 42, 827–833 (2003). [CrossRef]
18. J. de Jong and H. Meng, “Digital holographic particle validation via complex wave,” Appl. Opt. 46, 7652–7661 (2007). [CrossRef]
19. F. Soulez, L. Denis, C. Fournier, É. Thiébaut, and C. Goepfert, “Inverse-problem approach for particle digital holography: accurate location based on local optimization,” J. Opt. Soc. Am. A 24, 1164–1171 (2007). [CrossRef]
20. J. Fung, K. E. Martin, R. W. Perry, D. M. Kaz, R. McGorty, and V. N. Manoharan, “Measuring translational, rotational, and vibrational dynamics in colloids with digital holographic microscopy,” Opt. Express 19, 8051–8065 (2011). [CrossRef]
21. J. Öhman and M. Sjödahl, “Off-axis digital holographic particle positioning based on polarization-sensitive wavefront curvature estimation,” Appl. Opt. 55, 7503–7510 (2016). [CrossRef]
22. F. C. Cheong, B. J. Krishnatreya, and D. G. Grier, “Strategies for three-dimensional particle tracking with holographic video microscopy,” Opt. Express 18, 13563–13573 (2010). [CrossRef]
23. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley, 1998).
24. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (Dover, 1972), Vol. 9.
25. Y.-S. Choi, K.-W. Seo, M.-H. Sohn, and S.-J. Lee, “Advances in digital holographic micro-PTV for analyzing microscale flows,” Opt. Lasers Eng. 50, 39–45 (2012). [CrossRef]
26. P. Memmolo, L. Miccio, M. Paturzo, G. D. Caprio, G. Coppola, P. A. Netti, and P. Ferraro, “Recent advances in holographic 3D particle tracking,” Adv. Opt. Photon. 7, 713–755 (2015). [CrossRef]
27. Y. Pu and H. Meng, “Intrinsic aberrations due to Mie scattering in particle holography,” J. Opt. Soc. Am. A 20, 1920–1932 (2003). [CrossRef]