Determining particle size distributions from a single projection image

Imaging techniques employed to measure the structure of granular, particulate and porous materials are limited by scale, temporal resolution and, for biological samples, radiation exposure. This paper describes a technique for determining the distribution of particle sizes in opaque samples, for particle volume fractions less than ten percent, using a single projection radiograph. The method is based on the derived property of the additivity of the particles’ spatial autocorrelation function in projection images. Simulations and experiments demonstrate the ability to use this property to determine the distribution of particle sizes in a material. ©2012 Optical Society of America OCIS codes: (100.2960) Image analysis; (110.7440) X-ray imaging; (050.5080) Phase shift; (350.4990) Particles. References and links 1. M. E. Davis, “Ordered porous materials for emerging applications,” Nature 417(6891), 813–821 (2002). 2. H. A. Makse, S. Havlin, P. R. King, and H. E. Stanley, “Spontaneous stratification in granular mixtures,” Nature 386(6623), 379–382 (1997). 3. C. N. Davies, Aerosol Science, First ed. (Academic Press, 1966). 4. R. A. Dobbins, L. Crocco, and I. Glassmans, “Measurement of Mean Particle Sizes of Sprays from Diffractively Scattered Light,” AIAA J. 1(8), 1882–1886 (1963). 5. B. P. Flannery, H. W. Deckman, W. G. Roberge, and K. L. D’Amico, “Three-Dimensional X-ray Microtomography,” Science 237(4821), 1439–1444 (1987). 6. T. Narayanan, O. Diat, and P. Bösecke, “SAXS and USAXS on the high brilliance beamline at the ESRF,” Nucl. Instrum. Meth. A 467–468, 1005–1009 (2001). 7. L. Rigon, H.-J. Besch, F. Arfelli, R.-H. Menk, G. Heitner, and H. Plothow-Besch, “A new DEI algorithm capable of investigating sub-pixel structures,” J. Phys. D Appl. Phys. 36(10A), A107–A112 (2003). 8. H. Suhonen, M. Fernández, A. Bravin, J. Keyriläinen, and P. Suortti, “Refraction and scattering of X-rays in analyzer-based imaging,” J. Synchrotron Radiat. 14(6), 512–521 (2007). 9. R. Cerbino, L. Peverini, M. A. C. Potenza, A. Robert, P. Bosecke, and M. Giglio, “X-ray-scattering information obtained from near-field speckle,” Nat. Phys. 4(3), 238–243 (2008). 10. A. Snigirev, I. Snigireva, V. Kohn, S. Kuznetsov, and I. Schelokov, “On the possibilities of x‐ray phase contrast microimaging by coherent high‐energy synchrotron radiation,” Rev. Sci. Instrum. 66(12), 5486–5492 (1995). 11. S. W. Wilkins, T. E. Gureyev, D. Gao, A. Pogany, and A. W. Stevenson, “Phase-contrast imaging using polychromatic hard X-rays,” Nature 384(6607), 335–338 (1996). 12. M. J. Kitchen, D. Paganin, R. A. Lewis, N. Yagi, K. Uesugi, and S. T. Mudie, “On the origin of speckle in x-ray phase contrast images of lung tissue,” Phys. Med. Biol. 49(18), 4335–4348 (2004). 13. M. D. Alaimo, D. Magatti, F. Ferri, and M. A. C. Potenza, “Heterodyne speckle velocimetry,” Appl. Phys. Lett. 88(19), 191101 (2006). 14. A. Fouras, J. Dusting, R. Lewis, and K. Hourigan, “Three-dimensional synchrotron x-ray particle image velocimetry,” J. Appl. Phys. 102(6), 064916 (2007). 15. A. Fouras, D. Lo Jacono, C. V. Nguyen, and K. Hourigan, “Volumetric correlation PIV: a new technique for 3D velocity vector field measurement,” Exp. Fluids 47(4-5), 569–577 (2009). 16. S. Dubsky, R. A. Jamison, S. C. Irvine, K. K. W. Siu, K. Hourigan, and A. Fouras, “Computed tomographic xray velocimetry,” Appl. Phys. Lett. 96(2), 023702 (2010). 17. C. V. Nguyen, J. Carberry, and A. Fouras, “Volumetric-correlation PIV to measure particle concentration and velocity of microflows,” Exp. Fluids (in-press). 18. D. Paganin, S. C. Mayo, T. E. Gureyev, P. R. Miller, and S. W. Wilkins, “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” J. Microsc. 206(1), 33–40 (2002). #164928 $15.00 USD Received 19 Mar 2012; revised 25 May 2012; accepted 28 May 2012; published 28 Jun 2012 (C) 2012 OSA 2 July 2012 / Vol. 20, No. 14 / OPTICS EXPRESS 15962 19. A. Lipson, S. G. Lipson, and H. Lipson, Optical Physics, 4th ed. (Cambridge University Press, 2010). 20. M. Nieto-Vesperinas, Scattering And Diffraction in Physical Optics, 2nd ed. (World Scientific Pub Co Inc, 2006). 21. A. Pogany, D. Gao, and S. W. Wilkins, “Contrast and resolution in imaging with a microfocus x-ray source,” Rev. Sci. Instrum. 68(7), 2774–2782 (1997). 22. E. L. Crow and K. Shimizu, Lognormal Distributions: Theory and Applications (M. Dekker, 1988). 23. S. Goto, K. Takeshita, Y. Suzuki, H. Ohashi, Y. Asano, H. Kimura, T. Matsushita, N. Yagi, M. Isshiki, H. Yamazaki, Y. Yoneda, K. Umetani, and T. Ishikawa, “Construction and commissioning of a 215-m-long beamline at SPring-8,” Nucl. Instrum. Meth. A 467–468, 682–685 (2001). 24. C. Kittel, Introduction to Solid State Physics, 8th ed. (Wiley, 2005).


