Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

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

Open Access Open Access

Abstract

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.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. 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 [810], 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 [1720].

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 [3032], 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.

2. Materials and methods

2.1 Complex refractive index

2.1.1 Formulation

If the incident beam is considered as a monochromatic electromagnetic wave of energy $E$, an object indexed by $\boldsymbol {p}=(x,y,z)$ can be fully described by its complex refractive index spatial distribution $n_E\left (\boldsymbol {p}\right )$. 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

$$n_E \left( \boldsymbol{p} \right) = 1 - \delta_E \left( \boldsymbol{p} \right) +i \beta_E \left( \boldsymbol{p} \right)$$
where $\beta$ is the attenuation index and $\delta$ 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]:
$$n_E \left( \boldsymbol{p} \right) = 1 - \frac{r_e}{2 \pi} \left(\frac{hc}{E}\right)^2 \sum_{k\in K} a_k\left( \boldsymbol{p} \right) \left( f_{k,1}\left( E \right) - i f_{k,2}\left( E \right) \right)$$
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 $\delta _E$ over the attenuation index $\beta _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].

 figure: Fig. 1.

Fig. 1. Ratio of the refractive index decrement $\delta _E$ over the attenuation index $\beta _E$ for cortical bone and soft tissue (see ICRP report [37]) as a function of the energy. Computed via xraylib [39].

Download Full Size | PDF

2.1.2 Implementation

The attenuation index $\beta$ is already available via the g4emcalculator class of Geant4 [23] for every material of the phantom but not the refractive index decrement $\delta$, which could be computed from $\beta$ via a tedious integration process over all energies [40]. We choose instead to directly retrieve $\delta$ 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 $\delta$ (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 $\delta$ ultimately.

2.2 Refraction and total reflection

2.2.1 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 $n_k=1-\delta _k + i \beta _k$ with $k\in \{1,2\}$. The angles of incidence $\theta _1$ and refraction $\theta _2$ of a monochromatic electromagnetic wave upon a smooth surface are related by Snell’s law,

$$\left(1-\delta_{1}\right)\sin\theta_1 = \left(1-\delta_{2}\right)\sin\theta_2.$$
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
$$R_F\left(\theta_1\right) = \frac{1}{2}\left[\left|\frac{\left(1-\delta_{1}\right)\cos\theta_1-\left(1-\delta_{2}\right)\cos\theta_2}{\left(1-\delta_{1}\right)\cos\theta_1+\left(1-\delta_{2}\right)\cos\theta_2}\right|^2 + \left|\frac{\left(1-\delta_{1}\right)\cos\theta_2-\left(1-\delta_{2}\right)\cos\theta_1}{\left(1-\delta_{1}\right)\cos\theta_2+\left(1-\delta_{2}\right)\cos\theta_1}\right|^2\right]$$
and above the critical angle
$$\theta_{1}^{\:\!c} = \mbox{asin}\left(\frac{1-\delta_{2}}{1-\delta_{1}}\right)$$
reflection is total as illustrated in Fig. 2 for a 10 keV X-ray beam impinging on air-to-silicon interface (i.e. $\delta _1 = 0$ and $\delta _2 = 5\times 10^{-6}$). Total reflection therefore only occurs for grazing incidence: in this example we have $\theta _{1}^{\:\!c} = 89.82^ {\circ }$.

 figure: Fig. 2.

Fig. 2. Air-to-silicon interface: reflectance and transmittance of an unpolarized $10$ keV X-ray wave (see Eq. (4) with $\delta _1 = 0$ and $\delta _2 = 5\times 10^{-6}$).

Download Full Size | PDF

2.2.2 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^\circ$). 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.

2.3 Wave propagation

2.3.1 Free-space propagation model

