Feasibility testing of a pre-clinical coded aperture phase contrast imaging configuration using a simple fast Monte Carlo simulator

A simple method of simulating possible coded aperture phase contrast X-ray imaging apparatus is presented. The method is based on ray tracing, with the rays treated ballistically within a voxelized sample and with the phase-shift-induced angular deviations and absorptions applied at a plane in the middle of the sample. For the particular case of a coded aperture phase contrast configuration suitable for small animal pre-clinical imaging we present results obtained using a high resolution voxel array representation of a mathematically-defined ‘digital’ mouse. At the end of the article a link to the software is supplied. ©2013 Optical Society of America OCIS codes: (340.7440) X-ray imaging; (170.3880) Medical and biological imaging. References and links 1. G. Meggitt, Taming the Rays (Lulu Press Inc, Raleigh, NC, 2010) 2. A. Momose, T. Takeda, Y. Itai, and K. Hirano, “Phase-contrast X-ray computed tomography for observing biological soft tissues,” Nat. Med. 2(4), 473–475 (1996). 3. A. Momose, “Demonstration of phase-contrast X-ray computed tomography using an X-ray interferometer,” Nucl. Instrum. Methods Phys. Res. A 352(3), 622–628 (1995). 4. C. David, B. Nohammer, H. H. Solak, and E. Ziegler, “Differential x-ray phase contrast imaging using a shearing interferometer,” Appl. Phys. Lett. 81(17), 3287 (2002). 5. P. Cloetens, J. P. Guigay, C. De Martino, J. Baruchel, and M. Schlenker, “Fractional Talbot imaging of phase gratings with hard x rays,” Opt. Lett. 22(14), 1059–1061 (1997). 6. C. Kottler, V. Revol, R. Kaufmann, and C. Urban, “Dual energy phase contrast x-ray imaging with Talbot-Lau interferometer,” J. Appl. Phys. 108(11), 114906 (2010). 7. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nat. Phys. 2(4), 258–261 (2006). 8. A. Momose, W. Yashiro, H. Kuwabara, and K. Kawabata, “Grating-Based X-ray Phase Imaging Using Multiline X-ray Source,” Jpn. J. Appl. Phys. 48(7), 076512 (2009). 9. T. Takeda, A. Momose, K. Hirano, S. Haraoka, T. Watanabe, and Y. Itai, “Human carcinoma: early experience with phase-contrast X-ray CT with synchrotron radiation--comparative specimen study with optical microscopy,” Radiology 214(1), 298–301 (2000). 10. D. Chapman, W. Thomlinson, R. E. Johnston, D. Washburn, E. Pisano, N. Gmür, Z. Zhong, R. Menk, F. Arfelli, and D. Sayers, “Diffraction enhanced x-ray imaging,” Phys. Med. Biol. 42(11), 2015–2025 (1997). 11. V. N. Ingal and E. A. Beliaevskaya, “X-ray plane-wave topography observation of the phase contrast from a non-crystalline object,” J. Phys. D. 28(11), 2314–2317 (1995). 12. T. J. Davis, D. Gao, T. E. Gureyev, A. W. Stevenson, and S. W. Wilkins, “Phase-contrast imaging of weakly absorbing materials using hard X-rays,” Nature 373(6515), 595–598 (1995). 13. A. Bravin, J. Keyriläinen, M. Fernández, S. Fiedler, C. Nemoz, M. L. Karjalainen-Lindsberg, M. Tenhunen, P. Virkkunen, M. Leidenius, K. von Smitten, P. Sipilä, and P. Suortti, “High-resolution CT by diffraction-enhanced x-ray imaging: mapping of breast tissue samples and comparison with their histo-pathology,” Phys. Med. Biol. 52(8), 2197–2211 (2007). 14. P. Coan, A. Bravin, and G. Tromba, “Phase-contrast x-ray imaging of the breast: recent developments towards clinics,” J. Phys. D Appl. Phys. 46(49), 494007 (2013). 15. Y. Zhao, E. Brun, P. Coan, Z. Huang, A. Sztrókay, P. C. Diemoz, S. Liebhardt, A. Mittone, S. Gasilov, J. Miao, and A. Bravin, “High-resolution, low-dose phase contrast X-ray tomography for 3D diagnosis of human breast cancers,” Proc. Natl. Acad. Sci. U.S.A. 109(45), 18290–18294 (2012). #199649 $15.00 USD Received 25 Oct 2013; revised 29 Nov 2013; accepted 2 Dec 2013; published 4 Dec 2013 (C) 2013 OSA1 January 2014 | Vol. 5, No. 1 | DOI:10.1364/BOE.5.000093 | BIOMEDICAL OPTICS EXPRESS 93 16. K. A. Nugent, T. E. Gureyev, D. J. Cookson, D. Paganin, and Z. Barnea, “Quantitative Phase Imaging Using Hard X Rays,” Phys. Rev. Lett. 77(14), 2961–2964 (1996). 17. 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). 18. S. C. Mayo, A. W. Stevenson, and S. W. Wilkins, “In-Line Phase-Contrast X-ray Imaging and Tomography for Materials Science,” Materials 5(12), 937–965 (2012). 19. Y. Liu, J. Nelson, C. Holzner, J. C. Andrews, and P. Pianetta, “Recent advances in synchrotron-based hard x-ray phase contrast imaging,” J. Phys. D. 46(49), 494001 (2013). 20. 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). 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. D. H. Larsson, U. Lundström, U. K. Westermark, M. Arsenian Henriksson, A. Burvall, and H. M. Hertz, “First application of liquid-metal-jet sources for small-animal imaging: High-resolution CT and phase-contrast tumor demarcation,” Med. Phys. 40(2), 021909 (2013). 23. D. Zhang, M. Donovan, L. L. Fajardo, A. Archer, X. Wu, and H. Liu, “Preliminary feasibility study of an in-line phase contrast X-ray imaging prototype,” IEEE Trans. Biomed. Eng. 55(9), 2249–2257 (2008). 24. Z. Di, L. Zheng, H. Zhi-Feng, Y. Ai-Min, and S. Wei, “In-line phase contrast for weakly absorbing materials with a microfocus x-ray source,” Chin. Phys. 15(8), 1731–1737 (2006). 25. A. Olivo, F. Arfelli, G. Cantatore, R. Longo, R. H. Menk, S. Pani, M. Prest, P. Poropat, L. Rigon, G. Tromba, E. Vallazza, and E. Castelli, “An innovative digital imaging set-up allowing a low-dose approach to phase contrast applications in the medical field,” Med. Phys. 28(8), 1610–1619 (2001). 26. A. Olivo and R. Speller, “A coded-aperture technique allowing x-ray phase contrast imaging with conventional sources,” Appl. Phys. Lett. 91(7), 074106 (2007). 27. K. Ignatyev, P. R. T. Munro, D. Chana, R. D. Speller, and A. Olivo, “A Coded apertures allow high-energy x-ray phase contrast imaging with laboratory sources,” J. Appl. Phys. 110(1), 014906 (2011). 28. A. Olivo and R. Speller, “Modelling of a novel x-ray phase contrast imaging technique based on coded apertures,” Phys. Med. Biol. 52(22), 6555–6573 (2007). 29. T. P. Millard, M. Endrizzi, K. Ignatyev, C. K. Hagen, P. R. T. Munro, R. D. Speller, and A. Olivo, “Method for automatization of the alignment of a laboratory based x-ray phase contrast edge illumination system,” Rev. Sci. Instrum. 84(8), 083702 (2013). 30. Y. Hwu, W. L. Tsai, A. Groso, G. Margaritondo, and J. H. Je, “Coherence-enhanced synchrotron radiology: simple theory and practical applications,” J. Phys. D. 35(13), R105–R120 (2002). 31. K. M. Pavlov, T. E. Gureyev, D. Paganin, Y. I. Nesterets, M. J. Morgan, and R. A. Lewis, “Linear systems with slowly varying transfer functions and their application to x-ray phase-contrast imaging,” J. Phys. D. 37(19), 2746–2750 (2004). 32. Y. I. Nesterets, S. W. Wilkins, T. E. Gureyev, A. Pogany, and A. W. Stevenson, “On the Optimization of Experimental Parameters for X-ray In-Line Phase-Contrast Imaging,” Rev. Sci. Instrum. 76(9), 093706 (2005). 33. A. Peterzol, J. Berthier, P. Duvauchelle, C. Ferrero, and D. Babo, “New developments in simulating x-ray phase contrast imaging” Proceedings of International Symposium on Digital industrial Radiology and Computed Tomography (Lyon, France, 2007) 34. A. Bravin, V. Mocella, P. Coan, A. Astolfo, and C. Ferrero, “A numerical wave-optical approach for the simulation of analyzer-based x-ray imaging,” Opt. Express 15(9), 5641–5648 (2007). 35. P. R. T. Munro, K. Ignatyev, R. D. Speller, and A. Olivo, “The relationship between wave and geometrical optics models of coded aperture type x-ray phase contrast imaging systems,” Opt. Express 18(5), 4103–4117 (2010). 36. A. Peterzol, A. Olivo, L. Rigon, S. Pani, and D. Dreossi, “The effects of the imaging system on the validity limits of the ray-optical approach to phase contrast imaging,” Med. Phys. 32(12), 3617–3627 (2005). 37. W. P. Segars, B. M. W. Tsui, E. C. Frey, G. A. Johnson, and S. S. Berr, “Development of a 4-D digital mouse phantom for molecular imaging research,” Mol. Imaging Biol. 6(3), 149–159 (2004). 38. E. Spiller, “X-ray optics,” Advances in X-ray Analysis 42, 297 (2000). 39. A. Olivo, S. E. Bohndiek, J. A. Griffiths, A. Konstantinidis, and R. D. Speller, “A non-free-space propagation xray phase contrast imaging method sensitive to phase effects in two directions simultaneously,” Appl. Phys. Lett. 94(4), 044108 (2009). 40. T. Donath, F. Pfeiffer, O. Bunk, C. Grünzweig, E. Hempel, S. Popescu, P. Vock, and C. David, “Toward Clinical X-ray Phase-Contrast CT: Demonstration of Enhanced Soft-Tissue Contrast in Human Specimen,” Invest. Radiol. 45(7), 445–452 (2010). 41. A. V. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” J. Opt. Soc. Am. A 19(3), 472–480 (2002). 42. P. Cloetens, “Contribution to Phase Contrast Imaging, Reconstruction and Tomography with Hard Synchrotron Radiation: Principles, Implementation and Applications” Ph.D. thesis Vrije Universiteit Brussel (1999) 43. Z. Wang, Z. Huang, L. Zhang, Z. Chen, K. Kang, H. Yin, Z. Wang, and S. Marco, “Low dose reconstruction algorithm for differential phase contrast imaging,” J. XRay Sci. Technol. 19(3), 403–415 (2011). 44. T. Köhler, B. Brendel, and E. Roessl, “Iterative reconstruction for differential phase contrast imaging using spherically symmetric basis functions,” Med. Phys. 38(8), 4542–4545 (2011). 45. J. Zambelli, N. Bevins, Z. Qi, and G. H. Chen, “Radiation dose efficiency comparison between differential phase contrast CT and conventional absorption CT,” Med. Phys. 37(6), 2473–2479 (2010). #199649 $15.00 USD Received 25 Oct 2013; revised 29 Nov 2013; accepted 2 Dec 2013; published 4 Dec 2013 (C) 2013 OSA1 January 2014 | Vol. 5, No


