Open Access
Add to BibSonomyAbstract
We report a novel method of diffraction imaging flow cytometry to measure and analyze size distribution of microspheres. An automated and robust image processing software based on the short-time-Fourier-transform algorithm has been developed to analyze the characteristic and spatially varying oscillations of side scatters recorded as a diffraction image. Our results demonstrate that the new method allows accurate and rapid determination of single microspheres’ diameters ranging from 1 to 100μm. The capacity for analysis of light scattering by two-sphere aggregates has been demonstrated but analytical tools for characterization of aggregates by multiple microspheres remain to be developed.
©2012 Optical Society of America
1. Introduction
Microspheres of diameters ranging from 1μm to 100μm are widely used in various fields such as high-Q resonators, drug delivery, instrument calibration and optical phantoms [1, 2]. A critical need for any application is to accurately and rapidly measure the size distribution of a given microsphere sample [3]. Current methods for sizing microspheres include microscopy imaging followed by manual counting, electronic measurement with Coulter counter and optical measurement based on light scattering with flow cytometer. Among these the angle-resolved measurement of light scattered from flowing microspheres attract intense research interests because of the high accuracy of size determination based on the Mie theory and potentially high speed of measurement afforded by these methods. Several approaches have been reported for angle-resolved measurement which include the use of detector array [4, 5], scanning of the polar scattering angle via flight time mapping [6, 7] and diffraction imaging of forward scatters [8, 9]. Different from biological cells of heterogeneous structures, single microspheres yield highly symmetric and characteristic oscillation patterns in the spatial distribution of scattered light intensity as predicted by the Mie theory. The oscillatory patterns, however, can vary significantly depending on the beam profile, position of the illuminated microsphere relative to the center of incident beam and number of microspheres present in the incident beam. Light scattering by multiple microspheres is of particular research interest [10] since it occurs quite often in suspension due to aggregation or crowded flow. For these cases inverse determination of size parameters may become error-prone if scattered light intensity is known as a function of polar angle only. It is thus highly desirable to develop a method adequately robust to solve these inverse problems.
We have developed a diffraction imaging flow cytometry method for acquisition of side scatters for rapid analysis of 3D morphology of particles [11–13]. Compared to other methods, the side scatter based diffraction imaging approach allows acquisition of scattered light intensity as a function of both polar and azimuthal scattering angles with high signal-to-noise ratio, fast frame rate if a high-speed camera is used and relative simple design of optics and detection by taking advantage of the rapidly advancing imaging sensor technology. To accurately and rapidly extract a microsphere’s diameter from its diffraction image we present here a short-time-Fourier-transform (STFT) based algorithm and demonstrate its utility by obtaining the distribution of microsphere diameters. Its extension to analyze two-sphere aggregates has also been explored with a Gabor transform algorithm.
2. Method
The details of the diffraction imaging flow cytometer design have been presented elsewhere [11, 12]. Briefly, a microsphere suspension sample is injected from a stainless steel nozzle as the core fluid into a square glass cuvette of 3mmx3mm inside sides. A separate sheath fluid enters simultaneously into the same cuvette under a higher pressure than that on the core fluid. The pressure difference is controlled by syringe pumps to form a laminar flow with the core fluid’s diameter narrowing down to about 20μm at the focal point of an incident laser beam of wavelength of 532nm and power of 100mW. The laser beam is expanded and collimated to a diameter of 20mm and focused to about 20μm in diameter at the core fluid with a lens of 175mm in focal length. Figure 1(a) presents a schematic of the experiment with a microsphere illuminated by the incident beam. The flow speed was adjusted to about 4mm/s to reduce blurring of images acquired with a 12-bit 640x480 pixel CCD camera (LU075M, Lumenera) and an exposure time of 1ms. The CCD camera was oriented perpendicular to the x-axis with its sensor collecting the side scatters from the microsphere after an infinity-corrected microscope objective of 0.55 in numerical aperture (NA) and 13mm in working distance (M Plan Apo 50x, Mitutoyo) and a tube lens of 75mm in focal length. The diffraction images as shown in Fig. 1(b) were taken by the imaging unit translated to an off-focus position of Δx = 150μm toward the flow chamber. The “origin” position of Δx = 0 was determined through alignment at which one can clearly observe the interception of the core fluid, brightened by a white light illumination, with the focus of the incident laser beam. Previously we have shown that placing the imaging unit at positive Δx produces diffraction images correlated strongly with the coherently illuminated particle’s morphology while negative Δx introduces fringe patterns by the imaging optics [12]. Variation of Δx can change the angular range of the detected scatters as can be seen by comparing the images in Fig. 5.
Fig. 1 (a) Schematic of the experimental setup, Obj = microscope objective, Tub = tube lens; (b) examples of diffraction image by single polystyrene microspheres with nominal diameter, d, labeled on the left lower corner and microsphere ID number N on the right lower corner.
Four types of polystyrene microspheres have been purchased for this study with different nominal diameter d and coefficient of variation (CV): d = 9.6μm and CV ≤ 20% (No. 7510A, Duke Scientific); d = 7.9μm and CV ≤ 20% (No. 7508A); d = 5.7μm and CV ≤ 7.5% (No. 6-1-0500, Tianjin Baseline ChromTech Research Centre); d = 2.5μm and CV ≤ 7.5% (No. 6-1-0200). Before measurement, a microsphere suspension was diluted to a concentration of 4.4x105/ml with deionized water before injection into the core fluid chamber. At an acquisition rate of about 4 frames/s, a total of about 10,000 diffraction images of microspheres were acquired for each group in about an hour. From these data we selected the first 1000 images by single microspheres per group for development of an image analysis algorithm to extract their diameters.
As can be seen from Fig. 1(b), the patterns in the diffraction images, I(Z, Y), of different microspheres are highly symmetric with the pixel intensity oscillating mainly along the z-axis, where Z and Y are the pixel coordinates as shown in Fig. 1(a). Previously, we have applied a fast-Fourier-transform (FFT) based algorithm to diffraction images of microspheres and found that a sideband of largest peak frequency outside the dominant DC component may be used to determine the microsphere diameters [11]. Further evaluation of the FFT based algorithm indicated, however, that this method of microsphere sizing has poor stability due to the presence of numerous sidebands of similar amplitudes caused by the change in the period of pixel intensity oscillation in both z and y directions. To correct this problem, we have investigated several algorithms capable of analyzing signal intensity variations in both real and inverse frequency spaces such as Gabor transform, wavelet and short-time-Fourier-transform (STFT). While the Gabor transform is a special case of STFT with a Gaussian window, the wavelet algorithm utilizes different windows of varied width to extract spatial and frequency features at multiple scales. The diffraction images of concern here, however, exhibit rather limited scales for their oscillatory patterns relative to the image resolution defined by the number of pixels. For this type of images the STFT or Gabor transform is very suitable for its algorithm simplicity and processing speed with at least two-fold improvement in comparison to the wavelet analysis.
3. Results and discussions
We applied 1D STFT on a selected row of the input diffraction image which allows Fourier transforms on pixels defined by a window function g(Z-Z’) of width w. A set of spectra S(f; Z, Y) can be obtained by translating the window centered at Z along the z-axis to account for the local variation of the spatial frequency f over the row of Y. With a Gaussian window function, the STFT spectrum can be written in continuous form as the following [14]
whereThe width parameter w determines the effective number of pixels within the window centered at Z contributing to the STFT spectrum. A larger w yields better frequency resolution but poor spatial resolution. As the microsphere diameter decreases, the spatial period of intensity variations increases (see Fig. 1(b)) and one needs to use larger w to yield sufficient frequency resolution for determination of diameter, which places the lower threshold of this method for sizing microspheres at about 1μm with the imaging parameters described here.The STFT based algorithm starts by applying the discrete form of Eq. (1) on the image row with the highest pixel intensity. The value of w was initially set as twice the shortest distance between the center vertical dark fringe and the first local maximum of pixel intensity at each row. Figure 2 presents the 2D STFT power spectra P(f; Z) = |S(f; Z, Y)| on such a row of Y near the center of I(Z, Y) for two microspheres shown in Fig. 1(b). Among the 640 pixels in each row only 512 center pixels were picked for STFT calculation to reduce computing time without data padding. The 2D power spectra of PSTFT(f, Z) shown in Fig. 2 demonstrate that a well-defined single sideband exists outside the dominant DC component whose peak may serve as an anchor to determine a frequency for microsphere sizing.
Fig. 2 The 2D STFT power spectra and comparison of the 1D STFT and FFT power spectra of two microspheres. All spectral data were obtained from single rows with Δ as the pixel separation: (a) d = 7.9μm, N = 118 and row number Y = 185; (b) d = 2.5μm, N = 400 and Y = 272. The red arrow for each PSTFT(f) curve indicates the frequency location of sideband.
For comparison of STFT with FFT results performed on the same row, we compress the 2D STFT spectra into an 1D curve of PSTFT(f) by removing the dependence of P(f, Z) on Z with the following maximum operation on P for a fixed value of f
The above conversion allows rapid determination of a unique sideband associated with the characteristic oscillation of light intensity near the row center. The power spectra of STFT and conventional FFT performed on a single row Y are compared in Fig. 2. These results demonstrate clearly that the use of STFT markedly enhances the sideband associated with the characteristic oscillation of light intensity near the row center and significantly suppresses the DC component and other sidebands induced by the spatial variation of oscillation period.To identify accurately a frequency for particle sizing, we developed an iteration scheme of STFT analysis to obtain appropriate window widths w by taking advantage of the 2D nature of diffraction image data. After the STFT operation on a center row selected on the basis of highest pixel intensity, the sideband peak in PSTFT(f) was obtained as fs by searching for the first peak with f decreasing from its maximum value in the range of Fig. 2(b). Then the sideband frequency was obtained for other rows as fs(Y) with the initial value of w picked for each row using the same definition. A careful observation of the diffraction images for single microspheres in Fig. 1(b) shows that the initial values of w increase gradually from the center toward peripheral rows due to the increasing separation of the dark fringes. Consequently fs(Y) is expected to decrease gradually from the center to peripheral rows. Based on this observation, the values of w were varied to iterate the STFT calculation for those rows with fs deviating significantly from a polynomial fitting curve to fs(Y) until the deviation diminishes. The final results are presented in Fig. 3 which were used to accurately determine the true center row and the associated maximum value of fs as fs,max for a given diffraction image.
Fig. 3 The sideband frequency fs versus the row number Y on selected rows of (a) d = 7.9μm, and N = 118; (b) d = 2.5μm and N = 400. The solid lines represents the maximum sideband frequency fs,max used to determine the microsphere diameter.
An automated image processing software as described above has been developed using the MATLAB platform (MathWorks) in which the STFT function of MATLAB was used. For a diffraction image as shown in Fig. 1(b) the calculation of the maximum sideband frequency fs,max takes on average about 0.5 minute to complete on a 2GHZ computer. Figure 4 presents the histograms of fs,max within each group of the microspheres and the mean value of fs,max plotted against the nominal diameter d for each group. With the mean values and standard deviations of fs,max, we found the CV of fs,max distributions to be 15%, 13%, 3.3% and 3.2% for the groups of microspheres of d = 9.6, 7.9, 5.7 and 2.5μm, respectively, which agree well with the nominal values of CV provided by the vendors of microspheres. From the excellent linear relation between fs,max and d as shown in Fig. 4 we can establish the following equation to obtain the diameter of a microsphere ds in μm from fs,max in 1/Δ
Fig. 4 Histograms of maximum sideband frequency fs,max for four groups of 1000 microspheres with nominal diameter d. Inset: the mean values (symbols) and standard deviations (error bars) of fs,max versus d, the solid line represents a linear regression with r2 as the coefficient of determination.
It is interesting to note that microspheres may aggregate in suspension and present a significant challenge for characterization. The diffraction imaging method reported here allows detailed mapping of the scattered light intensity against both of the polar and azimuthal angles as image data. As a result, it affords the possibility for distinguishing and characterizing aggregated microspheres. To demonstrate such capability, we have performed numerical simulations of light scattering by two identical microspheres attached to each other. Simulated diffraction images of different diameters and orientations were obtained and analyzed to gain insight on the relation between parameters of aggregates and fringe patterns in the images. The two-sphere aggregate’s orientation is marked by (θ0, ϕ0) as the angles of the line connecting the spherical centers using the coordinate system shown in Fig. 1(a). The angle-resolved S11 element of the Mueller matrix was first obtained with a discrete dipole approximation (DDA) based software as the unpolarized intensity of light scatter [15, 16]. The S11 element was then projected to an area in the Z-Y place representing the CCD sensor within an angular range of Δθs with θs as the polar scattering angle with no consideration of the imaging optics. Figure 5 shows examples of both measured and simulated diffraction images. These results indicate that the two-sphere aggregates yield distinct fringe patterns in comparison to the three-sphere aggregates with the axis of rotational symmetry for the oscillatory patterns given by their connecting lines and characteristic oscillation periods determined by both of their individual and combined sizes.
Fig. 5 The diffraction images by two and three microspheres of d = 2.5μm measured at different off-focus positions or simulated with specified Δθ at wavelength λ = 532nm: (a) measured, two-sphere, Δx = 50μm; (b) measured, three-sphere, Δx = 50μm; (c) measured, two-sphere, Δx = 150μm; (d) simulated, two-sphere, Δθs = 70°, θ0 = 3π/4, ϕ0 = 7π/4, refractive indices of sphere nsp = 1.588 and host nh = 1.334 [17].
To extend our analysis to the cases of two-sphere aggregates, we adopted the Gabor transform proposed in [18] as the 2D extension to the 1D STFT defined in Eq. (1). The final transform employs a circular window and its continuous form can be written as
where w is the same window width as the one defined in Eq. (2), u0 = fcosγ and v0 = fsinγ are the Cartesian components of a sampling frequency vector in the (u, v) plane with f and γ as its magnitude and angle. By fixing f = 1/w and varying w, one can apply the 2D Gabor transform to search for specific oscillatory patterns in the input image. We implemented the above transform in MATLAB to obtain the Gabor power spectra, PG(u, v) = |G(u, v;1/w, γ, w)|, at fixed w and γ. By comparing the spectra of different γ, certain sidebands in the (u, v) plane can be clearly identified where their peak values maximize as γ reaches a value related to the orientation of the pattern in the input diffraction image I(Z’, Y’). Figure 6 below provides two examples of the spectral images for each of the two diffraction image shown in Figs. 5(c) and 5(d). One may notice that the DC peaks at u = v = 0 in Fig. 6(b) are much smaller than those in Fig. 6(a) for lack of the noise in the simulated diffraction image. The sidebands correlate with the oscillatory patterns at different angles from the horizontal direction in the input images I, which can serve as the keys to determine aggregate parameters such as the sphere diameter and orientations. The above findings demonstrate that the diffraction imaging method could be extended to analyze aggregates of multiple microspheres. A systematic study is underway on the measurement and analysis of different aggregates.Fig. 6 The Gabor power spectral images PG(u, v) with fixed window width w and different sampling vector angle γ: (a) PG obtained from the measured diffraction image in Fig. 5(c) with w = 38Δ; (b) PG obtained from the simulated diffraction image in Fig. 5(d) with w = 76Δ. Among each pair of PG images, the left one shows the sidebands of maximum peaks and the right one shows the sidebands of minimal peaks for comparison.
4. Summary
We have developed a novel method for accurate determination of diameters of single microspheres using a diffraction imaging flow cytometer and an STFT based image processing software. The presented method complements the existing ones with the advantages of reasonable speed (up to 1000 microspheres per second with a fast camera) and high accuracy. The current design of imaging optics allows determination of the microsphere diameters in the range from about 1 to 100μm by varying the off-focus position Δx. The above range can be further extended into submicron sizes with an incident beam of shorter wavelength and microscope objective of larger NA. Initial results of 2D Gabor transform based image analysis on two-sphere aggregates demonstrate its capacity for analysis of light scattering by the aggregates but more in-depth studies are needed to yield a satisfactory solution.
Acknowledgment
Y.F. acknowledges the grant support by the National Science Foundation of China (Grant #81171342).
References and links
1. K. K. Kim, D. W. Pack, M. Ferrari, A. P. Lee, and L. J. Lee, “Microspheres for drug delivery,” in Biological and Biomedical Nanotechnology, A. P. Lee, J. Lee, and M. Ferrari, eds. (Springer, 2006), pp. 19–50.
2. D. Svoboda, M. Kozubek, and S. Stejskal, “Generation of digital phantoms of cell nuclei and simulation of image formation in 3D image cytometry,” Cytometry A 75A(6), 494–509 (2009). [CrossRef] [PubMed]
3. B. P. Hanley, L. Xing, and R. H. Cheng, “Variance in multiplex suspension array assays: microsphere size variation impact,” Theor. Biol. Med. Model. 4(1), 31–38 (2007). [CrossRef] [PubMed]
4. I. K. Ludlow and P. H. Kaye, “A scanning diffractometer for the rapid analysis of microparticles and biological cells,” J. Colloid Interface Sci. 69(3), 571–589 (1979). [CrossRef]
5. S. L. Min and A. Gomez, “High-resolution size measurement of single spherical particles with a fast Fourier transform of the angular scattering intensity,” Appl. Opt. 35(24), 4919–4926 (1996). [CrossRef] [PubMed]
6. J. T. Soini, A. V. Chernyshev, P. E. Hänninen, E. Soini, and V. P. Maltsev, “A new design of the flow cuvette and optical set-up for the scanning flow cytometer,” Cytometry 31(2), 78–84 (1998). [CrossRef] [PubMed]
7. K. A. Semyanov, P. A. Tarasov, A. E. Zharinov, A. V. Chernyshev, A. G. Hoekstra, and V. P. Maltsev, “Single-particle sizing from light scattering by spectral decomposition,” Appl. Opt. 43(26), 5110–5115 (2004). [CrossRef] [PubMed]
8. C. Godefroy and M. Adjouadi, “Particle sizing in a flow environment using light scattering patterns,” Part. Part. Syst. Charact. 17(2), 47–55 (2000). [CrossRef]
9. B. Ovryn, “Three-dimensional forward scattering particle image velocimetry applied to a microscopic field-of-view,” Exp. Fluids 29(7), S175–S184 (2000). [CrossRef]
10. Y. L. Xu and B. A. S. Gustafson, “Comparison between Multisphere Light-scattering Calculations: Rigorous Solution and Discrete-Dipole Approximation,” Astrophys. J. 513(2), 894–909 (1999). [CrossRef]
11. K. M. Jacobs, J. Q. Lu, and X. H. Hu, “Development of a diffraction imaging flow cytometer,” Opt. Lett. 34(19), 2985–2987 (2009). [CrossRef] [PubMed]
12. K. M. Jacobs, L. V. Yang, J. Ding, A. E. Ekpenyong, R. Castellone, J. Q. Lu, and X. H. Hu, “Diffraction imaging of spheres and melanoma cells with a microscope objective,” J Biophotonics 2(8-9), 521–527 (2009). [CrossRef] [PubMed]
13. K. Dong, Y. Feng, K. M. Jacobs, J. Q. Lu, R. S. Brock, L. V. Yang, F. E. Bertrand, M. A. Farwell, and X. H. Hu, “Label-free classification of cultured cells through diffraction imaging,” Biomed. Opt. Express 2(6), 1717–1726 (2011). [CrossRef] [PubMed]
14. M. Portnoff, “Time-frequency representation of digital signals and systems based on short-time Fourier analysis,” IEEE Trans. Acoust., Speech Signal Process. 28(1), 55–69 (1980). [CrossRef]
15. J. Q. Lu, P. Yang, and X.-H. Hu, “Simulations of light scattering from a biconcave red blood cell using the finite-difference time-domain method,” J. Biomed. Opt. 10(2), 024022 (2005). [CrossRef] [PubMed]
16. M. A. Yurkin, A. G. Hoekstra, R. S. Brock, and J. Q. Lu, “Systematic comparison of the discrete dipole approximation and the finite difference time domain method for large dielectric scatterers,” Opt. Express 15(26), 17902–17911 (2007). [CrossRef] [PubMed]
17. X. Ma, J. Q. Lu, R. S. Brock, K. M. Jacobs, P. Yang, and X. H. Hu, “Determination of complex refractive index of polystyrene microspheres from 370 to 1610 nm,” Phys. Med. Biol. 48(24), 4165–4172 (2003). [CrossRef] [PubMed]
18. J. G. Daugman, “Complete discrete 2-D Gabor transforms by neural networks for image analysis and compression,” IEEE Trans. Acoust., Speech, Signal Process., 36(7), 1169–1179 (1988). [CrossRef]
