Simulated grating-based x-ray phase contrast images of CFRP-like objects

Because of their microstructure, carbon fibre reinforced polymers outperform conventional materials in many aspects. These materials are, given the low x-ray absorption contrast they provide, ideally suited for phase contrast x-ray computed tomography. For system prototyping and for benchmarking and testing acquisition and reconstruction procedures, computer simulations are a valuable tool. In this paper, we propose the application of an alternative simulation strategy that unites Monte Carlo simulations with numerical wave optics and involves detailed modelling of the internal carbon fibre reinforced polymer structure down to the individual fibre scale.


Introduction
Carbon fibre reinforced polymers (CFRPs) outperform conventional materials in terms of high strength, elasticity, durability, and weight [1].Therefore, they are very promising in the search for new, cost-effective, function-oriented, highly integrated, and lightweight components.The macroscopic properties are defined by the inner microstructure of the fibre bundles, making detailed investigations at the microscale vital.As these materials provide in general very limited x-ray absorption contrast (AC), phase contrast (PC) x-ray computed tomography (CT) is a particularly suitable technique for investigating such samples [2].Grating-based interferometers (GBIs), based on the Talbot effect, are a powerful tool for acquiring such phase images [3,4].This method yields, in addition to the differential phase contrast image (DPC), also a conventional absorption contrast image and a dark field contrast (DFC) image.The latter relates to small-angle x-ray scatter caused by the sample's microstructure [5], and it has been shown that this type of contrast can be used to characterise the laminate structure of woven CFRPs, even if the individual fibres cannot be separated or the absorption contrast between the carbon fibres and the epoxy matrix is very low [6].When investigating new GBI-PCCT acquisition and reconstruction procedures, computational x-ray simulations are a valuable tool for benchmarking and testing.Previous efforts have already been put into x-ray absorption simulations specifically designed for CFRP-samples [7], but using a ray-tracing approach becomes insufficient when adequate physics modelling is of interest, as is the case for GBI-PCCT.Furthermore, assessing the small-angle x-ray scattering caused by the internal fibre structure requires a much more detailed modelling of both the phase and absorption related interaction physics and the internal structure of the CFRP, taking into account individual fibres rather than just the larger fibre bundles.Several approaches for phase contrast simulations, based on either numerical wave propagation [8,9], Monte Carlo (MC) simulations [10,11], or a combination of both [12,13] have been proposed during the past years.More recently, the combination of analytical and empirical input data has been explored to reduce simulation times [14].The latter approach, however, currently does not include dark field imaging.Here, we propose to adopt a phase contrast simulation strategy for CFRPs that unites MC simulations with wave optics [13] and involves detailed modelling of the internal CFRP structure down to the individual fibre scale.We review the principles of GBI-PC imaging and discuss the simulation framework, as well as our proposed approach for modelling the elementary properties of CFRPs.Subsequently, our first simulation results are shown and discussed.

Methods
In this section, we will first briefly address the imaging principles that apply to a grating-based interferometer in Section 2.1, followed by a discussion concerning the general simulation framework in Section 2.2.Finally, our approach for defining CFRPlike phantoms is discussed in Section 2.3.

Grating-based interferometer
In general, a grating-based interferometer consists out of a coherent x-ray source that emits x-rays towards a detector.On their path to the detector (see also Figure 1), the x-rays encounter two gratings [3].The first grating (G 1 ) causes interference to occur in the x-ray wavefront, generating an interference pattern at regular distances which has the same period as the grating (or an integer fraction thereof).This effect is known as the Talbot effect, while the distances at which it occurs are accordingly known as Talbot distances.Placing a phase-shifting object in the beam will distort the wavefront, causing the interference pattern to change.
9th Conference on Industrial Computed Tomography, Padova, Italy (iCT 2019) As the fringes of the interference pattern are in general not resolvable by the detector, these changes cannot be measured directly.Therefore, the second grating (G 2 ) is placed directly in front of the detector, which has a pitch that matches the periodicity of the undistorted interference pattern at that position.By shifting this grating perpendicular to both the optical axis and the grating lines, a procedure known as phase stepping, the portion of the interference pattern intensity that is being transmitted varies, hereby effectively translating the undetectable phase modulation into a detectable intensity modulation.When incoherent sources are used for imaging, an additional source grating (G 0 ) is required to allow the realization of a sufficiently strong interference pattern [4].

Simulation approach
The simulation pipeline can be split up in consecutive parts.The first part is based on MC calculations, where the x-ray interactions in the phantom are simulated (Section 2.2.1).After these interactions, the x-rays are collected in a discretised wavefront at a virtual plane behind the phantom (Section 2.2.2).In the last step, the constructed wavefront interacts with the gratings and is propagated towards the detector (Section 2.2.3).