Given an incident monochromatic electromagnetic plane wave characterized by its spatial wave function $\Psi _E(\boldsymbol {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 ${\Psi }_E(\boldsymbol {p})$ at the exit surface of the object $z=z_0$:

$${\Psi}_E\left(x,y,z=z_0\right) \approx \exp\left\{ -i \frac{\pi E}{h c} \int_{z=0}^{z=z_0} \left[1 - n_E^2\left(\boldsymbol{p}\right)\right] \mbox{d}z\right\} {\Psi}_E\left(x,y,z=0\right).$$
Using the approximation
$$1 - n_E^2\left(\boldsymbol{p}\right) \simeq 2 \left[1 - n_E\left(\boldsymbol{p}\right)\right]$$
the wave envelope ${\Psi _E}\left (x,y,z=z_0\right )$ can be modeled as the multiplication of the incident wave envelope ${\Psi _E}\left (x,y,z=0\right )$ with the transmittance function $T_E(x,y)$ of the object
$$T_E \left(x,y \right) = \exp \left[ -B_E \left( \boldsymbol{p} \right)\right] \exp\left[i \Phi_E \left( \boldsymbol{p} \right) \right]$$
where the attenuation $B_E\left ( \boldsymbol {p} \right )$ and the phase shift $\Phi _E \left ( \boldsymbol {p} \right )$ are projections through the refractive index distributions, with
$$B_E \left( x,y,z=z_0 \right) =\frac{2 \pi E}{h c} \int _{z=0}^{z=z_0} \beta_E \left( \boldsymbol{p} \right) \mbox{d}z$$
and
$$\Phi_E \left( x,y,z=z_0 \right) = - \frac{2 \pi E}{h c} \int _{z=0}^{z=z_0} \delta_E \left( \boldsymbol{p} \right) \mbox{d}z.$$
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(\boldsymbol {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:
$${I_E}\left(x,y,z=z_0\right) = \left|{\Psi}\left(x,y,z=z_0\right)\right|^2 \approx \exp \left[ -2 B_E \left( \boldsymbol{p} \right)\right] {I_E}\left(x,y,z=0\right)$$
where the linear attenuation coefficient of the material is
$$\mu_E \left( \boldsymbol{p} \right) = \frac{4 \pi E}{h c} \,\beta_E \left( \boldsymbol{p} \right).$$

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]:

$$\widehat{P}_{E,D} \left( f_x, f_y \right) = \exp \left[ -i \frac{\pi h c}{E} \,D \left(f_x^2 + f_y^2\right) \right]$$
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.

 figure: Fig. 3.

Fig. 3. Fresnel scaling theorem.

Download Full Size | PDF

In a divergent geometry, given the magnification factor $M$, the Fourier transform of the Fresnel propagator now becomes

$$\widehat{P}_{E,D} \left( f_x, f_y \right) = \exp \left[ -i \frac{\pi h c}{E} \, {D} {M}\left(f_x^2 + f_y^2\right) \right]$$
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].

2.3.2 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.

2.4 Test cases

2.4.1 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.

 figure: Fig. 4.

Fig. 4. View of the X-ray test bench.

Download Full Size | PDF

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.

 figure: Fig. 5.

Fig. 5. Sketch of test cases. Figures are not to scale. The PMMA cylinder diameter is 14 mm. The silicon edge is 0.7 mm thick. The silicon edge in test case #1 is shifted by 0.18 mm and tilted by 89.9 degrees.

Download Full Size | PDF

 figure: Fig. 6.

Fig. 6. Simulated source spectrum (dashed red) and detector energy response (solid blue).

Download Full Size | PDF

2.4.2 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.

2.4.3 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.

3. Results

3.1 Test case #1: Refraction

Figure 7(a) shows the experimental image of the refraction protocol. The silicon edge is on the left hand side and the PMMA cylinder in the bottom-right part of the image. 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.

 figure: Fig. 7.

Fig. 7. Test case #1. (a) Acquired image where total reflection on the silicon edge can be seen on the left and refraction on the the PMMA cylinder surface on the right, the sampled profiles are shown as overlaid yellow lines. (b) Simulated standard attenuation image without X-ray refraction implementation. (c) Simulated attenuation image with the new X-ray refraction boundary process.

Download Full Size | PDF

 figure: Fig. 8.

Fig. 8. Attenuation profiles sample from images in Fig. 7. (a) Profiles normal to the silicon edge in the air. (b) Profiles normal to the silicon edge in the cylinder. (c) Profiles normal to the cylinder top surface.

Download Full Size | PDF

3.2 Test case #2: Fresnel diffraction

Figure 9(a) shows the experimental image of the diffraction protocol. Figures 9(b) and 9(c) respectively show the corresponding simulated images without and with the Fresnel diffraction module. The latter image is a mixture of Fresnel diffraction patterns of different photon energy 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.

 figure: Fig. 9.

Fig. 9. Test case #2. (a) Acquired image. (b) Attenuation simulation and (c) Fresnel diffraction simulation using the experimental conditions (X-ray tube and detector).

Download Full Size | PDF

 figure: Fig. 10.

Fig. 10. Fresnel diffraction simulation with a monochromatic source and an ideal detector.

Download Full Size | PDF

 figure: Fig. 11.

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).

Download Full Size | PDF

4. 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 laboratory-based 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.

5. 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.

Funding

ANR (ANR-11-IDEX-0007, ANR-11-LABX-0063).

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.

References

1. 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). [CrossRef]  

