Computational simulations and assessment of two approaches for x-ray phase contrast imaging

X-ray phase-contrast imaging is a high-resolution imaging that permits an increase of the perceptibility of the details in three-dimensional objects, such as human tissues compared to conventional absorption imaging. There are different approaches for implementing phase-contrast imaging and their introduction into clinical practice requires advanced computational tools. A long-term goal of our research is the development of computational models of breast phase-contrast imaging. The aim of this study is to develop a software module for implementing grating-based phase-contrast imaging. For this purpose, an existing in-house software application for x-ray imaging with a function to model and simulate propagation-based phase-contrast x-ray images has been extended to include a model of grating-based imaging. To test the new functionality, four computational phantoms reflecting features, which can be screened in the real breast tissue and which differ in their complexity, were designed. Planar x-ray images in absorption, propagation-based and grating-based modes were generated and compared. Results showed improved visual appearance of the simulated objects in images obtained by simulating grating-based imaging setup. The developed subroutine is planned to be experimentally validated at synchrotron facility. The new software functionality will be exploited in studies related to new x-ray imaging techniques for breast screening and diagnosing.


Introduction
Breast cancer is the most common incident and mortality form of cancer for women. From the estimated 3.91 million new cases of cancer in Europe in 2018, 28.2% are due to breast cancer, while from the estimated 1.93 million deaths from cancer, 16.2% are caused by breast cancer [1]. Nowadays, mammography continues to be the 'gold technique' in screening and diagnosing the breast. Mammography images present a 2D distribution of the intensity of the x-rays, which are transmitted through the woman breast and then reach the detector. The technique is based on the x-ray absorption properties of the tissues in the breast. In case of dense breasts, however, tumour tissues could be hidden by the dense breast tissue; thus the mammography technique in this particular situation may result in superposition of large areas of the glandular tissue and the cancerous ones on the 2D detector, which will influence the diagnose of these structures [2]. Therefore, new and more efficient x-ray imaging techniques dedicated to screening and diagnosing the breast tissue, such as phase-contrast mammography [3,4], breast computed tomography [5,6] and contrast-enhanced mammography [7], are continuously under development, extensive testing and optimisation.
Phase-contrast mammography is an emerging alternative approach to the absorption based mammography. The base of the x-ray phase-contrast imaging is both, (a) the attenuation of the x-rays, and (b) the change in the phase of the x-rays, which is related to diffraction and refraction effects when x-rays propagate in the breast tissue [8,9]. The amount of phase change in a tissue is associated with the local value of its wavelength-dependent refractive index n, expressed as IOP Publishing doi:10.1088/1742-6596/2162/1/012013 2 n(λ) = 1-δ + ιβ (1) where n(λ) is the complex refractive index of the tissue at wavelength λ, β is the absorption index, δ is the decrement (with respect to unity) of the real part of the refractive index n and ι 2 = -1 [8]. The coefficient β reflects the x-ray attenuation in a tissue thickness s, while the coefficient δ > 0 determines the total phase change integrated along the direction of propagation, z, of the x-ray beam.
During the years, different techniques have been proposed for detecting phase-contrast effects: propagation-based approach, diffraction-enhanced imaging, crystal interferometry, grating interferometry, the edge-illumination approach, and speckle-tracking [9,10]. Amongst the popular approaches that have been implemented for clinical work, is the so-called in-line or propagation-based phase-contrast imaging approach that already have clinical implementation for breast cancer imaging [11]. In this approach, there are no optical elements, such as gratings and apertures, placed in the beam path during acquisition of x-ray images. In this mode, intensity changes in the Fresnel near field at the surface of a fine-pitch area detector positioned at a distance from the object (as shown in figure 1a), are related to the transverse Laplacian of the phase changes in the object. With this arrangement, the visual appearance of phase-contrast effects in the final x-ray image presents a characteristic contrast enhancement at the interfaces between structures with different x-ray refractive indices [12]. This mode has been computationally modelled and used in successful simulations of phase-contrast mammography and breast tomosynthesis by our team [12][13][14][15][16]. In addition, this approach was implemented at ELETTRA synchrotron facility via the mammography SYRMEP project [17,18]. One of the limitations of the propagation based approach is the clinical implementation, which is done only at difficult to access synchrotron facilities. An encouraging approach is developing of grating interferometry (with a typical design shown in figure 1b), which may turn to be applicable for imaging of the whole breast. This approach is implemented by using phase and analyser gratings. The most flexible approach to carry out development and optimisation work is the approach of modelling the components of the imaging system and subsequently imaging simulation.
There is a limited number of computational studies in the field of phase-contrast grating imaging reported in the literature [19][20][21][22][23][24]. All these studies are based on the use of simple objects, such as spheres and cylinders, and aim at proof of concept and optimisation of dose and imaging parameters. To the best of our knowledge, studies in the field of mammography have not been reported. To construct a model of a breast imaging technique based on gratings and advance in this field, a realistic three-dimensional (3D) model of the female breast with accurately modelled anatomical details and a computational model of the grating-based imaging itself, are needed.
A long-term goal of our research is the development of computational models of phase-contrast imaging for screening and diagnosing the breast. The research is entirely original, aiming at a complete computational platform with anthropomorphic breast models and modelled different approaches for phase-contrast imaging allowing comparison of techniques and further optimisation. Our team has the expertise in creating realistic computational breast phantoms [25,26], which turns out to be the main tool in studying novel imaging techniques, amongst which is also the phase-contrast imaging. The propagation-based phase-contrast approach was already successfully modelled for mammography 2D and 3D imaging [12][13][14][15][16]. The model of the phase-contrast based on gratings is still to be developed. The first step in achieving this goal is the development of a dedicated software, capable of modelling the phase-contrast grating imaging. The aim of this study is to develop a software module for implementing grating-based x-ray phase-contrast imaging and to assess the results of applying this technique in comparison to both, absorption and propagation-based x-ray imaging approaches.  Figure 1 shows the simulated phase-contrast setups for generation of x-ray images. The propagationbased mode, shown in figure 1a, is modelled based on the Fresnel-Kirchhoff diffraction theory [8] and is described in details in [12]. The model considers pure coherent radiation of wavelength λ. The xrays are emerging from a point-like monochromatic source, located at a distance SID (source-toisocenter distance) from the object plane. The object is modelled in respect to the isocenter point, while the detector is located at a distance SDD (source-to-detector distance) from the source. The propagation-based model was validated to correctly reproduce phase-contrast images against images of simple and complex physical phantoms obtained from experimental work carried out at the European Synchrotron Radiation Facility (ESRF), Grenoble, France. In propagation mode, the detector is placed away from the object in order to register the phase shift, while no optics are used. The modelled grating-based phase-contrast setup is demonstrated in figure 1b. It is based on a coherent x-ray source that emits an x-ray beam towards the detector [27,28]. The first grating, G1, is the phase grating with a phase shift of π that generates an interference pattern, known as Talbot carpet registered at a distance (3)

