1. 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-metal-jet 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 x-ray 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.

2. Materials and methods

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 x-ray 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].

 

Fig. 1 Schematic illustration of the partial illumination x-ray phase-contrast imaging setup with the edge illumination concept.

Download Full Size | PDF

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.

Download Full Size | PDF

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–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 x-rays for every single simulation. The unit of total photon weight is the number of x-rays per second. The source flux $\Phi ,$ is defined as the number of x-rays emitted per second from a 1-cm² area on the source surface, with direction within a steradian solid angle, and with wavelength within 1 Å interval. For a constant $\Phi ,$the source area – S, the solid angle of the aperture as seen from the source surface –$\Delta \Omega ,$and the width of the wavelength interval in which photons are emitted, the total intensity of x-rays emitted towards a given aperture [30]

To simulate x-ray tube source an electron beam of the uniform intensity of transverse cross section and energy E₀ 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 $\overrightarrow{r}(x,y,z),$wave vector $\overrightarrow{k}$ (strictly speaking, its projections on the axes – $k_x,k_y,k_z$), phase $\varphi ,$ 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]

The refractive index data was calculated theoretically [36]. The $\delta$and$\beta$are defined in terms of real and imaginary parts of the complex atomic scattering factor $f^1(\omega )=f_1^1(\omega )-i{f}_2^1(\omega )$

The phase-contrast as a local distribution of the refractive index decrement $\delta$is related to the electron density${\rho }_e$and the energy of the x-ray photons for energies far away from absorption edges

The attenuation contrast image is calculated as a result of photoelectric absorption$\tau ,$ incoherent (Compton) scattering${\mu }_{\text{incoh}},$ and coherent (Rayleigh) scattering ${\mu }_{\text{coh}}$ [39]. Thus, linear attenuation coefficient μ

The incoherent scattering cross section${\sigma }_{\text{incoh}}$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,$F_{\text{at}}^{\text{W}}(x),$ have been calculated with the independent atomic scattering approximation from the oxygen,$F_{\text{O}}(x),$and hydrogen,$F_{\text{H}}(x),$relativistic atomic form factors given by Hubbell and Overbo [41]

The integrated coherent (Rayleigh) scattering cross sections ${\sigma }^{\text{w}}{}_{\text{s,mol}}$– per free water molecule, ${\sigma }^{\text{w}}{}_{\text{s,liq}}$– per liquid water molecule were obtained for $x>1.25{\AA }^{-1}$and the value $F_{\text{at}}^{\text{W}}(x),$has been taken as an approximation for $F_{\text{mol}}^{\text{W}}(x)$and$F_{\text{liq}}^{\text{W}}(x).$

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 ${\sigma }_a$and a scattering cross section is${\sigma }_s.$ For a given volume of a unit cell V the inverse penetration lengths are defined by ${\mu }_s={\sigma }_s/V$and ${\mu }_a={\sigma }_a/V.$ Hence, the attenuation function along this path