Introduction
Granular, particulate and porous materials and systems are of interest across a spread of disciplines including material, chemical, geological and biological science and engineering. Many existing and new porous materials have special properties that are a consequence of their unique structure [1]. To geologists, the behavior of many natural systems, such as those comprised of sand, soil and snow, can be understood by studying the behavior of granular media [2]. As such, simple methods to characterize these materials and systems are in demand.
Characterization techniques based on optical scattering [3,4] are well established but are limited to non-opaque samples. X-ray computed tomography (CT) and microtomography [5] are capable of imaging the structure of opaque samples, but suffer comparatively poor temporal resolution; furthermore, ionizing radiation dose can be problematic for biological samples. Diffraction based X-ray techniques [6][7][8][9] have to date been limited to studying particulates no larger than a few micrometres or have only been qualitative in nature.
In this paper we show that for low volume fraction suspensions of particles, the spatial autocorrelation function (SAF) of a projection image of the particles is additive, and use this property to determine the distribution of particle sizes from a single image.

Theory
Many systems of interest consist of spheroidal particles dispersed in a pseudo-random arrangement within some medium. Planar X-ray imaging of such a system can produce a discernible speckle pattern. However, when the difference in attenuation between the particles and medium is small, the speckle contrast is poor. Propagation-based phase contrast X-ray imaging (PCXI) [10,11], using a partially coherent source, can produce significantly higher speckle contrast [9,12]. A typical speckle pattern can be seen in Fig. 1. Fig. 1. Schematic of propagation-based phase contrast X-ray imaging setup and image processing steps. The sample and detector are in-line and z is the propagation distance for phase contrast imaging. A region of a typical X-ray image has been magnified to highlight the speckled appearance produced by particles. A spatial autocorrelation is performed; which is azimuthally averaged to produce the final one-dimensional autocorrelation function.
To obtain useful information from these speckled images, a statistical measure such as the SAF can be used. In the context of heterodyne imaging of particles, Alaimo et al. derive the property of linearity: that the autocorrelation function is a weighted sum of contributions from all the particles, which they use to recover the velocity distribution from a seeded flow [13]. For a side scattering (non-heterodyne) particle image velocimetry configuration, Fouras et al.
proposed that the total cross-correlation function is equal to the sum of the cross-correlation functions at each depth, and used this notion to perform three-dimensional velocity measurements from two-dimensional images [14][15][16]. Nguyen et al. extended this concept and used the SAF to infer particle concentration as a function of depth [17]. In what follows, the linearity property of the SAF: that the SAF of a projection image is equal to the sum of the SAF of each particle, will be justified in the context of X-ray imaging, with or without phase contrast, and be used to determine the distribution of particle sizes.
For a sample comprised of randomly oriented particles or pores, the fluctuations in X-ray absorption across the image are generally small. So long as the net attenuation itself is not too large we can approximate the exponential decay in X-ray intensity to first order terms: where I′(x,y) is the measured intensity, I in (x,y) is the incident intensity, µ is the material linear attenuation coefficient assuming a single homogenous material and T(x,y) is the projected thickness of the particles / pores. In phase contrast X-ray imaging, intensity is again well approximated as a linear function of particle thickness for sufficiently short propagation distances (refer to the treatment of the transport of intensity equation by Paganin et al. [18]). This approximation is particularly accurate when high spatial frequencies are attenuated, which can be a result of the limited spatial coherence of the beam, the point spread function of the detection system and/or image post processing.
In both cases, because the change in intensity can be treated as a linear function of the projected thickness, the intensity function for a distribution of particle sizes, in which the position of each particle is given by a randomly positioned delta function, can be expressed as a sum of convolutions: where I m denotes that the image mean has been subtracted. Here we assume spherical particles of radius r, but note that the theory is valid for particles of any shape if their orientation are either identical or completely random (see [19]). We define the minimum and maximum radii as R min and R max , respectively. Here N r is the number of particles of a particular radius. The power spectrum of this image is given by: where F denotes a Fourier transform. By application of the Wiener-Khinchin theorem, the autocorrelation function can be expressed as: Theoretically, the positions of particles in a suspension of non-penetrating particles cannot be truly random. However, for low volume fractions of particles, the positions of particles are only weakly correlated and the autocorrelation function can be approximated by the sum of the autocorrelation function of each particle: Since the expected power spectrum of N r randomly positioned delta functions can be shown to be N r if the mean intensity is removed [19]: For a distribution of particle sizes, if the SAF of all potential particle sizes is known, Eq. (6) is a linear inverse problem and the distribution, N r , can be determined numerically.