Monte Carlo simulation
The MC tool that is used in this work is GATE, an open source application layer on top of the Geant4 simulation software, originally developped for biomedical purposes [15,16].Hence, the physics related to the absorption effects is already present in the existing framework.However, as described in [13], taking phase effects into account requires the photons to have an additional property (a phase or optical path) and the inclusion of a new physics process.In order to define these, each material in the simulation should be associated with its complex index of refraction, n = 1 − δ + iβ .Then, a phase φ can be associated with each photon p.When travelling through a certain material, the phase added to the photon's original phase equals the optical path length, where k p denotes the wavenumber and the integral of the decrement of the real part of the refractive index δ is calculated over the path followed by the photon.
At the boundary between one material (with refractive index decrement δ 1 ) and another (with refractive index decrement δ 2 ), the photon undergoes a deterministic refraction process, governed by Snell's law [10] (1 Here, the direction of the incoming photon in the first material with respect to the surface normal is given by θ 1 , and the direction of the refracted photon in de second material by θ 2 .
Tracking of the optical path length was done using the processGate library [17] and the refraction process was implemented in the GATE source code, hereby extending the existing physics library [18].

Wavefront construction
The MC simulation described in Section 2.2.1 results in a large amount of information about each photon's position, energy and phase.This information is subsequently used for coherently adding the photons in order to construct a discrete wavefront.This resulting wavefront is characterised by its energy and phase in each pixel.Since complex amplitudes are added rather than intensities, as one would do for standard absorption imaging, interference effects can be taken into account in the final part of the simulation, where the wave propagation and the interaction with the gratings is implemented (Section 2.2.3).
In practice, this is done by discretising the two-dimensional virtual plane in pixels of size ∆x × ∆y, creating a N x by N y grid.To coherently add the photons arriving at the same pixel, this information needs to be formulated in such a way that we can write this as an addition of complex amplitudes Ψ, rather than a sum of intensities, as is done in a conventional absorption-only approach.
Given that every photon p represents a position in the plane r, an energy E and a phase φ , we write p → p{r, E, φ }.Using this notation, we transform the particle representation to a local wave representation at position r: where The position r is used to determine to which pixel of the wavefront a certain photon will contribute.Making use of the wave representation of the photons, we can write the coherent addition of complex amplitudes Ψ n as: In this expression, D(x, y) denotes the complex wavefront amplitude resulting from the contributions of the N(x, y) photons arriving at pixel (x, y).If we write r = (r x , r y ) and define the continuous coordinates of the central point of pixel (x, y) as (x c , y c ), the condition for arriving at a certain pixel is given by where P(x, y) is the set of photons arriving at wavefront pixel (x, y).
All photon-tracking data generated with GATE is saved as a tree-structure in ROOT-format [19].To combine the data holding information about position, energy and phase of each photon as described earlier, the ROOT environment can be used.A ROOT-program written in the C++ language is applied to the tree-structure in order to perform the coherent addition of photon contributions.In order to achieve this in a practical way, we make use of the fact that the complex amplitudes can be written as Hence, the sum can also be formulated as: As can be seen from this expression, it is possible to split the sum over amplitudes into two separate sums.We thus store the real and imaginary part of the wavefront in two data vectors with the ROOT-program, rather than its energy and phase.

Wavefront propagation
The imaging system can be described as consisting out of three major parts.The first part is the phase grating G 1 , which introduces a phase modulation in the wavefront.Subsequently, the resulting wavefront propagates through free space towards the absorption grating G 2 , which forms the final part of the imaging system, before detection.This three-step imaging protocol is modelled in MATLAB.First, the binary files containing the real and imaginary vectors of the wavefront are loaded in MATLAB and reshaped to the orignal size of the wavefront.If A denotes the real part of the wavefront and B the imaginary part, the wavefront can be retrieved as D(x, y) = A(x, y) + iB(x, y).
An ideal G 1 phase grating can be constructed by using a superposition of two identical rectangular pulse train functions, of which one is laterally shifted by the grating bar width and multiplied by an exponential phase factor e iφ [20].This yields a grating that has 100% transmission over its full range, but introduces a phase shift of φ in a regular pattern.For Talbot-interferometry, this phase shift is either π or π/2.Mathematically, a phase grating with bar width b and pitch g 1 can be expressed as where * denotes the convolution operator, and The effect of G 1 on the wavefront can subsequently be simulated through point-by-point multiplication with the 2D extension of the G 1 (x) function.Straightforward implementation can be achieved through the built-in MATLAB function pulstran.
Propagation of a wavefront can, according to scalar diffraction theory, be described with the Huygens-Fresnel principle [21].
Assuming the Fresnel-approximation holds, the propagation of the wavefront can be described using a propagator kernel h as follows: This expression yields the complex amplitude of the wavefront at position (x 1 , y 1 ) after propagating a distance z.
The convolution kernel is given by where λ and k are the wavelenght and wavenumber, respectively, of the x-rays.For numerical applications, this analytical expression cannot be used directly and requires discretisation.This discretisation procedure involves sampling the convolution kernel h, and as implied by the convolution theorem [21], this can in principle either be performed in real space or in Fourier space.
Assuming sampling is performed in Fourier space, we introduce the discrete formulation of the propagation [22]: where DFT denotes the discrete Fourier transform and k x and k y are discrete sampling positions in Fourier space, with grid spacings ∆k x and ∆k y .The latter are defined as Computing this expression yields the new wavefront at propagation distance z.
For our numerical simulation, the effect of G 2 can be simulated in a very similar way as was done for G 1 .Moreover, since G 2 is an absorption grating, only one pulse train function is needed: The effect of this grating on the wavefront can subsequently be applied through point-by-point multiplication with the 2D extension of G 2 (x).Finally, for simulating detection of the intensity, a detector pixel size p d should be defined.Since in general p d ∆x, ∆y, the wavefront intensity |D(x, y)| 2 should be integrated over every detector pixel.As a result, phase stepping images can be generated by shifting G 2 .

