1. Introduction

Interferometric particle imaging (IPI), based on Mie scattering theory, is often used for instantaneously measuring the size and spatial distributions of spherical particles by means of the fringe pattern formed in the out-of-focus image plane. This technique was first proposed by König et al. [1], and has subsequently received considerable attention as a popular method for application to sprays [1–6], flow fields [7], and other situations [8–11]. The particle diameter can be determined by counting the number of fringes, and the measurement accuracy relies mainly on the image processing technique used to estimate the fringe spacing/fringe counts. Generally speaking, the image processing of IPI involves two steps: locating each particle and obtaining the interference pattern, and then performing a Fourier transformation of the interference pattern of each particle to determine the particle size. It is important to accurately locate the particle fringe image in subsequent processing, especially for the subsequent determination of movement and velocity using the particle tracking velocimetry technique. Difficulty also occurs in evaluating the individual particle image of IPI from overlapping out-of-focus image, especially for the circular image with fringes when there is a high degree of overlapping of the fringe pattern in a high particle density region.

To date, among several processing algorithms suggested [3–6, 11–13], the convolution method is most commonly used to find the center of the particle interferogram [4, 12–16]. This method improves positioning accuracy by suppressing the fringes within each particle image. Glover et al. [5] first used the Gaussian blur technique to remove fringes as much as possible, then utilized Canny edge detection and a Hough transform to locate individual droplets in the image field. The fringe spacing of the particle fringe pattern was then extracted by least-squares fitting to a Chip function, which was applied to spatially sparse sprays.

More recently, our research teams [14–16] have developed two new IPI methods and evaluated their performance by simulation and experiment. In this paper, we continue our investigation of IPI algorithms, and present a novel algorithm that simultaneously extracts the location and the fringe count/fringe spacing based on a unidirectional gradient-matched algorithm, a Fourier transform, and the improved Rife method. This algorithm can remove the fringe from within the particle fringe pattern, leaving only the circular border of the particle fringe pattern, similar to the blurred fringe image created by Gaussian blur in [5]. Our algorithm performs very well, and as a result, the particle is located with high accuracy. Additionally, an improved Rife method significantly decreases the probability of the direction of frequency shift being judged incorrectly, so that high-accuracy frequency estimation is achieved. The method is described in detail in this paper, and its performance is validated by simulation and experiment. The method is also compared with that of two previously proposed algorithms for both location determination and frequency extraction. The experimental results match well with those from the simulation.

2. Interferometric particle imaging (IPI)

Figure 1 schematically shows a typical IPI experimental setup. When a transparent spherical particle is illuminated by a laser sheet, the reflected and first-order refracted rays scattered from the particle form two glare points on the focus image plane, and interference fringe patterns on the out-of-focus image plane (The intensity of a pair of scattered rays, the reflected and first-order refracted rays, is much stronger than that of the higher order refraction using homogeneous illumination). The diameter of the particle is related to the separation of the glare-point pair or the fringe spacing of the interferogram. The relationship between particle diameter d and the number of fringes N can be expressed by the following [17–19]

 

Fig. 1 Schematic of interferometric imaging of particle scattering light.

Download Full Size | PDF

3. Algorithm

In the present study, we consider a novel algorithm for analyzing IPI fringe images. Figure 2 shows the algorithm flowchart, which is similar to previously proposed methods [14–16], but the method of extraction of the location and the frequency are different. The edge images of the interferogram and the mask image are extracted respectively by the gradient along the vertical orientation of the interference fringes (here in the x direction), and the geometric center coordinate (x,y) of the particle fringe images can be detected using a 2D correlation operation for the two obtained edge images. The interferogram of each particle is extracted from raw images according to the center coordinate (x,y) and the shape and size of the particle fringe image. A Fourier transform is then performed for each extracted fringe pattern to evaluate the fringe angle spacing and the number of fringes. The sub-pixel accuracy of the extracted frequency is then acquired by the improved Rife algorithm, and the particle diameter can be subsequently calculated by Eq. (1). The particle center coordinate is given through the geometric center coordinate (x,y) of the fringe pattern using geometric relationship.

 

Fig. 2 Flowchart of particle measurement.

Download Full Size | PDF

3.1 Locating the center of particle fringe image

