Decoding the structure of granular and porous materials from speckled phase contrast X-ray images

Imaging techniques for studying the structure of opaque, granular and porous materials are limited by temporal resolution and radiation dose. We present a technique for characterising the structure of such materials by decoding three dimensional structural information from single, propagation based phase contrast X-ray images. We demonstrate the technique by measuring the distribution of diameters of glass microspheres in packed samples. We also present synthetic data, which shows that our inverse method is stable and that accuracy is improved by phase contrast Xray imaging. Compared to computed tomography, our technique has superior temporal resolution and lower radiation dose. ©2013 Optical Society of America OCIS codes: (100.2960) Image analysis; (110.7440) X-ray imaging; (120.5050) Phase measurement; (350.4990) Particles. References and links 1. D. P. Schmidt and M. L. Corradini, “The internal flow of diesel fuel injector nozzles: A review,” Int. J. Engine Res. 2(1), 1–22 (2001). 2. F. J. Leij, T. A. Ghezzehei, and D. Or, “Modeling the dynamics of the soil pore-size distribution,” Soil Tillage Res. 64(1-2), 61–78 (2002). 3. D. E. Carney, C. E. Bredenberg, H. J. Schiller, A. L. Picone, U. G. McCann II, L. A. Gatto, G. Bailey, M. Fillinger, and G. F. Nieman, “The mechanism of lung volume change during mechanical ventilation,” Am. J. Resp. Crit. Care 160(5), 1697–1702 (1999). 4. A. J. Hajari, D. A. Yablonskiy, A. L. Sukstanskii, J. D. Quirk, M. S. Conradi, and J. C. Woods, “Morphometric changes in the human pulmonary acinus during inflation,” J. Appl. Physiol. 112(6), 937–943 (2012). 5. R. R. Mercer, M. L. Russell, and J. D. Crapo, “Alveolar septal structure in different species,” J. Appl. Physiol. 77(3), 1060–1066 (1994). 6. J. Knust, M. Ochs, H. J. G. Gundersen, and J. R. Nyengaard, “Stereological estimates of alveolar number and size and capillary length and surface area in mice lungs,” Anat. Rec. (Hoboken) 292(1), 113–122 (2009). 7. A. Fouras, B. J. Allison, M. J. Kitchen, S. Dubsky, J. Nguyen, K. Hourigan, K. K. W. Siu, R. A. Lewis, M. J. Wallace, and S. B. Hooper, “Altered lung motion is a sensitive indicator of regional lung disease,” Ann. Biomed. Eng. 40(5), 1160–1169 (2012). 8. C. N. Davies, Aerosol Science (Academic Press, 1966). 9. C. E. Mandt, Y. Kuga, L. Tsang, and A. Ishimaru, “Microwave propagation and scattering in a dense distribution of non-tenuous spheres: experiment and theory,” Waves Random Media 2(3), 225–234 (1992). 10. J. R. Fletcher, G. P. Swift, D. C. Dai, J. A. Levitt, and J. M. Chamberlain, “Propagation of terahertz radiation through random structures: An alternative theoretical approach and experimental validation,” J. Appl. Phys. 101(1), 013102 (2007). 11. C. M. Sayers and R. L. Smith, “The propagation of ultrasound in porous media,” Ultrasonics 20(5), 201–205 (1982). 12. The ImPACT Group, Comparative Specifications: 128 to 320 Slice CT Scanner Technical Specifications (Centre for Evidence-based Purchasing, 2009). 13. R. Mokso, P. Cloetens, E. Maire, W. Ludwig, and J.-Y. Buffière, “Nanoscale zoom tomography with hard x rays using Kirkpatrick-Baez optics,” Appl. Phys. Lett. 90(14), 144104 (2007). 14. M. Langer, A. Pacureanu, H. Suhonen, Q. Grimal, P. Cloetens, and F. Peyrin, “X-ray phase nanotomography resolves the 3D human bone ultrastructure,” PLoS ONE 7(8), e35691 (2012). #191926 $15.00 USD Received 7 Jun 2013; revised 8 Jul 2013; accepted 22 Jul 2013; published 5 Aug 2013 (C) 2013 OSA 12 August 2013 | Vol. 21, No. 16 | DOI:10.1364/OE.21.019153 | OPTICS EXPRESS 19153 15. R. Mokso, F. Marone, D. Haberthür, J. C. Schittny, G. Mikuljan, A. Isenegger, and M. Stampanoni, “Following dynamic processes by X‐ray tomographic microscopy with sub‐second temporal resolution,” in AIP Conference Proceedings (2011), Vol. 1365, pp. 38–41. 16. J. Moosmann, A. Ershov, V. Altapova, T. Baumbach, M. S. Prasad, C. LaBonne, X. Xiao, J. Kashef, and R. Hofmann, “X-ray phase-contrast in vivo microtomography probes new aspects of Xenopus gastrulation,” Nature 497(7449), 374–377 (2013). 17. S. Fain, M. L. Schiebler, D. G. McCormack, and G. Parraga, “Imaging of lung function using hyperpolarized helium-3 magnetic resonance imaging: Review of current and emerging translational methods and applications,” J. Magn. Reson. Imaging 32(6), 1398–1408 (2010). 18. O. Glatter and O. Kratky, Small Angle X-ray Scattering (Academic Press, 1982). 19. U. Bonse and M. Hart, “Small angle X-ray scattering by spherical particles of polystyrene and polyvinyltoluene,” Z. Phys. 189(2), 151–162 (1966). 20. 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). 21. R. P. Carnibella, M. J. Kitchen, and A. Fouras, “Determining particle size distributions from a single projection image,” Opt. Express 20(14), 15962–15968 (2012). 22. J. H. Holland, Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence (MIT Press, 1992). 23. L. Randy, Haupt, Practical Genetic Algorithms, 2nd ed. (Wiley Interscience, 2004). 24. 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). 25. 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). 26. 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). 27. S. Torquato, T. M. Truskett, and P. G. Debenedetti, “Is random close packing of spheres well defined?” Phys. Rev. Lett. 84(10), 2064–2067 (2000). 28. S. B. Yuste and A. Santos, “Radial distribution function for hard spheres,” Phys. Rev. A 43(10), 5418–5423 (1991). 29. C. James, Spall, “Stochastic optimization,” in Handbook of Computational Statistics: Concepts and Methods, Wolfgang Härdle, James E. Gentle, and Yuichi Mori, eds. (Springer, 2004), pp. 169–197. 30. V. Alan, Oppenheim, Discrete-time Signal Processing, 2nd ed. (Prentice Hall, 1999). 31. A. W. Klein, R. F. Becker, and M. R. Bryson, “A method for estimating the distribution of alveolar sizes from histological lung sections,” Trans. Am. Microsc. Soc. 91(2), 195–208 (1972). 32. L. Knudsen, E. R. Weibel, H. J. G. Gundersen, F. V. Weinstein, and M. Ochs, “Assessment of air space size characteristics by intercept (chord) measurement: an accurate and efficient stereological approach,” J. Appl. Physiol. 108(2), 412–421 (2010). 33. C. Kloss, C. Goniva, A. Hager, S. Amberger, and S. Pirker, “Models, algorithms and validation for opensource DEM and CFD–DEM,” Prog. Comput. Fluid Dy. 12, 140–152 (2012). 34. 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). 35. M. Nieto-Vesperinas, Scattering And Diffraction in Physical Optics, 2nd ed. (World Scientific, 2006). 36. L. Eshelman and J. Schaffer, “Real-coded genetic algorithms and interval-schemata.,” in Foundations of Genetic Algorithms, D. Whitley, ed. (Morgan Kaufmann., 1993), pp. 187–202. 37. 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). 38. M. J. Berger, J. H. Hubbell, S. M. Seltzer, J. Chang, J. S. Coursey, R. Sukumar, D. S. Zucker, and K. Olsen, “XCOM: Photon Cross Sections Database,” http://www.nist.gov/pml/data/xcom.index.cfm. 39. M. J. Kitchen, R. A. Lewis, M. J. Morgan, M. J. Wallace, M. L. Siew, K. K. W. Siu, A. Habib, A. Fouras, N. Yagi, K. Uesugi, and S. B. Hooper, “Dynamic measures of regional lung air volume using phase contrast x-ray imaging,” Phys. Med. Biol. 53(21), 6065–6077 (2008). 40. R. A. Jamison, J. A. Armitage, J. Carberry, M. J. Kitchen, S. B. Hooper, and A. Fouras, “Functional imaging to understand biomechanics: a critical tool for the study of biology, pathology and the development of pharmacological solutions,” Curr. Pharm. Biotechnol. 13(11), 2128–2140 (2012). 41. J. Thurgood, S. Hooper, M. Siew, M. Wallace, S. Dubsky, M. Kitchen, R. A. Jamison, R. Carnibella, and A. Fouras, “Functional lung imaging during HFV in preterm rabbits,” PLoS ONE 7(10), e48122 (2012).