Incorporating wavefront construction in GATE
The advantage of separating the calculations in Sections 2.2.1 and 2.2.2, is that the sampling parameters of the wavefront can be chosen independently of the MC simulation.This means that after performing the MC simulation, there is still free choice in wavefront resolution and number of contributing photons.The output size of the GATE simulations then depends on the number of photons, however, and thus increases rapidly when more photons are generated.
To overcome this limitation, at the cost of the aforementioned flexibility, we incorporated Section 2.2.2 in Section 2.2.1, by further extending the GATE source code [18].Furthermore, this reduces the overall simulation time, as we replace the three-step process by a two-step process, where we avoid looping over all photons for constructing the wavefront, as they are immediately added to the front in GATE upon arrival at the virtual plane.
In the latter case, the wavefront resolution and number of contributing photons are fixed for a certain MC simulation.Changing these parameters thus requires a new simulation.
Figure 2: High-resolution synchrotron scan of a CFRP sample.It can be seen that the fibres in a given bundle have a nearly identical orientation and are to a good approximation cylindrical in shape.Scans were taken at Elettra Sincrotrone, Trieste.

CFRP-like phantom model
In general, CFRPs consist out of many thin (<10 µm) carbon fibres, organised in bundles that are embedded in a resin of a different material, often epoxy.In order to take the scattering effects caused by the microstructure of the fibre bundles into account during the simulation, the phantom should be defined at the level of the individual fibres.In this section, the proposed approach for designing CFRP-like phantoms is discussed.The primary goal of our work is not to develop large-scale phantoms with realistic overall dimensions, but rather to build simple phantoms that still exhibit the most elementary properties of CFRPs.
The phantoms should therefore consist out of relevant materials, e.g.carbon for the fibres and epoxy for the resin.Within GATE, phantoms can be defined using basic geometric shapes, such as spheres and rectangular boxes.A fibre bundle can be modelled as a collection of solid cylinders, as can be seen from high-resolution synchrotron scans (Figure 2).

The number of fibres in a bundle
The question now rises how many fibres there should be in a given fibre bundle.Starting from the simplified model of a cylindrical fibre bundle, we derive the following formula for the number of fibres in a given bundle, based on geometrical considerations: Here, R is the radius of the cylindrical fibre bundle, r f the radius of a single fibre, ρ f the density of the fibre material, ρ r the resin density, and c f the carbon fibre weight fraction.In the more general case of an arbitrary bundle shape, we can still write the number density ρ N of fibres in the cross section as as this value should be independent of the shape of the cross section.Using equations ( 16) and ( 17), we can determine the number of fibres that need to be defined in GATE.

The arrangement of fibres in a bundle
Arranging a certain number of parallel fibres in a cylindrical bundle can be done homogeneously by applying the so-called sunflower seed arrangement [23].An example of such an arrangement for 330 fibres is shown in Figure 3.
A more irregular, but uniformly distributed arrangement of fibres can be achieved through Poisson disk sampling [24], which produces points that are tightly-packed, but no closer to each other than a specified minimum distance (the fibre diameter, in this case), resulting in a more natural pattern.Such a pattern is also shown in Figure 3.
In case of Poisson disk sampling, it is more straightforward to work with the fibre number density rather than an explicit number of fibres.The number density can be used to determine the number of fibres in a rectangular area.By definition, the distribution of fibres in this area is uniform, hence each (sufficiently large) arbitrary portion of this area should have the same number density.Fibre bundles with a custom cross sectional size can thus be extracted from this randomly generated rectangular field.
The choice for either of these approaches depends on whether more homogeneity or more randomness is desired.Defining thousands, or even hundreds, of fibres one by one to create a phantom is an extremely tedious and undesirable task.Therefore, the fibres are generated in GATE with built-in, continuous geometrical building blocks using the repeater function.This function allows generating an arbitrary number of copies of a given object, at positions and orientations provided in a placements list.Following this approach, only one explicitly defined fibre is required to generate thousands of copies, resulting in a bundle.In addition, as only the original object is stored in memory, generating bundles from fibres comes at no extra memory cost.

