Towards Monte Carlo simulation of X-ray phase contrast using GATE.

We describe the first developments towards a Monte Carlo X-ray phase contrast imaging simulator for the medical imaging and radiotherapy simulation software GATE. Phase contrast imaging is an imaging modality taking advantage of the phase shift of X-rays. This modality produces images with a higher sensitivity than conventional, attenuation based imaging. As the first developments towards Monte Carlo phase contrast simulation, we implemented a Monte Carlo process for the refraction and total reflection of X-rays, as well as an analytical wave optics approach for generating Fresnel diffraction patterns. The implementation is validated against data acquired using a laboratory X-ray tomography system. The overall agreement between the simulations and the data is encouraging, which motivates further development of Monte Carlo based simulation of X-ray phase contrast imaging. These developments have been released in GATE version 8.2.


Introduction
X-ray phase contrast imaging has gained increasing attention over the last decades. Since X-ray phase is not directly measurable, several techniques have been developed for phase contrast, for example propagation-based imaging (or in-line holography [1]), analyzer-based imaging [2], Talbot interferometry [3], active pixel sensors [4] and speckle-based imaging [5]. The main interest in X-ray phase imaging is that the real part of the refractive index can be over three orders of magnitude larger than the imaginary part [6], thus potentially yielding a corresponding increase in sensitivity compared to attenuation-based imaging.
The spatial evolution of a monochromatic electromagnetic wave function is governed by the Helmholtz partial differential equation [7], whose solution can only be obtained in the general case using approximations, e.g. the paraxial case. Several wave-propagation mathematical models have been developed, in particular for phase retrieval imaging [8][9][10], but most methods are based on a linearization of the problem. The differences in their derivations come from the various assumptions that are made to derive the filter expression. A popular model is the description of the diffracted X-ray wave field according to the paraxial Fresnel diffraction theory [7,11,12]. Simpler approximations have also been proposed when the detected diffraction has a Laplacian signature, i.e. for small propagation distances or for large detector point spread functions with respect to the diffraction patterns, using analytical approaches based on geometrical optics [13] or simplification of the transport of intensity equation [14,15]. Intra-object scattering and diffraction may also be taken care of via a multi-slice approach, which does not require the whole object to be valid under the projected-object approximation [16]. The Wigner distribution formalism has also been proposed to model phase effects more accurately in order to take into account the changes in spatial coherence and wave-front curvature of X-rays during the radiation propagation [17][18][19][20].
The interest in accurate simulation of wave-object processes is manifold. First, access to synchrotron beamtime is limited. Accurate off-line simulation and optimization of imaging conditions would permit a reduction in the time required to find an adequate set-up, as well as foresee and minimize radiation dose. Further, realistic simulation would provide cheap data for testing of phase retrieval methods and training of machine learning based approaches. Many artifacts in phase retrieval, such as low frequency noise, are still poorly understood and may not be correctly modeled by wave based simulation. For example, scattering might explain the low frequency artifacts observed in phase retrieval, hence the interest in simultaneously taking into account phase contrast and scattering.
An interesting python toolkit has been recently implemented to provide a fully deterministic model [21]. The image formation is based on the computation of transmission functions and free-space propagation. Although photon interactions are not taken into account, polychromaticity and realistic noise are considered. A fast alternative to full-wave models has been proposed to simulate transmission and differential phase images [22] based on an a priori of empirical blurring, which is supported in practice by the source extent and the detector spatial response.
Accurate particle-based simulations for the modeling of photon-matter interactions (scattering, photoelectric effect, pair creation) are easily accessible via Monte Carlo (MC) engines such as Geant4 [23]. The use of MC techniques for realistic wave-based simulations of the propagation of an electromagnetic wave in a medium remains a challenge. Both the diffraction and refraction processes have to be accounted for to faithfully reproduce coherent X-ray imaging experiments by simulation [24]. Some works have investigated the simulation of wave propagation using MC techniques, for example for the interference modelling of elastic scattering [25,26], for phase contrast imaging using ray-tracing [27], and for phase contrast imaging using refraction [28], but these codes are application-specific. More generic Monte Carlo frameworks have been proposed to combine wave and particle interactions [29,30]. Scattered photons are for example handled separately in order to add them as incoherent effects to the intensity of the wave front calculated by the Fresnel propagation [29]. To accomodate wave and particle properties, the simulation framework is usually split in several stages as in [30][31][32], where the MC stage is used to generate a phase space which is then converted into a complex wave amplitude for propagation in the optics. Diffraction has also been tackled using the Huygens-Fresnel principle [33] using a post-processing step at the detector stage. These implementations based on Monte Carlo engines do not clearly consider polychromatic beams in their modeling and validation of the wave-effects.
The aim of this work is to contribute towards a reconciliation of the particle and wave perspectives in the simulation of phase contrast images of polychromatic beams to achieve a more realistic representation of the imaging process and to make these developments widely available through the Monte Carlo platform GATE [34]. The implementation of refraction and Fresnel diffraction processes are first presented. Then, validation results using experimental data acquired with a commercial 3D X-ray tomography system are presented on two test cases. Finally, the limits of the current implementation and its possible extension are discussed.