Introduction
The study of granular and porous materials is relevant across a range of disciplines including material science, mechanical and chemical engineering, geophysics and biology. The development of cavitation bubbles in diesel injectors [1], the change in soil pore size during wetting and drying [2] and the mechanics of alveoli in the breathing lung [3] are examples of dynamic systems, which existing imaging techniques are not well suited to studying.
Expanding briefly on the example of the lung, despite much research, respiratory physiologists' understanding of alveolar mechanics during respiration remains incomplete. The degree to which alveolar expansion and recruitment each contribute to inflation continues to be a topic of debate [3,4]. Dynamic imaging of lung alveoli is a particularly difficult problem because of their small size compared to that of the entire lung and the velocity at which they move. A free breathing mouse, for example, has millions of alveoli, each of the order of 100 um in diameter [5,6], packed inside a chest several centimetres wide, which respire 3-4 times a second. In addition to physiological studies, there is potential for disease detection applications using dynamic lung imaging techniques [7].
To non-destructively study the kinds of samples described, optical and other low energy imaging modalities, from visible light to terahertz frequencies [8][9][10], suffer from scattering and limited penetration depth. Ultrasound imaging has similar drawbacks [11].
X-ray computed tomography (CT) is a powerful and widely used technique for visualizing the three-dimensional structure of a sample. Despite not being inherently suited to studying dynamic systems because of the requirement for multiple projections at any time point, modern scanners are still capable of quite respectable temporal resolution. Third generation clinical CT scanners are capable of scanning at around 3-4 frames per second with sufficient spatial resolution to differentiate structures down to about half a millimetre in size [12]. However, this is still slower than many dynamic processes and the spatial resolution is too coarse for micro-imaging. Micro-tomography and even nano-tomography [13,14] are possible using an ultra-bright synchrotron source. A state of the art micro-tomography scanner, with an ultra-bright synchrotron radiation source, has been demonstrated capable of scanning an infant rat's chest in half a second [15]. Another interesting example is 4D tomography of Xenopus gastrulation [16]. However, these scanners are incapable of the frame rates required to dynamically image rapidly moving systems such as breathing lungs. In addition to temporal resolution, a significant problem with CT is the high radiation dose, which is particularly relevant for biological samples.
In lung imaging, hyperpolarized helium diffusion MRI has recently been used to try and measure alveolar dynamics [4]. There are a number of problems with this method, which include the difficulty and expense of obtaining hyperpolarized helium, the fact that diffusion is a relative measurement, and that airspace dimensions can only indirectly be obtained from these diffusion measurements [17].
At submicron scales, X-ray scattering techniques can be used to obtain structural information from single projection images [18][19][20]. Since scattering angles are inversely proportional to the size of the scattering particles, the problem of measuring increasingly smaller scattering angles limits the maximum feature size which can be studied.
For studying particles tens to hundreds of microns in size, Carnibella et al. [21] presented a method for determining basic morphological parameters of randomly packed particles at low packing fractions from a single projection X-ray image. This method was based on the encoding of three-dimensional structural information in speckled two-dimensional X-ray images of the samples. Specifically, this method exploited the linearity of the spatial autocorrelation function (SAF) of these images. For densely packed systems of particles, the positions of individual particles are more strongly correlated and the SAFs are no longer linear. Therefore, a non-linear inverse method is necessary to decode their structure. In this paper, we present a technique using a genetic algorithm [22,23] (GA) to recover morphological parameters of randomly packed particles, without restriction on the packing fraction.
We begin by describing the theoretical and technical background of the technique. We demonstrate, experimentally, its application in determining the distribution of diameters in a packed sample of glass microspheres. Finally, we use synthetic X-ray images of microspheres to characterise the performance of our GA and the efficacy of propagation based phase contrast.

