Analytical reconstructions of intensity modulated x-ray phase-contrast imaging of human scale phantoms

: This paper presents analytical approach to modeling of a full planar and volumetric acquisition system with image reconstructions originated from partial illumination x-ray phase-contrast imaging at a human scale using graphics processor units. The model is based on x-ray tracing and wave optics methods to develop a numerical framework for predicting the performance of a preclinical phase-contrast imaging system of a human-scaled phantom. In this study, experimental images of simple numerical phantoms and high resolution anthropomorphic phantoms of head and thorax based on non-uniform rational b-spline shapes (NURBS) prove the correctness of the model. Presented results can be used to simulate the performance of partial illumination x-ray phase-contrast imaging system on various preclinical applications.


Introduction
Conventional medical x-ray imaging system is based on the changes in linear attenuation coefficients between tissues that produce differences in photon fluence incident upon the detector [1]. However, these differences between soft and bone tissue are significant, but are small among the different types of soft tissue, which results in low contrast (or signal-to-noise ratio). Thus, the x-ray attenuation contrast is relatively poor and cannot achieve satisfactory sensitivity and specificity [2,3]. Low x-ray energies can improve the x-ray attenuation contrast, however the photon absorption causes damage to the living cells of the body and limits the range of the energies that may be used. Recent developments in imaging techniques provide the potential to reduce the dose with simultaneous improvements in signal-to-noise ratio. One of these methods is a phase-contrast imaging, which enables the detection of differences in the refractive index, especially for tissues with low absorption [4,5]. Inner features of the examined objects have been observed successfully using Bonse-Hart geometry or diffraction enhanced imaging (DEI), which offers the possibility to distinguish structures through diffraction from perfect crystals ("analyzer based" imaging) [6,7] or free space propagation [8,9].
Although direct phase measurement is impossible, interferometric methods can be used to determine the phase-shift [10,11] with the cross section of the order of 100 to 1000 times greater than the absorption cross section in the biological soft tissue over the 10-100 keV range [12]. However, in grating interferometry the gratings have small pitches of the order of a few microns and require very precise alignments (to a few tens of nm [13]). Contrary to the Talbot (grating-based) method, the partial illumination phase-contrast technique (PIXPCi) is an incoherent one and is based on another physical principle, the so-called edge illumination instead of Talbot's self-imaging, while grating interferometry requires at least spatial coherence [14]. In particular, for grating-based phase contrast imaging when the source dimensions are made larger, phase effects disappear very fast, due to the loss of flux at the gratings; hence, microfocal sources can be used instead of synchrotrons, but this results in excessive long exposure times even running up to hours [6,15]. However, using liquid-metaljet sources [16] a shorter exposure time of a few minutes can be obtained.
In general, the x-ray photon path can be measured on the detector and makes image formation possible. Measuring such images as projections from different angles or using sample rotation allows us to obtain 3D reconstruction of the object [17]. Images projected by divergent x-ray beams at long sample-detector distances must be acquired with detectors of reasonable dimensions, which can be challenging to use with typical flat panel detectors, especially in a clinical routine. But, phase-contrast imaging technique has the potential to expose detailed structural variations between different tissues, thereby providing a high contrast and resolution to distinguish healthy and malignant tissues. However, simulation of the PIXPCi at a human scale poses a challenge in the imaging model development [14].
In this study, the authors have proposed an analytical model of low dose partial illumination x-ray phase contrast technique to reconstruct the planar and volumetric phase contrast maps of the adult human scale objects as a two-dimensional images and tomographic set of projections respectively, obtained from the simulation of phase-shift detection and measurements. Projections of the phase-shift induced by the sample serve to reconstruct a 3D refractive index distribution using a tomographic reconstruction algorithm. The volumetric xray images have the potential to reveal subtle differences with high resolution, even in tissues of low absorption properties. Currently, there is no hardware x-ray phase contrast imaging system that can image large structures of human body and simulation to predict its potential is needed.