2. 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). [CrossRef]  

3. 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). [CrossRef]  

4. A. Olivo, C. Arvanitis, S. Bohndiek, A. Clark, M. Prydderch, R. Turchetta, and R. Speller, “First evidence of phase-contrast imaging with laboratory sources and active pixel sensors,” Nucl. Instrum. Methods Phys. Res., Sect. A 581(3), 776–782 (2007). [CrossRef]  

5. S. Berujon, E. Ziegler, R. Cerbino, and L. Peverini, “Two-dimensional X-ray beam phase sensing,” Phys. Rev. Lett. 108(15), 158102 (2012). [CrossRef]  

6. D. Attwood, Soft X-rays and extreme ultraviolet radiation (Cambridge University Press, 1999).

7. D. Paganin, Coherent X-Ray Optics, Oxford Series on Synchrotron Radiation (Oxford University Press, USA, 2006).

8. A. V. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” J. Opt. Soc. Am. A 19(3), 472–480 (2002). [CrossRef]  

9. M. Langer, P. Cloetens, J.-P. Guigay, and F. Peyrin, “Quantitative comparison of direct phase retrieval algorithms in in-line phase tomography,” Med. Phys. 35(10), 4556–4566 (2008). [CrossRef]  

10. A. Burvall, U. Lundström, P. A. C. Takman, D. H. Larsson, and H. M. Hertz, “Phase retrieval in X-ray phase-contrast imaging suitable for tomography,” Opt. Express 19(11), 10359–10376 (2011). [CrossRef]  

11. J. W. Goodman, Introduction to Fourier optics (Roberts and Co., 2005), 3rd ed.

12. K. Bliznakova, P. Russo, G. Mettivier, H. Requardt, P. Popov, A. Bravin, and I. Buliev, “A software platform for phase contrast X-ray breast imaging research,” Comput. Biol. Med. 61, 62–74 (2015). [CrossRef]  

13. P. Monnin, S. Bulling, J. Hoszowska, J.-F. Valley, R. Meuli, and F. R. Verdun, “Quantitative characterization of edge enhancement in phase contrast X-ray imaging,” Med. Phys. 31(6), 1372–1383 (2004). [CrossRef]  

14. A. Groso, R. Abela, and M. Stampanoni, “Implementation of a fast method for high resolution phase contrast tomography,” Opt. Express 14(18), 8103–8110 (2006). [CrossRef]  

15. T. E. Gureyev, Y. I. Nesterets, A. W. Stevenson, P. R. Miller, A. Pogany, and S. W. Wilkins, “Some simple rules for contrast, signal-to-noise and resolution in in-line X-ray phase-contrast imaging,” Opt. Express 16(5), 3223–3241 (2008). [CrossRef]  

16. E. R. Shanblatt, Y. Sung, R. Gupta, B. J. Nelson, S. Leng, W. S. Graves, and C. H. McCollough, “Forward model for propagation-based X-ray phase contrast imaging in parallel- and cone-beam geometry,” Opt. Express 27(4), 4504 (2019). [CrossRef]  

17. X. Wu and H. Liu, “A new theory of phase-contrast X-ray imaging based on Wigner distributions,” Med. Phys. 31(9), 2378–2384 (2004). [CrossRef]  

18. E. F. Donnelly, R. R. Price, and D. R. Pickens, “Experimental validation of the Wigner distributions theory of phase-contrast imaging,” Med. Phys. 32(4), 928–931 (2005). [CrossRef]  

19. I. V. Bazarov, “Synchrotron radiation representation in phase space,” Phys. Rev. Spec. Top.--Accel. Beams 15(5), 050703 (2012). [CrossRef]  

20. A. Yan, X. Wu, and H. Liu, “A general theory of interference fringes in X-ray phase grating imaging,” Med. Phys. 42(6Part1), 3036–3047 (2015). [CrossRef]  