Description of technique
When a sample of randomly packed particles is imaged using absorption or phase contrast Xray imaging [24,25] (PCXI), the resulting image has a speckled appearance [26]. Examples of these images can be found in Figs. 1 and 3. We make note that this speckle is not necessarily so called near field speckle [20] and that, in fact, it is also produced in absorption based imaging. Speckle contrast, however, is markedly improved in PCXI. Carnibella et al. [21] showed that at low packing fractions, the SAF of a projection image of a sample is equal to the sum of the SAF of the image of each particle in that sample. This implies that the speckled images contain information about the morphology of the particles, since the SAF of the image of a single particle is directly related to its morphology. In the same paper, it was also found that at higher packing fractions this behaviour began to break down, which corresponded with the development of oscillations in the previously single peaked SAFs. It was hypothesised that these oscillations were related to both the shape (form) and organisation (structure) of the particles and that it should be possible to recover details of these properties. Fig. 1. An overview of the complete process of solving for particle statistics. On the left, a single phase contrast image of the sample is taken at a propagation distance, z, and its spatial autocorrelation function (SAF) calculated. On the right the genetic algorithm (GA) iterates over simulated particle SAFs to find the particle parameters that produce the closest match between the two SAFs.
Ideally, we seek an analytical expression for the SAF as a function of some parameters related to the morphology and/or structure of the particles. Then, an appropriate optimisation method could be applied to solve for these parameters. To our knowledge, no such analytical expression exists describing the distribution of particles in samples approaching random close packing [27]. We note that for lower packing densities (suspensions), the Percus-Yevick approximation can be used to obtain the radial distribution function analytically [28], and that this could potentially be used as the basis of an inverse method.
Our method is based on stochastic modelling of particle packing, combined with physical modelling of the imaging process, using an iterative approach to solve the forward imaging problem until simulated results match experimental data. The particles are selected randomly from a distribution described by some parameters, and a GA attempts to find the set of these parameters that produces the best match between simulated and experimental SAFs. A simple outline of this procedure can be found in Fig. 1. GAs are particularly suited to this task because they tend to locate global minima and are robust in the presence of variability related to the stochastic nature of the simulations [29].
The experimental SAF is easily obtained after imaging a sample. Propagation based phase contrast X-ray images can be obtained by imaging with a coherent X-ray source such as a synchrotron or micro-focus lab source. The images are pre-processed by applying flat and dark field corrections: raw dark flat dark .
The SAF is obtained by dividing the image into N smaller windows, averaging their normalised 2D SAFs, and radially averaging: where M is the number of pixels that have a centre within a circular band defined by the radii: l-0.5 and l + 0.5. Ideally, the window size should be such that it can capture the full SAF. The point at which the SAF settles near zero depends on the size of the particles. We found choosing a window size that captures at least the first three peaks of the SAF to be sufficient. Windowing is necessary because a single autocorrelation is not a consistent estimator of the true SAF: that is, the variance of the autocorrelation function does not decay to zero as the window size is increased [30].
We now outline the process of obtaining a simulated SAF, produced by a given set of parameter estimates. In the context of a GA, a set of parameter estimates is referred to as a chromosome. The simulation process begins with a description of the distribution of particle diameters. We use a volume weighted, lognormal probability density function (PDF), which can be characterised by two parameters: its geometric mean (GM) and geometric standard deviation (GSD). Many collections of particles including alveoli are well approximated by a lognormal distribution [31,32], however, any other PDF could be similarly used. The pouring of particles is simulated using molecular dynamics software (LIGGGHTS [33]). We assume the particles are spherical and non-penetrating. The mean packing fraction achieved was 60.0 per cent.
Next, we simulate X-ray images of these particles. We consider the case of propagation based PCXI, the modality used in the experiments presented later. Following the method outlined by Kitchen et al. [26], the sample volume is first reduced to two dimensions by making a projection approximation [34]. Then, assuming a plane wave source, the exit wave field is calculated: where µ is the linear attenuation coefficient, δ is the refractive index decrement and t is the projected thickness of the sample. This wave field is then propagated a distance, z, to the plane of the detector using the angular spectrum method [35]. The intensity function of this wave field is convolved with a Gaussian kernel, described by a single parameter: its standard deviation, σ blur , which accounts for the combined effects of partial coherence, the point spread function of the detector and penumbral blurring: where N(µ,σ) is a normal distribution with mean, µ and standard deviation, σ. This last operation is simply the application of a Gaussian low pass spatial filter. The resulting image is binned to match the effective pixel size of the actual detector. The normalised two dimensional SAFs of smaller sample windows within the image are calculated. The SAF of multiple windows are averaged and that result radially averaged to produce the final simulated SAF. The GA requires an initial randomly generated population of N pop chromosomes. Each chromosome consists of the following parameters: the GM and the GSD of the distribution of particle diameters and the standard deviation of the Gaussian low pass filter kernel. Before running the algorithm it is necessary to set upper and lower limits for each parameter. Choosing limits which bracket smaller ranges will decrease convergence time. However, care must be taken to ensure that the range of each parameter encompasses the true, but unknown value. The cost of each chromosome is determined by calculating the sum of squared errors between the experimental SAF and a simulated SAF generated from each chromosome's parameters.
The N pop /2 chromosomes with the highest costs are discarded. Pairs of chromosomes are randomly selected, from a rank weighted population [23], to be parents. Uniform crossover using the BLX-0.5 operator [36] is performed to produce two offspring from each pair of parents. At the end of this process the population size remains unchanged. Mutations, according to a Gaussian distribution (standard deviation, σ mutation ) are randomly applied to a fraction, f mutation , of the parameters from any chromosomes other than that with the lowest cost (elitism): where x is the value of the parameter being mutated, N is a normal distribution and x lo and x hi are the lower and upper limits of that parameter. This process is repeated until the cost function shows no signs of improvement, at which point we assume the algorithm has reached the vicinity of the global minimum. Mutations are then disabled and the algorithm is allowed to continue until all the parameters have converged. The specific stopping criteria used in the following sections are described in each section.