Experimental setup
Partial illumination x-ray phase-contrast imaging (PIXPCi) technique is based on the edge illumination method [18] and its setup is illustrated schematically in Fig. 1. A conventional xray source which generates a diverging beam may be applied when the dimensions of the sample aperture are scaled down to provide highly sensitive phase response with fully polychromatic beams and uncollimated, unapertured source size of up to 100 μm i.e. current mammography sources [18]. Pixel edge illumination is employed by the x-ray coded apertures phase-contrast imaging system, which works by projecting small x-ray beams onto the edges of sensitive detector pixels. The details of the PIXPCi technique can be found in [19] and the partial pixel illumination condition is presented in [20,21] as shown schematically in Fig. 2. Fig. 2. Illustration of the partial pixel illumination condition, where a set of pre-sample slits creates an aperture and subdivides primary beam into small beamlets.

Simulation overview of PIXPC imaging system
Modeling of the partial illumination x-ray phase contrast was proposed to reconstruct a series of x-ray phase images of an object based on the geometrical ray-optics approach with comparable results obtained using paraxial approximation to the Fresnel-Kirchhoff diffraction theory [22,23]. The main disadvantage of the wave-optics mathematical approach [24-26] is that there is no straightforward solution for geometrically complicated objects that cannot be represented by simple geometrical forms.
Ray-optics approach is based on using ray-tracing methods from x-ray source to detector in relaxed coherence conditions [27] and multislice propagation model [27][28][29] is used for coherent source. Simulation calculates x-ray deflections inside the object due to refraction effects at the different tissue boundaries (represented numerically by different refractive index values), in particular the gradient of the real part of the refractive index inside the specimen was used to get the scatter angle of the photon incident upon the detector [6]. Calculated deviation angles should be sampled with a resolution corresponding to the beam dimensions, but higher resolution of the sampled object demands longer calculation times and superior amount of computational resources.
Simulation considers the physical processes of scattering and the overall system geometry (mainly parameters of the x-ray source, apertures, samples and detector) besides calculating the path of the photons that pass through any component of the system, absorption, angle deviation of the particle and map of refractive index, where each element modifies the characteristics of the beam.

Modeling of the x-ray tube source
The laboratory x-ray tube sources were introduced to the model through sampling of the full x-ray spectrum. The source was described by using bremsstrahlung energy spectra, depending on target atomic number, incident kinetic energy and fraction of energy radiated. Main function, which describes a source, is to determine a set of parameters for each x-ray i.e. position, velocity and emission time. The initial position is found from uniform distribution over the circular source surface. The initial velocity is selected within an interval of the corresponding energy. The model allows the choice of the emission time of each x-ray, which is determined using a Monte Carlo method.
The total intensity of the photon flux is defined as the sum of the weights of all emitted xrays for every single simulation. The unit of total photon weight is the number of x-rays per second. The source flux , Φ is defined as the number of x-rays emitted per second from a 1cm 2 area on the source surface, with direction within a steradian solid angle, and with wavelength within 1 Å interval. For a constant , Φ the source area -S, the solid angle of the aperture as seen from the source surface -, ΔΩ and the width of the wavelength interval in which photons are emitted, the total intensity of x-rays emitted towards a given aperture [30] .
For a given number N of simulated x-ray histories, the initial photon weight w 0 which statistically describe a given x-ray in the beam along its trajectory through the components of optical system.
To simulate x-ray tube source an electron beam of the uniform intensity of transverse cross section and energy E 0 impinges on the infinitely thick target material. Monte Carlo methods are used to determine photon generation, bremsstrahlung or the x-ray emission lines of the material. Using the Monte Carlo method, it is possible to transport electrons and photons inside the target and filter to obtain detailed information about the weight factors contributing to the production of the x-ray spectrum [31].
The produced bremsstrahlung is described by the spectral characteristics proposed by Kramers [32]. The set of Cauchy functions [33] allow the sampling of characteristic emission with spectral widths [34].