The first and most important step in the interference pattern analysis of IPI is to accurately identify the location of the particle fringe center. The interferogram of particle$I(x,y)$ and the mask image$P(x,y)$are operated on by a gradient algorithm in the x direction, and the x-direction edge images${\nabla }_x\{I(x,y)\}$and${\nabla }_x\{P(x,y)\}$are then extracted. The cross-correlation function between the two obtained edge images is calculated as

 

Fig. 3 Simulation results of the interferogram center-detecting using different algorithm: (a) the interferogram image; (b) the mask image; the results of (c) WMF algorithm, (d) erosion match algorithm and (e) unidirectional gradient-matched algorithm. Shown from left to right are the extracted edge image of (a) and (b), 3D and 2D correlation distribution.

Download Full Size | PDF

3.2 Estimation of fringe frequency

In [14–16], a Fourier transform and the modified Rife algorithm are employed to obtain the sub-pixel accuracy of the extracted frequency from the interference pattern of each particle. The modified Rife algorithm [20] is very much like the interpolation iterative algorithm in [21]. However, when the signal-noise ratio (SNR) is low, the direction of a shift in the frequency domain is likely to be incorrectly judged using the m-Rife method, eventually causing a large error in the frequency reading. Therefore, to avoid this problem, the improved Rife algorithm [22] is proposed by the use of improved judgment criteria, and a higher sub-pixel accuracy of the estimated frequency is acquired. The improved Rife algorithm can be described as follows.

The discrete Fourier transform of signal$s(n)$is defined as

Similarly, for the m-Rife method in [14–16], the frequency f_m is also estimated using Eq. (9) after frequency shift, but the r = ± 1 condition is different, such that $r=1$when $\left|S(k_1+1)\right|>\left|S(k_1-1)\right|$, and $r=-1$ when $\left|S(k_1+1)\right|\le \left|S(k_1-1)\right|$. However, when k₁ is very close to the true peak, $\left|S(k_1+1)\right|$and $\left|S(k_1-1)\right|$ are very close to the minima, and thus noise may reverse their amplitudes. If the amplitudes of $\left|S(k_1+1)\right|$and $\left|S(k_1-1)\right|$ are reversed, Eq. (6) will move the peak to the wrong direction and cause large errors in the frequency reading while the value of$\left|S(k_1\pm 0.5)\right|$is very close to the peak it will not be greatly affected by noise. Thus, the improved Rife algorithm is less noise-sensitive compared with the m-Rife method, which evidently increases the precision of frequency estimation.

4. Numerical results and analysis

4.1 Interference fringe pattern positioning

Here we evaluate the performance of the algorithm by simulation. The IPI numerical simulator was developed by Shen [9], in which the scattering light by the particles is calculated using a Debye series expansion of the Mie theory. The analysis of each particle’s image is very time consuming by this approach, especially for thousands of particles in field of view. Considering that the size d_i of the particle fringe pattern on the out-of-focus plane is independent of particle size, in order that the whole simulator will not be time consuming, we use a simplified model [13]. The IPI fringe is created by randomly distributing disks with fixed diameter d_i which is modulated by the cosine term, and the intensity profile is

To verify the validity of the gradient-matched algorithm we proposed in this paper, comparisons between the gradient-matched algorithm and the erosion matching algorithm, WMF algorithm are presented in this section. First, we examine the performance of the three algorithms: the erosion matching algorithm, WMF algorithm, and the gradient-matched algorithm. Figure 3 shows an example of numerical simulation using Eq. (9). Figure 3(a) is the simulated interferogram with 10 randomly located spherical particles and Fig. 3(b) is the mask image. Figure 3(c) shows sequentially from left to right the edge images of 3(a) and 3(b) extracted using the Mexican-hat wavelet, and 3D and 2D cross-correlation distributions, respectively. Figures 3(d) and 3(e) respectively display the corresponding results of the erosion matching algorithm and the gradient-matched algorithm. One can see from Fig. 3(e) that the fringe within each particle image dies away, leaving just the circular border of the particle fringe pattern as expected, which is significant for center determination. Figures 3(c)–3(e) show that the edge image extracted by the gradient-matched algorithm is finer and sharper than those from the erosion matching method and the WMF algorithm. The peak value of the cross-correlation distribution in Fig. 3(e) is also sharper and its bright center areas are smaller than those of Figs. 3(c) and 3(d). The positions of the interference disk in Fig. 3(a) are completely and accurately identified by the three methods, and Fig. 3(e) is very useful for the detection of particle centers in higher particle density regions.