Introduction
Medical X-ray imaging has been in use since 1896 [1].The image is conventionally formed by measuring the variations in absorption of X-ray photons along various ray paths through a sample.Improved detection sensitivity and the advent of computer tomography have greatly improved the range of applications of this modality but, until recently, the fundamental imaging principle has remained the same.The basis of this imaging modality is that contrast arises only from photon absorption.Since X-ray photons are ionizing, every absorption event in living subjects is associated with damage to cells.For a given contrast-to-noise level, the minimum number of interacting photons can be defined, resulting in a corresponding minimum level of damage.If the sample is composed of components with very similar attenuation properties, the resulting damage (i.e.dose) required to image them may be significant.Naturally, this is not a desirable property in any imaging modality, and X-ray absorption places a fundamental limitation of the application of X-ray imaging.
However, X-rays also exhibit wave-like behavior and this too can be exploited for imaging purposes, potentially reducing the imaging dose.One method which exploits X-ray wave-like properties is phase contrast imaging.When an X-ray photon traverses a sample, it undergoes a phase shift relative to a photon traversing free space.Measuring this phase-shift for several ray paths allows the formation of an image of the sample.Although direct measurement of the phase of an X-ray photon is currently impossible, interferometric methods can be used to infer the phase shift [2][3][4][5].Recent refinements have greatly speeded up the acquisition and reduced the X-ray source requirements [6][7][8].
When an X-ray traverses an interface within the sample, the induced X-ray phase-shift for non-normal angles of incidence with the gradient of the refractive index is also accompanied by a shift in the vector of propagation of the X-ray.The difference in the real part of the refractive index between tissues at X-ray energies is very small, in the range 10 −7 to 10 −8 [9], which makes the angular deflections relatively small and thus difficult to measure.Nevertheless, such small angular deflections can be analyzed by utilizing the extreme angular sensitivity in the reflectivity, in the X-ray ray wavelength range, of certain crystalline materials [10][11][12][13][14][15].Imaging by exploiting such geometries is often called analyzer-based imaging (ABI).However, this requires precise alignment and when used with a standard (polychromatic) X-ray source is quite inefficient as it relies on using the high wavelength selectivity in the reflectivity of the crystal.
Alternatively one can place a detector a long distance away from the sample (up to several meters) magnifying the lateral deflection caused by the small phase shifts at the sample interfaces.The latter technique is called free space propagation (FSP) or in-line holography.FSP is commonly used with synchrotron radiation [16][17][18][19], although it can also be achieved with micro-focus X-ray tubes [20][21][22][23][24].When using standard micro-focus X-ray sources, the available flux is low, consequently exposure times of up to hours may be necessary [21] however by using modern novel sources [22] an exposure time of a few minutes is possible.When using a micro-focus source the emitted beam will be divergent and thus the sample field of view is restricted since the projected image, at long sample-detector distances (magnifying the linear lateral deflection), must be acquired with practically-dimensioned detectors.This is especially true if the imaging system is to use standard flat panel digital detectors.The large distances required to make this technique work with standard detectors make it difficult to use this geometry in a reasonably sized pre-clinical imaging system.However, a related method to the FSP technique, the edge illumination method in its laboratory-based 'coded aperture' implementation, goes some way towards helping lift this requirement.