Simulation of phase effects and contrast formation
The simulated photon is projected onto the sample at a particular angle and later the position of the intersection point is calculated. The auxiliary normal vector to the surface is built at the intersection point and computation of the incidence angle of the photon serves to determine the transformed photon in the object. Refractive index is also calculated in the point of intersection (see Eq. (6)) and refraction of x-ray is determined for the photon traversing through the sample until it reaches the second interface (at the boundary between heterogeneities characterized by refractive index).
Subsequently, similar calculations are performed at the photon's exit point, which is used to compute the new direction of the outgoing photon. Transmission of the object, which describes the influence of the variation of the sample's structure on the incident photon, is calculated according to the absorption coefficient; then, the weight of the photon is adjusted.
Directional focusing method was applied and calculations were done only for the x-rays in the directions that have a chance of being detected and weight factors were accordingly corrected, since Monte Carlo sampling of the rays having their density proportional to the beam intensity can be extremely inefficient for sharply peaked distributions [30]. Weight factors w i : 0<w i <1 (because total emitted flux should decrease through the system) for projected rays allow to speed up the calculation process, since only a small amount of primary x-rays reach the detector and if the absorption of any element in the system is equal to the threshold, only rays which pass through are considered.
Each simulated photon is characterized by its position ( , , ), r x y z  wave vector k  (strictly speaking, its projections on the axes -, , x y z k k k ), phase , φ polarization E and weight w. The transformation of the incident photon during the interaction with the object was determined in terms of Snell's law [35] 1 1 where 1 θ is the incident x-ray angle with respect to the plane of the object interface, 2 θ is the outgoing photon ray angle and 1 , j δ − j = 1,2 denote the real part of the complex refractive index. For a complex angle, we assume .
The refractive index data was calculated theoretically [36]. The δ and β are defined in terms of real and imaginary parts of the complex atomic scattering factor where 1 ( ) f ω is the ratio of the electric field strength scattered by an atom to that of a single free electron, ω is the incident wave frequency, a n is the atomic density in atoms per unit volume, e r is classical electron radius, λ is the wavelength in vacuum and E is the x-ray photon energy. Hence, the refractive index can be defined as and where 1 1 f is written as having first order term Z, the number of electrons per atom, and a departure therefrom due to the degree of binding. C P indicates taking only the non-divergent Cauchy principal part of the integral. Thus, this provides a technique for numerically determining values of . σ depends only on the x-ray energy and at high photon energies the linear attenuation coefficient is nearly proportional to the electron density [39]. The bigger contribution of the Compton effect to the attenuation process, the more the attenuation-based image is proportional to the electron density, similarly to the phase contrast. Therefore, for low atomic number Z materials the contrast differences for phase contrast and attenuation images at high x-ray energies are similar [39].
The coherent scattering is modeled in terms of molecular effects, especially for biological imaging of tissues with atomic number Z~10. This effect is responsible for changes in the direction of motion according to a distribution that is determined by interference among the electrons in the atom. The forward direction is preferred, the peaking in this direction increasing with increasing energy and decreasing atomic number. As a fraction of the total mass attenuation coefficient, σ coh /ρ is maximal at atomic numbers around Z = 10 and photon energies in the interval 30-50 keV. At higher atomic numbers, the relative fraction decreases due to the strong increase of photoelectric absorption with increasing atomic number. In a molecule, coherent scattering is determined by the interference of all the electrons in the molecule. Unfortunately, the effect is not very well known for most molecules. However, water serves as a biological tissue equivalent in numerous modeling applications, i.e. dose distribution calculations in diagnostics or radiation therapy; thus for water, molecular effect has been determined experimentally and expressed in terms of scattering probabilities [40].
Water molecular form factors, W at ( ), F x have been calculated with the independent atomic scattering approximation from the oxygen, O ( ), F x and hydrogen, H ( ), F x relativistic atomic form factors given by Hubbell and Overbo [41] [ ] [ ] where   [42]. For liquid water and the energy 15 keV, E > the numerical integration procedure used logarithmically spaced points over the range 12
The beam attenuation by the scattering and absorption in the sample is taken into account in Monte Carlo simulations. Thus, the absorption cross section per unit cell of a sample is a σ and a scattering cross section is .
where l is the full path length in the sample. The probability for a given x-ray to be scattered from within the interval 1 1 [ , ] l l dl + and the probability for an x-ray to be scattered from within the interval into the solid angle Ω but not being absorbed or scattered further on the exit point of the sample where 1 l is the x-ray path length in the sample before scattering, 2 l is the x-ray path length in the sample after the scattering event and ( ) g Ω denotes the directional distribution of the scattered x-rays [30].
The pre-and post-sample apertures are introduced to the simulation as a two-dimensional meshes with slits in two perpendicular directions and allow 2D phase change registration. A given x-ray incident on the mask is traced to find an intersection point between photon and aperture. Assuming photon intersects the mask, the path length dl is calculated within the mask. Apertures are assumed to be made of gold material with the thickness of 130 μm, however it is possible to use custom material and thickness. The nominal density μ of the aperture as a function of wavelength is calculated by interpolating data extracted from the NIST x-ray database [44] and the photon weight is adjusted to 0 exp( ). w w dl μ = − × Similarly, custom x-ray source filters can be generated and introduced to the PIXPCi simulation using NIST data.
Signal diffusion effect [19] was also introduced to the simulation in the adjacent detector pixels. Hence, image contrast as a peak-to-peak ratio is reduced by the influence of opposite low and high peaks in the neighboring pixels. Signal induced in the pixels is given as a percentage of counts in the detector.