Furthermore, we also evaluate the center detecting errors for three methods. Figure 4(a) shows the simulated interferogram of 100 spherical particles randomly distributed in a field of view of 2048 × 2048 pixels, where the out-of-focus image size of a particle is 100 pixels. Figure 4(b) shows the detected disk center using the WMF algorithm, erosion matching algorithm, and gradient-matched algorithm, denoted with light green circles, pink squares, and red dots, respectively. Figures 4(c) and 4(d) show the discrepancies of particle image position, $(x_i-x_0)$ and $(y_i-y_0)$, which are the absolute errors of the extracted center coordinates in the x and y directions, where$(x_0,y_0)$is the preset coordinate and ${(x_i,y_i)}_{i=1,2,3}$ are the position coordinates obtained by the erosion matching algorithm, WMF algorithm, and gradient matched algorithm, respectively. One notices that the positions of the fringe patterns of particles in Fig. 4(a) are accurately presented in Fig. 4(b), with a horizontal and vertical precision of ± 3 pixels and ± 0.7 pixels, respectively. As expected, the gradient-matched algorithm has a higher positioning accuracy than the other two algorithms.

 

Fig. 4 (a) Simulation interferogram images; (b) the results of the center detecting using different methods; error distribution in (c) x-direction and (d) y-direction, where the pink “□” indicates erosion matching algorithm, the light green “ο” is WMF algorithm, and the red dot “•” is gradient-matched algorithm, respectively.

Download Full Size | PDF

Finally, we study the recognition ratio R_p for three methods with the different number interferograms from 50 to 10000, the results of which are shown as in Fig. 5. The data points of Fig. 5 are curve fitted by using the Boltzmann model, as depicted in solid line of Fig. 5. For each interferogram, the locations and the spatial fringe frequencies are randomly selected. The overall result of the interferogram detection is similar to the one presented in Fig. 4(b), and the interferogram location coordinate (x, y) can be obtained with high accuracy using the algorithm presented. However the number of identified particles decreased with the particle number density increasing per image, the main source of which overlaps between the interferograms. When the particle number density is low, the recognition ratio is high for fringe patterns with relatively low overlapping. Along with the increasing particle number density, the overlap coefficient goes up and the recognition ratio is further decreased. For R_p = 90%, the number of extracted particles is 208, 308, and 497 for WMF algorithm, erosion matching algorithm and gradient-matched algorithm.

 

Fig. 5 Comparison of recognition ratio R_p of three algorithms.

Download Full Size | PDF

From the above analysis, the gradient-matched algorithm has a higher positioning accuracy and recognition ratio R_p than the other two algorithms because of removing the fringe within each particle image, promising to apply to a high-density spray field. The recognition ratio R_p is related to the overlapped coefficient. Clearly, the higher overlap coefficient reduces the recognition ratio R_p, the detailed analysis of the relations between overlapped coefficient, the recognition ratio R_p and fringe frequency extraction will be reported in other paper.

4.2 Frequency estimation of particle fringe pattern