Formulation
If the incident beam is considered as a monochromatic electromagnetic wave of energy E, an object indexed by p = (x, y, z) can be fully described by its complex refractive index spatial distribution n E (p). For X-rays, the real part of n E is very close to, but smaller than, 1, with exceptions in the vicinity of absorption resonance frequencies of soft X-rays [35]. Therefore, the complex refractive index is usually expressed as a decrement where β is the attenuation index and δ is the refractive index decrement. For a compound K composed of several elements, the complex refractive index can be calculated from the atomic scattering factors f k,1 (E) and f k,2 (E) of element k using [36]: where a k is the atom number density of element k, r e is the classical electron radius, h the Planck constant, and c the speed of light in vacuum. Figure 1 shows the variations of the ratio of the refractive index decrement δ E over the attenuation index β E in terms of the photon energy for cortical bone and soft tissue (from tables 105 and 108 of ICRP [37]), also referred as 'Bone, Cortical (ICRP)' and 'Tissue, Soft (ICRP)' in NIST Standard Reference Database [38]. Ratio of the refractive index decrement δ E over the attenuation index β E for cortical bone and soft tissue (see ICRP report [37]) as a function of the energy. Computed via xraylib [39].

Implementation
The attenuation index β is already available via the g4emcalculator class of Geant4 [23] for every material of the phantom but not the refractive index decrement δ, which could be computed from β via a tedious integration process over all energies [40]. We choose instead to directly retrieve δ from xraylib [39], integrated to this end in the GATE platform via a software dependency. It is only in the vicinity of the discontinuities of δ (i.e. absorption edges) that possible differences might occur. This issue is not crucial in this early development stage but a calculation based on the attenuation cross sections of Geant4 will be more rigorous to estimate δ ultimately.

Snell's law
Refraction is a deterministic process due to changes in phase velocity in the traversed media. Assume we have an interface between two materials characterized by their refractive indices The angles of incidence θ 1 and refraction θ 2 of a monochromatic electromagnetic wave upon a smooth surface are related by Snell's law, When a wave enters a medium with lower refractive index, the wave will be refracted away from the surface normal and, as the incident angle approaches grazing incidence, the intensity of the reflected light increases and the intensity of the refracted light decreases. The Fresnel reflectivity for unpolarized waves is given by and above the critical angle reflection is total as illustrated in Fig. 2 for a 10 keV X-ray beam impinging on air-to-silicon interface (i.e. δ 1 = 0 and δ 2 = 5 × 10 −6 ). Total reflection therefore only occurs for grazing incidence: in this example we have θ c 1 = 89.82 • .

Implementation
Refraction was implemented in GATE as a discrete process (g4xrayboundaryprocess) which is called at the end of each step of particle propagation in the Monte Carlo simulation. When the process is triggered, it examines the previous (prestep) and the current (poststep) positions of a particle. If the two positions indicate a transition into another medium, identified by a change of refractive index, the process will be triggered to handle the refraction event.
When the control is passed to the refraction process handler, it retrieves the particle's movement information and two pointers to volumes on both sides of the boundary. Snell's law (Eq. (3)) is then applied to compute the particle direction after refraction. According to the calculated angle of refraction, the algorithm decides whether the interaction should be treated as refraction (if the photon incidence is below the critical angle) or as total reflection. In the current implementation, reflection and refraction are considered exclusive and deterministic, so that the type of interaction is selected based on the incident angle only. Fresnel reflectivity curves shown in Fig. 2 do not fully support this assumption but the binary approximation seems reasonable since the transition zone is very steep (90% drop within 0.03 • ). It is worth noting that Fresnel reflectivity is only valid for sharp, flat surfaces and more complex models would be required for rough surfaces [41]. This implementation is a direct extension to X-rays of g4opboundaryprocess, Geant4 optical photons processes, where the surface normal is given by g4transportationmanager.
It should also be noted that this discrete process for implementing refraction is deterministic at the moment but the reflectivity together with a surface roughness could be included as a probability term to improve on the binary approximation. The launching of photons and their other interactions with matter are still stochastic and handled by the Geant4 Monte Carlo engine.