Multislice simulations
The method is based on a discretization of the wave propagation along the direction of the incident wave from one plane to another. The electromagnetic wave field propagating through the tick sample is modeled by a multislice approximation [27] and beam propagation methods [27,45] that divide the 3D sample into a series of N thin slices, each of thickness z Δ with complex refractive index ( , ) 1 ( , ) ( , ). n x y x y i x y δ β = − + We assume, that the incident plane-wave radiation propagates in the positive z direction. The optical constants of the sample are estimated from the atomic scattering factors of its constituents (see Eqs. (4), (5) and (6)). The approximation for small-angle scattering (paraxial approximation) [46] is applied to the propagated x-rays. The phase and amplitude change is caused by propagation of photons through a given slice (single plane) followed by propagation through distance z Δ to the next slice. The distribution of complex transmission of the slice is denoted by is the wave number at wavelength . where the Fresnel diffraction integral [47] is The convolution form of the Fresnel integral is calculated as a product in Fourier space [48] where  and 1 −  are the Fourier transformation and its inverse respectively. The transmitted wave function ' Q ψ ψ = is computed for a given slice, propagated to the next slice (see Eq. (22)) and then becomes the incident wave function for that slice. The iteration proceeds until the wave function in the exit plane is obtained.
Hence, the intensity distribution of the wave is 2 ( , ) ( , ) , I x y x y ψ = as a phase contrast image. Thus, phase-shift can be reconstructed from the phase-shift images.
The multislice method for a polychromatic (incoherent or partially coherent) case of object, whose optical properties may vary among the slices and as a function of position within each slice consists of making repeated m monochromatic calculations at closely spaced wavelength intervals, and accumulating the recorded intensity from each wavelength after propagation from the exit plane of the sample to the detector plane. Thus, for partially coherent source the intensity distribution is The multislice method is expected to be valid if a sufficient number of closely spaced slices are taken such that spread of the wave function caused by Fresnel diffraction within the slice thickness does not exceed spatial resolution required in the calculation [49].
For propagation through a distance : z Δ 2 2 ( ) ( ) . x y zλ Δ + Δ << Δ The geometrical constraints should be small enough to represent the 3D distribution of the optical constants inside the specimen. To avoid artifacts at the ends of the sample array, the simulated object is stored in extended arrays that the infinite results assumed by the FFT are sufficiently separated and not interact. Such a procedure leads to a cumulative loss of high spatial frequency data in the Fourier convolutions, however the initial oversampling is set to compensate this effect.