The probability for a given x-ray to be scattered from within the interval $[l_1,l_1+dl]$

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 $w=w_0\exp (-\mu \times 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$\Delta z$with complex refractive index $n(x,y)=1-\delta (x,y)+i\beta (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$\Delta z$to the next slice.

The distribution of complex transmission of the slice is denoted by

The phase-shift $\phi$ and the attenuation $\mu$ induced by the sample can be described as projections through the absorption and refractive index distributions respectively

Assuming the paraboloidal approximation, the illumination function (defined in real space) to the next slice in the z direction is given by Fresnel propagator

Thus, the wave function ${\psi }_N$at the Nth slice is computed recursively according to the information in the previous slice by convolution (denoted by star) as

The convolution form of the Fresnel integral is calculated as a product in Fourier space [48]

The transmitted wave function $\psi \text{'}=\psi 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$I(x,y)={\left|\psi (x,y)\right|}^2,$ 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$I(x,y)={{\displaystyle {\sum }_m\left[{\left|\psi (x,y)\right|}^2\right]}}_m.$

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$\Delta z:$${(\Delta x)}^2+{(\Delta y)}^2<<\Delta z\lambda .$

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$\theta ,$which propagates along the positive z direction and incident upon a mixed phase and amplitude sample can be approximated as [46]

Volumetric image registration mode requires the specimen to be imaged at different projection angles – $\theta \in [0,\pi ]$ to obtain mixed phase and amplitude reconstructions

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.

 

Fig. 3 (a) Axial cross section of numerical sphere phantom with inserts imitating different types of tissue in water: 1-lung, 2-adipose, 3-muscles, 4-bones, 5-reference insert: dense bone, (b) 3D model of sphere phantom, and (c) numerical wire phantom with inhomogeneities.

Download Full Size | PDF

 

Fig. 4 Illustration of the 4D XCAT 2 anthropomorphic phantom with inhomogeneities generated throughout changes in attenuation coefficients in the volume of organs of interest recorded within binary reproduction: (a) axial scan with simulated 10 mm diameter spherical lesion (arrow), (b) 4D XCAT 2 overview of thorax and head region, and (c) sagittal cross section with 10 mm lesion (arrow).

Download Full Size | PDF

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.

3. Results

X-ray tube source simulations

An estimated spectra for 80 kVp and 100 kVp x-ray tube voltage emitted from a$12°$pure tungsten target with 1.2 mm Al filter for the beam hardening are plotted in Fig. 5.

 

Fig. 5 Plots of simulated x-ray spectra for 80 kVp and 100 kVp tube voltage for pure tungsten target with 1.2 mm Al beam hardening filter.

Download Full Size | PDF

The calculated spectrum is normalized to the total number of photons in the spectrum and shows characteristic tungsten x-rays at $K_{\alpha 1}(58\text{keV}),$ $K_{\alpha 2}(59.5\text{keV}),$ $K_{\beta 1}(67.5\text{keV})$ and $K_{\beta 2}(69\text{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.

Table 1. Comparison between measured (Fewell 1981, Bhat 1998), calculated (IPEM report 78) and simulated (presented model) HVLs and mean spectrum energy for spectra produced using 12° target angle and different tube voltages. The percentage difference between IPEM and our model estimates is also shown.

View Table

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.

 

Fig. 6 Diffraction patterns using a coherent source (a) and partially coherent source (b) in the simulation of Leitenberger et al. experiment.

Download Full Size | PDF

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 source-detector 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.

 

Fig. 7 Simulation images of sphere and wire numerical phantoms for 40 keV photons. For all simulations the number of tracked particles was 1.5 × 10⁹, the computational time per particle was ~0.5 × 10⁻⁶ s. Settings of pixel fraction illumination used to image both phantoms as a post-sample mask offset: 0% - (a, d); 66% - (b, e); 100% - (c, f) and percentage of attenuation image: 100% - (a, d); 60% - (b, e); 0.1% - (c, f).

Download Full Size | PDF

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 $({\beta }_{\text{w}}=8.97\times {10}^{-10},{\delta }_{\text{w}}=2.46\times {10}^{-6})$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.

 

Fig. 8 A series of plot profiles extracted from the reconstructed maps of the wires formed as the images shown in Fig. 7(a)-7(c) corresponding to a percentage of the attenuation image: 100% - (a), 60% - (b), 0.1% - (c).

Download Full Size | PDF

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

 

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.

Download Full Size | PDF

 

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

Download Full Size | PDF

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.

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

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 x-rays 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.

 

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.

Download Full Size | PDF

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.

Presence of scattered radiation degrades the contrast and the scatter-to-primary ratio (the ratio of the contributions to the signal from scattered and primary photons) and may take values as high as 5–10 in x-ray diagnostic examinations of thick body parts using large field areas, thus diminishing image quality.

Following the possibility of realization of the three-dimensional set-up and volumetric reconstructions through modeling, it is hoped that the usability of the system presented in this study should be achieved in diagnostic imaging of the human scale objects. Especially, in cases where low doses during on-line imaging are required and hard constraints for the imaging are imposed to reduce overall treatment time, i.e. image guided radiotherapy during stereotactic treatments.

5. Conclusions

Fast model of planar and volumetric intensity modulated partial illumination x-ray phase-contrast effect at a human scale has been presented. This simulation can be completely controlled by the user. Its design was tailored to PIXPCi technique that can produce phase-contrast maps at an adult human scale and was not based on any x-ray propagation framework, thus making it robust to any further modifications, i.e. dose or x-ray source modeling, custom geometry changes, detector response simulations, image restoration employment or stress tests in extreme cases and quantitative analysis.

Developing the model by using GPU with Nvidia CUDA platform significantly speeds up the calculations and shortens acquisition and reconstruction time, which makes it possible to further reduce time on the x-ray-tracing level. This step will be optimized in the next version of the model. The dose reduction obtained with concurrent higher peak-signal-to-noise ratio can be achieved by using sparse regularization reconstruction transforms i.e. compressed sensing [56], which have been developed to use increasingly less raw data. However, even with dose modeling capability within the imaged objects and simulation of the system response to the dose, there is still a shortage of full clinical implementation of such a technique to realize its advantages in clinical practice [57,58]. Furthermore, image quality improvements in terms of contrast enhancement with simultaneous noise reduction should be considered toward building a fully clinical prototype of a human scale phase-contrast imaging system and the framework presented in this study allows to predict the performance of reconstruction algorithms and such system.

Here presented model of partial illumination x-ray phase contrast technique may be obtained from http://www.fuw.edu.pl/~jkpie/medsoft/