Demonstration and discussion
As an example of the proposed approach, we present a phantom representing the most fundamental building block that characterises a CFRP sample, namely a pair of crossing carbon fibre bundles embedded in an epoxy resin (see Figure 4).A solid carbon cylinder the same size as the bundles has also been included as a reference.The total field of view is 3.6 mm x 2.0 mm, and the fibre bundles were defined as having a radius of 250 µm and lengths of 2.0 mm (horizontally oriented fibre) and 1.5 mm (vertically oriented fibre), respectively.Given a single carbon fibre radius of 3.5 µm and a weight fraction of 60 percent, this resulted in 2126 simulated fibres per bundle, assuming a sunflower seed arrangement.We simulated parallel beam x-rays at 25 keV, for a total amount of approximately 8 • 10 9 photons in GATE [25].With respect to the wavefront propagation, we assumed a grating pitch of 4 µm for the phase grating and a phase-shift of π, yielding a grating pitch of 2 µm for the absorption grating [26].To further mimic the experimental practice, phase-stepping at a 20 µm pixel size detector was performed, followed by a phase-retrieval algorithm [27] used to extract the AC, DPC, and DFC images.Since the grating lines were vertically oriented, the simulated setup is only sensitive to phase effects occurring in the horizontal direction, causing the horizontal bundle to vanish in the DPC and DFC image.Indeed, no phase gradient is present in the horizontal direction and the scatter directions are perpendicular to the main axis of the fibre bundle.Only the edges of this bundle are visible in the DFC image, as the local reflection of x-rays causes a drop in visibility.Furthermore, the visibility drops in the vertically oriented bundle because of the x-ray scatter that occurs at the fibre bundles, yielding problems with retrieving the correct phase (black and white lines), which has also been observed in experiments [27].
In the fibre bundle regions, vertical and horizontal lines are visible in the AC image.The vertical lines are also visible in the DFC image (recall that the horizontal bundle cannot be detected).These lines are most likely due to the internal structure of the fibre bundle showing through.Indeed, the fibres, having a diameter of 7 µm, are parallel to the 20 µm pixels of the detector, resulting in the overlay of different patterns.This overlay effectively leads to the occurrence of a moiré-like effect.
In Figure 5, experimental images of a CFRP sample acquired with a grating-based interferometer are shown.From the grainy structure of the DPC image, it is clear that the strong dark field signal, originating from scattering, is an indication for errors in the determination of the phase.Often, these artefacts can be reduced by preprocessing the images [28].
These observations indicate that there is a qualitative agreement between simulation and experiment and that, although additional polishing is needed, some of the key properties of CFRP samples can already be included in simulations with the proposed approach.
Figure 5: Actual DFC (left) and DPC (right) images of a CFRP.The DFC image has been log-corrected, such that a stronger dark field signal corresponds to a higher grey value.Images were acquired with a Bruker SkyScan 1294 phase contrast micro-CT system.

Conclusion and outlook
In this work, we proposed a simulation approach for grating-based phase contrast imaging of CFRP-like phantoms using GATE.
We have shown the general ability of the approach to reproduce the most important effects that play a role in phase contrast imaging of CFRPs.
For future work, it would be of interest to investigate different fibre configurations, where the fibres are for example no longer perfectly parallel.This, in combination with a larger detector pixel size, could furthermore reduce the observed moiré effect.
In addition, the simulations should be extended beyond the model of the cylindrical fibre bundle.This can for example be done using Poisson disk sampling in a rectangular box, from which an arbitrarily shaped area with fibre positions can be extracted.

Figure 1 :
Figure 1: Schematic overview of a grating-based interferometer (not to scale).

Figure 3 :
Figure 3: Schematic overview of 330 fibre positions, arranged in a sunflower seed configuration (left) and a configuration resulting from Poisson disk sampling (right).

Figure 4 :
Figure 4: Top left: the simulated phantom, the epoxy resin is made transparent to show the carbon fibre bundles.Top middle: zoomed-in image.Top right: actual orientation.Bottom: final images, from left to right: absorption, dark field and (differential) phase contrast.