Simulated setup
In this equation, g1 is the grating pitch, while m is a positive integer value. The second grating, G2, is the analyser grating, placed just in front of the detector at a distance dm from the phase grating G1. The pitch g2 of the analyser grating matches the periodicity of the undistorted interference pattern at the corresponding detector position.
In the computational implementation of the grating-based phase-contrast imaging, G2 is always placed at a distance at which the Talbot carpet is obtained. The main steps in implementing this imaging are shown in figure 1c. The plane wave is propagated through the object and the first grating G1 by applying the x-ray tracing algorithm [12] followed by propagation of the field, implemented by using the Fresnel diffraction integral. Finally, the propagation of the wave through the G2 is performed. For the simulations performed in this study, the phase grid G1 was designed to have a pitch of 2.5 μm, while the analyser grid G2 had a pitch of 1.25 μm.

Computational phantoms
The computational phantoms were designed with the in-house developed software application XRayImagingSimulator [12]. This simulator has three main modules: a module for creation of solidbased 3D phantoms, a module that creates x-ray images, when simulated x-rays with a given energy (in keV), following transmission through the computational phantoms, reach the detector, and a module for visualization of both, created computational phantoms and generated x-ray images.

Figure 2.
Computational phantoms used in the simulations: a) a nylon wire with a diameter of 2 mm; b) two water spheres located at different z heights and with diameters of 2 mm and 1 mm, respectively; c) a PMMA container with a thickness of 40 mm and four air cylinders with a diameter of 1 mm, d) a paraffin container with a thickness of 5 mm and an area 10 mm x 10 mm and 473 water spheres with a diameter of 1 mm.
Four computational phantoms, shown in figure 2, were created in order to test and validate the software code. The first phantom, shown in figure 2a, consists of a nylon wire 2 mm in a diameter placed in the air. The second phantom, shown in figure 2b, is represented by two water spheres with diameters of 2 mm and 1 mm placed at different positions in air, while the third phantom contains four cylinders filled with air, having diameters of 1 mm, arranged at different planes, as shown in figure 2c. Finally, the complex phantom, phantom 4 shown in figure 2d, is represented by a paraffin container filled with 473 water spheres with a diameter of 1 mm. The phantoms in figure 2 are shown together with the coordinate system of the modelled x-ray imaging geometry, depicted in figure 1. The center of the x-ray imaging system is the isocenter of the system. In +z direction, a point-like x-ray source is modelled. In -z direction, the flat panel detector is modelled as a parallelepiped, absorbing all the xrays reaching its surface.
Three acquisition setups were simulated by implementing planar acquisition x-ray geometry, i.e. one x-ray image is generated. The acquisition geometry is described with the two distances SID and  figure 1), which were set to 1000 mm and 1300 mm, respectively. The high-resolution projection images are of a pixel size of 5 µm and an image size between 1500 pixels x 1500 pixels and 12000 pixels x 12000 pixels, depending on the computational phantom size. Incident x-ray beam energy was set to 20 keV. Detector noise was not considered in this simulation. Registered x-ray images present line integrals. The propagation of the x-rays until the G1 are based on the Fresnel-Kirchhoff diffraction theory, computationally described in [12]. The propagation through the gratings G1 and G2 is implemented by the approach discussed in section 2.1. Figure 3 shows side by side the simulated x-ray projection images of the first three computational phantoms for the simulated absorption, propagation-based and phase-contrast x-ray imaging acquisition geometries. Although simple in design, these phantoms are of imaging interest, i.e. optimisation of acquisition distances, x-ray energy, tissue composition of objects, gratings material. These computational objects, modelled in with shapes and sizes, as shown in figure 2, could approximate well normal and abnormal tissues in real breast. Rows 1 and 2 correspond to the first computational phantom (the nylon wire, shown in figure 2a), rows 3 and 4 correspond to the second computational phantom (the two water spheres, shown in figure 2b), while rows 5 and 6 correspond to the third computational phantom (the PMMA container with air channels, shown in figure 2c). For each phantom, there are three x-ray projection images simulated in absorption mode (figure 3a), propagation mode (figure 3b) and grating-based approach (figure 3c). The absorption mode is implemented by applying the geometry, shown in figure 1a with the detector, placed behind the object.