21. T. Faragó, P. Mikulík, A. Ershov, M. Vogelgesang, D. Hänschke, and T. Baumbach, “syris: a flexible and efficient framework for X-ray imaging experiments simulation,” J. Synchrotron Radiat. 24(6), 1283–1295 (2017). [CrossRef]  

22. J. Vignero, N. W. Marshall, K. Bliznakova, and H. Bosmans, “A hybrid simulation framework for computer simulation and modelling studies of grating-based X-ray phase-contrast images,” Phys. Med. Biol. 63(14), 14NT03 (2018). [CrossRef]  

23. J. Allison, K. Amako, J. Apostolakis, P. Arce, M. Asai, T. Aso, E. Bagli, A. Bagulya, S. Banerjee, G. Barrand, B. Beck, A. Bogdanov, D. Brandt, J. Brown, H. Burkhardt, P. Canal, D. Cano-Ott, S. Chauvie, K. Cho, G. Cirrone, G. Cooperman, M. Cortés-Giraldo, G. Cosmo, G. Cuttone, G. Depaola, L. Desorgher, X. Dong, A. Dotti, V. Elvira, G. Folger, Z. Francis, A. Galoyan, L. Garnier, M. Gayer, K. Genser, V. Grichine, S. Guatelli, P. Guèye, P. Gumplinger, A. Howard, I. Hřivnáčová, S. Hwang, S. Incerti, A. Ivanchenko, V. Ivanchenko, F. Jones, S. Jun, P. Kaitaniemi, N. Karakatsanis, M. Karamitros, M. Kelsey, A. Kimura, T. Koi, H. Kurashige, A. Lechner, S. Lee, F. Longo, M. Maire, D. Mancusi, A. Mantero, E. Mendoza, B. Morgan, K. Murakami, T. Nikitina, L. Pandola, P. Paprocki, J. Perl, I. Petrović, M. Pia, W. Pokorski, J. Quesada, M. Raine, M. Reis, A. Ribon, A. R. Fira, F. Romano, G. Russo, G. Santin, T. Sasaki, D. Sawkey, J. Shin, I. Strakovsky, A. Taborda, S. Tanaka, B. Tomé, T. Toshito, H. Tran, P. Truscott, L. Urban, V. Uzhinsky, J. Verbeke, M. Verderi, B. Wendt, H. Wenzel, D. Wright, D. Wright, T. Yamashita, J. Yarba, and H. Yoshida, “Recent developments in Geant4,” Nucl. Instrum. Methods Phys. Res., Sect. A 835, 186–225 (2016). [CrossRef]  

24. A. G. Tur’yanskii, O. V. Konovalov, S. S. Gizha, and N. D. Beilin, “Edge diffraction effect at the refraction of x rays in a diamond prism,” JETP Lett. 100(8), 540–542 (2014). [CrossRef]  

25. D. M. Cunha, A. Tomal, and M. E. Poletti, “Evaluation of scatter-to-primary ratio, grid performance and normalized average glandular dose in mammography by Monte Carlo simulation including interference and energy broadening effects,” Phys. Med. Biol. 55(15), 4335–4359 (2010). [CrossRef]  

26. J. Ulseth, Z. Zhu, Y. Sun, and S. Pang, “Accelerated X-ray diffraction (tensor) tomography simulation using OptiX GPU ray-tracing engine,” IEEE Trans. Nucl. Sci. 66(12), 2347–2354 (2019). [CrossRef]  

27. 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). [CrossRef]  

28. Z. Wang, Z. Huang, L. Zhang, Z. Chen, and K. Kang, “Implement X-ray refraction effect in Geant4 for phase contrast imaging,” in Proc. IEEE Nuclear Science Symp. Conf. Record (NSS/MIC), (2009), pp. 2395–2398.

29. P. Bartl, J. Durst, W. Haas, T. Michel, A. Ritter, T. Weber, and G. Anton, “Simulation of X-ray phase-contrast imaging using grating-interferometry,” in 2009 IEEE Nuclear Science Symposium Conference Record (NSS/MIC), (IEEE, 2009).

30. S. Peter, P. Modregger, M. K. Fix, W. Volken, D. Frei, P. Manser, and M. Stampanoni, “Combining Monte Carlo methods with coherent wave optics for the simulation of phase-sensitive X-ray imaging,” J. Synchrotron Radiat. 21(3), 613–622 (2014). [CrossRef]  