Coded Aperture Phase Contrast Imaging
Coded aperture imaging attempts not to measure intensity variations due to the deflection of photons onto adjacent pixels (as in FSP), but rather attempts to exploit the phase-shift by measuring detector pixel variations due to deflections of small beams (beamlets) onto non detecting areas on the same pixel.Figure 1 shows a diagram of such a system.The technique works both with synchrotron sources [25] and conventional micro-focus X-ray tubes [26].The beamlets are produced by passing a broad X-ray beam through an attenuating mask with a series of well-defined apertures: the coded aperture mask.The beamlets emerge from the apertures, pass through the sample and then through a second coded aperture mask of similar construction to the first.The second mask is parallel to the first and is aligned in such a way as to pass a given percentage of the undeflected beamlet through and onto an X-ray intensity detector.Small deflections of each beamlet by the sample vary the proportion of the beamlet passing through the second mask.The pitch of both mask apertures is set so as to match the pitch of pixels (or groups of pixels) in the X-ray detector.If there were no attenuation in the sample the variation in signal across the detector would be entirely due to phase-shifts induced by the sample.
In a coded aperture imaging configuration it is no longer required that photons cross over into different pixels in order to detect phase-shift induced angular deviations; it is merely necessary to measure the variations in the signal intensity collected by individual detector pixels.If the mask apertures have sharp edges and strongly attenuating materials are used in the construction of the mask the intensity variation will be a sensitive metric of deflection.Construction of suitable masks is well within the scope of standard lithographic techniques, as described in [27].The alignment of the coded aperture masks to each other, and to the detector pixels, is crucial as this defines the response of the system to phase shifts [28].However, the tolerated misalignment of the masks is related to the beamlet dimensions rather than to the much lower X-ray wavelength.The system is thus robust to small misalignments (1-2 μm) [29].Another advantage of such a system is that only the dimensions of the detector and those of the coded aperture masks limit the field of view.
Such a system seems a possible contender for incorporation in a reasonably sized preclinical or clinical imaging machine.However, the critical dimensions for a specific application (i.e.coded aperture pitch, source-to-first aperture distance, first aperture-tosample, and sample-to-detector distances) for such a system are not obvious.The detector pitch in particular is not a completely free parameter but is, in practice, a property of available detectors.It is thus desirable to have a method that fixes the parameter ranges that are applicable to a practical system.Some form of simulation is therefore highly desirable to optimize the design and to justify the design parameters.
The purpose of this article is thus to create and use simple simulations to test the feasibility of a practical low dose focused pre-clinical coded aperture imaging system, and to help decide if it is worth further investigation and the undertaking of more detailed and time consuming simulations.