Free-space propagation model
Given an incident monochromatic electromagnetic plane wave characterized by its spatial wave function Ψ E (p) and propagating in the z direction, the inhomogeneous Helmholtz partial differential equation together with the paraxial approximation [42] lead to the following approximate expression of the wave field envelope Ψ E (p) at the exit surface of the object z = z 0 : Using the approximation the wave envelope Ψ E (x, y, z = z 0 ) can be modeled as the multiplication of the incident wave envelope Ψ E (x, y, z = 0) with the transmittance function T E (x, y) of the object where the attenuation B E (p) and the phase shift Φ E (p) are projections through the refractive index distributions, with and This projection approximation holds if the scattering is sufficiently weak in the object. If we take the squared modulus of the monochromatic wave envelope we get the intensity I E (p) of the wave-field, in other words the number of transmitted photons of energy E, which is expressed by the Beer-Lambert attenuation law: where the linear attenuation coefficient of the material is Propagation over a relatively small distance D downstream the object (i.e to z = z 0 + D) can be modeled as a linear system with respect to the wave, and hence as a product in Fourier space with the Fourier transform of the bi-dimensional Fresnel propagator [11]: where f x and f y are the frequency variables conjugate to (x, y). All previous developments were derived for a plane incident wave. The generalization to a divergent geometry is possible via the Fresnel scaling theorem [7]. This is shown schematically in Fig. 3. In a divergent geometry, given the magnification factor M, the Fourier transform of the Fresnel propagator now becomes where the divergence remains small enough for the paraxial approximation to hold. We finally need to ensure that the propagator is always sampled appropriately above the Nyquist limit. Insufficient sampling rate can give rise to spurious features in simulations [43].

Implementation
Fresnel diffraction was integrated in an existing variance-reduction actor in GATE for fixed-force detection, namely the GateFixedForcedDetectionActor class. This technique has been proposed for scatter estimation in X-ray CT [44]. The real and imaginary parts of the transmittance function are independently calculated by a ray-casting algorithm, which is discretized and performed using the reconstruction toolkit RTK [45]. The ray-casting is independent for all pixels and is multi-threaded on CPU. We created a functor called Transmittance which takes the two previously computed integrals as input and returns the complex transmittance array. The propagator is then applied in Fourier space. This diffraction actor is deterministic and currently only handles mono-energetic beams. The Fresnel intensity image was added as an additional output of the fixed-force detection actor, alongside primary photons and scattered contributions.