31. J. Sanctorum, J. De Beenhouwer, and J. Sijbers, “X-ray phase-contrast simulations of fibrous phantoms using GATE,” in 2018 IEEE Nuclear Science Symposium and Medical Imaging Conference Proceedings (NSS/MIC), (IEEE, 2018).

32. J. Sanctorum, J. De Beenhouwer, J. Weissenböck, C. Heinzl, J. Kastner, and J. Sijbers, “Simulated grating-based X-ray phase contrast images of CFRP-like objects,” in 9th Conference on Industrial Computed Tomography, Padova, Italy (2019).

33. S. Cipiccia, F. A. Vittoria, M. Weikum, A. Olivo, and D. A. Jaroszynski, “Inclusion of coherence in Monte Carlo models for simulation of X-ray phase contrast imaging,” Opt. Express 22(19), 23480 (2014). [CrossRef]  

34. D. Sarrut, M. Bardiès, N. Boussion, N. Freud, S. Jan, J.-M. Létang, G. Loudos, L. Maigne, S. Marcatili, T. Mauxion, P. Papadimitroulas, Y. Perrot, U. Pietrzyk, C. Robert, D. Schaart, D. Visvikis, and I. Buvat, “A review of the use and potential of the GATE Monte Carlo simulation code for radiation therapy and dosimetry applications,” Med. Phys. 41(6Part1), 064301 (2014). [CrossRef]  

35. M. Magnuson and C. F. Hague, “Refractive index of vanadium determined by resonant diffraction of soft X-rays,” J. Electron Spectrosc. Relat. Phenom. 137-140, 519–522 (2004). [CrossRef]  

36. E. M. Gullikson, “Optical properties of materials,” in Experimental Methods in the Physical Sciences, (Elsevier, 1998), chap. 13, pp. 257–270.

37. ICRP, Report No 23: Report on the Task Group on Reference Man (Pergamon Press, 1974).

38. M. J. Berger, J. S. Coursey, M. A. Zucker, and J. Chang, “Stopping-powers and range tables for electrons, protons, and helium ions,” http://physics.nist.gov/Star [2020, April 9], National Institute of Standards and Technology, Standard Reference Database 124 (2017).

39. T. Schoonjans, A. Brunetti, G. B M. Sanchez del Rio, V. A. Solé, C. Ferrero, and L. Vincze, “The xraylib library for X-ray – matter interactions. recent developments,” Spectrochim. Acta, Part B 66(11-12), 776–784 (2011). [CrossRef]  

40. B. L. Henke, E. M. Gullikson, and J. C. Davis, “X-ray interactions: Photoabsorption, scattering, transmission, and reflection at E = 50-30,000 eV, Z = 1-92,” At. Data Nucl. Data Tables 54(2), 181–342 (1993). [CrossRef]  

41. J. Als-Nielsen and D. Mcmorrow, Elements of Modern X-ray Physics (John Wiley & Sons, 2011).

42. D. M. Paganin and D. Pelliccia, “Tutorials on X-ray phase contrast imaging: Some fundamentals and some conjectures on future developments,” arXiv preprint arXiv:1902.00364 (2019).

43. N. Siegel, J. Rosen, and G. Brooker, “Faithful reconstruction of digital holograms captured by FINCH using a Hamming window function in the Fresnel propagation,” Opt. Lett. 38(19), 3922 (2013). [CrossRef]  

44. G. Poludniowski, P. M. Evans, V. N. Hansen, and S. Webb, “An efficient Monte Carlo-based algorithm for scatter correction in keV cone-beam CT,” Phys. Med. Biol. 54(12), 3847–3864 (2009). [CrossRef]  

45. S. Rit, M. V. Oliva, S. Brousmiche, R. Labarbe, D. Sarrut, and G. C. Sharp, “The reconstruction toolkit (rtk), an open-source cone-beam CT reconstruction toolkit based on the insight toolkit (ITK),” J. Phys.: Conf. Ser. 489, 012079 (2014). [CrossRef]  

46. F. H. Attix, Introduction to Radiological Physics and Radiation Dosimetry (Wiley-VCH, 2004).

47. A. Peterzol, J. Berthier, P. Duvauchelle, and F. C D. Babot, “X-ray phase contrast image simulation,” Nucl. Instrum. Methods Phys. Res., Sect. B 254(2), 307–318 (2007). [CrossRef]  