Simulation
Simulation of a phase contrast imaging system can take two fundamental forms.One possibility is to use the precise mathematical approach based on the application of the Kirchoff formula in the Fresnel diffraction case, as explained in [30][31][32][33] with an alternative method presented in [34].However, solving this is not trivial when using complex samples, particularly samples that are unlikely to be described well by simple geometric shapes.The other possibility is to use a simpler ray tracing approach, similar to that used by Olivo and Speller in [28] which then allows the simplicity of using voxelized samples.While this is simpler for many practical cases, for cases of limited coherence this approach produces the same results that can be obtained using diffraction theory, as shown in [35,36].
The ray optic approach involves tracing rays from a source to the detector, calculating any deflection due to refractions at material interfaces traversed by the rays.Since the resultant deviation is strongly dependent on the interface angle, the sample interfaces must be sampled at a resolution of the order of the size of the beamlet at the sample.However, a high resolution voxelized sample requires a significant amount of memory and is not optimal for use in a simulation.
For X-rays of energies appropriate to imaging organic samples, the deflection angles will be very small (100s of nanoradians to several microradians).This allows us to assume that, for any ray path, the deflections within the sample can be merely added together.In other words, the X-rays will travel ballistically through the sample.The maximum combined linear deflection within practical samples is of the order of 10s of nanometres, so ignoring this seems entirely reasonable.It should thus be reasonable to project the effects of photons traversing the sample into one plane.This "collapsing" of the sample array to 3 slices, known as "thin sample approximation", provides a simple and fast way to deal with large complex voxelized samples.
While this method will ultimately provide a somewhat less accurate result than that produced using the precise mathematical approaches, it is fast and applicable to any voxelized sample.It can thus provide useful and rapidly-obtained representative images that demonstrate the possible performance of a particular imaging configuration.