Simulations
Computer simulations were undertaken to validate this theory and to demonstrate the ability to perform particle sizing. Synthetic samples were generated consisting of non-penetrating spherical glass particles randomly suspended in rectangular volumes (10 mm thick) with a total volume fraction of 2%. The position of each particle was generated randomly: a position was accepted if the particle to be inserted did not collide with any particles already in the volume. If a collision did occur, a new position was generated until an unoccupied location was found. At a volume fraction of 2%, the number of collisions is small, implying that the particles' positions are only weakly correlated PCXI was simulated by reducing the volume to a projected thickness, calculating the transmitted wave function and propagating using the angular spectrum method [20]. The imaging source was a coherent monochromatic beam at 33 keV and the propagation distance between the sample and detector, z, was 50 cm (Fig. 1). The glass particles had a refractive index decrement, Re(1-n), of 4.863 × 10 −7 and an absorption coefficient of 127.36 m −1 [12]. The propagated intensity function was filtered with a Gaussian kernel to simulate the point spread function of the detector and binned to achieve an effective pixel size of 6 µm. This low pass filtering also accounts for the limited spatial coherence of an actual synchrotron beam [21]. The image mean was subtracted before a two-dimensional spatial autocorrelation was performed. The SAF was averaged over multiple windows (128 × 128 pixels) and azimuthally averaged (Fig. 1).
For a simple mixture with an equal number of particles of two discrete diameters, the SAF should be an equally weighted sum of the SAFs of single spheres of each diameter, as is shown to be the case in Fig. 2.
To model a more realistic sample, a mixture was generated in which particle diameters were sampled from a log-normal distribution: a skewed distribution whose logarithm is normally distributed, which is commonly used for the characterization of particle sizes [22]. The mean and standard deviation of the diameter's logarithm were µ = 4.72 and σ = 0.231, respectively. The SAF was averaged over 1000 windows and 50 calibration SAFs were generated covering the range of particle sizes in the sample, each generated from the image of a single sphere. Solving for the distribution of particle diameters, N r in Eq. (6), is a linear inverse problem. A least squares solver with a non-negative constraint and a second order finite-difference regularizer to penalize non-smooth solutions was used. The calculated distribution of particle diameters is presented in Fig. 3. The solution is shifted slightly away from the skew, which is the result of smoothing applied by the regularization term and depends on the choice of the regularization parameter.