48. N. Freud, J.-M. Letang, and D. Babot, “A hybrid approach to simulate X-ray imaging techniques, combining Monte Carlo and deterministic algorithms,” IEEE Trans. Nucl. Sci. 52(5), 1329–1334 (2005). [CrossRef]  

49. S. Rit, “Fresnel_FFD & XRayRefraction,” https://github.com/OpenGATE/GateContrib/tree/master/imaging/CT [2020, April 11] OpenGATE / GateContrib user-oriented public repository (2018).

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (11)

Fig. 1.
Fig. 1. Ratio of the refractive index decrement $\delta _E$ over the attenuation index $\beta _E$ for cortical bone and soft tissue (see ICRP report [37]) as a function of the energy. Computed via xraylib [39].
Fig. 2.
Fig. 2. Air-to-silicon interface: reflectance and transmittance of an unpolarized $10$ keV X-ray wave (see Eq. (4) with $\delta _1 = 0$ and $\delta _2 = 5\times 10^{-6}$ ).
Fig. 3.
Fig. 3. Fresnel scaling theorem.
Fig. 4.
Fig. 4. View of the X-ray test bench.
Fig. 5.
Fig. 5. Sketch of test cases. Figures are not to scale. The PMMA cylinder diameter is 14 mm. The silicon edge is 0.7 mm thick. The silicon edge in test case #1 is shifted by 0.18 mm and tilted by 89.9 degrees.
Fig. 6.
Fig. 6. Simulated source spectrum (dashed red) and detector energy response (solid blue).
Fig. 7.
Fig. 7. Test case #1. (a) Acquired image where total reflection on the silicon edge can be seen on the left and refraction on the the PMMA cylinder surface on the right, the sampled profiles are shown as overlaid yellow lines. (b) Simulated standard attenuation image without X-ray refraction implementation. (c) Simulated attenuation image with the new X-ray refraction boundary process.
Fig. 8.
Fig. 8. Attenuation profiles sample from images in Fig. 7. (a) Profiles normal to the silicon edge in the air. (b) Profiles normal to the silicon edge in the cylinder. (c) Profiles normal to the cylinder top surface.
Fig. 9.
Fig. 9. Test case #2. (a) Acquired image. (b) Attenuation simulation and (c) Fresnel diffraction simulation using the experimental conditions (X-ray tube and detector).
Fig. 10.
Fig. 10. Fresnel diffraction simulation with a monochromatic source and an ideal detector.
Fig. 11.
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).

Equations (14)

Equations on this page are rendered with MathJax. Learn more.

n E ( p ) = 1 δ E ( p ) + i β E ( p )
n E ( p ) = 1 r e 2 π ( h c E ) 2 k K a k ( p ) ( f k , 1 ( E ) i f k , 2 ( E ) )
( 1 δ 1 ) sin θ 1 = ( 1 δ 2 ) sin θ 2 .
R F ( θ 1 ) = 1 2 [ | ( 1 δ 1 ) cos θ 1 ( 1 δ 2 ) cos θ 2 ( 1 δ 1 ) cos θ 1 + ( 1 δ 2 ) cos θ 2 | 2 + | ( 1 δ 1 ) cos θ 2 ( 1 δ 2 ) cos θ 1 ( 1 δ 1 ) cos θ 2 + ( 1 δ 2 ) cos θ 1 | 2 ]
θ 1 c = asin ( 1 δ 2 1 δ 1 )
Ψ E ( x , y , z = z 0 ) exp { i π E h c z = 0 z = z 0 [ 1 n E 2 ( p ) ] d z } Ψ E ( x , y , z = 0 ) .
1 n E 2 ( p ) 2 [ 1 n E ( p ) ]
T E ( x , y ) = exp [ B E ( p ) ] exp [ i Φ E ( p ) ]
B E ( x , y , z = z 0 ) = 2 π E h c z = 0 z = z 0 β E ( p ) d z
Φ E ( x , y , z = z 0 ) = 2 π E h c z = 0 z = z 0 δ E ( p ) d z .
I E ( x , y , z = z 0 ) = | Ψ ( x , y , z = z 0 ) | 2 exp [ 2 B E ( p ) ] I E ( x , y , z = 0 )
μ E ( p ) = 4 π E h c β E ( p ) .
P ^ E , D ( f x , f y ) = exp [ i π h c E D ( f x 2 + f y 2 ) ]
P ^ E , D ( f x , f y ) = exp [ i π h c E D M ( f x 2 + f y 2 ) ]
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.