Method
The simulation was written in C++ and ran on one CPU of a standard desktop computer (Dell precision T3500).As a preliminary test, water droplets were simulated in both FSP and coded aperture geometries.Following this, an appropriate sample resolution was selected and a digital mouse phantom (Moby [37]) was used to generate a flat binary file describing a cubic voxelized digital phantom of a mouse's head at a resolution of 50 μm.Figure 2 shows a render of the whole Moby phantom and two orthogonal views of the head region.Before the simulation begins some sample pre-processing is completed.First the voxel array was searched to find voxels that correspond to sample interfaces in the direction of the X-ray beam.At the identified interface voxels locations, edge-finding routines were used to evaluate the interface angle for that particular voxel in both left-right (LR) and up-down (UD) directions.Next, the whole array was processed to determine a second array specifying the refractive index of each voxel.This was performed by first converting the sample linear attenuation values to electron density by scaling from the known water value, then the electron densities were used to calculate the real part of the refractive index values using Eq.(1) as described in [38]: where E d is the electron density, r 0 the classical electron radius and λ the wavelength of the Xray photons.The result of the processing is four files containing arrays of refractive index, linear attenuation, interface angles in the LR direction and interface angles in the UD direction respectively.The linear attenuation array was then projected through, along the beam direction, simply by adding all the corresponding values in each slice.Next, the refractive index array was projected through and, whenever an interface occurred, the values in the corresponding interface and refraction arrays were used to calculate the expected angular deflection in the two directions using Snell's law.
The final result of the sample pre-processing is then three images that correspond to the attenuation and angular deflections in the LR and UD directions that would be expected when traversing the sample.Each pixel value in the attenuation slice can be replaced by e −1 × value to give the best possible attenuation image that could be achieved (in the limit of infinite numbers of photons or dose) at the sample resolution.The angular deflection images represent the best possible images that could theoretically be achieved by measurement of angular deflection at the sample voxel resolution within the ballistic assumption.Figure 3 shows an example of three pre-processed image files.The generation of coded aperture phase contrast images requires using the pre-processed images in a Monte Carlo simulation along with the following extra elements: a source of photons, a pre-sample coded aperture, a coded aperture pixel mask and finally a pixelated detector, see Fig. 4. FSP images can also be generated by removing the coded aperture elements.For each new sample, the pre-processing to obtain the "collapsed" image files takes a few minutes.After this initial operation subsequent Monte Carlo runs can be very fast, processing 10 6 photons in around 3.5 seconds on a 2.5 GHz CPU.As Monte Carlo modeling is an inherently parallel computing method, further optimization would be possible but was not deemed necessary for this particular work.