Application to experiments
To test our technique, we conducted experiments on beamline BL20B2 at the Spring-8 synchrotron, Japan [37]. A detuned Si(111) double crystal monochromator was used to generate a monochromatic beam, with an energy of 34 keV to image solid glass microspheres. The sieved microspheres were nominally of 4 different ranges of diameters (63-75 µm, 75-90 µm, 90-106 µm, 106-125 µm). Plastic cuvettes (10 × 50 × 10 mm 3 ) filled with these microspheres were placed at a propagation distance of 0.5 m from the detector. The detector was a Hamamatsu C4880-41S, with an effective pixel size of 5.9 um. Exposure times were 2 seconds. The SAF was averaged over approximately 200 windows (96 × 96 pixels, 1132.8 × 1132.8 µm 2 ).
To simulate SAFs for the GA, the sample volume had the imaging area of a single experimental window and was 1 mm deep (1132.8 × 1132.8 × 1000 µm 3 ). Spheres were modelled as soda lime glass, which when imaged at 34 keV have a linear absorption coefficient of 198.0 m −1 and a refractive index decrement of 4.67 × 10 −7 [38]. The SAF was averaged over 16 windows (chosen as an acceptable balance between accuracy and computational effort). Other relevant parameters were the same as for the experimental setup.
The GA was run with a population size of 12 and with the following limits on each parameter's range: 20 µm to 120 µm for the GM, 1.01 to 1.3 for the GSD and 6 µm to 30 µm for the standard deviation of the Gaussian kernel. These limits were chosen to encompass all reasonable solutions, given we had prior information on the microspheres' diameters. The mutation distribution, σ mutation , was 0.2 of each parameter range and f mutation was 0.2. The GA was deemed to have reached the global minimum when the mean slope of the cost function for the last 10 iterations was equal to or greater than zero. With mutations disabled, the algorithm was allowed to continue until the difference between parameters across the population was less than 1 per cent. At this point the chromosome with the lowest cost was chosen as the solution.  Computational effort is almost solely attributable to the simulation of pouring spheres. Running on a single 3 GHz processor core, for the window size described, the pouring simulations took on average around 5 minutes. For a population size of 12, when averaging over 16 windows, this means 192 pouring simulations per iteration of the GA. Fortunately, since every simulation window is completely independent, if enough cores are available, an entire iteration could be completed in a few minutes. Typically, less than 40 iterations were required for the GA to converge.
Our results are displayed in Fig. 2 and quantified in Table 1. As can be seen, the mean microsphere diameters agree well with both the nominal particle sizes and those obtained by a commercial measurement system (Mastersizer 2000). Interestingly, the spread of diameters is less than that measured by the commercial system, and more closely agrees with the nominal sieve limits.
Two second exposure times were used to obtain a high signal to noise ratio for testing this technique on a strongly attenuating sample and does not represent the best temporal resolution of the technique. For example, at this beamline much shorter exposure times (tens of milliseconds) are possible within the lungs [7,39,40], with frame rates as high as three hundred frames per second even achievable [41]. In comparison to the technique presented here, a high resolution CT scan may require of the order of 1000 projections. Hence we could obtain the same morphological information at 1000 times the frame rate and with 1/1000 the radiation dose of CT.

