- Email: [email protected]
High-accuracy particle sizing by interferometric particle imaging
Optics Communications 312 (2014) 312–318
Contents lists available at ScienceDirect
Optics Communications journal homepage: www.elsevier.com/locate/optcom
High-accuracy particle sizing by interferometric particle imaging Lü Qieni n, Jin Wenhua, Lü Tong, Wang Xiang, Zhang Yimo College of Optoelectronics & Precision Instrument Engineering, Tianjin University, Key Laboratory of Opto-electronics Information Technology, Tianjin University, Ministry of Education, Tianjin 300072, China
art ic l e i nf o
a b s t r a c t
Article history: Received 10 June 2013 Received in revised form 4 September 2013 Accepted 22 September 2013 Available online 6 October 2013
A method of high-accuracy estimation of fringes number/fringes frequency of interferogram based on erosion match and the Fourier transform technique is proposed. The edge images of the interference pattern of particles and the particle mask image are detected respectively by erosion operating firstly and then subtracted with the respective original image, and the center coordinate of particles can be extracted through the 2D correlation operation for the two edge images obtained. The interference pattern of each particle can then be achieved using the center coordinate, the shape and size of the particle image. The number of fringes/fringe spacing of the interferogram of the particle is extracted by Fourier transform and the modified Rife algorithm, and sub-pixel accuracy of the extracted frequency is acquired. Its performance is demonstrated by numerical simulation and experimental measurement. The measurement uncertainty is 7 0.91 μm and the relative error 1.13% for the standard particle of diameter 45 μm. The research results show that the algorithm presented boasts high accuracy for particle sizing as well as location measurement. & 2013 Elsevier B.V. All rights reserved.
Keywords: Interferometric particle imaging (IPI) Erosion match Fourier transform Modified Rife algorithm
1. Introduction Particle sizing plays an important role in many fields of engineering and research [1–12]. Interferometric particle imaging (IPI), based on Mie scattering theory, is a popular method for measuring instantaneous size and spatial distribution of the spherical particle, utilizing the interference fringes generated by reflection and first-order refraction ray scattered by a particle on a defocused image plane. And the particle diameter is proportional to the angular frequency/fringe counts of this interference pattern, and its measurement accuracy depends on the image processing technique for estimating the number of fringes/fringe spacing. Owing to out-of-focus image overlapping, especially, the highdegree overlapping of the fringe pattern in a higher particle number density region, it is very difficult to evaluate the individual particle image of IPI. No doubt, automatically extracting the location and fringe spacing/fringe counts of a particle from interferometric fringe pattern of IPI is both hard but crucial. Several methods have been proposed for IPI image processing so far [7–12,14–16]. For example, Glover et al. [9] first used the Gaussian blur, Canny edge detection and Hough transform to locate the individual droplet in the image field and then detected fringe spacing by least-squares fitting to a Chip function, which was applied to spatially sparse sprays. Niwa et al. [10] detected the
n
Corresponding author. Tel.: þ 86 22 27401903. E-mail address: [email protected] (L. Qieni).
0030-4018/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.optcom.2013.09.049
position of bubble by using the template matching method with circular pattern in low number density of bubbles, and the fringe frequency of sub-pixel accuracy was achieved by Fourier transform and Gaussian fitting. The method includes 2D convolution and FFT was proposed to extract the special coordinate and fringe angle frequency of fringe pattern from interferometric particle image in Ref [14]. Quérel group [11] developed a global algorithm based on 1D Fourier transform to calculate the real-time droplet size distribution from IPI interferogram, and such method demands high speed to perform real time processing during a flight through clouds. Lacagnina et al. [12] carried out firstly the cross correlation between a Fourier spectrogram of interference fringe of a particle in the wavenumber domain and a reference image, and the area around the cross-correlation local maxima as the center was selected by a square window, and then the FFT and auto correlation for this area was performed to evaluate the fringe wave number and then the bubble diameter. Recently, we [16] developed a method of extracting the number of fringes/fringe spacing of circular interferogram based on wavelet matched filter and the Fourier transform technique, and sub-pixel accuracy of the extracted frequency was attained with the modified Rife algorithm and we are currently studying a new algorithm of extracting fringe counts/fringe spacing of interference fringe pattern. The present algorithm is based on erosion match and the Fourier transform technique, and frequency is subdivided by the modified Rife method. The paper recommends this algorithm for estimating the fringe spacing, whose performance is demonstrated through simulation and experimental measurements. Comparing this new
L. Qieni et al. / Optics Communications 312 (2014) 312–318
313
Fig. 1. Schematic of interferometric imaging of particle scattering light.
method with that of the wavelet matched filter, for the same algorithm recognition rate and measurement accuracy, the measurable particle number density doubles, then higher accurate position measurement can be realized, promising to apply to the high-density spray field.
2. Interferomemtic particle imaging A diagram of interferometric particle imaging is shown in Fig. 1. When a transparent spherical particle is illuminated by a laser sheet, the reflected and first-order refracted ray scattered from the particle form two glare points on the focus image plane, one glare point representing the reflected ray and the other the refracted ray. And the interference fringe pattern is produced in an out-offocus image plane. The shape of the defocused particle image relies on that of the receiving aperture, and size di of the defocused particle image is given by [7] 1 1 ð1Þ di ¼ da 1 zr f z0 where da is the diameter of the aperture, f is the focal length of the imaging lens, zo is the object distance, zr is the distance between the defocused image plane and the lens, as shown in Fig. 1. The size of the particle fringe pattern on the out-of-focus plane is independent of particle size. The diameter of particle is related to the fringe spacing of interferogram. Using geometrical optics, particle diameter d can be expressed as the following [13,17]: !1 2λN θ m sin ðθ=2Þ d¼ cos þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð2Þ α 2 m2 2m cos ðθ=2Þ þ 1 Where λ is the wavelength of the laser light source, m is the refractive index of the particle, θ and α are the scattering angle and collecting angle, respectively, N ¼α/Δφ is the number of fringes and Δφ the angular spacing of the interference fringe. From Eq. (2), for a given experimental system, that is, θ, m, λ and α given, particle diameter d is proportional to the number of fringes N, and can be calculated by measuring N or Δφ. Consequently, the accuracy of particle sizing using IPI depends upon the interferogram image analysis, and the automatic extraction of fringe spacing or fringe counts is one of the key aspects of the technique.
3. Algorithm Fig. 2 shows the algorithm flowchart of automatically processing fringe image of IPI we present. The interference pattern of particles and the particle mask image (to construct a particle mask image according to Eq. (1)) are processed by eroding and subtracting, respectively; then, the edge images are obtained, and the position coordinate (x,y) of particle images can be extracted through
I (x, y)
P (x, y)
interferogram
particle mask image
erosion edge
I (x, y)
P (x, y)
erosion edge
particle (x,y) coordinate each particle interferogram
FT
the modified Rife algorithm
number of fringes/fringe spacing particle diameter and (x,y) coordinate Fig. 2. Flowchart of particle measurement.
the 2D correlation operation for the two edge images obtained. The interference fringe image of each particle is extracted from raw images according to the center coordinate (x,y), the shape and size of the particle fringe image; then Fourier transform is performed to evaluate fringe angle spacing or the number of fringes, and subpixel accuracy of the extracted frequency is acquired by the modified Rife algorithm and then particle diameter.
3.1. Detection of particle center The accurate identification of the location of particle center is crucial in interference fringe image processing of IPI. At first, the interference pattern of particles and the particle mask image are eroded using a disc with one pixel radius, the edge images are then extracted by subtracting the erosion image from the respective original image, respectively; the cross-correlation function between the two edge images obtained is calculated. Suppose that the interference pattern of particles I(x,y) is the input image, and the particle mask image P(x,y) the target image to be identified, after eroded and edge detected, are denoted by I′(x,y) and P′(x,y), respectively; the 2D correlation operation for the two edge images commanded is computed I′ðx; yÞ P′ðx; yÞ ¼ ∬ I ′n ðx′; y′ÞP′ðx′ þ x; y′þ yÞdx′dy′
ð3Þ
where denotes correlation operation, and the peak values of the cross-correlation function of Eq. (3) is the center position of particle image. We designate this algorithm as erosion match algorithm or erosion correlated algorithm.
314
L. Qieni et al. / Optics Communications 312 (2014) 312–318
3.2. Estimation of fringe spacing
amplitude of power spectrum is A2. Introduce a parameter δ,
After detecting particle location by the erosion match algorithm, the interference pattern of each particle can be achieved using the position coordinate, the shape and size of the particle fringe image, and then Fourier transform is performed. As discrete Fourier power spectrum, the accuracy of frequency estimating is limited for the resolution of FFT, the modified Rife algorithm [18] is then proposed, which is very much like the interpolation iterative algorithm in Ref. [19], and sub-pixel accuracy of the extracted frequency is acquired. Let the discrete Fourier transform of signal s(n) be S(k), which can be shown to be
δ¼
N1 2πnk SðkÞ ¼ ∑ sðnÞexp i N n¼0
ð4Þ
where k ¼ 0; 1; 2; ⋯; N 1. Suppose that k1 is the discrete frequency with the maximum power of the discrete spectrum S(k), and the corresponding amplitude is A1, the discrete frequency with the secondary maximum power is k2 ¼k1 71, and the corresponding
A2 A1 þ A2
ð5Þ
If δ satisfies 1 2 o jδjo 3 3
ð6Þ
then, the frequency estimation f f¼
k1 þrδ NΔx
ð7Þ
where Δx is the sampling interval, r is a coefficient. When k2 ¼ k1 þ 1, r ¼ 1; when k2 ¼k1 1, r¼ 1. If Eq. (6) is not satisfied, then s(n) will undergo frequency shift, and 1 Δδ ¼ jδj 2
ð8Þ
2πn rΔδ s1 ðnÞ ¼ sðnÞexp i N
ð9Þ
For Eq. (9), the frequency can be estimated by completing Eqs. (4)–(7).
Fig. 3. Simulation results of the interferogram center detecting: (a) interferometric image; (b) the mask image; (c) the extracted edge image of the mask image; (d) the extracted edge image of (a); (e) distribution of cross-correlation value and (f) the results of the center detecting. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
L. Qieni et al. / Optics Communications 312 (2014) 312–318
4. Simulation and analysis 4.1. Detection of interference fringe pattern center Computer simulation has been first carried out to evaluate the performance of the algorithm we present. Fig. 3 shows the results of simulation: the field of view consists of 2048 2048 pixels, and the out-of-focus image size of a particle is 100 pixels. Fig. 3(a) is
315
the simulated interferograms with 50 spherical particles with random placement. Fig. 3(b) is the particle mask image with the diameter of 100 pixels. Fig. 3(c) and (d) shows the edge images of Fig. 3(a) and (b) obtained by eroding and subtracting, respectively, and the cross-correlation function between Fig. 3(c) and (d) is calculated and the result is displayed in Fig. 3(e). Fig. 3(f) is the results of detected particle center denoted with the red spot, and the detected location coordinates absolutely agree with the preset
Fig. 4. The results of the center detecting for different number of particles: (a) 100; (b) 400; (c) 700 and (d) 1000.
Fig. 5. Particle position detected versus the preset values: (a) X coordinate and (b) Y coordinate.
316
L. Qieni et al. / Optics Communications 312 (2014) 312–318
values. Comparing Fig. 3(a) with Fig. 3(f), it can be seen that the position of fringe pattern of particles in Fig. 3(a) is completely and accurately identical in Fig. 3(f).
1.0
Recognition ratio Rp
0.8
4.2. Recognition ratio analysis
0.6 0.4 0.2 0.0 0.0
0.2
0.4
0.6
0.8
Overlap coefficient
The reliability of the extracted position method highly matters for particle sizing from the fringe pattern of IPI. For a given IPI system, the interferometric images are overlapped, which depends on particle number density. The higher the particle number density, the higher probability of particle image overlap, the recognition is being affected. Let two parameters recognition rate Rp and overlap coefficient γ be defined. Recognition rate Rp is defined as the ratio of the number of extracted particles to the total number of the particles as follows [20]: Rp ¼
Fig. 6. Variation of the recognition ratio Rp with an overlap coefficient γ.
Recognition ratio Rp
Fig. 4 shows the results of located particle center with the erosion match algorithm for the number of particles of 100, 400, 700 and 1000, respectively, the correlation value between the detected coordinate (x,y) and the preset value (x0,y0) is shown in Fig. 5, where x/x0 ¼ 1, 1, 0.99998, 0.99999 and 1.00007, y/y0 ¼ 1, 1, 0.99996, 1.00004 and 1.00002, respectively. As shown in Figs. 4 and 5, the high accuracy of particle position coordinate(x,y) can be possessed with algorithm presented in this paper.
The number of extracted particles The number of seeded particles
ð10Þ
1.0
An overlap coefficient γ is defined as the ratio of the area of particle image overlap to the overall area of the particle images [21], and can be expressed as
0.8
γ¼
0.6
0.4
0.2
0.0
2D erosion 2D wavelet 100
1000 Number of particles
Fig. 7. Recognition ratio Rp comparison of two algorithms.
10000
The area of overlap The total area
ð11Þ
Indeed, recognition ratio Rp depends on the number of extracted particles and the number of extracted particles is related to overlap coefficient γ. Fig. 6 illustrates the dependence of Rp on overlap coefficient γ. The figure demonstrates distinctly that the number of extracted particles decreases with increased overlap coefficient γ. For Rp 40.8, γo0.2345, and Rp 40.9, γo0.1618. Fig. 7 shows the estimated recognition ratio Rp versus particle count with the erosion match algorithm and the wavelet-matched filter algorithm and those fitting curves obtained by using the Boltzmann model. It can be seen from Fig. 7 that the recognition rate Rp based on the erosionmatched algorithm is higher than that of the wavelet-matched filter algorithm; for Rp 40.9, the number of extracted particles is 532 and 250 for the erosion-matched algorithm and the wavelet-matched filter algorithm, respectively.
Fig. 8. Simulation results of the frequency estimation: (a) the interferometric image; (b) 2D frequency spectrum of Fig.(a); (c)fine peak detection and frequency estimation.
L. Qieni et al. / Optics Communications 312 (2014) 312–318
317
4.3. Frequency extraction of particle fringe pattern
5. Experimental results
After extracting fringe pattern of each particle and calculating 1D Fourier transform for each fringe pattern, then the fringe counts/fringe spacing of the interferogram of a particle is calculated with the modified Rife algorithm. For a laser interference pattern, the frequency shift will occur only along X-axis to evaluate frequency with modified Rife algorithm, which can be written as
To validate the method presented of particle size measurement experimentally, the measurement system setup is assembled. A 532 nm CW semiconductor laser with the maximum power of 1.5 W is used as a light source, and a laser beam, expanded, spatially filtered, collimated, compressed by a pair of cylindrical lenses, becomes a vertical sheet with a thickness of 1.25 mm. The CCD is a 10-bit digital CCD sensor (DS-21-04M15camera) with 2048 2048 pixels and pixel size of 7.4 μm 7.4 μm. The imaging lens is a AF50 mm f/1.8D Nikon lens. The particle used in the experiment is GBW(E)120045 polystyrene particle with diameter 45 μm, and the refraction index of 1.59. In experiment, the particles are immersed in deionized water, the scattering angle is set at θ¼451, the measuring magnification M¼0.174, zr ¼60.20 mm and the collecting angle α¼4.721 in air. Fig. 9(a) is a captured interference fringe pattern, the number of which in each interferogram is the same, about 5 fringes. Using particle sizing method shown as in Fig. 2, the edge image of each fringe pattern is acquired by eroding and subtracting and the correlation function is calculated and then the center of each particle is achieved, the results of correlation value and location determination are exhibited as in Fig. 9(b) and (c), respectively. The fringe angular spacing/the number of fringes is obtained by means of Fourier transform and the modified Rife algorithm for each fringe extracted. Fig. 9(d) and (e) represent a fringe pattern and its frequency spectrum of a particle in Fig. 9(a), the 1D Fourier spectrum is shown as in Fig. 9(f) and (g) is 1 order discrete frequency spectrum of Fig. 9(f). The distance between the zero-order spectrum and first-order spectrum of the frequency spectrum is the oscillation frequency of interference fringe, the estimated frequency f¼0.064lp/pix with the m-Rife method, the number of fringe N¼5.16, and the particle diameter d¼ 43.28 μm. Each particle in Fig. 9(a) is analyzed, the result of which is that the particle diameter d¼44.4970.91 μm, the absolute error is 0.51 μm and the relative error of 1.13%.
2πx I 1 ðx; yÞ ¼ Iðx; yÞexp i rΔδ N
ð12Þ
Fig. 8(a) displays a simulated circular interferogram of a particle with the diameter of 201 pixels and the frequency of 0.0501lp/pix and Δx¼ 1. The signal is corrupted by Gaussian noise with variance of 0.15. The 2D Fourier spectrum is shown in Fig. 8(b), the discrete frequency with the maximum power and the secondary maximum of the discrete spectrum k1 ¼101 91¼10, k2 ¼101 90¼ 11(see Fig. 8(c)). Using Eq. (5), we get δ¼0.2027 and this value does not satisfy Eq. (6); therefore, the frequency shift is needed and completed, currently δ¼0.3862, which meets the Eq. (6), the estimated frequency fM ¼0.0502lp/pix using Eq. (7), as shown in Fig. 8(c), the absolute error and relative error are 0.0001 and 0.1996%, respectively. The extracted frequency fG ¼0.0503lp/pix with Gaussian fitting interpolation, as also shown in Fig. 8(c), the absolute error and relative error are 0.0002 and 0.3992%, respectively. Compared with Gaussian fitting interpolation, the m-Rife interpolation has high accuracy in frequency estimation. At the same time we have also investigated the RMSE(root mean square error) and the MAE (mean absolute error) of the estimated frequency by employing Gaussian, m-Gaussian, Rife, m-Rife and cubic B-spline interpolation. The simulation results show that the RMSE of these interpolation methods are almost the same; however, the MAE of m-Rife method holds superiority to other several algorithms in a wide range [16].
Fig. 9. The measurement results of the standard particles: (a) the interferometric image; (b)the cross-correlation value; (c)the result of the center detecting; (d)–(g) an interferometric image of a particle in (a), 2D and 1D frequency spectrum of (d) and 1 order discrete spectrum of (f).
318
L. Qieni et al. / Optics Communications 312 (2014) 312–318
6. Conclusions
References
The automatic fringe spacing/the number of fringes extracting assumes critical significance for the accuracy of particle size measurement in IPI technique. We present an algorithm of estimating the number of fringes/fringes frequency based on erosion match and the Fourier transform technique, and sub-pixel accuracy of the extracted frequency is obtained by the modified Rife algorithm, and simulation and experiment are completed. The simulation results are that the recognition ratio is higher than 0.9 and the extracted frequency error is less than 0.0185% for the overlap coefficient γo0.1618, by which IPI image processing is effectively addressed. The experiments are conducted for standard particles of diameter 45 μm. The measurement results are that the particle diameter is 44.4970.91 μm, the relative error 1.13%, respectively. What matters is that this method can simultaneously detect particle center location with higher accuracy and the velocity can then be obtained with the particle tracking velocimetry technique. The research results demonstrate the feasibility of the method presented in this paper, which applies to interferometric out-of-focus imaging of bubbles and droplet particle, such as spray filed.
[1] Lin Jia, E.L Thomas, Physical Review B 84 (2011) 125128. [2] Lin Jia, E.L. Thomas, Journal of the Optical Society of America B 26 (2009) 1882. [3] H. Ahn, P. Thiyagarajan, L. Jia, S. Kim, J. Yoon, E.L. Thomas, J. Jang, Nanoscale 5 (2013) 1836. [4] B.W.S. Kim, L. Jia, E.L. Thomas, Adv. Mater. 21 (2009) 1921. [5] C. Lacour, D. Durox, S. Ducruix, M. Massot, Experiments in Fluids 51 (2011) 295. [6] C.F Hess, D.L' Esperance, Experiments in Fluids 47 (2009) 171. [7] N. Damaschke, H. Nobach, T.I Nonn, N. Semidetnov, C. Tropea, Experiments in Fluids 39 (2005) 336. [8] M. Maeda, Y. Akasaka, T. Kawaguchi, Experiments in Fluids 33 (2002) 125. [9] A.R Glover, S.M. Skippon, R.D. boyle, Applied Optics 34 (1995) 8409. [10] Y. Niwa, Y. Kamiya, T. Kawaguchi, M. Maeda, Proceedings of 10th International Symposium on Application of laser techniques to fluid mechanics, Lisbon, vol. 1, 2000, 〈http://in3.dem.ist.utl.pt/downloads/lxlaser2000/pdf/38_1.pdf〉. [11] A. Quérel, P. Lemaitre, M. Brunel, E. Porcheron, G. Gréhan, Measurement Science and Technology 21 (2010). (015306-1). [12] G. Lacagnina, S. Grizzi, M. Falchi, F. Di Felice, G.P Romano, Experiments in Fluids 50 (2011) 1153. [13] K.H Hesselbacher, K. Anders, A. Frohn, Applied Optics 30 (1991) 4930. [14] Pu Shiliang, Pu Xiangguo, Yuan Zhenfu, Cen Kefa, L. Mees, D. Lebrun, Proceedings of the Chinese Society for Electrical Engineering 24 (2005) 201. (in Chinese). [15] K. Jorabchi, R.G. Brennan, J.A. Levine, A. Montaser, Journal of Analytical Atomic Spectrometry 21 (2006) 839. [16] Lü Qieni, Ge Baozhen, Chen Yiliang, Zhang Yimo, Acta Optica Sinica 31 (2011) 0412009-1 (in Chinese). [17] N. Semidetnov, C. Tropea, Measurement Science and Technology 15 (2004) 112. [18] Deng Zhenmiao, Liu Yu, Wang Zhizhong, Journal of Data Acquisition and Processing 21 (2006) 473. (in Chinese). [19] J.B. Tsui, Digital Techniques for Wideband Receivers, 2nd ed., SciTech Publishing Inc., Raleigh, NC (2004) 119–126. [20] S. Kim, S.J. Lee, Experiments in Fluids 44 (2008) 623. [21] N. Damaschke, H. Nobach, C. Tropea, Experiments in Fluids 32 (2002) 143.
Acknowledgments The work is supported by National Natural Science Foundation of China (No. 61275019) and 2011 Project of State Key laboratory of Engines (No. K2011-10).
Contents lists available at ScienceDirect
Optics Communications journal homepage: www.elsevier.com/locate/optcom
High-accuracy particle sizing by interferometric particle imaging Lü Qieni n, Jin Wenhua, Lü Tong, Wang Xiang, Zhang Yimo College of Optoelectronics & Precision Instrument Engineering, Tianjin University, Key Laboratory of Opto-electronics Information Technology, Tianjin University, Ministry of Education, Tianjin 300072, China
art ic l e i nf o
a b s t r a c t
Article history: Received 10 June 2013 Received in revised form 4 September 2013 Accepted 22 September 2013 Available online 6 October 2013
A method of high-accuracy estimation of fringes number/fringes frequency of interferogram based on erosion match and the Fourier transform technique is proposed. The edge images of the interference pattern of particles and the particle mask image are detected respectively by erosion operating firstly and then subtracted with the respective original image, and the center coordinate of particles can be extracted through the 2D correlation operation for the two edge images obtained. The interference pattern of each particle can then be achieved using the center coordinate, the shape and size of the particle image. The number of fringes/fringe spacing of the interferogram of the particle is extracted by Fourier transform and the modified Rife algorithm, and sub-pixel accuracy of the extracted frequency is acquired. Its performance is demonstrated by numerical simulation and experimental measurement. The measurement uncertainty is 7 0.91 μm and the relative error 1.13% for the standard particle of diameter 45 μm. The research results show that the algorithm presented boasts high accuracy for particle sizing as well as location measurement. & 2013 Elsevier B.V. All rights reserved.
Keywords: Interferometric particle imaging (IPI) Erosion match Fourier transform Modified Rife algorithm
1. Introduction Particle sizing plays an important role in many fields of engineering and research [1–12]. Interferometric particle imaging (IPI), based on Mie scattering theory, is a popular method for measuring instantaneous size and spatial distribution of the spherical particle, utilizing the interference fringes generated by reflection and first-order refraction ray scattered by a particle on a defocused image plane. And the particle diameter is proportional to the angular frequency/fringe counts of this interference pattern, and its measurement accuracy depends on the image processing technique for estimating the number of fringes/fringe spacing. Owing to out-of-focus image overlapping, especially, the highdegree overlapping of the fringe pattern in a higher particle number density region, it is very difficult to evaluate the individual particle image of IPI. No doubt, automatically extracting the location and fringe spacing/fringe counts of a particle from interferometric fringe pattern of IPI is both hard but crucial. Several methods have been proposed for IPI image processing so far [7–12,14–16]. For example, Glover et al. [9] first used the Gaussian blur, Canny edge detection and Hough transform to locate the individual droplet in the image field and then detected fringe spacing by least-squares fitting to a Chip function, which was applied to spatially sparse sprays. Niwa et al. [10] detected the
n
Corresponding author. Tel.: þ 86 22 27401903. E-mail address: [email protected] (L. Qieni).
0030-4018/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.optcom.2013.09.049
position of bubble by using the template matching method with circular pattern in low number density of bubbles, and the fringe frequency of sub-pixel accuracy was achieved by Fourier transform and Gaussian fitting. The method includes 2D convolution and FFT was proposed to extract the special coordinate and fringe angle frequency of fringe pattern from interferometric particle image in Ref [14]. Quérel group [11] developed a global algorithm based on 1D Fourier transform to calculate the real-time droplet size distribution from IPI interferogram, and such method demands high speed to perform real time processing during a flight through clouds. Lacagnina et al. [12] carried out firstly the cross correlation between a Fourier spectrogram of interference fringe of a particle in the wavenumber domain and a reference image, and the area around the cross-correlation local maxima as the center was selected by a square window, and then the FFT and auto correlation for this area was performed to evaluate the fringe wave number and then the bubble diameter. Recently, we [16] developed a method of extracting the number of fringes/fringe spacing of circular interferogram based on wavelet matched filter and the Fourier transform technique, and sub-pixel accuracy of the extracted frequency was attained with the modified Rife algorithm and we are currently studying a new algorithm of extracting fringe counts/fringe spacing of interference fringe pattern. The present algorithm is based on erosion match and the Fourier transform technique, and frequency is subdivided by the modified Rife method. The paper recommends this algorithm for estimating the fringe spacing, whose performance is demonstrated through simulation and experimental measurements. Comparing this new
L. Qieni et al. / Optics Communications 312 (2014) 312–318
313
Fig. 1. Schematic of interferometric imaging of particle scattering light.
method with that of the wavelet matched filter, for the same algorithm recognition rate and measurement accuracy, the measurable particle number density doubles, then higher accurate position measurement can be realized, promising to apply to the high-density spray field.
2. Interferomemtic particle imaging A diagram of interferometric particle imaging is shown in Fig. 1. When a transparent spherical particle is illuminated by a laser sheet, the reflected and first-order refracted ray scattered from the particle form two glare points on the focus image plane, one glare point representing the reflected ray and the other the refracted ray. And the interference fringe pattern is produced in an out-offocus image plane. The shape of the defocused particle image relies on that of the receiving aperture, and size di of the defocused particle image is given by [7] 1 1 ð1Þ di ¼ da 1 zr f z0 where da is the diameter of the aperture, f is the focal length of the imaging lens, zo is the object distance, zr is the distance between the defocused image plane and the lens, as shown in Fig. 1. The size of the particle fringe pattern on the out-of-focus plane is independent of particle size. The diameter of particle is related to the fringe spacing of interferogram. Using geometrical optics, particle diameter d can be expressed as the following [13,17]: !1 2λN θ m sin ðθ=2Þ d¼ cos þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð2Þ α 2 m2 2m cos ðθ=2Þ þ 1 Where λ is the wavelength of the laser light source, m is the refractive index of the particle, θ and α are the scattering angle and collecting angle, respectively, N ¼α/Δφ is the number of fringes and Δφ the angular spacing of the interference fringe. From Eq. (2), for a given experimental system, that is, θ, m, λ and α given, particle diameter d is proportional to the number of fringes N, and can be calculated by measuring N or Δφ. Consequently, the accuracy of particle sizing using IPI depends upon the interferogram image analysis, and the automatic extraction of fringe spacing or fringe counts is one of the key aspects of the technique.
3. Algorithm Fig. 2 shows the algorithm flowchart of automatically processing fringe image of IPI we present. The interference pattern of particles and the particle mask image (to construct a particle mask image according to Eq. (1)) are processed by eroding and subtracting, respectively; then, the edge images are obtained, and the position coordinate (x,y) of particle images can be extracted through
I (x, y)
P (x, y)
interferogram
particle mask image
erosion edge
I (x, y)
P (x, y)
erosion edge
particle (x,y) coordinate each particle interferogram
FT
the modified Rife algorithm
number of fringes/fringe spacing particle diameter and (x,y) coordinate Fig. 2. Flowchart of particle measurement.
the 2D correlation operation for the two edge images obtained. The interference fringe image of each particle is extracted from raw images according to the center coordinate (x,y), the shape and size of the particle fringe image; then Fourier transform is performed to evaluate fringe angle spacing or the number of fringes, and subpixel accuracy of the extracted frequency is acquired by the modified Rife algorithm and then particle diameter.
3.1. Detection of particle center The accurate identification of the location of particle center is crucial in interference fringe image processing of IPI. At first, the interference pattern of particles and the particle mask image are eroded using a disc with one pixel radius, the edge images are then extracted by subtracting the erosion image from the respective original image, respectively; the cross-correlation function between the two edge images obtained is calculated. Suppose that the interference pattern of particles I(x,y) is the input image, and the particle mask image P(x,y) the target image to be identified, after eroded and edge detected, are denoted by I′(x,y) and P′(x,y), respectively; the 2D correlation operation for the two edge images commanded is computed I′ðx; yÞ P′ðx; yÞ ¼ ∬ I ′n ðx′; y′ÞP′ðx′ þ x; y′þ yÞdx′dy′
ð3Þ
where denotes correlation operation, and the peak values of the cross-correlation function of Eq. (3) is the center position of particle image. We designate this algorithm as erosion match algorithm or erosion correlated algorithm.
314
L. Qieni et al. / Optics Communications 312 (2014) 312–318
3.2. Estimation of fringe spacing
amplitude of power spectrum is A2. Introduce a parameter δ,
After detecting particle location by the erosion match algorithm, the interference pattern of each particle can be achieved using the position coordinate, the shape and size of the particle fringe image, and then Fourier transform is performed. As discrete Fourier power spectrum, the accuracy of frequency estimating is limited for the resolution of FFT, the modified Rife algorithm [18] is then proposed, which is very much like the interpolation iterative algorithm in Ref. [19], and sub-pixel accuracy of the extracted frequency is acquired. Let the discrete Fourier transform of signal s(n) be S(k), which can be shown to be
δ¼
N1 2πnk SðkÞ ¼ ∑ sðnÞexp i N n¼0
ð4Þ
where k ¼ 0; 1; 2; ⋯; N 1. Suppose that k1 is the discrete frequency with the maximum power of the discrete spectrum S(k), and the corresponding amplitude is A1, the discrete frequency with the secondary maximum power is k2 ¼k1 71, and the corresponding
A2 A1 þ A2
ð5Þ
If δ satisfies 1 2 o jδjo 3 3
ð6Þ
then, the frequency estimation f f¼
k1 þrδ NΔx
ð7Þ
where Δx is the sampling interval, r is a coefficient. When k2 ¼ k1 þ 1, r ¼ 1; when k2 ¼k1 1, r¼ 1. If Eq. (6) is not satisfied, then s(n) will undergo frequency shift, and 1 Δδ ¼ jδj 2
ð8Þ
2πn rΔδ s1 ðnÞ ¼ sðnÞexp i N
ð9Þ
For Eq. (9), the frequency can be estimated by completing Eqs. (4)–(7).
Fig. 3. Simulation results of the interferogram center detecting: (a) interferometric image; (b) the mask image; (c) the extracted edge image of the mask image; (d) the extracted edge image of (a); (e) distribution of cross-correlation value and (f) the results of the center detecting. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
L. Qieni et al. / Optics Communications 312 (2014) 312–318
4. Simulation and analysis 4.1. Detection of interference fringe pattern center Computer simulation has been first carried out to evaluate the performance of the algorithm we present. Fig. 3 shows the results of simulation: the field of view consists of 2048 2048 pixels, and the out-of-focus image size of a particle is 100 pixels. Fig. 3(a) is
315
the simulated interferograms with 50 spherical particles with random placement. Fig. 3(b) is the particle mask image with the diameter of 100 pixels. Fig. 3(c) and (d) shows the edge images of Fig. 3(a) and (b) obtained by eroding and subtracting, respectively, and the cross-correlation function between Fig. 3(c) and (d) is calculated and the result is displayed in Fig. 3(e). Fig. 3(f) is the results of detected particle center denoted with the red spot, and the detected location coordinates absolutely agree with the preset
Fig. 4. The results of the center detecting for different number of particles: (a) 100; (b) 400; (c) 700 and (d) 1000.
Fig. 5. Particle position detected versus the preset values: (a) X coordinate and (b) Y coordinate.
316
L. Qieni et al. / Optics Communications 312 (2014) 312–318
values. Comparing Fig. 3(a) with Fig. 3(f), it can be seen that the position of fringe pattern of particles in Fig. 3(a) is completely and accurately identical in Fig. 3(f).
1.0
Recognition ratio Rp
0.8
4.2. Recognition ratio analysis
0.6 0.4 0.2 0.0 0.0
0.2
0.4
0.6
0.8
Overlap coefficient
The reliability of the extracted position method highly matters for particle sizing from the fringe pattern of IPI. For a given IPI system, the interferometric images are overlapped, which depends on particle number density. The higher the particle number density, the higher probability of particle image overlap, the recognition is being affected. Let two parameters recognition rate Rp and overlap coefficient γ be defined. Recognition rate Rp is defined as the ratio of the number of extracted particles to the total number of the particles as follows [20]: Rp ¼
Fig. 6. Variation of the recognition ratio Rp with an overlap coefficient γ.
Recognition ratio Rp
Fig. 4 shows the results of located particle center with the erosion match algorithm for the number of particles of 100, 400, 700 and 1000, respectively, the correlation value between the detected coordinate (x,y) and the preset value (x0,y0) is shown in Fig. 5, where x/x0 ¼ 1, 1, 0.99998, 0.99999 and 1.00007, y/y0 ¼ 1, 1, 0.99996, 1.00004 and 1.00002, respectively. As shown in Figs. 4 and 5, the high accuracy of particle position coordinate(x,y) can be possessed with algorithm presented in this paper.
The number of extracted particles The number of seeded particles
ð10Þ
1.0
An overlap coefficient γ is defined as the ratio of the area of particle image overlap to the overall area of the particle images [21], and can be expressed as
0.8
γ¼
0.6
0.4
0.2
0.0
2D erosion 2D wavelet 100
1000 Number of particles
Fig. 7. Recognition ratio Rp comparison of two algorithms.
10000
The area of overlap The total area
ð11Þ
Indeed, recognition ratio Rp depends on the number of extracted particles and the number of extracted particles is related to overlap coefficient γ. Fig. 6 illustrates the dependence of Rp on overlap coefficient γ. The figure demonstrates distinctly that the number of extracted particles decreases with increased overlap coefficient γ. For Rp 40.8, γo0.2345, and Rp 40.9, γo0.1618. Fig. 7 shows the estimated recognition ratio Rp versus particle count with the erosion match algorithm and the wavelet-matched filter algorithm and those fitting curves obtained by using the Boltzmann model. It can be seen from Fig. 7 that the recognition rate Rp based on the erosionmatched algorithm is higher than that of the wavelet-matched filter algorithm; for Rp 40.9, the number of extracted particles is 532 and 250 for the erosion-matched algorithm and the wavelet-matched filter algorithm, respectively.
Fig. 8. Simulation results of the frequency estimation: (a) the interferometric image; (b) 2D frequency spectrum of Fig.(a); (c)fine peak detection and frequency estimation.
L. Qieni et al. / Optics Communications 312 (2014) 312–318
317
4.3. Frequency extraction of particle fringe pattern
5. Experimental results
After extracting fringe pattern of each particle and calculating 1D Fourier transform for each fringe pattern, then the fringe counts/fringe spacing of the interferogram of a particle is calculated with the modified Rife algorithm. For a laser interference pattern, the frequency shift will occur only along X-axis to evaluate frequency with modified Rife algorithm, which can be written as
To validate the method presented of particle size measurement experimentally, the measurement system setup is assembled. A 532 nm CW semiconductor laser with the maximum power of 1.5 W is used as a light source, and a laser beam, expanded, spatially filtered, collimated, compressed by a pair of cylindrical lenses, becomes a vertical sheet with a thickness of 1.25 mm. The CCD is a 10-bit digital CCD sensor (DS-21-04M15camera) with 2048 2048 pixels and pixel size of 7.4 μm 7.4 μm. The imaging lens is a AF50 mm f/1.8D Nikon lens. The particle used in the experiment is GBW(E)120045 polystyrene particle with diameter 45 μm, and the refraction index of 1.59. In experiment, the particles are immersed in deionized water, the scattering angle is set at θ¼451, the measuring magnification M¼0.174, zr ¼60.20 mm and the collecting angle α¼4.721 in air. Fig. 9(a) is a captured interference fringe pattern, the number of which in each interferogram is the same, about 5 fringes. Using particle sizing method shown as in Fig. 2, the edge image of each fringe pattern is acquired by eroding and subtracting and the correlation function is calculated and then the center of each particle is achieved, the results of correlation value and location determination are exhibited as in Fig. 9(b) and (c), respectively. The fringe angular spacing/the number of fringes is obtained by means of Fourier transform and the modified Rife algorithm for each fringe extracted. Fig. 9(d) and (e) represent a fringe pattern and its frequency spectrum of a particle in Fig. 9(a), the 1D Fourier spectrum is shown as in Fig. 9(f) and (g) is 1 order discrete frequency spectrum of Fig. 9(f). The distance between the zero-order spectrum and first-order spectrum of the frequency spectrum is the oscillation frequency of interference fringe, the estimated frequency f¼0.064lp/pix with the m-Rife method, the number of fringe N¼5.16, and the particle diameter d¼ 43.28 μm. Each particle in Fig. 9(a) is analyzed, the result of which is that the particle diameter d¼44.4970.91 μm, the absolute error is 0.51 μm and the relative error of 1.13%.
2πx I 1 ðx; yÞ ¼ Iðx; yÞexp i rΔδ N
ð12Þ
Fig. 8(a) displays a simulated circular interferogram of a particle with the diameter of 201 pixels and the frequency of 0.0501lp/pix and Δx¼ 1. The signal is corrupted by Gaussian noise with variance of 0.15. The 2D Fourier spectrum is shown in Fig. 8(b), the discrete frequency with the maximum power and the secondary maximum of the discrete spectrum k1 ¼101 91¼10, k2 ¼101 90¼ 11(see Fig. 8(c)). Using Eq. (5), we get δ¼0.2027 and this value does not satisfy Eq. (6); therefore, the frequency shift is needed and completed, currently δ¼0.3862, which meets the Eq. (6), the estimated frequency fM ¼0.0502lp/pix using Eq. (7), as shown in Fig. 8(c), the absolute error and relative error are 0.0001 and 0.1996%, respectively. The extracted frequency fG ¼0.0503lp/pix with Gaussian fitting interpolation, as also shown in Fig. 8(c), the absolute error and relative error are 0.0002 and 0.3992%, respectively. Compared with Gaussian fitting interpolation, the m-Rife interpolation has high accuracy in frequency estimation. At the same time we have also investigated the RMSE(root mean square error) and the MAE (mean absolute error) of the estimated frequency by employing Gaussian, m-Gaussian, Rife, m-Rife and cubic B-spline interpolation. The simulation results show that the RMSE of these interpolation methods are almost the same; however, the MAE of m-Rife method holds superiority to other several algorithms in a wide range [16].
Fig. 9. The measurement results of the standard particles: (a) the interferometric image; (b)the cross-correlation value; (c)the result of the center detecting; (d)–(g) an interferometric image of a particle in (a), 2D and 1D frequency spectrum of (d) and 1 order discrete spectrum of (f).
318
L. Qieni et al. / Optics Communications 312 (2014) 312–318
6. Conclusions
References
The automatic fringe spacing/the number of fringes extracting assumes critical significance for the accuracy of particle size measurement in IPI technique. We present an algorithm of estimating the number of fringes/fringes frequency based on erosion match and the Fourier transform technique, and sub-pixel accuracy of the extracted frequency is obtained by the modified Rife algorithm, and simulation and experiment are completed. The simulation results are that the recognition ratio is higher than 0.9 and the extracted frequency error is less than 0.0185% for the overlap coefficient γo0.1618, by which IPI image processing is effectively addressed. The experiments are conducted for standard particles of diameter 45 μm. The measurement results are that the particle diameter is 44.4970.91 μm, the relative error 1.13%, respectively. What matters is that this method can simultaneously detect particle center location with higher accuracy and the velocity can then be obtained with the particle tracking velocimetry technique. The research results demonstrate the feasibility of the method presented in this paper, which applies to interferometric out-of-focus imaging of bubbles and droplet particle, such as spray filed.
[1] Lin Jia, E.L Thomas, Physical Review B 84 (2011) 125128. [2] Lin Jia, E.L. Thomas, Journal of the Optical Society of America B 26 (2009) 1882. [3] H. Ahn, P. Thiyagarajan, L. Jia, S. Kim, J. Yoon, E.L. Thomas, J. Jang, Nanoscale 5 (2013) 1836. [4] B.W.S. Kim, L. Jia, E.L. Thomas, Adv. Mater. 21 (2009) 1921. [5] C. Lacour, D. Durox, S. Ducruix, M. Massot, Experiments in Fluids 51 (2011) 295. [6] C.F Hess, D.L' Esperance, Experiments in Fluids 47 (2009) 171. [7] N. Damaschke, H. Nobach, T.I Nonn, N. Semidetnov, C. Tropea, Experiments in Fluids 39 (2005) 336. [8] M. Maeda, Y. Akasaka, T. Kawaguchi, Experiments in Fluids 33 (2002) 125. [9] A.R Glover, S.M. Skippon, R.D. boyle, Applied Optics 34 (1995) 8409. [10] Y. Niwa, Y. Kamiya, T. Kawaguchi, M. Maeda, Proceedings of 10th International Symposium on Application of laser techniques to fluid mechanics, Lisbon, vol. 1, 2000, 〈http://in3.dem.ist.utl.pt/downloads/lxlaser2000/pdf/38_1.pdf〉. [11] A. Quérel, P. Lemaitre, M. Brunel, E. Porcheron, G. Gréhan, Measurement Science and Technology 21 (2010). (015306-1). [12] G. Lacagnina, S. Grizzi, M. Falchi, F. Di Felice, G.P Romano, Experiments in Fluids 50 (2011) 1153. [13] K.H Hesselbacher, K. Anders, A. Frohn, Applied Optics 30 (1991) 4930. [14] Pu Shiliang, Pu Xiangguo, Yuan Zhenfu, Cen Kefa, L. Mees, D. Lebrun, Proceedings of the Chinese Society for Electrical Engineering 24 (2005) 201. (in Chinese). [15] K. Jorabchi, R.G. Brennan, J.A. Levine, A. Montaser, Journal of Analytical Atomic Spectrometry 21 (2006) 839. [16] Lü Qieni, Ge Baozhen, Chen Yiliang, Zhang Yimo, Acta Optica Sinica 31 (2011) 0412009-1 (in Chinese). [17] N. Semidetnov, C. Tropea, Measurement Science and Technology 15 (2004) 112. [18] Deng Zhenmiao, Liu Yu, Wang Zhizhong, Journal of Data Acquisition and Processing 21 (2006) 473. (in Chinese). [19] J.B. Tsui, Digital Techniques for Wideband Receivers, 2nd ed., SciTech Publishing Inc., Raleigh, NC (2004) 119–126. [20] S. Kim, S.J. Lee, Experiments in Fluids 44 (2008) 623. [21] N. Damaschke, H. Nobach, C. Tropea, Experiments in Fluids 32 (2002) 143.
Acknowledgments The work is supported by National Natural Science Foundation of China (No. 61275019) and 2011 Project of State Key laboratory of Engines (No. K2011-10).