Results
A simple FSP setup was modeled first.The sample was represented by 10 μm cubic voxels describing a 1 mm diameter water-equivalent sphere with two off-center bubbles (100 μm and 300 μm).This sphere was positioned 1 meter from a monoenergetic 30 keV point X-ray source and was 'imaged' by a detector 3 meter away.A 3000 × 3000 pixel detector was used, with pixel dimensions of 2.5 × 2.5 μm.This resulted in the image shown in Fig. 5. Fig. 5.A simulated result of a 1 mm water sphere imaged using a FSP setup described in the text. 2 × 10 9 30 keV photons were tracked in the simulation.The attenuation of the water (~4% max) is not obvious due to the contrast being stretched by the edge enhancement but it is present.
The simulation was then changed to a more realistic configuration by moving the detector to 1.5 m away from the sample, enlarging the detector pixels to 50 × 50 μm, reducing the sensor size to 1500 × 1500 pixels and enlarging the water sphere and the two off axis air bubbles to 10 mm, 3 mm and 1 mm respectively.These changes render imaging the edge enhancement due to the induced phase-shifts by FSP impossible as the larger and now closer pixels integrate and thus average out the enhancement resulting in attenuation only images.Hence the geometry was switched to a coded aperture configuration.Coded aperture imaging results depend very strongly on the percentage of the undeflected beamlet passing through the second coded aperture mask [28].Data were thus acquired with various settings, as given in Table 1.The simulations resulted in images shown in Fig. 6 with plot profiles through the center 10 rows of these images shown in Fig. 7.The next stage was to attempt virtual imaging of a more realistic sample.Keeping the same geometry as for the 10 mm water sphere, the sample was replaced with the head section of the digital Moby mouse phantom.The head region provides the most easily identified anatomical features in projection images.Using the same settings as those shown in Table 1 another selection of images was acquired.These are shown in Fig. 8.