Synthetic studies
To characterise the performance of our solver and study the efficacy of X-ray phase contrast for sizing particles, we applied the technique to measure the distribution of sphere diameters from synthetically generated X-ray images. The synthetic images were generated using the same method as previously described. The distribution of sphere diameters had a GM of 56.7 µm and a GSD of 1.14. These parameters are those of the largest sieved microspheres (106-125 µm), which we had previously measured (Fig. 2). The width of the Gaussian kernel was also obtained from the previous experiment, having imaged samples at a number of propagation distances. Photon and readout noise were added according to the camera's specifications and assuming photon counts similar to those recorded experimentally. The synthetic images were scaled to match the mean intensity and contrast of the 10 mm deep samples we measured experimentally (in our simulations other speckle characteristics were independent of sample depth). Synthetic images were generated at propagation distances between zero and three metres. Figure 3 shows that, to the eye, these synthetic speckle images closely match the corresponding experimental images. Synthetic experimental SAFs were produced by averaging 400 windows and radially averaging. The GA was run with the following limits on each parameter's range: 48 µm to 68 µm for the GM, 1.05 to 1.2 for the GSD and 6 µm to 35 µm for the standard deviation of the Gaussian kernel. To produce an estimate of uncertainty in the solution, upon convergence the algorithm was restarted (by keeping the best and randomly generating 11 new chromosomes) and run until the variables converged again. In this way the solver was allowed to converge a total of 10 times, which produced 10 solutions. In Figs. 4(a) and 4(b) the parameters belonging to each of these solutions is plotted. The spread of the GAs solutions (the points in 4(a) and (b)) at each propagation distance isn't obviously correlated with the propagation distance or image contrast (Fig. 4(e)). This is consistent with this variability being the result of the stochastic nature of the simulations. Therefore, we expect that this variance can be reduced by averaging the SAF over more windows, at the expense of increased computation time.
Overall, the best solutions are quite close to the true distributions of the microspheres' diameters. Comparing the extremes of these solutions (at z = 0 m and z = 1 m) in Fig. 4(c) highlights that even in the worst case (z = 0 m), the solution is still quite accurate: the GM and GSD are in error by 1.5% and 2.9%, respectively. In Fig. 4(d) it can be seen that accuracy is worst at very short propagation distances, where phase contrast is weak, and at the largest propagation distances, where the size of the Gaussian kernel is greatest (not shown). The inverse relationship between the solution error and the speckle contrast is illustrated by comparing Figs. 4(d) and 4(e). Since speckle contrast is proportional to the signal to noise ratio, this relationship is not surprising. We note that our experimental images were of static samples, with a high intensity synchrotron source and significant exposure times. We expect that the benefit afforded by phase contrast would be of greater significance when imaging conditions are less ideal.
We encountered an unexpected result in the behaviour of the Gaussian low pass filter kernel over the range of propagation distances. We found that this parameter increased almost linearly over the range of propagation distances, at a rate greater than expected by penumbral blurring for the known source size of 150 µm (horizontal) by 10 µm (vertical). As mentioned, this accounts for the decrease in contrast at the largest propagation distances. Our hypothesis is that this may be the result of scattering/refraction within the volume of the sample, which we assume is negligible when we make the projection approximation. Alternatively, or additionally, it may be the result of scattering by optical elements and/or air along the beam path which hasn't been accounted for in our model. These hypotheses warrant further investigation.

Conclusions
In this paper we have outlined a new technique for the quantification of useful parameters related to random granular and porous systems without restriction on the mode of packing. The basis of the technique lies in retrieving structural information encoded in speckled phase contrast X-ray images of such systems. The experimental and synthetic results presented demonstrate the accuracy and robustness of this technique and that propagation based phase contrast significantly improves the accuracy. We have also shown that there is a limit to the gains that can be obtained by increasing the propagation distance due to an increase in blurring. The advantages of our technique are high temporal resolution and low radiation dose, which suggest potential applications for imaging dynamic and biological systems. Specifically, we envisage the technique being applied to the measurement of dynamic lung morphology, for physiological studies or disease detection.