Phase contrast volumetric reconstruction
Biological objects exhibit both strong amplitude and phase modulation, even in the hard x-ray region. The intensity distribution of a monochromatic plane wave at distance d and angle of rotation , θ which propagates along the positive z direction and incident upon a mixed phase and amplitude sample can be approximated as [46] 2 , , 0 ( , ) ( , )exp ( , ) , The two terms are considered independently since the volumetric (tomographic) reconstruction is a linear operation. The first term denotes the projection through the absorption index and the second term is the 2D Laplacian of the refractive index decrement projection in the image plane. Hence, mixed phase and attenuation reconstruction image is ( ) as a superposition of 3D Laplacian of the refractive index decrement and linear attenuation coefficient. However, such approximation is valid for small refraction angles. The phase contrast images are used directly as input to a tomographic reconstruction algorithm. Volumetric image registration mode requires the specimen to be imaged at different projection angles - [0, ] θ π ∈ to obtain mixed phase and amplitude reconstructions x y x y z x y z θ π μ β λ =  (31) where ' cos sin , ' , ' sin cos x x z y y z x z θ θ θ θ = − = = + and the integration is done for a given x-ray path.

Data acquisition system and analysis
Here, the proposed model was based on quasi-Monte Carlo and wave optics methods, where propagation of a single ray (with weight factors) emitted from the source and the reconstruction of three-dimensional volume of the object were calculated in parallel using graphics processing unit (GPU) multithreading, implemented with Nvidia CUDA. The synergy of GPU combined with Nvidia CUDA technology was applied to optimize the level of simulation processing and it was a strategy useful in making significant advances in speed of data treatment. The calculated multiple values were grouped in the process within the GPU, reducing the number of memory reads. To model PIXPCi system, the x-ray source, numerical phantoms, pre-and post-sample masks and detector had to be designed numerically as a C++ code, tested on the cluster computer system at our department and specialized Dell T5600 workstation with dual Intel Xeon E5 2 GHz processors (12 cores), 32 GB of RAM memory, Nvidia Quadro 600 and Nvidia GTX770 GPU graphics.
Multiple sets of sphere and wire phantoms with different refractive index values and easily defined geometry were produced as a numerical code for simplicity to generate clear refraction and absorption maps. Special inserts allowed for taking into account system response to different types of tissues in numerical phantom material with cross sections as illustrated in Fig. 3. The phantom has water as a medium and soft tissue inhomogeneities, inter alia, lung, adipose, muscle tissues, bones and main insert as a reference-dense bone. Furthermore, two copies of the anthropomorphic phantom were used. The NURBS-based extended human body-head and cardiac-torso phantom (4D XCAT 2) [50] as presented in Fig. 4 with well-defined structures-served as an input for testing the model's usability in diagnostic and preclinical applications.  The cross sections of imaged sample were taken into account by simulation and the map of refractive index was calculated within the object. After performing an image sequence at a given angular sample position or by using sample rotation, projections from arbitrary angles allowed to produce a three-dimensional reconstruction of the specimen. Such an approach makes it possible to reveal volumetric structural variations inside the object or test the influence of misalignment of the apertures on the quality of images, in terms of contrast reduction. The phase and attenuation images were calculated by Fourier analysis of the intensity curves in each detector pixel.

X-ray tube source simulations
An estimated spectra for 80 kVp and 100 kVp x-ray tube voltage emitted from a12°pure tungsten target with 1.2 mm Al filter for the beam hardening are plotted in Fig. 5. The calculated spectrum is normalized to the total number of photons in the spectrum and shows characteristic tungsten x-rays at 1 (58 keV), K α 2 (59.5 keV), K α 1 (67.5 keV) K β and 2 (69 keV).
The mean spectrum energies are 40.6 keV and 47.1 keV respectively and first HVLs are 1.93 mm and 2.51 mm Al. The results as shown in Table 1 are in good agreement with the IPEM Report [51], Bhat et al [52] and Fewell et al [53] however, binning the data into 0.5 keV energy intervals causes small shift in characteristic x-ray energy. Our model allows using a tungsten source in 10-100 keV mean spectrum energy range.
In addition, the anode heel effect and off-axis x-ray spectra were assessed for different anode angles, but these details are not within the scope of this paper. Following the experiment in Leitenberger et al [54], we simulated a source (120 μm, 40 keV energy) projected on a plane with two slits of 2 μm diameter separated by 10 μm. The intensity profiles as a diffraction patterns calculated for a coherent and partially coherent beam are shown in Figs. 6a and 6b respectively. The result of the simulation is in very good agreement with the experimental results shown in [54].

