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 [1] and the measurement of atmospheric particle size distributions [2]. Conventional methods employ forward propagation through Fourier methods and image processing methods such as pixel counting to obtain the particle size information [3–5]. 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 [6–10]. For example, Yang et al. [7] 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 [8,9] 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. [7] and the optimization method of Soulez et al. [8]: 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 [11]. 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)|², 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 $I(\eta ,\xi )=\delta (\eta ,\xi )-P(\eta ,\xi )H_{z_0}(\eta ,\xi )-P^{*}(\eta ,\xi )H_{z_0}^{*}(\eta ,\xi )$, where P(η,ξ) is the Fourier transform of the particle aperture function p(x,y) in the particle object plane (x,y) at position z₀ along the optical axis, and H_z(η,ξ) = exp[−iπλz(η² +ξ²)] [12]. 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 $J_1\left(a\sqrt{{\eta }^2+{\xi }^2}\right)/\left(a\sqrt{{\eta }^2+{\xi }^2}\right)$, the angular spectrum of the circular aperture with radius a, where J₁ is the Bessel function of the first kind. Omitting the effect of the twin image, and denoting F⁻¹{} 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

Since both the filter function and the spectrum of a detected particle aperture function are real and circularly symmetric, the modified aperture function p_m,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 = z₀, is sharply peaked as illustrated in Fig. 1. This feature can be understood through Fourier theory [13]: 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 [6, 7]. (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 = z₀) 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.

Download Full Size | PDF

It should be noted that the form and properties of the sign matched filter are similar to those of the “phase-only” matched filter [14] well known in pattern recognition. The angular spectrum of a circular aperture is a Jinc function (i.e., J₁(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 [14].

This filter approach is analogous to that described by Yang et al. [7], 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 = z₀. 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 [7]. 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.

Download Full Size | PDF

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 p_m,a(u, v) = F⁻¹{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.

Download Full Size | PDF

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 [15] 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.

Download Full Size | PDF

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 z₀ = 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.

Download Full Size | PDF

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 [5]. 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.

Download Full Size | PDF

5. Summary and conclusions

In this paper, we have proposed and verified using simulated and experimental data, that particle size estimation from in-line holograms can be greatly improved through a sign matched filtering method, compared to typical reconstruction and pixel counting methods. The Fourier-domain filter is equivalent to the sign operator on the angular spectrum of a preselected aperture function. The approach is easily implemented within the conventional hologram reconstruction method, applying a sign matched filter in the Fourier domain. When the aperture function has a diameter equal to that of a measured particle there is a strong peak in the reconstructed field. A range of filters for different particle sizes is applied, and the maximum in the peak-intensity versus filter particle size gives the correct particle size. As implemented here, conventional reconstruction and pixel counting is used to obtain coarse particle size information, and then the size estimate is refined through the filter method. Furthermore, when the accurate size information is obtained we simultaneously have obtained an improved estimate of the particle position along the optical axis (improved depth resolution). It is plausible that the filter method discussed here could be used directly as a particle detection method, although computational expense would need to be considered. Finally, we note that the method is robust to additive noise and can be applied to a wide range of particle sizes, even if those particles are in the same hologram.