1. Introduction
Digital in-line holography of dispersed particles is able to extract the 3D position of particles and provide information on the particle size and shape. These abilities have led to a wide range of applications, including 3D particle image velocimetry in fluid mechanics [
1J. Katz and J. Sheng, “Applications of holography in fluid mechanics and particle dynamics,” Annu. Rev. Fluid Mech. 42, 531–555 (2010). [CrossRef]
] and the measurement of atmospheric particle size distributions [
2J. P. Fugal and R. A. Shaw, “Cloud particle size distributions measured with an airborne digital in-line holographic instrument,” Atmos. Meas. Tech. 2, 259–271 (2009). [CrossRef]
]. Conventional methods employ forward propagation through Fourier methods and image processing methods such as pixel counting to obtain the particle size information [
3S. Soontaranon, J. Widjaja, and T. Asakura, “Improved holographic particle sizing by using absolute values of the wavelet transform,” Opt. Commun. 240, 253–260 (2004). [CrossRef]
–
5J. Lu, J. P. Fugal, H. Nordsiek, E. W. Saw, R. A. Shaw, and W. Yang, “Lagrangian particle tracking in three dimensions via single-camera in-line digital holography,” New J. Phys. 10, 125013 (2008). [CrossRef]
]. These methods lead to limited depth resolution in the particle position estimation and significant uncertainty in size estimation due to the finite CCD pixel size and ambiguity in threshold setting. Considerable efforts have been devoted to overcoming the depth resolution problem [
6G. Pan and H. Meng, “Digital holography of particle fields: reconstruction by use of complex amplitude,” Appl. Opt. 42, 827–833 (2003). [CrossRef] [PubMed]
–
10L. Dixon, F. C. Cheong, and D. G. Grier, “Holographic deconvolution microscopy for high-resolution particle tracking,” Opt. Express 19, 16410–16417 (2011). [CrossRef] [PubMed]
]. For example, Yang et al. [
7W. Yang, A. B. Kostinski, and R. A. Shaw, “Depth-of-focus reduction for digital in-line holography of particle fields,” Opt. Lett. 30, 1303–1305 (2005). [CrossRef] [PubMed]
] developed a Fourier-domain inverse filter, obtained with a knowledge of the particle shape and size, that resulted in a narrow ‘intensity peak,’ which in turn is used for position estimation along the optical axis (depth resolution). The inverse method of the Lyon group [
8F. Soulez, L. Denis, C. Fournier, E. Thiébaut, and C. Goepfert, “Inverse problem approach for particle digital holography: accurate location based on local optimisation,” J. Opt. Soc. Am. A 24, 1164–1171 (2007). [CrossRef]
,
9J. Gire, L. Denis, C. Fournier, E. Thiébaut, F. Soulez, and C. Ducottet, “Digital holography of particles: benefits of the ‘inverse problem’ approach,” Meas. Sci. Technol. 19, 074005 (2008). [CrossRef]
] represents a fundamental advance because it allows for simultaneous improvement of size and position information through an iterative optimization algorithm.
The aim of this work is to outline a method for improved particle size estimation from digital in-line holograms, relative to conventional reconstruction and pixel counting methods. The approach is inspired by the filtering method of Yang et al. [
7W. Yang, A. B. Kostinski, and R. A. Shaw, “Depth-of-focus reduction for digital in-line holography of particle fields,” Opt. Lett. 30, 1303–1305 (2005). [CrossRef] [PubMed]
] and the optimization method of Soulez et al. [
8F. Soulez, L. Denis, C. Fournier, E. Thiébaut, and C. Goepfert, “Inverse problem approach for particle digital holography: accurate location based on local optimisation,” J. Opt. Soc. Am. A 24, 1164–1171 (2007). [CrossRef]
]: we employ a matched filtering method for simultaneous determination of particle position and size, on a particle-by-particle basis. The new matched filtering method leads to improved accuracy compared to typical pixel counting methods, and can be implemented as a relatively simple extension of the conventional holographic reconstruction method. We first present an overview of the principle of operation, then a numerical exploration of the method, and finally an experimental demonstration of its actual implementation. The contribution of this paper is to introduce the proposed method and its underlying theoretical basis, and to make an initial characterization of its inherent precision and performance in the presence of noise. The method is compared to a standard reconstruction and pixel counting approach, and detailed comparison with other filtering or inverse methods is the focus of future work.
2. Principle of operation
In the reconstruction of a digital hologram we can consider the profile of intensity along an axis centered on a particle and parallel to the optical axis, here referred to as the ‘axial profile.’ The position of the particle along the optical axis can be determined from the resulting axial profile, but not necessarily from the maximum along the profile since it may consist of several peaks depending on the particle size, the wavelength, and the detector pixel size [
11C. Fournier, C. Ducottet, and T. Fournel, “Digital in-line holography: influence of the reconstruction function on the axial profile of a reconstructed particle image,” Meas. Sci. Technol. 15, 686–693 (2004). [CrossRef]
]. The method proposed here is relevant when the particle size is not known. In that case, we can apply a Fourier filter for a range of different particle sizes and check the resulting profile peak intensity; when the filter size matches the true particle size, the peak intensity is highest. This filter, therefore, provides improved estimates of both particle size and depth position.
This section starts with a brief description of filter implementation in the Fourier domain and proceeds to introduce the proposed matched filter, and to analyze some of the basic properties of the reconstructed field at the particle plane and along the optical axis. The typical digital in-line holographic configuration contains only collimating optics and a CCD that records the interference intensity pattern between the collimated on-axis reference light and the forward-scattered light from the particles. The light intensity in the hologram plane (
s,t) is
i(
s,t) = 1 −
o(
s,t) −
o*(
s,t) +
|o(
s,t)
|2, where
* denote the complex conjugate. After omitting the contribution from the higher-order (diffraction) term, the Fourier transform of
i(
s,t) can be approximated as
, where
P(
η,
ξ) is the Fourier transform of the particle aperture function
p(
x,y) in the particle object plane (
x,y) at position
z0 along the optical axis, and
Hz(
η,
ξ) = exp[−
iπλz(
η2 +
ξ2)] [
12J. Goodman, Introduction to Fourier Optics , 2nd ed. (McGraw-Hill, Boston, Mass., 1996).
]. For a spherical particle
p(
x,y) can be assumed to be a circular aperture with the diameter of the particle. Then
P(
η,
ξ) is proportional to
, the angular spectrum of the circular aperture with radius
a, where
J1 is the Bessel function of the first kind. Omitting the effect of the twin image, and denoting
F−1{} as the inverse Fourier transform, the reconstructed field in an arbitrary plane (
u,v) at distance −
z from the hologram is
We can apply an arbitrary filter by multiplying
P′(
η,
ξ) with the Fourier spectrum of the recorded hologram, such that the reconstructed field becomes
When
P′ is the
z-dependent Gaussian function or the inverse of the angular spectrum of detected particles, the axial profile of the reconstructed field has a maximum intensity value at
z =
z0, as described by Fournier et al. [
11C. Fournier, C. Ducottet, and T. Fournel, “Digital in-line holography: influence of the reconstruction function on the axial profile of a reconstructed particle image,” Meas. Sci. Technol. 15, 686–693 (2004). [CrossRef]
] and Yang et al. [
7W. Yang, A. B. Kostinski, and R. A. Shaw, “Depth-of-focus reduction for digital in-line holography of particle fields,” Opt. Lett. 30, 1303–1305 (2005). [CrossRef] [PubMed]
], respectively. Here we propose defining a filter
P′(
η,
ξ) as the sign function of the angular spectrum of a circular aperture with radius
a′:
This binary filter is characterized by radius
a′ of a preselected circular aperture function, and the resulting reconstructed field is equivalent to the diffraction field from a modified aperture
pm,a(
u, v) =
F−1{
P′(
η,
ξ) ×
P(
η,
ξ)} (i.e., at
z =
z0 the propagation kernel is
Hz0−z(
η,
ξ) = 1 and therefore the reconstructed field has the simple form
ez0 (
u, v) =
pm,a(
u, v), which is the form of a matched correlation function.)
Since both the filter function and the spectrum of a detected particle aperture function are real and circularly symmetric, the modified aperture function
pm,a will be real and circularly symmetric. Most importantly, when radius
a′ matches the detected particle radius
a, the modified aperture function, or the reconstructed field at the focus
z =
z0, is sharply peaked as illustrated in
Fig. 1. This feature can be understood through Fourier theory [
13J. D. Gaskill, Linear Systems, Fourier Transforms, and Optics (Wiley, New York, 1978).
]: when the filter size is matched, all spectrum values in the Fourier domain are positive and the center value of the modified aperture is the sum of all these positive values; in contrast, off-center values of the modified aperture are the sum of these positive values, but they are modulated by different harmonic waveforms that contain positive and negative values. Therefore we can expect that the modified aperture is a circular symmetric function that is center peaked. As seen in the figure, the modified aperture function allows some negative values and does not have a clear aperture edge; however, it has a significantly intensified field value at the center. The intensified center values effectively reduce the size of aperture, thereby resulting in a greatly reduced depth of focus [
6G. Pan and H. Meng, “Digital holography of particle fields: reconstruction by use of complex amplitude,” Appl. Opt. 42, 827–833 (2003). [CrossRef] [PubMed]
,
7W. Yang, A. B. Kostinski, and R. A. Shaw, “Depth-of-focus reduction for digital in-line holography of particle fields,” Opt. Lett. 30, 1303–1305 (2005). [CrossRef] [PubMed]
]. (Details of the axial profile of the reconstructed field will be discussed in the next section.) When
a′ of the filter does not match with the size of the detected particle we expect the modified aperture still to be a circularly symmetric real function, but with much lower center values and an overall wider aperture compared to the matching case. This is because when not matching, the sign filter would not be able to turn all the spectrum values to positive values, thus the center value of the modified aperture would be smaller than the the center value when size matches. Therefore a much wider depth of focus and lowered axial field peak can also be expected.
Fig. 1 The filter-modified aperture (on the reconstructed field at focus z = z0) when a′ matches with the detected particle size a. The black curve is the modified aperture function, exhibiting greatly intensified center values. In contrast, the aperture function before applying the filter is shown as the blue curve. This computation assumes a circular aperture of 250 pixels in radius on a computational array of 2048 × 2048.
It should be noted that the form and properties of the sign matched filter are similar to those of the “phase-only” matched filter [
14J. L. Horner and P. D. Gianino, “Phase-only matched filtering,” Appl. Opt. 23, 812–816 (1984). [CrossRef] [PubMed]
] well known in pattern recognition. The angular spectrum of a circular aperture is a Jinc function (i.e.,
J1(
x)/
x), which contains no explicit phase angle variables, however if we confine the amplitude to be positive only, the positive or negative sign of the Jinc function can be substituted by exp(
i ·
π) or exp(−
i ·
π). In effect, the sign filter captures the binary phase of the aperture angular spectrum. As such, we can expect that the sign matched filter should also have the properties of low sensitivity to noise and high sensitivity to features of the detected object as exist for the phase-only matched filter [
14J. L. Horner and P. D. Gianino, “Phase-only matched filtering,” Appl. Opt. 23, 812–816 (1984). [CrossRef] [PubMed]
].
This filter approach is analogous to that described by Yang et al. [
7W. Yang, A. B. Kostinski, and R. A. Shaw, “Depth-of-focus reduction for digital in-line holography of particle fields,” Opt. Lett. 30, 1303–1305 (2005). [CrossRef] [PubMed]
], in which the transformed intensity is divided by
P(
η,
ξ). In that case the result is identical to the field generated by a point source, having an infinite intensity peak at
z =
z0. In practice, however, the multiplication by
P′(
η,
ξ) as described in this paper avoids the noise problems that arise from division by very small values or zeros in
P′(
η,
ξ). As a consequence this method is much more robust to additive noise and can be applied for
P′(
η,
ξ) over a range of particle sizes.
It is now possible to apply this ‘sign matched filtering’ method for a specified range of particle sizes near the initially estimated particle size. That means, we assume different particle sizes and for each size we apply the corresponding matched filter
P′(
η,
ξ) and then reconstruct the hologram as in
Eq. (2). We can then select the maximum value from each axial intensity profile (i.e., for each
P′(
η,
ξ)) and use the maximum values to construct a plot of maximum intensity versus particle size
a′ for the applied filter. The highest intensity corresponds to the sharpest axial profile, that for which
a′ =
a, the radius of the actual particle.
3. Numerical results
We first explore the improvement in particle size estimation using simulated in-line holograms. The wavelength of illuminating light is taken to be
λ = 0.527
μm and the hologram is a 4008 × 2672 array with a pixel size of 9
μm. The particle diameter is
d = 56
μm. Panel (a) of
Fig. 2 shows the axial (i.e., along the center of the particle) intensity profiles without imposing the matched filtering, analogous to
Eq. (1). Panel (b) shows the axial intensity profiles for a matched filter corresponding to the actual particle size,
P′(
η,
ξ) = sgn[
P(
η,
ξ)]. This demonstrates how the peak intensity profile changes to a single peak through the multiplication of
P′(
η,
ξ), which greatly improves the depth of focus [
7W. Yang, A. B. Kostinski, and R. A. Shaw, “Depth-of-focus reduction for digital in-line holography of particle fields,” Opt. Lett. 30, 1303–1305 (2005). [CrossRef] [PubMed]
]. We note that here we have generated
P′(
η,
ξ) directly from the Bessel function, as in
Eq. (3) instead of a direct pixel-generated method. Finally, panels (c) and (d) show axial intensity profiles for matched filters with diameters 10% less than and 10% greater than the actual particle diameter, respectively. The intensity scale for panels (b), (c), and (d) are all normalized to the peak intensity for
P′(
η,
ξ) = sgn[
P(
η,
ξ)], so it can be seen that the single, central peak in the axial profile is maintained, but that it becomes less pronounced as the filter diameter changes from the actual particle diameter. An additional check, the results of which are not shown in detail here, was to vary the lateral position of the particle: up to the maximum tested distance of 1/40 of the hologram width from the hologram edge, there was no significant degradation in the peak quality.
Fig. 2 Axial intensity profiles verses depth along the optical axis. Panel a shows the conventional, unfiltered result and panel b shows the axial profile after filtering with P′(η,ξ) = sgn[P(η,ξ)], i.e., for filter diameter equal to the particle diameter of 56 μm. The bottom two images show the axial intensity profiles for matched filters corresponding to the actual diameter minus and plus ten percent (50 and 62 μm respectively). The intensity scales for panels b, c, and d are normalized by the same factor, such that the best matched filter results in a maximum intensity of unity. The plots are made using 10 μm steps along the z-axis.
When selecting different filter size parameters a single peak may remain, but the peak intensity of the profile changes. If selecting the exact or true size parameter, the profile becomes most narrow, and the peak intensity reaches the highest value. In this way the maximum value for each different filtering size parameter is obtained, as shown in
Fig. 3. It can be seen that the maximum value in the maximum-intensity versus filter-size plot corresponds to the correct particle size. The origin of this peaked curve can be better understood by noting that it is identical to the function
pm,a(
u, v) =
F−1{
P′(
η,
ξ) ×
P(
η,
ξ)}, with
P′ varying with the filter diameter
a′ and the peak appearing when
a′ =
a.
Fig. 3 The peak-intensity versus filter-size, showing a maximum when filter size is equal to the true particle size.
An important consideration for any method to be applied in practical systems is its robustness in the presence of noise. An initial investigation of sensitivity of the sign filter method to noise has been made in two ways. First, the dependence of the diameter estimation on additive gaussian white noise has been checked. The additive noise has a signal-to-noise ratio (SNR) defined as
μ/
σ, where
μ is the mean of the hologram ‘signal’ and
σ is the gaussian standard deviation. Results for two representative values, SNR of 1 and 1/10, are shown in
Fig. 4: plots in the left column are profiles of intensity along the optical axis, and plots in the right column are profiles of peak height versus filter diameter. Even for the very low SNR case the peak is quite evident. It is noted that the width of the curve in the peak height versus filter diameter profiles does not undergo any significant change with increasing noise. By taking 30 independent realizations of single particle holograms with SNR of 1 a diameter mean of 55.94
μm and standard deviation of 0.11
μm are obtained. In practice each size estimate is accomplished via a parabolic fitting to obtain the position of the peak. The second noise test considers the influence of particle density in the hologram. The third row in
Fig. 4 shows results of depth and diameter profiles for a realization of a hologram containing 1000 particles. The 1000 particles are randomly placed within a 4008 × 9
μm by 2672 × 9
μm by 2672 × 9
μm box. By taking 30 independent realizations of particle positions a diameter mean of 55.97
μm and a standard deviation of 0.10
μm are obtained. These initial characterizations of additive and multiple-particle noise indicate that the proposed filter method is robust under a variety of realistic conditions for implementation. This is further borne out in the following section with experimental results. A more fundamental analysis of the Cramer-Rao bound [
15C. Fournier, L. Denis, and T. Fournel, “On the single point resolution of on-axis digital holography,” J. Opt. Soc. Am. A 27, 1856–1862 (2010). [CrossRef]
] for diameter estimation with the sign filter is a topic of ongoing work.
Fig. 4 Axial intensity profiles verses depth along the optical axis (left) and peak-intensity versus filter-size (right) for three different noise levels. The top panels are for additive Gaussian white noise with SNR of 1; the middle panels are for the same but with SNR of 1/10. The bottom panels are for ‘noise’ resulting from 1000 particles within the roughly cube-shaped volume.
4. Experimental results
To verify the effectiveness of this method, we also obtained analyzed hologram recorded in the laboratory. The optical system consists of a frequency doubled Nd:YLF laser (
λ = 0.527
μm). The beam is passed through a beam expander with the addition of a pinhole for spatial filtering. The camera is an Illunis XMV-11000 with 4008 × 2672 with 9
μm pixels and 12-bit output. Water droplets are injected by a piezoelectric injector that generates monodisperse droplets. The mean droplet diameter is
d̄ = 56.1
μm with an upper bound on the size variability of Δ
d = 0.5
μm, as determined from a high resolution, telecentric microscope with a spatial resolution of 0.9
μm (the diameter can be determined to better than the spatial resolution by considering the full area of the image). The experimental data consists of 136 holograms taken at
z0 = 531 mm, with one particle taken from each hologram. We note that the particle is taken at the same location in the hologram to avoid any possible evaporation effects. First, we use conventional reconstruction and pixel counting to obtain the coarse particle size and location. Then each particle in the hologram was processed using the sign matched filtering method described here. We then repeat the reconstruction with filters
P′(
η,
ξ) corresponding to particle diameters in 0.2
μm increments, obtaining the peak axial-profile intensity as a function of the particle filter size, as shown in
Fig. 5. The peak is around 56
μm as expected.
Fig. 5 The peak-intensity versus filter-size for a single hologram recorded in the lab. The filter size varies from 44 μm to 68 μm in 0.2 μm steps. The maximum agrees with the actual particle diameter to within experimental uncertainty.
The red histogram in
Fig. 6 shows the statistical results after reconstructing the 136 holograms and using the filter method to obtain the best particle size estimate. For comparison, the blue histogram shows the results from a direct reconstruction and pixel counting method (pixel counting based on a global intensity threshold). When using pixel counting the mean particle diameter is
d̄ = 55.7
μm and the standard deviation is
σd = 2.8
μm. This standard deviation is consistent with the uncertainty scaling with square root of pixel size [
5J. Lu, J. P. Fugal, H. Nordsiek, E. W. Saw, R. A. Shaw, and W. Yang, “Lagrangian particle tracking in three dimensions via single-camera in-line digital holography,” New J. Phys. 10, 125013 (2008). [CrossRef]
]. For the matched filtering method the mean size is
d̄= 56.0
μm and the standard deviation is
σd = 0.23
μm. This standard deviation is consistent with the expected variability from the piezoelectric droplet generator. This clearly demonstrates the ability of the new method to greatly improve the precision of particle size estimation.
Fig. 6 Distributions of estimated particle diameters from 136 holograms. The blue histogram, with d̄ = 55.7 μm and σd = 2.8 μm, is for the the conventional pixel counting method based on a global intensity threshold. The red histogram, with d̄ = 56.0 μm and σd = 0.23 μm, is for the same holograms analyzed with the sign matched filter method.