A hybrid simulation framework for computer simulation and modelling studies of grating-based x-ray phase-contrast images

Clinical studies performed using computer simulation are inexpensive, flexible methods that can be used to study aspects of a proposed imaging technique prior to a full clinical study. Typically, lesions are simulated into (experimental) data to assess the clinical potential of new methods or algorithms. In grating-based phase-contrast imaging (GB-PCI), full wave simulations are, however, computationally expensive due to the high periodicity of the gratings and therefore not practically applicable when large data sets are required. This work describes the development of a hybrid modelling platform that combines analytical and empirical input data for a more rapid simulation of GB-PCI images with little loss of accuracy. Instead of an explicit implementation of grating details, measured summary metrics (i.e. visibility, flux, noise power spectra, presampling modulation transfer function) are applied in order to generate transmission and differential phase images with large fields of view. Realistic transmission and differential phase images were obtained with good quantitative accuracy. The different steps of the simulation framework, as well as the methods to measure the summary metrics, are discussed in detail such that the technique can be easily customized for a given system. The platform offers a fast, accurate alternative to full wave simulations when the focus switches from grating/system design and set up to the generation of GB-PCI images for an established system.

is exploited frequently (Young et al 2013, Solomon andSamei 2014), however, in GB-PCI, as far as we know, no such studies have yet been performed. A potential hurdle might be that the simulation framework required, has to produce large data sets of realistic images that have large fields of view, all to be done relatively quickly. Simulation frameworks of GB-PCI images described in the literature aim to simulate the imaging setup in high detail. Most frameworks are based on wave-optical numerical analysis (Weitkamp 2004, Malecki, Potdevin and while some combine wave simulations with Monte Carlo frameworks in order to include x-ray scattering (Bartl et al 2009, Peter et al 2014. Others use computationally expensive Monte Carlo methods, where refraction is included via ray tracing (Wang et al 2009) or where diffraction effects based on the Huygens-Fresnel principle are included (Cipiccia et al 2014). While these methods offer an accurate and detailed simulation of signal transfer, the complexity of these frameworks make them computationally expensive with limited applicability for simulation of simulation studies of clinical tasks.
This work is based on the hypothesis that full system modelling may not always be required. When many technical details of the system have been established and/or when a prototype system is being produced, the focus may shift to the study of specific (imaging) tasks. The topic of this study is, therefore, not detailed modelling of system components and their influence on image quality, but rather the application of summary data acquired on a prototype GB-PCI system to generate realistic images. To this end, a hybrid image modelling technique is introduced that combines theoretical equations with empirical data to simulate transmission and differential phase images. Dark field images are not included, since the lack of data on the microscale structure of biological tissues responsible for the dark field signals is limited. In this work, three imaging tasks of differing complexities are simulated with a hybrid, fast simulation framework. The results are then compared to experimental data in terms of structural similarity and contrast.

Materials and methods
When x-rays propagate through matter, both the phase and amplitude are influenced by the material specific coefficients δ(λ) and β(λ) respectively, which are related to the complex refractive index n = 1 − δ + iβ. For a given wavelength, the x-ray wave function ψ 0 after transmission through a material with thickness t is given by ψ: in which k = 2π/λ is the wavenumber. A changing phase component, perpendicular to the beam propagation direction ( x), causes the beam to refract under a specific angle α (Bravin 2013): In the case of x-rays, typical values for α lie in the µrad range and in order to extract such small effects dedicated methods have been developed. In GB-PCI, the object to be imaged is placed in a region containing a high frequency intensity pattern. Comparison of the object-disturbed pattern with the original pattern enables the construction of transmission (Tr), differential phase (dP) and dark field (DF) images (David et al 2002, Momose et al 2003, Pfeiffer et al 2006. Three gratings are used to create and measure the intensity pattern (see figure 1). The first grating G 0 creates individually coherent beams, which constructively interfere at the detector plane. The G 1 grating induces a periodic phase shift within the x-ray wave front, creating a high intensity fringe pattern at a fractional Talbot distance d. The pitch ( p 2 ) of the third grating G 2 matches the intensity pattern period and by stepping this grating, the pattern can be sampled and reconstructed. The amplitude to mean intensity of the fringe pattern is the so-called 'visibility' of the system and determines the noise in the dP and DF image. A detailed description of the technology can be found in Pfeiffer et al (2006).
As the technology relies on wave interference effects, most simulation platforms are based on numerical wave propagation methods. Such methods are flexible and suitable for exploring GB-PCI system geometry and set up, but are demanding in terms of memory usage and simulation time as the high grating pitch, typically between 2 and 5 µm, requires high spatial sampling of the wave matrix which is computationally expensive for realistic fields of view. Instead of modelling the wave interactions, it is also possible to compute the expected signal based on theory and this approach forms the basis of our method. For transmission imaging, the expected signal (A) can be predicted from the Beer-Lambert law: where t is the thickness of the object. For the differential phase signal (ϕ), a similar calculation can be made. The intensity shift at G 2 due to refraction angle α is tan (α) · d, using equation (2): where S is the signal loss due to the offset distance (d 0 ) between the object and G 1 ; S = 1 − d0 L (Engelhardt et al 2007, Chabior et al 2012.

The simulation framework
The hybrid framework developed here combines theoretical equations with experimentally measured data to generate realistic images. First, the phase and attenuation properties of the projected object (βt and δt ) are calculated. Each material present in the object is modelled as a grid of thicknesses (T i (x, y)) with a maximum object thickness T z . The parameters x and y denote the pixel coordinates. A summation over all material constituents results in a phase matrix (∆) and an attenuation matrix (B), for example: where bg indicates the background material. A similar expression can be derived for B using β and the same thickness function T i . Secondly, the finite source size is accounted for by convolving the object with a Gaussian filter (G(u, v)) with a standard deviation equal to the scaled focal spot blur. Blurring due to the finite source size and the x-ray detector is included via a 2D presampling modulation transfer function (MTF) obtained from x-ray detector characterization data (see section 2.2.1). A standard Fast Fourier Transform (FFT) technique is used to apply the MTF. In combination with equations (3) and (4): Lastly, noise ( N i , where i indexes the image type) of the correct texture is added, using the approach of Saunders and Samei (2003): In this method, first, a matrix of normally distributed random values (R) with zero mean and unit variance is generated and then filtered by the square root of the noise power spectrum (NPS): The NPS was measured as described in section 2.3.1. We assume that the x-ray detector remains quantum noise limited and therefore the shape of the NPS remains invariant as exposure level at the x-ray detector changes. The correct magnitude of noise in the final noise image ( N i (x, y)) is obtained by scaling the noise image N 0 i . To ensure a unit variance in N 0 i , the noise image is first scaled with a factor 1/σ Ni , with σ Ni the average standard deviation The variance of the noise is calculated for given system parameters and entrance air kerma applied using the method supplied by Revol et al (2010) or by Engel et al (2011). The level of quantum noise (σ i ) depends on the image type and is given by: where A is the transmission signal in each pixel, therefore the flux is modulated depending on attenuation in the object. The dark field induced coherence loss (V) is generally not known and therefore set to one. N G2 is the number of steps used for an acquisition, v is the visibility, which was set to the average experimentally measured value in a central region of interest. The coefficient f r is a detector specific constant that corrects for deviations in the Poisson nature of the detector response (Revol et al 2010). The dose factor a 0 is the average pixel value in all projections and can be calculated by a 0 = a K · K, where a K is a calibration factor and K the entrance air kerma applied. The calibration factor a K is the average pixel value generated per unit air kerma and is experimentally determined as described in section 2.2.1. In experimental setups there will always be some beam divergence, magnifying the object. The magnification of the object is then: where d 0 is some additional offset distance due to the finite object to G 1 distance. This beam divergence is included by a simple approach of rescaling the pixel spacing by a factor 1/M 0 in the x and y directions. As the flux was scaled to experimental data, this does not influence the total signal in the pixel or the variance of the noise.

The experimental setup
Experimental data for the model development were acquired using a grating-based phase contrast system designed by Carestream Health and installed in the molecular Small Animal Imaging Centre of the University of Leuven (moSAIC). The source is a Varian x-ray tube (Varian Medical Systems, California, USA) operated at 40 kVp for a system design energy of 27.7 keV. The gratings (G 0 ,G 1 ,G 2 ) have nominal periods of respectively 73 µm, 3.901 µm and 2 µm and duty cycles of 0.3, 0.5 and 0.5. The distance L between G 0 and G 1 is 1.7 m and the distance between G 1 and G 2 (d) is 4.35 cm. Objects are centred 5.7 cm in front of G 1 . This corresponds to a magnification factor of 1.06 and a sensitivity correction S of 0.966. The CsI/CMOS detector (Teledyne-Dalsa Xineos 1313, Eindhoven, Netherlands) has 100 × 100 µm 2 pixels and f r calculated to be 1.127. In the central region of interest, the system has an average visibility of 0.22. The nominal focal spot size of the x-ray tube was 0.6 mm, corresponding to a scaled focal spot blur of 23 µm.

System characterization metrics to tune the framework
To construct Tr and dP images with the proposed framework, in addition to system dimension, the following data from the experimental setup should be obtained: the visibility, the MTF, the NPS and the calibration factor a K . Methods to derive the data accurately are elaborated next.

2.2.1.1.The MTF
The detector presampling MTF was determined separately for the transmission and differential phase images. The 1D detector presampling MTF was measured for the grating stepping direction in the image using an implementation of the well-established oversampling technique (Fujita et al 1992, Samei et al 1998, Buhr and Neitzel 2003, Marshall 2006) via an edge test object. The 2D distribution was then formed assuming an isotropic MTF, which is a good approximation following (Wells and Dobbins 2012). A test object with a sharp, straight edge was positioned in front of G 1 and the edge rotated slightly (between 1° and 4°) to the x axis. For the Tr image a 0.8 mm thick steel plate with a polished straight edge was used while for the dP image the test object was a sharp PMMA cube rotated 45°. From each image, a super-sampled edge spread function (ESF) was generated, which was then differentiated to give the line spread function (LSF). Applying an FFT to the LSF followed by normalization to the zero bin value gave the MTF; the toolkit developed by Zhang et al (2017) was used for this. Instead of exporting the MTF directly, the MTF was fitted with a polynomial of degree 6 and the fitting parameters were exported.

2.2.1.2.The NPS
The NPS was calculated from uniformly exposed (flood) Tr and dP images acquired at an entrance air kerma of 0.608 mGy using a standard algorithm as described in Marshall (2006) and executed using MATLAB. The ensemble consisted of 27 images for which a central region of interest (512 × 512) was selected. The images were subdivided in overlaying regions of interest (ROI) of 128 by 128 pixels. Any trend in the ROIs was removed by subtracting a fitted 2D surface to the analysis region. The NPS ROI in the ROI is given by |F {ROI}| 2 and the overall NPS is found from the ensemble average of all NPS ROI data. Still, this NPS contains noise and is not entirely isotropic. Therefore the logarithm of NPS as a function of radial frequency was fitted with a polynomial of degree 2 and the fit parameters exported. Only the shape of the NPS was used in the simulations as this ensures the final image has the correct texture (Saunders and Samei 2003). The absolute magnitude of the noise is calculated in a different way (see equations (12) and (13)). In the image reconstruction of GB-PCI, pixel-to-pixel comparisons between reference and object scans eliminate the contribution of structural noise, which is therefore neglected in this analysis. Quantum noise is the only source of noise considered in this work.

2.2.1.3.The calibration factor a K
Instead of modelling all the successive steps of the flux-from tube output, intensity loss in the gratings, up to the detector efficiency-the signal level in each pixel was scaled with an experimentally derived pixel value as a function of entrance air kerma. The entrance air kerma (K) was measured with a calibrated R100 dose probe of RTI Piranha (RTI Electronics AB, Molndal, Sweden) positioned at d 0 , the location of the object (see figure 1). The calibration coefficient a K was calculated as the mean pixel value per K and for our setup a K equalled 2.45 × 10 5 per mGy at d 0 .

Comparison of simulation framework to experimental data
Accurate sample modelling is also an important aspect of image simulation studies. Analytical models can be implemented with precision, however relevant models of sufficient complexity and realism have not been described. Therefore in addition to a simple analytical model (section 2.3.1), two voxel models based on µCT data were considered when comparing simulation data to experimental results. For the analytical model, simulated images were compared to experimental data using the structure similarity index (SSIM) described in Wang et al (2004). Pixel values for each image of the same type were rescaled to range between 0 and 4000, while window and the Kw constant were set to the default value (Wang et al 2004). For the voxel models, where perfect registration cannot be achieved, this method is not applicable and instead maximum contrasts in regions of interest were compared as the difference in maximum and minimum value of averaged ROIs in the y-direction.

Analytical object model
The first object was a PMMA sphere with a radius of 0.805 cm Assuming PMMA (C 5 H 8 O 2 ) with a density of 1.19 g cm −3 and an x-ray energy of 27.7 keV, the parameters δ and β were found to be respectively 2.94 ×10 −7 and 1.18 ×10 −10 . The entrance air kerma used was 3 mGy, while for the modelling of the object in the simulation framework a perfect sphere was assumed.

Voxel models
Two more complex models were studied in addition: a 3D printed heart and a living mouse. Both objects were first scanned on a µCT (SkyScan 1278, Kontich, Belgium) and subsequently imaged using the GB-PCI system described above. The µCT had a reconstructed voxel size of 51.17 × 51.17 × 51.17 µm 3 reconstructed, with pixel values ranging between 0-255. A voxel model was constructed for each case by applying thresholds to the reconstructed µCT slices. For the 3D printed heart phantom, a threshold was set at 26 to divide the material into air and the 3D printed material. The object (Materialise, Leuven, Belgium) was made of 'Vero magenta' (Stratasys, Minnesota, USA) which has a density between 1.17-1.18 g/cm 3 . The exact chemical composition is, however, unknown and the β and δ coefficients were estimated from PMMA with a nominal density of 1.175 g/ cm 3 . Therefore, at 27.7 keV, β is calculated to be 1.49 ×10 −10 using NIST (Hubbell 1982, Seltzer 1993) and δ equals 3.43 ×10 −7 via Henke et al (1993). The mouse model consisted of five materials: air, compressed lung tissue, adipose tissue, muscle and bone with corresponding thresholds set at 45, 78, 100 and 200. For the latter, material compositions were taken from ICRP110 (ICRP 2009) and also evaluated using Henke et al (1993) and NIST data at an energy of 27.7 keV. The bony tissue was assumed to be 10% mineral bone and 90% spongiosa.

Analytical object model
For a simple object like the PMMA sphere, close agreement between simulation and experiment is found as demonstrated in figure 2. Quantitatively a SSIM between simulation and experiment of 0.99 and 0.81 is found for Tr and dP respectively. Note that the higher magnitude of noise in the differential phase image compared to the transmission image leads to a reduction in SSIM value.

3D printed heart
In figure 3 the result for the printed heart demonstrates that similar images are produced for more complex shapes. After registration, the simulated data (a,d) matches the experimental data (b,e) with similar contrasts as shown in figure 3.
An offset in Tr values can be seen in figure 3(c) which we attribute to inaccurate material coefficient modelling-we would expect a reduction in this offset if the exact chemical composition were known. However, it can be seen than similar contrasts are achieved. The maximum contrast in each of the ROIs (separated with blue vertical lines in figures 3(c) and (f), for Tr are respectively 0.21 0.07 0.08 and 0.12 in simulation and in experiment 0.19, 0.06, 0.08 and 0.10. Similarly in dP the maximum contrasts in experiment were respectively 0.34, 0.17, 0.14 and 0.15, while in simulation 0.34, 0.16, 0.13 and 0.12.

Chest RX mouse
Differential phase and transmission images of a mouse chest using the simulation framework match the appearance of the experimentally acquired datasets (figures 4(a), (b), (d) and (e)). Imperfect registration of the µCT and GB-PCI data causes some geometric variations, however, importantly the bony, soft tissue and lung field structures within these images are closely matched. This is also confirmed quantitatively where the maximum contrast in the three regions indicated in figures 4((c) and (f)) is respectively 0.16 0.33 and 0.21 for the experimental Tr data and 0.14 0.35 and 0.20 for the simulated data. For the dP data similarly the contrast is 0.25 0.20 and 0.27 for the experimental data and 0.17, 0.27 and 0.26 for the simulated data. The contrast is both over-and underestimate, which indicates that this is due to geometrical variations. As the images in the GB-PCI setup are taken horizontally, while in the µCT vertically, gravity acts differently on the animal and the anatomy in the experimental GB-PCI images is more compressed. It is also apparent that the simulated images have a slightly higher level of blurring, due to the modelling based on the µCT data with pixel sizes not an inverse multiple of the GB-PCI pixel sizes.
Bone was modelled as 90% rib spongiosa and 10% mineral bone, however, in the experimental data a higher degree of attenuation is seen. The simulated dP images underestimate the level of noise at the level of the lungs, since the dark field coherence loss was not accounted for.

Computational cost
A distinct advantage of this type of modelling over numerical wave propagation simulation platforms is the relaxation of the requirement for high spatial sampling. As a consequence, computational cost is greatly reduced. For a 200 × 200 pixel FOV (pixel size is 100 µm, p 2 is 2 µm and 8 G 2 steps) and a sampling step of 0.25 µm in the numerical wave propagation framework, the largest array in WIM consists of complex numbers and requires a memory of 8192 megabytes. Efficient variable use and memory handling is therefore required if memory allocation is not to become a limiting factor for this type of simulation platform. In the proposed hybrid framework maximum, estimated memory requirement for one variable is only 5.12 megabytes and CPU time to generate 100 images with a field of view of 200 by 200 pixels is 2.5 min. This allows the rapid generation of large data sets.

Discussion
The proposed framework differs from published frameworks as the fundamental idea is not to produce an extremely detailed simulation at the grating level that can assist in system design and development, but rather use established summary metrics as a means of generating realistic images quickly. This speed-up is paid for with a number of limitations. First and most straightforward, it is not possible to examine hypothetical cases, for example studying what would be the effect of increased G 1 -G 2 distance, as the simulation is fixed for the system for which the summary data were measured. The implications on visibility, flux, magnification, etc are therefore not implemented within the simulation and detailed grating modelling using for example a wave propagation framework should be used to first estimate the impact of this change. Simulating dark field images using the hybrid method requires an analytical relationship between the object properties and the measured signal. In Bech et al (2010)    et al 2016). The focal spot size, pixel size, magnification, beam hardening and strong phase gradients will induce pseudo-dark field contrast. All these effects should be included in order to simulate the dark field images with sufficient accuracy. Moreover, no tabulated values of are currently available in the literature for the relevant materials, further restricting the application of this type of modelling. For this reason, dark field modelling was not included in this work.
One of the consequences of the omission of the dark field component is an underestimation of the noise in regions of the dP images with a high level of small angle scatter. In the case of prior knowledge of the dark field signal (e.g. via previous experiments), this coherence loss of the signal can be included empirically in the noise calculation of the dP signal. The simulation currently assumes a monochromatic energy. Extending the framework to use a polychromatic source is possible but would require knowledge of visibility and a K response versus energy. As this information is difficult to access and the accuracy gained is expected to be limited, this was not considered. As a consequence, mean energy values for δ and β were used and beam hardening was not included. Lastly, large angle scatter was not accounted for. Most applications of interest are limited in size and the contribution of scatter to the detector is assumed to be minimal (Vedantham et al 2014, Vignero 2018. There is a movement towards the use of simulation and modelling to assist in the regulatory evaluation of medical devices in general (Morrison and Dreher 2017) and imaging products in particular (VICTRE, FDA (2018)). These types of studies can be referred to as computational modelling and simulation (CM&S) (Morrison and Dreher 2017) or as virtual clinical studies (Bakic et al 2016(Bakic et al , 2018 and are likely to become increasingly important in the development and justification of new medical imaging modalities. In addition to the system characterization required for the modelling (Mackenzie et al 2017) these methods require simulation of the anatomy (Kiarashi et al 2014, Segars et al 2018 and pathology of interest (Shaheen et al 2014). The presented framework allows rapid modelling of dP and Tr images for a well-established GB-PCI setup. After the user has characterized system setup as described in section 2.2.1, images suitable for use in clinical simulation study can be produced quickly, with a CPU time shorter than 3 min to generate 100 images with a 200 by 200 FOV. The framework described here is therefore ideally suited for these types of computational modelling and simulation studies where large numbers of images are needed. For example, nodules can be simulated into experimental or simulated lung images to study the effect of location and size on the detectability performance and to benchmark against the performance of current, standard imaging techniques. This should enable in-depth, quantitative studies into the added benefit of using dP images for specified tissue groups or for other new or even theoretical materials such as new contrast agents. Exploration of specific applications via computational modelling and simulation helps to justify animal and human experiments and can be a crucial step for GB-PCI to evolve into a clinical tool.
Of course there is the question of the required level of accuracy in the simulation, as it is crucial that accuracy is sacrificed for simulation speed. While there are some differences between the experimental and simulated data presented here, the accuracy is believed to be sufficient for these types of studies. In fact the required accuracy . The transmission and differential phase images resulting from simulation ((a) and (d)) and from experiment (b) and (e). In plot ((c) and (f)) the average cross section in the red box is plotted, demonstrating close quantitative agreement between simulation and experiment. The blue vertical lines indicate the different ROIs where the maximum contrast is calculated (c) and (f). is likely to be task dependent and is difficult to generalise. As for any simulation study, realistic modelling of the objects is crucial: here experimentally acquired µCT data was used for the modelling. The results suffered slightly from the influence of the µCT resolution and ideally a higher resolution µCT scanner should be used. Further improvement of the simulated images could also come from manual segmentation of the images.

Conclusion
A simulation framework for producing GB-PCI differential phase and transmission images has been introduced with the goal of generating images for use in computational modelling and simulation studies. The framework uses 3D voxel models and typical metrics measured/estimated for the system in question, namely the detector presampling MTF, NPS, visibility and a dose-to-noise calibration factor. The major advantage of the proposed simulation framework is the gain in simulation speed and reduced memory requirements. With only a small loss of accuracy, transmission and differential phase images with large fields of view can be generated quickly without the requirement of high spatial sampling.