Results and discussions
The plotted profiles, shown below the corresponding x-ray projection images, are taken through the center of the projection images. The first column, figure 3a, shows the simulated projection images when the detector was placed just behind the phantoms (i.e. no phase-contrast effects are observed), while second and third columns, figure 3b,c, show the simulated projection images when the detector was placed at 300 mm away from the phantoms. For all phantoms, the visual comparison of the projection images shows that both, the propagation-based and the grating-based acquisition geometries result in improved object edges compared to the absorption image (first column in figure 3). This observation is well supported by the comparison of the line profiles, which show quantitatively the edge effect due to the phase effect. Specifically, the objective evaluation, which included calculation of the edge enhancement in phase-contrast images (figure 3b,c) by using edge-enhancement index, EEI [29]: where P and T are the peak and trough intensity values at the edge in the plot graphs in figure 3b,c, while H and L are values taken in the high-and low-intensity regions next to the edge [29]. The results for the EEI shows that for the nylon wire the edge enhancement due to the phase effect is approximately 5 times higher for grating-based acquisition geometry than for propagation-based imaging. Similarly, the edge enhancement for the second and the third phantoms, i.e. the water spheres and the PMMA container with four air channels, in case of grating-based acquisition imaging were approximately 4 and 3 times, respectively compared to the propagation-based approach, demonstrated in figure 3b. Finally, the resulting phase-contrast images with gratings visually match with these reported in the study of Fu et al [23], based on numerical calculations for a phantom composed of spheres and one cylinder. To speed the process of producing the images, the software code for producing x-ray projection images based on grating approach was rewritten for parallel programming using the omp library. While for the third phantom, the time needed to obtain an image with resolution of 5 μm and size of 1500 pixels x 1500 pixels was 5 hours, the image obtained with the parallel code on a PC with 24 threads was only few seconds.
Further, this edge-enhancement is well demonstrated in the simulated phase-contrast images from the texture phantom, shown in figure 2d. This complex phantom is initially used to test the newly developed software module to produce grating-based phase-contrast images and to evaluate the visibility of the texture on these images.  The visual comparison of x-ray projection images showed well differentiated projected spheres on the detector plane for the grating-based acquisition geometry in respect to the other two simulated xray imaging acquisitions. This improvement is well demonstrated by the peaks in the plotted profiles, shown below the corresponding x-ray projection images, the latter taken though the center of the projection images. The calculated EEI for the spherical edges are 2 to 3 times higher in grating-based images compared to propagation-based images. The result from this experiment showed that texture is well visualised with apparent edge enhancement.
This simulation study showed that both, phase-contrast techniques result in improved edge enhancement effect, as the images obtained through the grating-based approach are superior in comparison to the images obtained in propagation-based setup. These simulation results are to be experimentally validated at synchrotron facility. The experimental procedure requires: (a) physical models of the computational phantoms, shown in figure 2, and (b) experimental acquisition of x-ray images of these physical model at synchrotron facility with the same acquisition setup, shown in figure  1. For the creation of physical models, the 3D printing facilities in our laboratory will be used [30,31]. They will be placed on a hardware platform at a synchrotron beamline, providing x-rays, suitable for mammography applications. Images are then acquired at different distances and with the addition of gratings. Once the results of the proposed computational approach are validated, it will be used in investigations with anthropomorphic computational breast models.

Conclusions
This computational study demonstrated the development of a software module for generating x-ray grating-based phase-contrast images and was tested to generate phase-contrast images in comparison to x-ray simulated images, obtained in absorption and propagation-based acquisition geometries. The simulation results showed that the module can produce phase-contrast images in grating-based acquisition geometry. When compared to absorption and propagation-based geometries, the gratingbased x-ray phase-contrast images demonstrate higher edge enhancement for all computational phantoms used in the study, even for the phantom containing highly structured background. The module will allow comparison of different phase-contrast techniques, applied for detection of various types of lesions, which can be found in the breast tissue. A necessary step to be undertaken is the experimental validation of this software module, which validation could be easily performed at a synchrotron facility. Once validated, it will be used in breast imaging studies with anthropomorphic physical and computational breast models.