Simulated partial illumination x-ray phase contrast radiography
First, our simulation system was set up to acquire planar phase-contrast images, which can be feasible in a diagnostic system to identify structures that might be invisible on attenuation projections. The PIXPCi simulation presented in this study calculates both the amplitude and phase of electric field at the detector plane. Image registration is performed for the object positioned at a source-sample distance (SSD) of 100 cm and the detector placed at a sourcedetector distance (SDD) of 260 cm from 40 keV-and 60 keV-point x-ray sources. The optics has been set to an effective pixel size of 50 μm. Phase-shift effect with edge enhancement was simulated for materials of varied structural and physical properties, which provided mixed strong and weak absorption and refraction signals at selected energies.
Planar images of simple cases as sphere and wire phantoms are presented in Fig. 7. Four spheres of 60 μm with a central sphere of 120 μm radius served as a sphere phantom. Set of three wires of 38 μm, 30 μm and 50 μm width served as a wire phantom. Both phantoms were imaged with appropriate second aperture shifts to assure pixel illumination conditions with 0%, 66% and 100% of individual beam to assess different set-ups from attenuation through phase contrast to dark field scans, as shown in Fig. 7. Imaging of low complexity objects provided proof of model's accuracy.
Dark field (ultra-small angle scatter) images were based on measuring the illumination curve obtained by scanning the pre-sample aperture with the fixed rest of the set-up [55]. Profiles measured in x direction for the wire phantom are presented in Fig. 8. These profiles provide the evidence of the edges detection capability of phase contrast imaging system for the different materials, which occur as negative and positive peaks in the signal. The highest signal values were obtained for bones, especially dense bone which has the largest mass density of all investigated materials in the phantom, as well as the highest atomic number Z. On the contrary, lung tissue material yields the lowest peak-to-peak signal and was harder to distinguish from the surrounding water-imitating substance 10  in attenuation contrast image. However, mixed phase and attenuation contrast images (see Figs. 7(b) and 7(e)) provide higher level of perceptibility of such material with concurrent better edge detection and structure differentiation.

Volumetric simulations of partial illumination x-ray phase contrast imaging
Simulation of much more complex examples of irregular tissue structures was performed. Heterogeneous numerical phantoms made of tissue-like substance and NURBS-based, anthropomorphic head and thorax 4D XCAT 2 phantoms with known anatomical features were imaged at the same settings for clarifying the results. Images of a head part, with a 10 mm spherical lesion in the brain region, and a thorax part of non-uniform rational b-splines phantom were acquired. The lesion was clearly visible as an abnormal structure. Furthermore, 4D XCAT 2 phantoms were scanned in axial, sagittal and antero-posterior direction for revealing phase-shift effects in clinically relevant plane directions. As done previously, geometrical and irradiation conditions were applied to the imaging of these phantoms to assure realistic conditions for the imaging of human body. The sample was rotated over 360 degrees with a 1 degree angular step.
As a result, maps of δ and β were reconstructed with the filtered backprojection algorithm (FBP). The volumetric projections of attenuation, edge enhancement due to phase contrast and dark field images were registered with the detector aperture positioned such that 0%, 50% and 100% of each beamlet falls onto the uncovered pixel area. These are shown in Figs. 9 and 10 as a three dimensional reconstruction of a recorded set of images taken during rotation of the sample (presented as cross sections in relevant directions). Reconstruction and visualization of volumetric structural variations in the specimen were obtained. Maps of the mixed signal were also reconstructed using FBP from the projections containing combination of phase and attenuation contrast. The high level of details is visible, as shown in Figs. 9(c) and 9(g).
It was also possible to freely move and rotate any component of the system, which enabled translation and rotation in three dimensions, to make the model scalable and resistant to any custom settings. These features allowed for the modeling of volumes in phase-contrast computed tomography (CT). However, such an approach requires further quantitative analysis of attenuation and phase contrast, especially in biological materials. Fig. 9. Images of head region of 4D XCAT 2 phantom, simulated with 60 keV x-ray source. Reconstructions in sagittal (upper row) and axial (bottom row) directions presented as phase image (a, e), attenuation (b, f), mixed phase and attenuation contrast (c, g) and dark field (d, h) corresponding to a pixel-illuminated fractions of 50% -(a, e, c, g); 100% -(b, f); 0% -(d, h). Lesion of 10 mm diameter is clearly visible (arrow) on (c) and (g) as an irregular structure in the occipital part of the brain. Fig. 10. Images of thorax region of 4D XCAT 2 phantom, simulated with 60 keV x-ray source. Reconstructions in antero-posterior (upper row) and axial (bottom row) directions presented as phase projections (a, e), attenuation (b, f), mixed phase and attenuation contrast (c, g) and dark field (d, h) corresponding to a pixel-illuminated fractions of 50% -(a, e, c, g); 100% -(b, f); 0% -(d, h).