Description of the test bench
To validate the implemented simulations, some tests were performed using an experimental bench for tomography. The X-ray tube has a transmission target consisting of a 20 µm diamond window coated with 1 µm tungsten (anode) and the source is 0.7 µm wide. The detector is a high resolution CCD camera coupled to a Gd 2 O 2 S scintillator, and the pixel pitch is 11.8 µm in both directions. Figure 4 shows a picture of the experimental bench. The detector is 47 cm away from the source and the high voltage was set to 40 kV. Two test cases were designed, one for refraction (#1) and one for Fresnel diffraction (#2). The corresponding set-ups are depicted in Fig. 5. The X-ray source spectrum used in the MC simulations is the same for both test cases. This X-ray source spectrum is the product of two contributions: • A linearly decreasing Bremsstrahlung X-ray radiant energy [46] which means that the number of photons drops in 1/E, • The transmission factor of the 100-micrometer diamond window.
The resulting X-ray source spectrum is plotted as a dashed-red curve in Fig. 6. The energy response of the Gd 2 O 2 S scintillator is modeled as the product of the mass attenuation coefficient by the photon energy (solid-blue curve in Fig. 6). Those two contributions should not be fused into a single effective spectrum in the simulation because secondary photons subsequent to radiative processes with energy change (i.e. Compton scattering and Fluorescence) would not be adequately handled.

Test case #1: Refraction
A 14 mm diameter PMMA cylinder is placed 48 mm away from the source. A 700 µm thick silicon edge is inserted between the X-ray source and the PMMA cylinder so that total reflection occurs. Figure 5(a) illustrates the protocol. For the GATE simulations, 10 9 incident particles were used and the two objects were modeled as GateBox and GateCylinder volumes respectively.

Test case #2: Fresnel diffraction
The silicon slab is now placed perpendicular to the beam so that no significant refraction should occur, and very close to the X-ray source (9.7 mm). Figure 5(b) illustrates the protocol. The magnification factor is around 48 which gives a geometric unsharpness (penumbral blur due the focal spot size) at the detector level of about 3 pixels. For the simulation, this unsharpness is modeled as four standard deviations of a Gaussian filtering. Only one particle per source position triggers this deterministic calculation. The deterministic GateFixedForcedDetectionActor was used in the GATE simulation and the slab was modeled as a GateImageNestedParametrisedVolume voxelized volume sampled at 1  µm to minimize the step discontinuities for the traversed X-ray paths which could cause large phase contrast fringes, like in CAD facetted models [47]. Each pixel value is the average of 10 different rays (pixel size 1.18 µm) to get a better sampling of the fringes.  Figures 7(b) and (c) show the corresponding simulation images without and with the proposed X-ray boundary processes for refraction. For better readability, the images have been normalized via a flat-field correction. Three profiles have been extracted to compare the simulated images to the experimental one. The location of these profiles is shown in Fig. 7(a). Figure 8(a) corresponds to the profile sampled horizontally across the silicon edge in the air region. Figure 8(b) corresponds to the profile sampled horizontally across the silicon edge in the PMMA region. Finally, Fig. 8(c) corresponds to the profile sampled vertically across the cylinder top surface.     Fig. 11. Attenuation profiles sampled per pixel for the experiment and for the GATE simulation with or without diffraction process (original images in Fig. 9). over the input spectrum (shown in Fig. 6): three such Fresnel diffraction components have been shown in Fig. 10. The detector point-spread-function has been computed from the geometrical unsharpness which is around 3 pixels. Figure 11 shows attenuation profiles sampled normally to the silicon edge.

Discussion
Considering Fig. 7, we can see that qualitatively, the simulation matches the refraction experiment quite well. On line profiles (Fig. 8), however, we can see that the refraction profile is always overestimated. This can be explained by several contributions. The most important is that currently, we consider the reflection a deterministic event, purely decided by the refractive indices of the interfacing materials and the incidence angle of the photon. This means that we are not taking into account any surface properties such as roughness, but consider all surfaces to be perfect like cleaved crystal surfaces. Further, the accurate estimation of the surface normal used in Snell-Descartes' law is straight-forward to determine in analytically defined objects like the ones used in test case #1. A remaining challenge is to calculate refraction in a voxelized geometry. In addition, the experimental setup might not have permitted a precise enough alignment of the edge to achieve the maximum reflection. The resulting image should be sufficient to validate the qualitative behavior of the code, however. Finally, no effort was made to precisely model the spectrum. A rough estimate was used to demonstrate the functioning of the code. Continuing with the diffraction case, as can be see from Fig. 9, the Fresnel diffraction simulation correctly reproduces the Laplacian-like phase-contrast image that was observed experimentally. Currently, the simulation is limited to the monochromatic case to compute the diffraction pattern since the propagator requires a complex wavefront of a single energy. The polychromatic beam is emulated by summing over diffraction images at different wavelengths (Fig. 10), i.e. via a deterministic loop over all energy bins of the X-ray spectrum. A remaining challenge is how to incorporate ballistic effects in this kind of simulation. The reconciliation of wave and particle effects in a single Monte Carlo simulation toolkit is only at the beginning. In the proposed models, diffraction and refraction processes differ in the way they are implemented in GATE: • The refraction mode is a boundary process which is defined in addition to standard physical processes in the Monte Carlo engine. This means that scattered photons are already taken into account by the simulation. No variance reduction technique is involved when the refraction process is activated.
• The Fresnel diffraction mode is based on fixed-forced detection for the simulation of the phase of the wave front, which means that this is a deterministic module, similar to a digital reconstructed radiograph. The contribution of the scattered photons is not considered in this module.
The two modules can be used together and linearly combined as in [29,48], in particular if the contribution of the scattering effects has to be considered at the same time as diffraction. Advances in simulation of phase contrast images are of great interest. Currently, high quality phase contrast images can only be obtained at synchrotron radiation facilities, although laboratorybased imaging systems are steadily showing increasing image quality (as can be attested from the experimental images in this work, for example). New developments in X-ray sources promise to push this even further.

Conclusion
We presented the implementation of X-ray refraction and Fresnel diffraction processes in GATE. This represents the first steps to an integrated X-ray phase contrast simulator taking into account both wave optical and particle effects. These new modules are part of GATE v8.2 (released in 02/2019) and macro examples can be found in the user-oriented public repository [49]. Preliminary validation tests with a laboratory X-ray source showed encouraging results. Future works include modeling of surface roughness, use of voxelized objects and closer integration of the two perspectives. This work paves the way towards fully taking into account coherent effects in a Monte Carlo model of phase contrast imaging.

Acknowledgments
Authors are very grateful to Jérôme Adrien and the MATEIS laboratory in Lyon for the experimental images.

Disclosures
The authors declare that there are no conflicts of interest related to this article.