Experimental results
A physical experiment, to further demonstrate the validity of the theory outlined here, was performed in experimental Hutch 3 of beamline 20B2 at the SPring-8 synchrotron in Japan [23]. The source-to-sample distance of ~210 m provided an incident wavefield of sufficient spatial coherence for propagation-based phase contrast imaging.
Because of the difficulty in obtaining particles for calibration purposes, a slightly different problem was tackled: to determine the ratio of mixtures of two particle sizes (where each nominal size is known to contain a distribution of diameters).
Nine mixtures of glass particles of two sizes (nominally 69 µm and 116 µm; by Master sizer particle size analyzer log-normal distributions: µ = 4.26, σ = 0.227 and µ = 4.72, σ = 0.231) were prepared in cuvettes (10x10x60 mm 3 ) by suspension in white petroleum jelly. The total volume fraction was again chosen to be 2%. Imaging was performed with a monochromatic beam at an energy of 33 keV, with a propagation distance of 50 cm between sample and detector. The detector's effective pixel size was 6 µm. The results for the fractions of each particle are presented in Fig. 4, and it can be seen that there is good agreement with the true values. Simulations of bimodal mixtures of particles with the same lognormal distributions used in the previous experiment were performed, and the results are presented in Fig. 5, which highlight the relationship between the error metric and the number of windows averaged over. The experimental error seen in Fig. 4 is more significant than that predicted by simulation. It is likely this discrepancy is largely due to inaccuracies in the preparation of the samples (e.g. uneven mixing). An apparent limitation evident in this experiment is the need for calibration samples. In principle it is possible to produce these samples, although a more practical approach would be to use simulated calibration samples in their place. However, the imaging system is complex and modeling it with the required accuracy is not trivial.
A theoretical limit on the maximum volume fraction for this method has already been touched on. Decreasing randomness at higher volume fractions corresponds to a breakdown of the assumption that the particles' positions are uncorrelated and the development of organization in their arrangement. This is analogous to the concept of the structure factor [24] in scattering techniques. Based on the SAF of synthetic samples, evidence of non-randomness was observed for volume fractions greater than around 2-3%. However, in simulated experiments, comparable to those used to produce Fig. 5, increasing the volume fraction of particles only had a significant impact on the accuracy of solutions above 10%.

Conclusions
The methodologies discussed in this paper are expected to be useful for studying granular and particulate materials and systems with less radiation and higher temporal resolution than CT imaging, and at a scale beyond that where diffraction based X-ray scattering techniques have been demonstrated to be applicable.
Simulations suggest the upper limit of packing fractions for these methods is around 10%. As the packing fraction increases, the assumption of randomness of the position of the particles begins to break down as the sample starts developing a degree of organization. By using the extra information present in the SAF of such systems, it is conceivable that these methodologies could be extended to more densely packed systems.