Discussion
In this paper we present a wide range of results obtained from imaging of the objects characterized by different properties with simple structure, complex pseudo-organic systems and even high complexity approximations of the human body using head and thorax numerical phantoms. Modeled samples were not completely transparent to x-rays, where the presence of the phase-shift effect was ambiguous. However, our preclinical research focused on the objects relevant to further investigations toward a fully clinical use of such a simulator and dose modeling, although obtained images were not as clear as for samples with low attenuation only.
Anthropomorphic human tissue phantoms based on NURBS allowed us to predict system behavior in case of use in clinical situations. However, these phantoms are still too simple to serve as a precise imitation of biological structures, where the boundaries between different types of tissues are very complex and cause higher contrast images. Fig. 11. Mixed phase and attenuation contrast reconstructions on sagittal cross sections of 4D XCAT 2 human scale phantom (head part) imaged at: (a) 40 keV, which served as an extreme case lower limit for hard x-ray energy, (b) 60 keV energy of conventional x-ray tube source, (c) 100 keV as a maximum hard x-ray energy. Definitely more details can be observed at (a), but by using higher x-ray energy, the dose to "body" is significantly lower and the differences between structures can still be distinguished with edge enhancement effect.
The energies chosen for the purpose of the present work were 40 keV, 60 keV and 100 keV, where at lower energies we observed larger differences in refractive index and higher contrast, as presented in Fig. 11. However, the lower energy used in the model served as a comparison to make sure that simulation works on a wide range of conditions, but low hard xrays are not clinically relevant due to increased absorbed dose. The 60 keV and 100 keV energies made the phase-shift effect less noticeable, but still possible to visualize and served as energies typical for conventional x-ray tube sources. The level of pixel illumination made it possible to assess the images as a result of various physical effects; specifically, by continuously switching the system settings to work in attenuation, phase contrast or dark field imaging regime. Edge illuminated images give an impression of a spatial and for diagnostic purposes it makes it easier to distinguish boundaries even between tissues of similar attenuation. However, ray tracing or Monte Carlo methods are unsuitable for image formation because they do not model coherent scattering and diffraction processes with sufficient accuracy [19,21,25]. In addition, such a simulation becomes more difficult as the size of phantoms increases.
Volumetric reconstruction done with a set of projections made it possible to reveal structural variations in any direction. However, the quality of these reconstructions could be improved by implementing iterative reconstruction algorithms. Also the quality of planar images could be highly improved at the software level to deal with noise of the imaging system and blurring caused by scattering, which is less expensive in comparison to hardware solutions. Furthermore, sampling causes noise due to aliasing of high-frequency signal components, and digitization produces quantization errors, hence the ideas of such an enhancement on the quality was introduced to the model based on edge preserving sparse regularization methods, which will be presented later. The sparse regularization transforms can be very useful as a volumetric phase contrast reconstruction techniques and dose reduction methods [56].
In general, objects characterized by low phase-contrast tend to exhibit additional interference fringes at the edges as a negative effects of finite source sizes and limited spatial resolution [6] but their presence can be minimized through image restoration improvements that are now under investigation.
In particular, the design of the presented system provides attenuation of the primary beam through the introduction of a pre-sample mask, which makes it possible to control the dose to the sample. However, our present approach does not incorporate dose distribution analysis inside the object; our intention being that it should be the focus of our future work that would also require using Monte Carlo methods for further optimization.