Discussion
The majority of previous phase contrast work has been performed with either very small samples or with larger samples that are almost transparent to X-rays.This is only natural as phase contrast imaging provides an obvious and high benefit in these situations.However, our research interests ultimately require imaging of larger volumes of relatively attenuating (organic) material.As a result, our images are less impressive than is common but this is due to the limitations on the accuracy of the current model and the use of conditions that are practical for our future research applications.We have also used similar brightness and contrast for sub-images in comparison (Figs. 6 and 8).Whilst this allows a fairer comparison it does prevent appreciation of the higher contrast leading to higher signal-to-noise in exposing only a small fraction of the pixels compared to standard attenuation based imaging.
To see this effect we advise examination of the plot profiles shown in Fig. 7.
The results can, in many ways, be thought of as a worst case since they contain structures which are smoothed models of real organs with simple interfaces, where the edge of each structure is normally only traversed twice.In practice organs interfaces may be complex or gradual.Many possible biological structures could thus cause far more deviation and make much better images than those resulting from this theoretical model approximating them.It would be possible to incorporate some of these complex edge features into the model to improve the contrast; however, this moves away from the worst-case scenarios and design decisions based on the improved results could be somewhat unjustified.
The energy of the X-ray photons is set at 30 keV for all the above simulations but we could easily use other energies or a range of X-ray energies to more accurately approximate an X-ray tube spectrum.Lower energies will naturally have larger differences in refractive index and 30 keV was chosen merely as this is a reasonable energy to serve as a worst-case.It is worth noting that an X-ray tube, even with 1 mm of Aluminum beam hardening, needs to have about 55 kVp applied to have a mean spectral energy of 30 keV which is relatively high for pre-clinical imaging.
All the simulations were completed with the system set up to measure phase-shifts that resulted in LR deviations.The simulations work just as well in the UD direction but the morphology of the mouse phantom in the head regions led to very slightly clearer results in the LR direction.Using more complex coded aperture mask shapes would allow phase contrast imaging in simultaneously in both LR and UD or a linear mixture of both as shown in [39].
When considering coded aperture style imaging the benefit is linked strongly to the proportion of the beamlet falling on the pixels and to the beamlet size.There is a smooth continuum of set-ups from pure attenuation based imaging through phase contrast imaging to dark field imaging as seen in Figs. 6 or 8. Assessing the optimal parameters is dependent not only on the actual sample but, more importantly, on how the information contained in the image will be used.For example, a diagnostic system may be set up to acquire planar phasecontrast projection images to identify structures that might be invisible on attenuation images [40].In this case, the delivered dose may be the same and the phase contrast used only to provide further information, as shown in Fig. 9, which shows a normal attenuation-only projection and a phase contrast enhanced image.Fig. 9. Images of a normal attenuation only projection (left) and a phase contrast enhanced coded aperture projection (right).These images were produced by tracking 10 9 30 keV photons.The phase contrast geometry was set so that the second coded aperture blocked 20% of the undeflected beamlet.
Despite both images being formed from the same number of photons the phase contrast image is sharper due to the edge enhancing effect of phase contrast.Due to the edge enhancement the simple interpretation of pixel values as a measure of attenuation is no longer valid at the edges.However for most diagnostic applications it is more important to detect the edge of a feature than to assess the actual attenuation.
It is however also possible that a system might be set up to form volumetric phase contrast reconstructions [41,42].Previous work has shown the possibility of volumetric reconstruction with phase contrast projections [43][44][45].As phase contrast imaging allows the acquisition of very low dose projections low dose volumetric reconstructions should thus also be possible.These low dose reconstructions would likely be ideal for image-guided radiotherapy (IGRT), as this only requires the resolution of the reconstruction to be about the same as the resolution of the treatment beam.It is the shape and position of nearby organs at risks or the actual target which are most important as opposed to the finer details.Hence, for most IGRT applications, the main limitation on the technique is the imaging dose, especially for cases where multiple volumetric images are required in a relatively short period of time.In the limit of very low dose imaging, phase contrast images can provide much more information, as demonstrated in Fig. 10.Coded apertures can also be adapted to work at higher energies as shown in [27] which may allow yet further dose reduction.Fig. 10.Images of a phase contrast enhanced projection (left) and a normal attenuation only projection (right).These images were both produced by tracking 1.25 × 10 6 30 keV photons corresponding to a sample dose of about 500 nGy.As both projections are composed of relatively few hits (~2 × 10 5 and ~10 6 photons respectively) they were re-binned to 500 × 500 pixels using bilinear interpolation in order to optimise displayed dynamic range.Brightness and contrast were further individually adjusted to make each projection more appropriate for reproduction.
Given the growing clinical push towards very low dose imaging, this may be a very important aspect of phase contrast imaging in the future.However, the current lack of facilities capable of routine phase contrast imaging makes it rather difficult to estimate how useful this imaging mode would be in clinical or pre-clinical practice.Some previous work [40,[46][47][48][49][50] has indicated promise, although the practical and current technical limitations may reduce its scope.

Conclusion
The results of these simulations confirm that investigation into the development of small preclinical phase-contrast imaging machines with coded aperture geometry should continue.A practical set-up might use beamlet dimensions of about three times the maximum deflection and have 67% of each beamlet falling inside the aperture with 33% falling on the absorbing septa of the post sample coded aperture.This would make the maximum phase contrast change in any pixel equal to 50% of the maximum possible signal from attenuation of the beamlet.For example, assuming similar deflections at both edges of the same detail, this would create a 100% peak-to-peak phase contrast for the detail causing maximum deflection, leaving a wide range of possible contrasts for all details creating smaller deflection.For a detector placed 1 m away, the phase-shift induced linear deviation at the detector would be roughly equal to the resulting angular deviation in radians.The maximum linear deflection we observed in the simulation of the Moby phantom was of the order of 15 microradians, suggesting that the ideal dimensions of the aperture in the mask should be about 45 μm which ties in well with current X-ray detector technology.In truth, these are larger than apertures used so far in most experiments [26][27][28][29], which leaves scope to significantly increase the phase sensitivity should the specific application require it.In the future, we hope to extend these simulations to use measured high definition experimental images from MRI or micro CT.
Ultimately this model, its future use and experience with existing coded aperture phase contrast imaging apparatus will allow us to construct a prototype pre-clinical imager and assess the practical applications of coded aperture phase contrast imaging.
An executable of the Monte Carlo simulation software used in this paper is freely available from http://users.ox.ac.uk/~atdgroup/downloads.shtml Source code is available on request.

Fig. 1 .
Fig. 1.The basic geometry of coded aperture phase contrast X-ray imaging.

Fig. 2 .
Fig. 2. Top panel: Volumetric ender of the Moby phantom.Transparency and organ color values have been set to maximize visibility.Bottom panels: Cross sections through the head area of the Moby phantom in the coronal (a) and sagittal (b) planes.

Fig. 3 .
Fig. 3. Example of the three images produced when the voxelized sample is collapsed.They represent the attenuation (a) and the angular deflection in the LR (b) and UD (c) direction that would be expected when photons traverse the sample

Fig. 4 .
Fig. 4. A flow chart showing the various elements in the Monte Carlo simulation when it is simulating a coded-aperture phase contrast imaging geometry.

Fig. 6 .
Fig.6.Series of simulation images of a 1 cm water drop with two air bubbles obtained using the settings shown in Table1.1.25 × 10 8 30 keV photons were tracked to create each image.

Fig. 7 .
Fig. 7.A Series of plot profiles formed from the center 10 rows of the images shown in Fig. 6.Individual traces are ordered by LR offset of the post sample coded aperture and have been vertically offset from each other by 450 for clarify.

Table 1 .Fig. 8 .
Fig.8.Series of simulation images of the head region of the Moby phantom using the settings shown in Table1.1.25 × 10 8 30 keV photons were tracked to create each image.