Open Access
Add to BibSonomy
Abstract
Edge illumination x-ray phase contrast imaging (XPCI) provides increased contrast for low absorbing materials compared to attenuation images and sheds light on the material microstructure through dark field contrast. To apply XPCI in areas such as non-destructive testing and inline inspection, where scanned samples are increasingly compared to simulated reference images, accurate and efficient simulation software is required. However, currently available simulators rely on expensive Monte Carlo techniques or wave-optics frameworks, resulting in long simulation times. Furthermore, these simulators are often not optimized to work with computer-aided design (CAD) models, a common and memory-efficient method to represent manufactured objects, hindering their integration in an inspection pipeline. In this work, we address these shortcomings by introducing an edge illumination XPCI simulation framework built upon the recently developed CAD-ASTRA toolbox. CAD-ASTRA allows for the efficient simulation of x-ray projections from CAD models through GPU-accelerated ray tracing and supports ray refraction in a geometric optics framework. The edge illumination implementation is validated and its performance is benchmarked against GATE, a state-of-the-art Monte Carlo simulator, revealing a simulation speed increase of up to three orders of magnitude, while maintaining high accuracy in the resulting images.
© 2024 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
X-ray phase contrast imaging (XPCI) setups can be used to acquire projections with higher contrast between low attenuating materials compared to conventional attenuation-based x-ray projections [1]. Moreover, most XPCI setups provide information on the microstructure distribution of the scanned sample through the dark field contrast [2–4]. Hence, XPCI can be effectively applied to industrial inspection [5,6], including integration into inline setups [7].
The XPCI setup considered in this work is the conventional two mask edge illumination (EI) setup [8], although the presented experiments can be easily adapted to accommodate a single mask setup [9]. A conventional EI setup relies on two absorbing masks with slit-shaped apertures, one to split the source x-ray beam into smaller beamlets, the other to cover the edges of every detector pixel row or column (depending on the mask aperture orientation). By combining acquisitions for different mask positions, the attenuation, phase and dark field contrasts can be obtained.
In non-destructive testing (NDT), the current standard x-ray inspection workflow is to extract a surface mesh from a computed tomography (CT) image and compare it with a sample computer-aided design (CAD) model [10]. As at each projection angle in EI multiple images are acquired at different mask positions, recording a complete EI-CT data is a time consuming process. Alternatively, sample projections can be compared directly to simulated reference images, reducing the required number of projections [11] and providing opportunities for high throughput EI industrial inspection. However, assessing the sample quality in projection space requires simulating reference projections from the sample CAD model, and the current lack of efficient EI simulation software makes this unfeasible in high throughput scenarios.
The EI setup is non-interferometric, relying on detecting beam refraction when a phase shift occurs instead of detecting wave interference effects. Hence, a geometric optics approach can be used to model it [12,13]. Advantages of a geometric optics over a wave optics approach are a reduced complexity of the physics involved, increased flexibility in the setup geometry and typically shorter simulation times. Previously, a geometric optics approach to modeling EI has been implemented as stochastic photon tracing algorithms in the form of Monte Carlo simulators in software packages such as Geant4 [14,15], GATE [16,17], and McXtrace [18,19]. A well-known downside to the stochastic nature of Monte Carlo simulations is the requirement to cast an abundant number of rays/photons to obtain workable statistics in the resulting images. Consequently, although they provide valuable results, the current Monte Carlo and CPU EI implementations often require a long simulation time, limiting the amount of data that can be generated. These long simulation times also hinder their application in iterative routines, such as the design and/or optimization of XPCI setups, and makes the simulation of extensive datasets, e.g. for the evaluation of reconstruction algorithms or the training of deep learning networks, practically unfeasible. A CPU-based geometric optics ray tracing algorithm to model the EI setup was presented in [12]. However, ray tracing lends itself to full parallelization over the rays, favoring a GPU implementation, as has been done previously when solely considering x-ray attenuation in [20–22].
In this paper, an x-ray EI phase contrast simulation framework built upon CAD-ASTRA, an open-source GPU ray tracing mesh projection toolbox [23], is presented. Preliminary results of a GPU EI implementation, based on the CAD-ASTRA toolbox, were presented in [24]. We expand upon those results here by introducing a new source model to better capture the effect of a finite focal spot and a new material parameter to account for microscopic sample surface roughness. We also describe an EI mask generation tool, perform several validation experiments and benchmark against state-of-the-art Monte Carlo reference images.
2. Methods
2.1 Edge illumination
A conventional EI setup consists of two absorbing masks with slit-shaped apertures. The first (sample) mask is placed in front of the sample and the second (detector) mask is placed right in front of the detector, as illustrated in Fig. 1. Mask dimensions are chosen such that they have the same projected aperture width and a pitch that matches the pixel size. Both masks are aligned to be parallel to the detector plane, with their apertures oriented along the detector columns. The detector mask is positioned such that the apertures are aligned with the center of a detector pixel column, while the absorbing bars block the pixel edges, creating insensitive regions on the detector. An acquisition consists of acquiring projections at different lateral steps of the sample mask, shifting the position of the sample mask aperture relative to that of the detector mask.
Fig. 1. Top view sketch of a typical EI setup (not to scale), where $z_{so}$ is the source - sample mask distance (also approximately the source - object distance), $z_{od}$ is the sample mask - detector distance (also approximately the object - detector distance) and $z_{sd}$ the source - detector distance, $t_i$ is the thickness, $p_i$ the pitch and $s_i$ the aperture width of mask $i$, with $i$ = 1 for the sample mask and $i$ = 2 for the detector mask. The mask movement direction to perform mask stepping is indicated by $f$.
The intensity modulation caused by mask stepping is called the illumination curve (IC). Due to a finite focal spot and detector blurring it has an approximately Gaussian shape. Introducing a sample and repeating the stepping process will result in a modified IC due to beam attenuation and refraction by the sample. Examples of a flatfield and sample IC are shown in Fig. 2.
Fig. 2. An example illustration of a flatfield and sample IC, which can be constructed by performing multiple acquisitions at different lateral sample mask positions, first without the sample present and then with the sample in field of view.
The attenuation, phase and dark field contrasts can be derived, per pixel, by comparing ICs for an acquisition without (flatfield) and with a sample present. The IC is defined as
with $A$ the amplitude, $\mu$ the mean position, $\sigma ^2$ the variance, and $b$ an intensity offset caused by the x-rays that leak through the absorbing masks. It is through the comparison of the flatfield and sample IC’s, for each pixel, that the different contrasts can be extracted [25]. Attenuation will cause a reduction in the area under the curve, while refraction will cause the mean position of the IC to shift. Finally, in case the sample has many micro structures that cannot be individually resolved, i.e. their projected dimensions are smaller than the aperture dimensions, the manyfold refraction by these sub-pixel structures will result in a broadening of the measured IC. This third source of contrast is termed the dark field contrast. The process of retrieving the different contrasts from the projection data is called phase retrieval and is performed through the following formulas: where $T$ stands for the transmission through the sample (which can be expressed in terms of attenuation through the Beer-Lambert law) and $\alpha$ is the beamlet refraction angle (under a small angle approximation). $M = z_{sd} / z_{so}$ is the geometric magnification of the sample mask, $z_{od}$ the object-detector distance and $\sigma ^2_{DF}$ is the angular broadening of the beamlet. The subscripts $s$ and $f$ indicate the sample and flatfield IC parameters, respectively. All parameters are also represented graphically in Fig. 2. Since Eq. (1) has four unknown parameters, at least four different mask steps are required in the acquisition process to be able to fit the ICs.2.2 EI using CAD-ASTRA
2.2.1 CAD-ASTRA
CAD-ASTRA is a recently developed toolbox for efficient x-ray projection simulations of mesh models [23]. The toolbox relies on NVIDIA OptiX to parallelize ray tracing on the GPU [26]. Its projection geometry can be defined as a standard ASTRA [20] cone beam or cone vector geometry, allowing for easy reconstruction of simulated sinograms in ASTRA. Complex samples can be modeled by nesting meshes inside each other, where each mesh can be given different material properties. The material properties include a linear attenuation coefficient $\mu$, the real part $\delta$ of the complex refractive index $n$, and a measure for the material surface roughness (see below). Ray intensity at the detector is calculated using $\mu _k$ and ray path length $l_k$ through each mesh $k$ in the Beer-Lambert law:
where $I_0$ is the flatfield intensity without a sample present. Refraction is modeled using Snell’s law, and thus depends on the incoming ray direction $\vec {i}$, the normal on the mesh surface $\vec {n}$, the angles $\theta _i$, $\theta _r$ between the incoming and the refracted ray and the surface normal, and the refractive indices $n_i$ and $n_r$ on either side of the interface:As material properties are energy-dependent, Eqs. (5) and (6) are valid only for a monochromatic x-ray source. The EI setup, however, can be used with a standard, polychromatic laboratory x-ray source [27]. To account for this, a polychromatic source spectrum was calculated using GATE Monte Carlo simulations [16]. The simulation included CAD models of the FleXCT source [28] and generated accelerated electrons hitting the source anode, hence generating x-rays. The EM standard physical model was used and an ideal detector was placed behind the source beryllium window and tungsten aperture to capture the x-rays, providing weights for the different x-ray energies, as shown in Fig. 3(a). Polychromatic projections in CAD-ASTRA are generated by calculating a weighted sum of monochromatic projections, using these energy weights.
Fig. 3. Source spectrum and spot shape estimates of the FleXCT [28] x-ray source used for some of the simulation experiments.
By default, CAD-ASTRA uses a point source from which a single ray is cast per detector pixel, with the rays aimed at the pixel centers. The number of rays can be independently increased along the detector row and/or column direction. If multiple rays are cast, they are evenly distributed across the detector pixel. A beamlet generated by an EI mask has a finite width, while a single ray is infinitesimally thin. Therefore, the effect of the number of rays on how well the beamlets are modeled will be investigated in the experiments. The EI setup measures beamlet refraction in the direction orthogonal to the mask apertures, so only the number of rays cast per pixel in this direction is expected to influence the simulation accuracy.
2.2.2 Source spot size sampling
The IC shape of an EI acquisition is influenced by the source spot size [2]. A point source results in triangular ICs, whereas (realistic) finite sources result in Gaussian ICs. This can be understood by noticing that mathematically the IC shape is a convolution between the source shape and the two mask functions, the latter two represented by step functions which are 1 inside the mask apertures and 0 elsewhere. Hence, to accurately simulate a laboratory EI setup, finite spot sizes have to be included. One method to include a finite source spot is by taking a weighted superposition of point sources, as described in [23]. Alternatively, the finite spot can be simulated through random sampling of the ray starting locations according to a Gaussian distribution representing the source shape. In this random sampling model, the ray starting positions for a projection are first set in one point, just like for a point source. Then, two Gaussian distributions are defined, one representing the source shape along the detector columns, the other along the detector rows. Next, each ray’s starting point is shifted from its initial location along the detector row and column directions. The magnitude of the translation along these directions is determined by randomly drawing a sample from the corresponding Gaussian distributions (cfr. Figure 3(b)). In this way, a single projection inherently includes the finite spot size, while the weighted superposition method requires looping over multiple point source projections. In this work, the random sampling method will be employed.
2.2.3 Mask generator
The EI masks are modeled as a collection of box-shaped surface meshes, each box representing a mask bar. The distance between consecutive bars, as well as their individual width and thickness, can be set to generate slit-shaped apertures. The most common case is a mask with a regular aperture pitch equal to the demagnified pixel size at the mask plane and an aperture width around 1/5th to 1/10th of the demagnified pixel size. A tool was created to generate the bars of EI masks based on the desired setup geometry. The bars of each mask are concatenated into a single object which can be given appropriate material properties and can be placed such that one mask is in front of the sample and the other mask is in front of the detector. For the experiments described below, this standard type of EI masks was modeled. Additional parameters can, however, be added to the mask generator tool to accommodate for more specialized mask designs, such as those described in [7,29].
2.2.4 Surface roughness
Surface meshes of CAD models (including the EI masks) are constructed out of perfectly flat triangles, which are typically too large to capture the microscopic surface roughness of a real object. In 3D printing, for example, objects exhibit a certain surface roughness resulting from the printing process, which is not present in the source CAD model. This can lead to simulation artefacts such as an unrealistic amount of total reflection at the sample edges. Theoretically, the surface roughness of an object can be modeled by repeatedly subdividing the mesh and slightly adjusting the orientation of the resulting small triangles. In practice, however, this is seldom a feasible solution. Instead, an artificial roughness can be introduced to the CAD model by applying a small update to the surface normal of each mesh triangle during the ray tracing process, before applying Snell’s Law, as described in [30]. The surface roughness implementation in CAD-ASTRA uses two orthogonal unit vectors, $\vec {u_1}$ and $\vec {u_2}$, in the triangle plane (i.e. also both orthogonal to the surface normal $\vec {n}$). Both unit vectors are scaled with random samples, $s_1$ and $s_2$, from a Gaussian distribution $s_{1,2} \in \mathcal {N}(0, \sigma _{\text {rough}})$, centered around zero and with the standard deviation $\sigma _{\text {rough}}$ determining the roughness of the mesh surface. Adding both scaled vectors to the surface normal gives the new rough surface normal:
Just like the linear attenuation coefficient and refractive index, the material roughness parameter $\sigma _{\text {rough}}$ can be set for each mesh individually.
3. Experiments
Simulation experiments were designed to validate the EI CAD-ASTRA implementation. The influence of the number of rays cast per pixel on attenuation and phase contrast images was investigated to determine the optimal number of rays required to generate qualitative images. Different source spot sizes were simulated to examine their influence on the shape of the IC, hence validating the random ray-origin sampling described in 2.2.2. The surface roughness parameter was tested by comparing projections of a rough mesh with those of a smooth mesh, where the roughness parameter was added to the smooth mesh. Then, using the optimal parameters found in these validation experiments, CAD-ASTRA simulations were compared to Monte Carlo GATE reference simulations for both a single projection and a CT reconstruction.
The CAD-ASTRA simulations were run on a desktop machine consisting of a Intel Core i7-5960X 3 GHz CPU processor and NVIDIA GeForce RTX 3070 GPU using CUDA version 12.0 and NVIDIA OptiX version 7.4.0. The GATE simulations were performed using an extended version of Gate 8.0 that includes ray refraction [17] on server hardware with AMD EPYC 7452 32-Core 2.35 GHz CPU processors.
3.1 Simulation setup
The acquisition geometry used for the EI simulations was based on the phase contrast setup that was recently installed in the FleXCT scanner at the University of Antwerp [28,31]. The total length of the system is 1.8 m, with $z_{so}$ = 1.2 m and $z_{od}$ = 0.6 m the source-sample mask and sample mask-detector distances, respectively, resulting in a magnification $M$ of 1.5. The detector had 150 $\times$ 150 µm2 pixels, and both masks were made of gold, with a projected aperture width of 30 µm and a thickness of 225 µm. The setup was aligned such that the 0 µm mask step corresponded with perfectly overlapping sample and detector mask apertures (i.e. peak flatfield IC intensity). More detailed setup specifications, or any deviations from these settings, are given for each experiment separately.
3.2 Reference images
A phase-sensitive version of GATE was used to generate reference images [17]. To optimize simulation time, the GATE simulations were parallelized over multiple CPU cores by using its job splitter [32]. The total simulation time given in Table 1 is the sum over all parallel jobs. For both simulation tools, the three contrasts are retrieved by applying the IC-fitting phase retrieval method detailed above, using the Gpufit curve fitting toolkit [33].
Table 1. Simulation run times for the CAD-ASTRA and GATE comparisons. For the GATE simulations, the total run time (i.e. sum over split jobs) is given.
3.3 Influence of number of rays
In the first experiment, the influence of the number of simulated rays on the retrieved phase contrast signal was investigated. Projections of a carbon cylinder were simulated for increasing numbers of rays, from 1 up to 5000 rays per pixel. Rays were only added in the phase-sensitive direction of the setup (the vertical direction in Fig. 1), as this is the only direction in which refraction is measurable. The cylinder had a diameter of 8 mm and a height of 16 mm, and the detector counted 128 by 128 pixels. The linear attenuation coefficient and refractive index were set for a monochromatic x-ray beam of 30 keV according to [34] and [35], respectively. The source spot dimensions were based on the FleXCT source [28], with a horizontal and vertical full-width at half-maximum (FWHM) of 30 µm and 26 µm, respectively, as illustrated in Fig. 3(b). Each simulation consisted of 11 steps of the sample mask, equally spaced in the [−25 µm, 25 µm] interval.
3.4 Influence of source spot size
For the second experiment, flatfield acquisitions for different source spot sizes were generated, illustrating the effect of a finite source spot on the ICs. The mask material properties were set for a monochromatic x-ray beam of 30 keV and the IC was sampled over 21 equally spaced steps in the [−40 µm, 40 µm] interval. As only effects along the (horizontal) phase sensitive direction were of interest, the vertical FWHM of the focal spot was kept constant at 26 µm, while the horizontal FWHM was increased from 0 (i.e. a point source) up to 30 µm.
3.5 Influence of surface roughness
The surface roughness parameter implementation was validated in the third experiment. An octagon with a diameter of 10 mm was extruded into a smooth 3D surface mesh. A rough version of this mesh was created by subdividing the original 32 face triangles into 32 768 triangles and then applying a small random push sampled from a Gaussian distribution to all mesh vertices. Both the smooth and rough versions of the mesh are shown in Fig. 4. The mesh material parameters and the source shape were identical to the first experiment, with 5000 rays being cast per pixel. Only the central (perfectly aligned) sample mask position was considered. A 128 pixels wide line detector was used and the sample was gradually rotated over 45 degrees to capture varying levels of total reflection on the perfectly smooth mesh surface. In total, 45 projections were acquired. First, the smooth and rough mesh projections were simulated, after which the smooth mesh simulation was repeated, but now with a nonzero surface roughness material parameter.
Fig. 4. The smooth (blue) and rough (red) octagonal surface meshes used in the surface roughness experiment.
3.6 Radiograph comparison with GATE
As a fourth experiment, retrieved transmission and refraction values for a CAD-ASTRA and a GATE simulation were compared to each other by overlaying line profiles of a water, polyethylene and carbon cylinder. All three cylinders had a diameter of 6 mm and their material properties were set for the polychromatic x-ray spectrum shown in Fig. 3(a). The detector had 64 by 256 pixels, and 5000 rays were cast per pixel for the CAD-ASTRA simulation. Mask stepping was performed over 11 equally spaced steps in the [−25 µm, 25 µm] interval, and the source spot dimensions were the same as in the first experiment. For the GATE reference data, the same geometric setup, samples, spectrum, and steps were used, with a perfect photon counting detector [17], while simulating 2$\times 10^8$ photons per mask step.
3.7 CT comparison with GATE
In the fifth and final experiment, a CT dataset was simulated. The sample was a polyethylene cylinder with a diameter of 20 mm, which contained both an area with approximately 80 000 thin carbon fibers of 5 µm diameter, and an area with 39 thicker 300 µm diameter carbon fibers. The sample also contained three larger 1 mm diameter cylinders, made of carbon, epoxy and nylon 66. A render of the nested cylinder meshes is shown in Fig. 5, alongside a zoom-in on the patch of carbon microfibers. Projections were simulated for 360 equally spaced projection angles in the [0°, 360°) interval, where each angle contained 7 equally spaced mask steps in the [−30 µm, 30 µm] interval. The detector consisted of a single row of 256 pixels. After phase retrieval, the CAD-ASTRA projection geometry could be reused in the reconstruction, which was performed using a gradient descent method with the Barzilai-Borwein step size [36]. The source shape was kept identical to that of the first experiment and the energy was monochromatic at 30 keV, with the attenuation coefficient and refractive index again set based on [34] and [35], respectively. Each mesh was given a surface roughness $\sigma _{\text {rough}}$ of 0.002. For the GATE reference data, the same setup geometry and source energy was used, while simulating 5$\times 10^6$ photons per projection angle and per mask step. However, built-in GATE geometries were used to model the sample, instead of providing surface mesh STL files. This change was made because the STL photon tracing that is available in GATE has poor performance, and hence slows down considerably with samples containing tens of thousands of microstructure meshes. Furthermore, the GATE sample was given a surface roughness using the implementation described in [30], with a variance of 0.05.
Fig. 5. The CT sample, with labeled and color-coded materials. The zoom-in shows the carbon microfibers. See the text for details.
4. Results
4.1 Influence of number of rays
The results of the first experiment are shown in Fig. 6. The columns show increasing numbers of rays cast in the phase sensitive direction. From the images, it can be seen that the attenuation contrast (apart from the case with 1 ray per pixel) remains largely unaffected, while the retrieved phase contrast does depend on the number of rays cast in the phase sensitive direction. More rays per pixel correspond to a smoother, more clearly defined phase contrast, while practically no phase contrast can be observed for the case with only one ray per pixel.
Fig. 6. CAD-ASTRA EI simulations of a carbon cylinder with different numbers of rays cast in the phase sensitive direction. Increasing the number of rays improves the shape of the phase signal.
4.2 Influence of source spot size
The influence of the source spot size on a flatfield IC is presented in Fig. 7. The red markers are the simulated mask steps, while the solid blue curve is the fitted Gaussian that is used for phase retrieval. IC curves simulated with a point source (in the phase sensitive direction) result in triangular ICs, as theoretically predicted [37]. ICs simulated with a finite width source, due to the blurring the source introduces in the forward projection model, gradually approximate a Gaussian shape for increasing source widths. Good agreement is found between the simulated mask steps and the Gaussian IC fit for a horizontal source FWHM of 30 µm.
4.3 Influence of surface roughness
Results of the surface roughness experiment are presented in Fig. 8. For the simulation of the smooth mesh with an added surface roughness parameter, the surface roughness parameter was estimated to be 0.001 by testing a range of values. In Fig. 8(a), sinograms over an angular range of 45 degrees are shown for the smooth mesh, rough mesh, and smooth mesh with surface roughness parameter. It can be seen that for a few angles (just past 20 degrees rotation) strong edge effects occur, which are exacerbated by the smooth mesh. The accompanying line profiles for one of these angles, shown in Fig. 8(b) and Fig. 8(c), indicate that adding a virtual roughness parameter to the smooth mesh can give similar results to the rough mesh without having to change the mesh geometry.
Fig. 8. Sinogram and corresponding line profiles for the angle indicated with the dashed lines for the smooth mesh, rough mesh, and smooth mesh with added surface roughness.
4.4 Radiograph comparison with GATE
The transmission and refraction profiles of three different material cylinders, for both a CAD-ASTRA and a Monte Carlo simulation using GATE, are shown in Fig. 9. Good agreement is found between the transmission and refraction angle line profiles, which were averaged over the detector rows to increase the GATE simulation statistics. The only noticeable difference is in the attenuation contrast of the polyethylene cylinder, where the CAD-ASTRA simulation shows a slight attenuation underestimation. The dark field contrast is not shown, as no microstructures were present in the sample. Simulation times are presented in Table 1, showing two orders of magnitude simulation speed increase for CAD-ASTRA compared to GATE.
Fig. 9. Transmission and refraction profiles generated with CAD-ASTRA and GATE for a water (W), polyethylene (PE) and carbon (C) cylinder. The line profiles are averaged over the detector rows.
4.5 CT comparison with GATE
Reconstructed slices of the CT sample (Fig. 5) for all three contrasts are shown in Fig. 10. The top row shows the reconstruction based on CAD-ASTRA projection data, while the bottom row shows the GATE projection data reconstruction. The phase contrast provides higher contrast than the attenuation contrast for the largest nested cylinders, consisting of different materials. Moreover, the thicker carbon fibers (in the bottom of the sample) are more clearly delineated from the polyethylene background in the phase contrast than in the attenuation contrast. Both the attenuation and phase contrast show the patch of carbon microfibers, picking up that the bulk material composition inside the patch is a mixture of polyethylene and carbon. The dark field contrast, however, clearly shows the presence of unresolvable small features in this area. Boxplots of several single material regions of interest (ROIs) are compared to their expected theoretical value in Fig. 11. For the attenuation contrast, good agreement is found between the mean reconstructed values and the theoretical linear attenuation coefficients for both CAD-ASTRA and GATE, with similar error bars for the two tools. For the phase reconstruction the error bars for both tools are noticeably smaller, but the difference with the theoretical delta values is larger. The values found for the two phase reconstructions are similar to each other, with the CAD-ASTRA values being slightly higher than the GATE values. Simulation times are found in Table 1, showing a three orders of magnitude simulation speed increase for CAD-ASTRA compared to GATE.
Fig. 10. Reconstructions of the CT sample for all three contrasts. The top row uses CAD-ASTRA projections, while the bottom row uses GATE projections. Note that the phase reconstruction shows the refractive index decrement $\delta$ and is dimensionless.
Fig. 11. Different material ROIs (drawn on the CAD-ASTRA phase reconstruction to illustrate) and corresponding box plots for all reconstructed voxels in the color-coded boxes. The color-coded dashed lines in the box plots represent the theoretical value.
5. Discussion
The experimental results in section 4 serve as a demonstration and illustration of the capability to generate qualitative EI simulations with the CAD-ASTRA toolbox. In the results of the first experiment, it is shown that the number of rays cast per pixel influences the smoothness of the resulting phase signal. This is caused by the fact that we use a number of zero-width rays to model finite-width x-ray beamlets. Consider a sample mask aperture beamlet being modeled by only a single ray. After being refracted by the sample, the ray can either go through the detector mask aperture or be absorbed by the mask septa. This results in a binary signal on the detector. Increasing the number of rays cast per pixel increases the number of beamlet rays (rays going through the sample mask apertures). Moreover, the beamlet rays are spread out across the sample mask apertures, hence modeling the finite beamlet width. The more rays are cast per pixel, the more densely packed the beamlet rays are, and the less relative weight each individual beamlet ray has. Thus, the more rays are cast per pixel, the lower the impact is of a single beamlet ray being refracted toward the detector mask septa on the pixel intensity. This effectively smooths the beamlet refraction signal, i.e. the phase contrast, as shown in Fig. 6. For an EI setup, which has a single phase sensitive direction, it is sufficient to increase the rays per pixel only in the phase sensitive direction. However, it should be noted that for other XPCI methods, such as speckle-based imaging [38], which have 2D phase sensitivity, the number of rays should be increased in both the row and column directions.
The IC model error that is introduced when simulating a point source is visualized in Fig. 7. Strictly speaking, not a point source but a line source was simulated. However, just like above for the number of rays cast, only the source size in the phase sensitive direction influences the IC. Thus, the line source can in this context be treated as a point source. Theoretically, the IC phase retrieval framework can also be used by employing a triangle fitting procedure for the ICs, instead of a Gaussian one. The attenuation, phase and dark field contrasts would be related to the triangle surface area, peak location, and width, respectively. Real-world laboratory x-ray sources, however, have a finite spot size and thus require the Gaussian fitting procedure. Additionally, a finite but small source spot size, such as the 10 µm FWHM of the middle IC in Fig. 7, results in a mixed shape between a triangle and a Gaussian. A dedicated study on the topic of how the source spot size influences EI image quality has recently been conducted [39], indicating the importance of including this parameter in an EI simulation framework.
Figure 8 shows how the surface roughness parameter model compares against an actually rough surface mesh, while the smooth mesh simulation in the figure serves as a demonstration why surface roughness is important. Even though a rough surface mesh was used in the experiments, this should not imply using rough meshes is straightforward. Generating the rough surface mesh required increasing the number of mesh triangles of the smooth mesh by three orders of magnitude to reach the high resolution upon which to add the roughness. Furthermore, this process would have to be repeated for each new mesh. Parameterizing the surface roughness and adding it to smooth meshes offers a much higher flexibility in tuning the desired material roughness and makes it easily transferable between different meshes, while providing simulation results which are very similar to an actually rough mesh.
On the subject of mesh detail, it should be pointed out that using an insufficient number of triangles to model a surface mesh (especially a quadric surface) can lead to phase-contrast artifacts. These artifacts are alleviated by more fine triangular meshing, which results in more memory intensive mesh models. Depending on the surface roughness, however, the roughness can smooth out these artifacts for a lower level of mesh detail than would be required for the case without mesh surface roughness. Actual required mesh detail, and associated memory demand, will thus depend on an interplay between the mesh shape and surface roughness. Furthermore, even memory intensive meshes, such as the roughly 500 MB mesh containing the microfibers used in the CT experiment, still exhibit good performance when used in CAD-ASTRA simulations.
A first comparison between CAD-ASTRA and state-of-the-art GATE simulations (Fig. 9) shows good agreement between the CAD-ASTRA and GATE refraction profiles of a water, polyethylene and carbon cylinder, indicating the Snell’s law ray refraction model in CAD-ASTRA is correctly implemented. Small differences are visible for the transmission profiles, which can be ascribed to the different particle attenuation models in the two simulation tools: GATE uses a complex list of physical processes to determine how the simulated photons interact with the sample, while CAD-ASTRA bundles all these processes in the linear attenuation coefficient and the Beer-Lambert law given in Eq. (5). This hypothesis is supported by existing literature [40,41], where similar attenuation differences were found when comparing between different x-ray databases. Note that, when enough compute resources are available, the GATE simulation can be parallelized over multiple CPU cores [32]. To remove this hardware dependency from the comparison, the total simulation times given in Table 1 are the summed simulation times of these parallel jobs. For this experiment, the GATE simulation was parallelized over 11 CPUs, which still resulted in an order of magnitude slower simulations than CAD-ASTRA. Furthermore, the profiles in Fig. 9 were averaged over the detector rows to reduce the noise in the GATE simulations. To remove this averaging without changing the photon statistics, e.g. to compare 2D images, the number of photons in the GATE simulation would have to be increased, hence further increasing the simulation time.
Reconstructions of the CT sample, shown in Fig. 10, show similar results between the CAD-ASTRA and GATE simulated data. The phase reconstructions exhibit a higher signal to noise (SNR) ratio than their attenuation counterpart, as indicated by the smaller error bars on the boxes in Fig. 11. However, the reconstruction error, compared to the theoretical value, is larger for the phase than for the attenuation reconstruction. This is likely caused by the reconstruction method, and not the projection data. The phase reconstruction involves a numerical integration of the refraction angle sinogram to convert the measured angle into a phase shift, which is prone to generate reconstruction artefacts. Note that, while a EI-CT acquisition requires many more projections than a EI radiograph, the CAD-ASTRA simulation time went down compared to the previous experiment. This is due to the specific implementation of CAD-ASTRA routines, where some mesh transformation have to be performed on the CPU, while others are performed on the GPU. As the monochromatic CT simulation could leverage more GPU operations, it performed better than the polychromatic simulation. The GATE simulation was split over 40 CPU cores, reducing the effective simulation time to just over 4 hours. This number of cores requires dedicated server hardware and is rarely available on a desktop machine.
6. Conclusion
We demonstrated the ability to efficiently and accurately generate EI simulations using the CAD-ASTRA toolbox. Several new additions were made to the toolbox to enable these simulations, including an improved source model, an EI mask generator, and a material surface roughness parameter. The forward projection model was successfully validated and a comparison was made with state-of-the-art GATE simulations, showing good agreement for all three contrasts between the two simulation tools. CAD-ASTRA outperformed GATE simulation times by two to three orders of magnitude, depending on the use case. This speed up could be increased even further by implementing current CAD-ASTRA CPU methods on the GPU.
Funding
Fonds Wetenschappelijk Onderzoek (FoodPhase (S003421N), G090020N, G094320N).
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. M. Endrizzi, “X-ray phase-contrast imaging,” Nucl. Instrum. Methods Phys. Res., Sect. A 878, 88–98 (2018). [CrossRef]
2. M. Endrizzi, P. C. Diemoz, T. P. Millard, et al., “Hard X-ray dark-field imaging with incoherent sample illumination,” Appl. Phys. Lett. 104(2), 024106 (2014). [CrossRef]
3. K. Taphorn, F. De Marco, J. Andrejewski, et al., “Grating-based spectral X-ray dark-field imaging for correlation with structural size measures,” Sci. Rep. 10(1), 13195 (2020). [CrossRef]
4. S. Meyer, S. Shi, N. Shapira, et al., “Quantitative analysis of speckle-based X-ray dark-field imaging using numerical wave-optics simulations,” Sci. Rep. 11(1), 16113 (2021). [CrossRef]
5. D. Shoukroun, L. Massimi, M. Endrizzi, et al., “Edge illumination X-ray phase contrast imaging for impact damage detection in CFRP,” Mater. Today Commun. 31, 103279 (2022). [CrossRef]
6. M. Endrizzi, B. Murat, P. Fromme, et al., “Edge-illumination X-ray dark-field imaging for visualising defects in composite structures,” Composite Structures 134, 895–899 (2015). [CrossRef]
7. M. Endrizzi, A. Astolfo, F. Vittoria, et al., “Asymmetric masks for laboratory-based X-ray phase-contrast imaging with edge illumination,” Sci. Rep. 6(1), 25466 (2016). [CrossRef]
8. A. Olivo, “Edge-illumination x-ray phase-contrast imaging,” J. Phys.: Condens. Matter 33(36), 363002 (2021). [CrossRef]
9. F. A. Vittoria, G. K. N. Kallon, D. Basta, et al., “Beam tracking approach for single-shot retrieval of absorption, refraction, and dark-field signals with laboratory x-ray sources,” Appl. Phys. Lett. 106(22), 224102 (2015). [CrossRef]
10. W. Dewulf, H. Bosse, S. Carmignato, et al., “Advances in the metrological traceability and performance of X-ray computed tomography,” CIRP Ann. 71(2), 693–716 (2022). [CrossRef]
11. A. Presenti, J. Sijbers, and J. De Beenhouwer, “Dynamic few-view X-ray imaging for inspection of CAD-based objects,” Expert Syst. with Appl. 180, 115012 (2021). [CrossRef]
12. A. Olivo and R. Speller, “Modelling of a novel x-ray phase contrast imaging technique based on coded apertures,” Phys. Med. Biol. 52(22), 6555–6573 (2007). [CrossRef]
13. P. R. Munro, K. Ignatyev, R. D. Speller, et al., “The relationship between wave and geometrical optics models of coded aperture type x-ray phase contrast imaging systems,” Opt. Express 18(5), 4103–4117 (2010). [CrossRef]
14. J. Allison, K. Amako, J. Apostolakis, et al., “Recent developments in Geant4,” Nucl. Instrum. Methods Phys. Res., Sect. A 835, 186–225 (2016). [CrossRef]
15. L. Brombal, F. Arfelli, F. Brun, et al., “X-ray differential phase-contrast imaging simulations with Geant4,” J. Phys. D: Appl. Phys. 55, 045102 (2021). [CrossRef]
16. D. Sarrut, M. Bardiès, N. Boussion, et al., “A review of the use and potential of the GATE Monte Carlo simulation code for radiation therapy and dosimetry applications,” Med. Phys. 41(6Part1), 064301 (2014). [CrossRef]
17. J. Sanctorum, J. De Beenhouwer, and J. Sijbers, “X-ray phase contrast simulation for grating-based interferometry using GATE,” Opt. Express 28(22), 33390–33412 (2020). [CrossRef]
18. E. Knudsen, A. Prodi, J. Baltser, et al., “McXtrace: A Monte Carlo software package for simulating X-ray optics, beamlines and experiments,” J. Appl. Crystallogr. 46(3), 679–696 (2013). [CrossRef]
19. T. P. Millard, M. Endrizzi, P. C. Diemoz, et al., “Monte Carlo model of a polychromatic laboratory based edge illumination x-ray phase contrast system,” Rev. Sci. Instrum. 85(5), 053702 (2014). [CrossRef]
20. W. van Aarle, W. J. Palenstijn, J. Cant, et al., “Fast and flexible X-ray tomography using the ASTRA toolbox,” Opt. Express 24(22), 25129–25147 (2016). [CrossRef]
21. Q. Gong, R.-I. Stoian, D. S. Coccarelli, et al., “Rapid simulation of X-ray transmission imaging for baggage inspection via GPU-based ray-tracing,” Nucl. Instrum. Methods Phys. Res., Sect. B 415, 100–109 (2018). [CrossRef]
22. Á. Marinovszki, J. De Beenhouwer, and J. Sijbers, “An efficient CAD projector for X-ray projection based 3D inspection with the ASTRA toolbox,” in 8th Conference on Industrial Computed Tomography, Wels, Austria, (2018).
23. P. Paramonov, N. Francken, J. Renders, et al., “CAD-ASTRA: a versatile and efficient mesh projector for x-ray tomography with the ASTRA-toolbox,” Opt. Express 32(3), 3425–3439 (2024). [CrossRef]
24. N. Francken, P. Paramonov, J. Sijbers, et al., “Enhancing industrial inspection with efficient edge illumination x-ray phase contrast simulations,” in IEEE EUROCON 2023 - 20th International Conference on Smart Technologies, (2023), pp. 723–727.
25. M. Endrizzi and A. Olivo, “Absorption, refraction and scattering retrieval with an edge-illumination-based imaging setup,” J. Phys. D: Appl. Phys. 47(50), 505102 (2014). [CrossRef]
26. S. G. Parker, H. Friedrich, D. Luebke, et al., “GPU ray tracing,” Commun. ACM 56(5), 93–101 (2013). [CrossRef]
27. P. C. Diemoz, C. K. Hagen, M. Endrizzi, et al., “Sensitivity of laboratory based implementations of edge illumination X-ray phase-contrast imaging,” Appl. Phys. Lett. 103(24), 244104 (2013). [CrossRef]
28. B. De Samber, J. Renders, T. Elberfeld, et al., “FleXCT: a flexible X-ray CT scanner with 10 degrees of freedom,” Opt. Express 29(3), 3438–3457 (2021). [CrossRef]
29. P.-J. Vanthienen, J. Sanctorum, B. Huyge, et al., “Grating designs for cone beam edge illumination X-ray phase contrast imaging: a simulation study,” Opt. Express 31(17), 28051–28064 (2023). [CrossRef]
30. J. Sanctorum, J. Sijbers, and J. De Beenhouwer, “Dark Field Sensitivity In Single Mask Edge Illumination Lung Imaging,” in 2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI), (2021), pp. 775–778.
31. J. Sanctorum, N. Six, J. Sijbers, et al., “Augmenting a conventional x-ray scanner with edge illumination-based phase contrast imaging: how to design the gratings,” in Developments in X-Ray Tomography XIV, vol. 12242 (SPIE, 2022), p. 1224218.
32. J. De Beenhouwer, S. Staelens, D. Kruecker, et al., “Cluster computing software for gate simulations,” Med. Phys. 34(6Part1), 1926–1933 (2007). [CrossRef]
33. A. Przybylski, B. Thiel, J. Keller, et al., “Gpufit: An open-source toolkit for gpu-accelerated curve fitting,” Sci. Rep. 7(1), 15722 (2017). [CrossRef]
34. C. T. Chantler, “Theoretical Form Factor, Attenuation, and Scattering Tabulation for Z=1-92 from E=1-10 eV to E=0.4-1.0 MeV,” J. Phys. Chem. Ref. Data 24(1), 71–643 (1995). [CrossRef]
35. S. Brennan and P. L. Cowan, “A suite of programs for calculating x-ray absorption, reflection, and diffraction performance for a variety of materials at arbitrary wavelengths,” Rev. Sci. Instrum. 63(1), 850–853 (1992). [CrossRef]
36. J. Barzilai and J. M. Borwein, “Two-Point Step Size Gradient Methods,” IMA J. Numer. Anal. 8(1), 141–148 (1988). [CrossRef]
37. P. Munro, K. Ignatyev, R. Speller, et al., “Source size and temporal coherence requirements of coded aperture type x-ray phase contrast imaging systems,” Opt. Express 18(19), 19681 (2010). [CrossRef]
38. K. S. Morgan, D. M. Paganin, and K. K. W. Siu, “X-ray phase imaging with a paper analyzer,” Appl. Phys. Lett. 100(12), 124102 (2012). [CrossRef]
39. A. Astolfo, I. Buchanan, T. Partridge, et al., “The effect of a variable focal spot size on the contrast channels retrieved in edge-illumination x-ray phase contrast imaging,” Sci. Rep. 12(1), 3354 (2022). [CrossRef]
40. K. Amako, S. Guatelli, V. Ivanchenko, et al., “Comparison of Geant4 electromagnetic physics models against the NIST reference data,” IEEE Trans. Nucl. Sci. 52(4), 910–918 (2005). [CrossRef]
41. G. Cirrone, G. Cuttone, F. Di Rosa, et al., “Validation of the Geant4 electromagnetic photon cross-sections for elements and compounds,” Nucl. Instrum. Methods Phys. Res., Sect. A 618(1-3), 315–322 (2010). [CrossRef]