In this section, we test the performance of the improved Rife algorithm. From Eq. (1), the particle diameter d is proportional to the particle spatial frequency f, the interferogram presented here is simulated using the numerical simulator of IPI developed by Shen [9], the simulation of which is completely different from Section 4.1. Figure 6(a) shows the simulated interferogram of a particle of diameter 40 μm according to the setup in Fig. 1, where Mie theory is used to calculate the scattered intensity of the reflected and first order refracted rays scattered by the particle. The imaging system consists of an aperture and an equivalent thin lens and is modeled by using Fourier optics. The parameters used in the numerical simulation are: λ = 532 nm, the imaging lens is a circular aperture with a focal length f = 50 mm and $F^{\#}$ = 1.8, the magnification M = 0.35, the distance between the CCD plane and the in-focus image plane, i.e. the defocusing distance, is 1.68 mm, then z_r = 69.2 mm, the field area is 201 × 201 pixels with a pixel size of 3.45μm × 3.45μm, the scattering angle θ = 76.4°, and the refraction index m = 1.426. The signal is corrupted by Gaussian noise of 5 dB.

 

Fig. 6 Simulation results of the frequency estimation: (a) interferometric image; (b) 2D frequency spectrum of (a); (c) fine peak detection and frequency estimation.

Download Full Size | PDF

The frequency shift of Fig. 6(a) is performed only along the x-axis. Figure 6(b) shows the 2D Fourier spectrum of Fig. 6(a) and the −1 order spectrum of the 1D frequency spectrum is shown in Fig. 6(c). The discrete frequencies with the maximum power and the secondary maximum of the discrete spectrum $k_1$ = 101–91 = 10, and$\left|S(k_1-0.5)\right|<\left|S(k_1+0.5)\right|$, hence, r = 1. Using Eqs. (4) and (5), Δδ = 0.1252, and then the frequency shift is completed along the x-axis toward the right using Eq. (6). The Fourier spectrum is calculated using Eq. (7), and the values $\left|S(k_1-r\Delta \delta )\right|$ = 6134 and $\left|S(k_1-r\Delta \delta +r)\right|$ = 5167 are obtained. The estimated frequency f_i = 0.0481 lp/pix is calculated using Eq. (9), and as shown in Fig. 6(c), the number of fringes N = 9.64, the particle diameter d = 40.18 μm, and the absolute error and relative error are 0.18 μm and 0.45%, respectively. The extracted frequency f_m = 0.0484 lp/pix is calculated using the m-Rife interpolation [16], and as also shown in Fig. 6(c), the number of fringes N = 9.70, the particle diameter d = 40.42 μm, and the absolute error and relative error are 0.42 μm and 1.1%, respectively. The result proves that a higher precision of frequency estimation is obtained using i-Rife interpolation.

Furthermore, we also numerically explore the performances of estimated frequency by means of Monte Carlo simulation for different frequencies and different SNR. Figure 7 shows the error curve of the estimated frequency obtained by 1000’s of Monte Carlo simulations for the m-Rife and the i-Rife interpolations as a function of SNR, where the frequency is 0.0795 lp/pix. Figure 7(a) shows the root mean square error (RMSE) and Fig. 7(b) shows the mean absolute error (MAE). Here, the SNR is defined as$\text{SNR}=A^2/2{\sigma }^2$, where A is the amplitude and ${\sigma }^2$is variance. The frequency estimation variance of the Cramer-Rao bound is $\mathrm{var}\left\{f\right\}=12{\sigma }^2/[A^2\Delta xN(N^2-1)]$ [23]. It can be seen from Fig. 7 that the i-Rife method is superior to the m-Rife method, especially at low SNR, revealing its insensitivity to noise and a higher precision of frequency estimation, which is highly useful in most practical situations, such as spray.

 

Fig. 7 Variation of RMSE (a) and MAE (b) with respect to SNR

Download Full Size | PDF

5. Experimental results

To validate the performance of the method presented, we have performed a simulated experiment to verify the accuracy of interference fringe pattern positioning and extraction frequency. Now we test the method experimentally, for both particle location detection and particle size measurement accuracy.

5.1 Detection of interferogram center

It is well-known that the exact particle positions are unknown for a real particle field. Consequently, we first need to design an experimental method for test algorithm of the particle centers detection, and then apply this algorithm to the particle field. Figure 8 shows the experimental set-up designed. Light from a 300 W Xenon bulb is compressed by an imaging lens (f = 100 mm) to form a point light source. The light, collimated by a lens (f = 450 mm) and filtered by a filter with a high transmission wavelength of 632.8 nm (Center wavelength is 632.8 ± 0.2nm, Full width half Maximum (FWHM) is 10nm and Peak transmission T>55%), acts as a plane wave with 30 mm diameter, then impinges onto the surface of a SLM (1024 × 768 pixels with a pitch of 10.8 μm), which is controlled via a computer. The modulated beam carrying a simulated interferogram (interference disk) is imaged in the CCD plane by an imaging lens (f = 350 mm), with the imaging system having a magnification of M = 1.797. The CCD used in the experiment is a 10-bit digital CCD sensor (DS-21-04M15 camera) with 2048 × 2048 pixels and pixel size of 7.4 μm × 7.4 μm. Note that the interferograms presented here are made using Eq. (10), and each disk takes as a particle image in the out-of-focus image plane.

 

Fig. 8 Experimental setup for particle center detection.

Download Full Size | PDF

Figure 9 shows the results produced by the processing algorithm. Figure 9(a) is an example of a simulated 20 interferogram randomly distributed in a field of view of 1024 × 768 pixels, in which all particle disk locations are known a priori. Figure 9(b) is the image of 9(a) in the CCD, which is the amplification of the inverted real image. The center of each fringe pattern in Fig. 9(b) is extracted by gradient-matched algorithm, denoted with red dots, as shown also in Fig. 9(b). Let the preset disk center coordinate in Fig. 9(a) be$(x_{0n},y_{0n})$, where n is the particle serial number (see Fig. 4), namely, the first disk coordinate is$(x_{01},y_{01})$, and the second of $(x_{02},y_{02})$and so on. Take any disk position in Fig. 9(a) as starting position, and the distance between the first disk and any other is given by$d_n=\sqrt{{(x_{0n}-x_{01})}^2+{(y_{0n}-y_{01})}^2}$. Similarly, the corresponding disk center coordinate extracted in Fig. 9(b) is$({x^{\prime }}_n,{y^{\prime }}_n)$, and the center coordinate of the corresponding starting position is$({x^{\prime }}_1,{y^{\prime }}_1)$, and thus we can obtain ${d^{\prime }}_n=\sqrt{{({x^{\prime }}_n-{x^{\prime }}_1)}^2+{({y^{\prime }}_n-{y^{\prime }}_1)}^2}$. For every particle pair, the ratio of the distance, ${d^{\prime }}_n/d_n$, is calculated, and theoretically should be equal. The value of ${d^{\prime }}_n/d_n$ = 19.41 is calculated for the experimental system after calibration. Figure 9(c) shows$d^{\prime }/d$ for particle counts of 10, 20, and 30, and the maximum deviation is 0.28. Figure 9(d) shows the results of the center detection using the three methods for 50 particles. The values$d^{\prime }/d$are 19.54 ± 0.07, 19.36 ± 0.02, and 19.45 ± 0.01 for the erosion matching algorithm, WMF algorithm, and gradient-matched algorithm, respectively. The gradient-matched algorithm clearly has higher center-detection accuracy, and the high errors when using the erosion matching algorithm are because of its sensitivity to noise in binarization. The results show good agreement between experiment and simulation.

 

Fig. 9 (a) A simulated interferogram; (b) the image of (a); (c) the result of the center location for different disk counts; (d) comparison of the location center location with the method in this work and the previously proposed algorithms.

Download Full Size | PDF

Although the experiment here is not identical to the real particle field experiment and particle is not a true particle, we believe that the result may have proved the validity of the location of extraction algorithm.

5.2 Particle size measurement

The detailed experimental configuration for particle size measurement has been described in [15,16]. The particle used in the experiment is GBW(E)120045 polystyrene particle with nominal diameter 45μm, and the refraction index of 1.59. In the experiment, the particles are immersed in deionized water, the scatting angle is set at θ = 45°, the measuring magnification M = 0.174 and z_r = 60.20 mm and the collecting angle α = 4.72° in air.

Figure 10(a) is a captured interference fringe pattern. Using the particle sizing method shown as in Fig. 2, the edge image of each fringe pattern is extracted by a unidirectional gradient-matched algorithm, and the correlation function is calculated to obtain the center of each particle. The result of center detection is shown in Fig. 10(b). The fringe angular spacing/the number of fringes is obtained using Fourier transformation and the improved Rife algorithm for each extracted fringe pattern. Figures 10(c) and 10(d) show, respectively, the fringe pattern of particle 3 in Fig. 10(a) and its 2D frequency spectrum. The spacing between the zero-order spectrum and first-order spectrum of the frequency spectrum is the frequency of the interference fringe, which is obtained by the i-Rife method and the m-Rife method. The estimated frequency f_i = 0.0501 lp/pix with the i-Rife method, as shown in Fig. 10(e), with N = 5.013, d = 41.58 μm, and the absolute error and relative error being 3.42 μm and 7.6%, respectively. The estimated frequency f_m = 0.0493 lp/pix with the m-Rife method is also shown in Fig. 10(e), with N = 4.931, d = 40.89 μm, and the absolute error and relative error being 4.11 μm and 9.13%, respectively. The reason for large discrepancy occurring between the two methods is that the shift direction of the peak is determined by$\left|S(k_1\pm 0.5)\right|$with the i-Rife method and by$\left|S(k_1\pm 1)\right|$with the m-Rife method. Because $\left|S(k_1-0.5)\right|>\left|S(k_1+0.5)\right|$, the signal is shifted in frequency domain towards the left using the i-Rife method, and $\left|S(k_1-1)\right|<\left|S(k_1+1)\right|$,the signal is shifted towards the right using the m-Rife method. Moreover, $\left|S(k_1\pm 0.5)\right|>\left|S(k_1\pm 1)\right|$ (see Fig. 10(e)), the value of $\left|S(k_1\pm 0.5)\right|$is very close to the peak and will not be strongly affected by noise, whereas the value of $\left|S(k_1\pm 1)\right|$is very close to the minima and noise reverses the amplitudes, and thus use of the m-Rife method can cause an error in the shift direction.

 

Fig. 10 Measurement results of the standard particles: (a) the interferometric image; (b) result of the center detection; (c) an interferometric image of particle 3 in (a); (d) 2D frequency spectrum of (c); (e) fine peak detection and frequency estimation; and (f) absolute errors of particle diameter measurement for the i-Rife method and the m-Rife method.

Download Full Size | PDF

Each particle in Fig. 10(a) is analyzed and the result is shown as in Fig. 10(f). It can be seen from Fig. 10(f) that the error of the extracted frequency from the i-Rife algorithm is less than that of from the m-Rife method, and that there are three particles (particle 3, 6 and 8) with large frequency errors using both methods. The maximum deviation between the nominal value and the measured value is for particle 8, and is 1.31 μm for the i-Rife method and −2.84 μm for the m-Rife method. Therefore, the experimental results indicate that the estimated frequency can be obtained with higher accuracy by the improved Rife algorithm.

Finally, we test the algorithm presented using a lot of experimental images, the overall result of the interferogram detection and frequency estimation are similar to the one presented in Fig. 10, and is thus very globally satisfactory, then the particle diameter is calculated using Eq. (1). Figure 11 shows particle size distribution. The particle diameter d = 45.46 ± 0.65 μm, the absolute error is 0.46 μm and the relative error is 1.02%. Note that the result in Fig. 11 is after the geometric correction. In Fig. 1, with the camera at 45°, the laser sheet is not perpendicular to the imaging axis, and thus the defocusing distance and magnification are not uniform across the image, and a correction is required [6]. Additionally, we do not take account of other parameters error (such as z_r, z₀,$\theta$, $\alpha$) and sources of noise present (diffraction effect, the laser sheet thickness, light scattering, particle non-sphericity, refractive index, etc.), the detailed error estimation of IPI has been analyzed in [4] and [24]. The final precision of the particle diameter measurement has been improved in our experiment, and the method presented is very powerful in particle sizing.

 

Fig. 11 Measured diameter distribution of monodisperse polystyrene particle of 45μm.

Download Full Size | PDF

6. Conclusion

In this paper, we propose an algorithm for automatically and simultaneously extracting particle location coordinates and the number of fringes/the fringe frequency based on a unidirectional gradient-matched algorithm and a Fourier transform technique. Sub-pixel accuracy of the extracted frequency is obtained with an improved Rife algorithm. The unidirectional gradient-matched algorithm has higher accuracy of particle center location and is able to tolerate greater particle number density than the WMF and erosion matching algorithms. The m-Rife algorithm is always inferior to the i-Rife algorithm in frequency extraction performance in the case of low SNR, showing that the m-Rife algorithm has large errors because of its unpredictable dependence on noise. The experimental results are in very good agreement with the theoretical predictions. In summary, the method presented in this paper is highly accurate both in particle location detection and particle size measurement, and is promising for application to the spray field, which is undertaken by us now.