Open Access
Add to BibSonomy
Abstract
We demonstrate two versatile, flexible, and accurate frameworks based on numerical and Monte Carlo approaches to simulate the X-ray speckle-based (SBI) technique for lab-based systems. The established tools can reproduce experimental setups in a cone-beam geometry and with polychromatic sources. Furthermore, they are computationally efficient to enable a fast virtual multi-modal tomography of digitized inhomogeneous phantoms. The proposed methods were evaluated and validated by analytical and experimental data for various samples. The Monte Carlo approach provides a realistic and accurate simulation, which is useful in diffuser design and dosimetry studies, while the numerical method is very efficient for parametric and tomographic studies. These approaches will be used for the optimization of lab-based X-ray SBI setups and generating sample images for enhancing phase retrieval algorithms.
Published by Optica Publishing Group under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
1. Introduction
A beneficial upgrade to well-known attenuation-based X-ray imaging are phase-contrast techniques, which can significantly improve the detection of material features with low attenuation contrast [1]. The advantage originates from the fact that these imaging methods are additionally sensitive to the phase shift of X-ray wave front upon refraction and diffraction from sample structures. Particularly, they are able to differentiate small electron density changes in low-contrast materials, which consequently reveals additional information about the microstructures [2].
X-ray speckle-based imaging (SBI) is a multimodal imaging technique, which can measure the first derivative of the phase in addition to the intensity of the X-ray wave front [3]. Moreover, it enables to extract an ultra-small-angle scattering signal, so called dark-field, to measure higher-order modulations in the wave front [4,5]. This method takes advantage of a simple setup with only an additional random mask, i.e. diffuser, between the source and detector. The diffuser (speckle mask) can be any suitable material with porous or granular texture, as depicted in Fig. 1. This mask acts as a wave modulator and creates near field speckles on the detector plane. Any change of the wave front induced by a sample can be detected through tracking the distortion of the speckle pattern, and then, the corresponding refraction angle or the related phase gradient can be calculated at each image point [6].
Fig. 1. Schematic of the X-ray SBI setup. A micro-focus X-ray tube and an aSi-CsI flat panel detector were set up with diffuser in between, positioned downstream the object.
Many interesting applications have been advised for this relatively new technique in material science, metrology, and biomedical sciences [7–9]. Since the SBI technique can be implemented efficiently in a single shot, also called X-ray speckle-tracking (XST), it is appealing for the tomographic imaging, combined with proper reconstruction algorithms [10–12]. As such, another potential application among others would be the tracking of dynamic processes in living systems [13].
A practical advantage of SBI technique is the ability to work with broad-spectrum laboratory sources with no demand for complicated and expensive optical elements, unlike other existing phase-sensitive X-ray tomography methods [14–16]. Since the diffusers are typically low absorbing, SBI also uses X-ray photons more efficiently compared to grating-based techniques or non-interferometric methods [14,17]. Regarding the two above practical merits, this technique is very promising specifically in a laboratory, which can pave the road for several medical and industrial applications.
A few studies have been conducted to investigate the effect of some SBI system parameters, however, real systematic and parametric assessments are still missing to fully optimize the technique for laboratory sources [18–21]. Therefore, a reliable and efficient simulation toolkit that enables one to reproduce experimental imaging systems as accurate as possible while keeping the computational time reasonably low would be highly beneficial [22]. Such tools can be employed at pilot studies before spending real-world resources for the design and optimization of an experimental setup in a synchrotron or laboratory [23]. It would be also very useful to employ it in comparing different tomography reconstruction algorithms or phase retrieval methods [24]. Moreover, one can utilize virtual imaging simulations to try new ideas on different operational modes, predict possible experimental challenges, and validate experimental data [25]. Furthermore, such simulated images can be used for improving iterative methods and training image processing algorithms faster and more efficiently [26]. A large database of simulated images can be produced for a specific application instead of making a significant effort to acquire experimental database, which is more expensive and less feasible.
There are some packages designed to simulate X-ray imaging experiments on synchrotrons, such as SRW [27], SHADOW [28], Phase [29], syris [22], McXtrace [30], OASYS [31], usually in propagation-based imaging (PBI) systems. These packages include complex modules to specifically simulate the generation of X-rays inside a synchrotron and all relevant components (like mirrors, monochromators, lenses, etc.) and create a parallel beam with expensive computations in ray tracing and/or full wave-optics approach. Therefore, those packages are not suitable and efficient to simulate divergent, partially coherent laboratory sources.
In this study, we present an efficient virtual imaging framework including two numerical and statistical approaches, specifically to simulate image formation in laboratory-based X-ray speckle-based imaging (SBI) systems. Both approaches are proven to be computationally efficient and have specific benefits, depending on the simulation purpose. They enable a large flexibility regarding the simulation of different aspects of SBI configurations for a broad range applications. To the best of our knowledge, Zdora et al. [18] presented the only simulation of SBI technique but limited only to parallel beams and 2D projections. They used near-field Fresnel diffraction approximation to numerically evaluate the effect of energy spectrum and beam hardening on their radiographs from a homogeneous single-material cuboid. Such evaluation of the diffraction integral can be challenging and numerically expensive for more realistic setups and objects, and not suitable for parametric optimization. Moreover, in a very recent article [32] that was published while the current paper was under review, a comparison of possible simulators for SBI setups was reported based on three different approaches (ray tracing, Fresnel wave propagation, and Monte Carlo method). Quénot et al [32]. discussed the advantages and disadvantages of these approaches in a projection imaging setup, such as ray optical validity in the ray tracing method, and aliasing artifact (sampling error) and computational cost in the Fresnel wave propagation approach. However, the authors don't suggest any solution to address these limitations and enhance the computation accuracy and throughput for various samples.
In this work, together with a Monte Carlo method, we demonstrate an accurate but more efficient numerical approach coupled with a multi-voxel-slice calculation to overcome some of these limitations, as discussed in the next paragraph. MC approach is used to benefit from its complementarity to address real-world effects like incoherent scatter that is not captured in the numerical approaches. We will go further and discuss how the MC method can provide more accurate results compared to the numerical approach. Nonetheless, we also provide a proof-of-concept for a virtual lab-based tomography from a digitized phantom with the proposed numerical method.
Most existing wave propagation forward models for x-ray imaging are limited to single-material low-absorbing objects [33]. Therefore, another problem that we address in this work is to propagate the cone-beam (CB) X-rays through inhomogeneous, thick objects with spatially strongly varying refractive index in a SBI system. Shanblatt et al. [34] suggested a multi-slice (MS) approach based on the Fresnel diffraction integral to consider the diffraction effects from inside the objects for the PBI technique with monochromatic beams. Simulations of the full Fresnel integral can be challenging in terms of choosing the appropriate sampling rate for the phase chirp function, which otherwise can lead to aliasing [35]. Moreover, as mentioned before, its computational cost makes it inefficient in many practical cases with polychromatic sources.
To avoid these problems, here, we propose a MS approach but based on transport-of-intensity equation (TIE) that gives a good approximation of the Fresnel integral for near-field propagation [36,37]. Using this numerical approach, we also adopt a higher resolution (for steps) in numerical calculations by choosing a voxel slice in our MS approach. Our proposed methods and validated simulation examples address more general aspects of a lab-based SBI configuration, including cone-beam geometry and polychromacity, and are applicable to inhomogeneous objects as well.
Here, we describe the simulation methods and their implementation. We will then evaluate the performance of two proposed approaches by comparison with analytical and experimental data using various objects to validate the tools for further optimization studies.
2. Simulation scope and implementation
In this section, we describe two methodologies and simulation implementation.
2.1 Monte Carlo simulation
Geant4 is an open-source MC simulation toolkit based on C++ programming language to transport various types of sub-particles and particles including but not limited to photon, electron, neutron, and ions inside the matter [38]. The program benefits from comprehensive interaction cross-section libraries for radiations and particles. This toolkit allows the user to define various types of geometries with desired material properties and to track the particles inside a phantom or detector. It also enables to design the source shape and radiation spectrum for a realistic simulation. In this study, MC simulation was performed through routines in the Gate platform (v.8.2) [39], a user-friendly toolkit based on Geant4.10.2 libraries. Gate has provided dedicated modules for radiation transport and scoring for medical applications from computed tomography to nuclear medicine, radiotherapy, and dosimetry.
In Gate routines, the source was defined by the general particle source (GPS) module according to the desired X-ray source properties. The monochromatic beam was made as an isotropic point source shape. For the polychromatic source, an elliptical focal spot of 2 × 2 µm was shaped. An isotropic and divergent beam was built up so that its solid angle covers the whole desired field of view (FoV) to efficiently sample emitted photons and decrease the computation time. The energy of photons was sampled based on a user-defined spectrum in 100 energy bins, calculated from the XrayTool simulator [40] for the relevant laboratory micro-focus X-ray tubes (see sections 3.3 and 3.4).
Detectors were simulated with different properties for each example (section 3). The detector was defined as a general scanner with a GateBox geometry in which all physical layers according to the detector type and characteristics were implemented. The scintillation layer was attached to a sensitive detector in appropriate voxel meshes (see section 3). The readout was performed in the layer of amorphous silicon via a digitizer module. This module enables one to introduce a spatial blurring or pixel cross-talk for a given detector and define energy resolution in the case of a counting detector. It also allows to set an energy threshold window for the detector. Nonetheless, it is possible to introduce various types of electronic and dark noise into the readout system in this module, which was not employed in our simulation. A DoseActor was formed to also measure the energy deposited in each detector voxel.
The Geometry of spheres and cylinders were modelled through GateSphere and GateCylinder volumes, respectively. The diffuser was designed as silicon carbide (SiC) spherical grits at predefined diameter ranges (section 3.2), uniformly distributed inside a cubic cellulose substrate through a random generator. Materials were set up and listed in a GateMaterials.db file. Through a Penelope physics model, all possible interactions of photons and electrons within the source energy range were considered in the physics list, including photoelectric, Rayleigh (coherent), and Compton (incoherent) scattering, Bremsstrahlung, multiple electron scattering, ionization, and fluorescence/Auger emission. The production threshold for photon and electrons were set as 250 eV and 1 keV, respectively, in all volumes.
Since Geant4 code considers the particle nature of X-rays to transport them in the matter, the wave nature of X-rays such as refraction and diffraction effects in the projection approximation should be included in the source codes [41,42]. To do so, a refraction process was created and added to the electromagnetic processes in the Geant4 physics modules. This new process indeed applies the Snell's law the geometric optics to calculate the refraction of photons through the matter [41,42]. The relevant dependencies were generated to correctly initiate the refraction process and change the momentum vectors of photons when each photon reaches a boundary of two materials.
The maximum computation steps in the object, diffuser, and detector volumes were set at 2 nm so that a sufficient precision (uncertainty < 10%) in all detector pixels is achieved while the computation is kept within a reasonable timeframe. The deposited energy inside each crystal (or voxel) element is recorded together with other parameters like the source and detected position of each particle, time of flight, energy, etc. in an ASCII file. This file is further processed by MATLAB scripts to retrieve the multi-modal contrast images. Several signal extraction methods have been demonstrated for the X-ray speckle-based imaging [10–12,43]. We used an implicit method to retrieve refraction and phase signals from the SBI simulations [12]. For calculation of the dark-field signal, a cross-correlation algorithm [10] with a 9 × 9 subset window was utilized.
2.2 Numerical simulation
We apply a numerical wave-propagation approach based on the transport-of-intensity equation (TIE) [36,44,45]. The TIE has been shown to be a good approximation to the Fresnel diffraction integral for paraxial waves in small propagation distances [36,37]. For propagation in matter, the so-called projection approximation is typically applied for a thin object or objects with slowly varying properties in space [36,46]. In other words, the projection approximation is valid as far as the diffraction inside the object can be neglected during the propagation through the matter.
To address this problem and make our approach more generic, we opted a multi-slice treatment [34,46,47] and modified it to propagate the intensity through inhomogeneous, thick objects. In our proposed method, we cut the sample along its thickness to thinnest possible slice, i.e. a voxel depth, and apply the TIE within each voxel slice to propagate the intensity through the whole volume.
At each point in the 3D matrix of the object and diffuser (diffuser) geometry, a complex refractive index $n\,=\,1 -\delta+ i\beta$ was assigned to create the spatial distribution of the material throughout the geometrical volumes. At the exit surface of each voxel slice, the intensity ($I$) is attenuated and the phase ($\Phi $) is shifted, as described by the projection approximation in (1) and (2), respectively:
For the cone-beam geometry, the propagation process is formulated based on the first-order Taylor expansion of TIE [48,49], scaled according to the Fresnel scaling theorem :
In a general view, the intensity is updated upon applying the propagation, followed by attenuation and phase-shift processes within each slice, in which the sub-voxel interpolated refractive indices at each slice are used (Fig. 2). This procedure is sequentially repeated for all voxel slices through the whole object/diffuser thickness.
Fig. 2. Schematic of the multi-voxel-slice calculations in a cone-beam geometry. The 3D refractive index matrix of the object/diffuser is resampled as a function of transverse position (x, y) for each slice, starting at the entrance surface of the geometry. The resampling is done by sub-voxel interpolation of β and δ between normal matrix points (black dots at corners) and the updated red points, upon magnification at each further slice.
To introduce the polychromacity as well as energy dependence of the detector response, the simulated X-ray spectrum was weighted by the calculated energy response function of that specific detector. This weighted spectrum was divided into 100 energy bins. A MATLAB function linked with Xraylib libraries [50] calculates the complex refractive index of every voxel in the sample/diffuser at each energy bin. The main TIE-MS simulation program uses these refractive indices to compute the intensity and phase of the X-rays propagating through each voxel of sample or diffuser (Eqs. (1) and (2)). Then the simulations carried out at different energy bins are integrated over the whole weighted energy spectrum. On the detector screen, the final image is convolved with a Gaussian point spread function (PSF) of the flat panel detector to account for the finite focal spot of the micro-focus tube.
The noise can be simulated through adding a normally distributed electronic noise, a uniformly distributed quantization noise, and a Poisson-distributed quantum noise [22]. Here, we only introduced a quantum noise for the polychromatic examples in our numerical simulations. The multi-contrast signals were retrieved by the same methods as used for the experiments and MC simulation (see section 2.1).
2.3 Analytical cone-beam simulation
To evaluate the accuracy of the CB geometry implementation through both proposed simulation methods, we conducted a comparison study with analytical calculations for a spherical object first in a PBI setting. The analytical solution for the magnified sphere's thickness can be obtained by:
3. Simulation examples
Four objects were simulated through the implemented TIE-MS and MC methods, as described below. For the last two examples, experiments were performed in a laboratory setup.
3.1 Monochromatic-PBI simulation of Al sphere
The first example is a homogeneous Aluminum sphere in a PBI setup for which the accuracy of the image formation and CB geometry implementation is evaluated by comparing to an analytical CB calculation (section 2.3). The Aluminum sphere with a diameter of 500 µm was placed 8 cm from a monochromatic point source at the energy of 17.5 keV. The detector was defined as a digital flat panel detector with a 500-µm CsI scintillator with a 25-µm pixel element. The detector was positioned 1.017 m away from the source and were meshed in 400 × 400 voxels.
3.2 Monochromatic-SBI simulation of POM sphere
The second example is a polyoxymethylene (POM) sphere with a diameter of 200 µm which was simulated in a SBI setup with a monochromatic point source at 10 keV. The sphere was placed 1 cm from the source and a diffuser was positioned 13 cm downstream the sample. The diffuser was modeled as a cellulose substrate containing 6000 random SiC spheres with diameters at the range of 34-136 µm. The same number and diameter range of spheres and distribution model was used to produce similar diffusers in both TIE-MS and MC methods. A Gadox scintillator detector with 20-µm thickness, coupled to a CCD camera was simulated at 1 m distant from the source. The detector pixel size was 100 µm in a 256 × 256 matrix in this example.
3.3 Polychromatic-SBI of PMMA cylinder
A PMMA cylinder with a diameter of 4 mm was simulated in a SBI setup with a laboratory micro-focus X-ray tube (VISCOM X9160-TXD). A Perkin Elmer (XRD 1621CN3ES) digital flat panel detector with a 500-µm-thick CsI-scintillator layer was modeled at 1.017 m away from the x-ray tube, as in the experimental setup. The image matrix was 360 × 360 with 200-µm for each pixel. The diffuser was modeled according to the sandpaper (SP120) through the method described in section 3.2, but with 12000 SiC spheres with an average diameter of 126 µm. For validation purposes, an experiment was performed with the aforementioned X-ray source and detector and similar geometry.
3.4 Virtual microtomography of an inhomogeneous phantom
Here, we demonstrate a virtual SBI micro-tomography of an inhomogeneous phantom. To do so, a POM cylinder (diameter = 4 mm), a hollow carbon-fiber cylinder (diameter = 2 mm), and a hollow high-density polyethylene (HD-PE) tube (diameter = 2.65 mm) were placed inside a polystyrene (PS) cylindrical container with the thickness of 1 mm. This phantom was scanned in the EasyTom XL Ultra 160 (RX Solution) microCT scanner in a normal tomography setup. The micro-focus source was a high power Hamamatsu L10801 X-ray tube operated at 50 kVp and 300 µA and placed 85 cm distant from a CsI flat panel detector with a pixel size of 127 µm. The phantom was positioned about 13.9 cm away from the source and 1088 projections were acquired in total. Then the diffuser (SP120) was placed in 19 cm distant from the source and the phantom was scanned in an XST tomography procedure with the same scan parameters.
For the simulation, the volume geometry of the phantom was reconstructed from the first tomography experiment and then was segmented to obtain a 3D digital phantom. The known material information was embedded into the virtual phantom, which was then fed into the simulation program. The energy spectrum for the micro-focus source was calculated by Xraytool software. The XST microtomography setup together with the digital phantom was simulated through our proposed TIE-MS method on 12 cores in parallel. Multi-contrast images were extracted via the same retrieval methods in section 2.1 by MATLAB scripts. The volume images were calculated via an in-house tomography reconstruction algorithm for both experimental and simulated XST microCT.
4. Results and discussion
4.1 PBI simulation and analytical evaluation
Projections of the Aluminum sphere simulated in a CB PBI setup are shown in Fig. 3(a-c). In this figure, the proposed approaches are compared with the TIE-analytical simulation of the sphere in the CB geometry. A number of 4 × 108 photons were transported in the MC simulation in this example.
Fig. 3. Propagation-based phase-contrast images of a 500-µm Al sphere simulated by the proposed (a) TIE-MS, (b) MC, and (c) the TIE-analytical methods. The propagation distance was 93.7 cm. (d) Mean intensity profiles (yellow dash region) of images (a), (b), and (c). (e) Retrieved phase line profiles of three simulations.
The intensity profile of the sphere was averaged over 60 vertical pixel columns (yellow dash rectangle) to give the mean intensity profile plotted in Fig. 3(d). This graph indicates a good conformity between the TIE-analytical results and both TIE-MS and MC simulations. As can be seen in Fig. 3(b,d), the edge signal has been relatively decreased in the MC image. It is mainly due to the incoherent scattering from the sample/diffuser, activated in the MC physical processes, which we believe makes it a more realistic imaging. This leads to some blurring in the formed image which is more evident at the edges. The MC image also consists of a quantum noise that was not added to the analytical and numerical images. The phase images were extracted by the Paganin's method for homogeneous, single-material objects [48] and their central horizontal line profile are plotted in Fig. 3(e). This plot verifies a very good agreement between the theoretical phase shift and the extracted phase images from three simulations based on the normalized root mean square error (NRMSE), shown in Table 1.
Table 1. Speckle visibility and NRMSE values for simulated objects in various setups.
The results of simulations for this homogeneous sample compared with the TIE-analytical calculation demonstrate that the CB geometry in the TIE-MS and MC methods have been implemented correctly. On the other hand, the comparison of phase images with the calculated theoretical phase proves the accuracy of the PBI implementation in both simulation approaches. Thus, the proposed methods can be applied for arbitrary shapes with inhomogeneous materials and unknown analytical solutions. Thus, the proposed methods can be applied for arbitrary shapes with inhomogeneous materials and unknown analytical solutions.
4.2 Monochromatic-SBI simulation
Figure 4 shows the results from the simulations of a POM sphere in a cone-beam SBI setup with a monochromatic beam. The multi-contrast images were retrieved through the Pavlov's method [12] as well as a cross-correlation algorithm [7] that are also depicted in Fig. 4. The speckle visibility was calculated from the speckle image by Eq. (6) [18], and reported in Table 1.
where ${\sigma _{speckle}}\; $and ${\bar{I}_{speckle}}\; $are the standard deviation and mean pixel value of an 150 ×150-pixel region of interest (ROI) on the speckle image. The number of photons run in the MC simulation was 4 × 108 that led to an uncertainty values less than 5% in all pixels. The speckle noise resulting from the statistical nature of the MC method can be reduced by using a higher number of initial photons and setting appropriate variance reduction techniques.Fig. 4. X-ray SBI images of a POM sphere simulated by TIE-MS (top) and MC (bottom) method and the retrieved multi-contrast images:(a) horizontal refraction angle, (b) vertical refraction angle, (c) integrated phase shift, (d) single-shot projection from the object, (e) visibility (dark-field), (f) transmission, (g) thickness, (h) speckle image without sample.
The profiles of horizontal refraction angle, phase, and thickness from both simulations are compared in Fig. 5. The retrieved thickness and phase from the projections show a good accordance between the simulations and theoretical values as implied by the NRMSE values shown in Table 1. These results manifest a correct implementation of SBI physics through the proposed methods. We also investigated how excluding Compton scattering from MC simulations can affect the simulation results. As we expected, for this monochromatic example at 10 keV, we didn't observe a significant effect on the result, because at this energy, the Compton scattering cross-section is not significant and therefore, it is not an important interaction here.
Fig. 5. Intensity profiles along the red dash line in Fig. 4-(a): (a-c) are corresponding to (a), (g), and (c), respectively, in Fig. 4.
4.3 Polychromatic-SBI simulation and experimental validation
As the next example, the simulated and experimental radiographs of a PMMA cylinder in a polychromatic CB SBI system are shown in Fig. 6. The retrieved transmission, horizontal and vertical refraction angle, absolute integrated phase shift, and the visibility signal are presented in Fig. 7. The horizontal profiles of some images in Fig. 7 are averaged vertically from the dashed yellow rectangle and plotted in Fig. 8.
Fig. 6. Single-shot XST images without (top) and with (bottom) a PMMA cylinder: (a) Monte Carlo simulation, (b) numerical TIE-MS approach, (c) experiment.
Fig. 7. Retrieved multi-contrast signals for the PMMA cylinder: (a) horizontal refraction, (b) vertical refraction, (c) visibility (dark-field), (d) transmission, (e) total absolute phase-shift.
Fig. 8. Horizontal mean profile (yellow dash contour) of retrieved signals from XST of the PMMA cylinder: (a) horizontal refraction, (b) vertical refraction, (c) transmission, and (d) total phase shift.
Note that in Fig. 8.b, there are residual peaks in the line profile of the vertical refraction angle for the experimental data that is expected to show a flat signal. This is simply due to the small inclination of the PMMA cylinder in the experiment. This tilt (about 1.5°) has created tiny peaks in the vertical refraction signal (±53 nrad), which has also decreased the horizontal refraction signal (Fig. 7.a). Following the experiment, we applied the same tilt to the simulated objects to account for this geometry, as evident in Figs. 6–8. According to Fig. 7 and Fig. 8, a good agreement is confirmed between simulations and the experiment, specifically for the MC method. The NRMSE of the transmission and phase signals have been reported in Table 1.
The main reason for some differences between the polychromatic experimental and simulation results is the very likely difference between the useful energy spectrum from the simulation and real x-ray tube, including all relevant attenuations and scattering from the tube window and other close components. Another source of small deviations from the experimental result could be that the detector energy-dependent response was simulated very approximately and can cause inaccuracies.
We inactivated the Compton scattering (CS) in the MC simulation to investigate if it contributes in the difference between two simulation methods. According to our results, the NRMSE of the MC results were increased substantially for all multi-contrast signals without the CS interaction (Table 1). This suggests that the Compton scattering should be considered, specifically at higher energies where the CS cross section increases, like in the case of this example. At a closer look, the refractive index calculated by Xraylib in TIE-MS considers the attenuation and scattering cross sections, however, the contribution of scattered photons that are absorbed inside the phantom or those not reaching the detector are not considered. Hence, this leads to about 12% lower signal in TIE-MS than MC simulation, in which these scattered photons are correctly transported.
However, considering all these uncertainties caused by a laboratory polychromatic source, both numerical and MC simulations provide a reasonably accurate simulation framework for qualitative and quantitative SBI.
Here, a total number of 4 × 108 photons were simulated in the MC calculation for each projection image, which lasted about 2 h on a single core (IntelR coreTM i7-9850H CPU 2.6 GHz processor) or 20 min in parallel mode on all cores of the same system. The time can be reduced proportionally by including more processing cores or on a GPU. The computation time for the numerical TIE-MS method with 10-point sub-pixel interpolation was 2.2 min for the whole spectrum for a 360 × 360 × 254 object matrix and 360 × 360 × 36 diffuser matrix, run on the same processing cores. The numerical method reveals a more efficient and fast performance for our example.
Generally speaking, the final computation time for the numerical implementation considerably depends on the volume image size, the number of slices, and the number of energy points. This dependency is much more relaxed in MC simulations so that it has no significant effect on the computational cost as far as it doesn't decrease the average number of photons per pixel to make the quantum noise intolerable. The run time in the MC simulation is mostly determined by the number of transported photons needed for an acceptable signal-to-noise ratio (SNR), the step size in stepping calculations, and variance reduction techniques for higher output efficiency and accuracy. It should be noted that no extra effort has been yet made to further optimize the memory and computation efficiency on either simulation methods.
4.4 Virtual microCT with inhomogeneous phantom
The simulation time for a 951 × 101 × 910 voxelized phantom was 180 s per projection for 10-bin energy spectrum. The multi-modal projections of the digitized phantom made of different polymer cylinders were retrieved from an XST microtomography simulation and actual experiment as well. The speckle visibility was 23% and 31% for the experiment and TIE-MS simulation. The tomographic reconstruction was performed on the retrieved images by the Feldkamp-Davis-Kress (FDK) implementation of the filtered back projection method, using the Xact (RX SolutionTM) and an in-house software. In Fig. 9(b-c), the reconstructed transmission and phase slices of simulated volumes (bottom row) are compared with those from the experimental XST microCT (top row). The histogram of an experimental and simulated projection (Fig. 9-(a)-top/down) is shown in Figure S1 in Supplement 1.
Fig. 9. Experimental (top) and virtual (bottom) microtomography of an inhomogeneous phantom. (a) Projection image of the phantom, (b) a reconstructed slice of the phase volume (δ is the real component of the refractive index), and (c) the transmission volume.
Figure 9 demonstrates the capability of the proposed numerical approach in reproducing the experimental X-ray speckle-based tomography setups within a reasonable computational time. Furthermore, this example proves that the SBI tomography of any arbitrary inhomogeneous 3D digital phantom can be simulated in this numerical framework. For the purpose of quantitative imaging, it is sufficient to know the material of the segmented regions in the real phantom or tissue to design a digitized phantom and simulate it in the desired X-ray SBI experiment.
5. Conclusion
Two model-driven simulation tools, one statistical and one numerical, have been developed for the simulation of lab-based X-ray SBI systems. The former method is based on an embedded geometric optics into a full Monte Carlo calculations and the latter method is established based on the TIE coupled with a multi-voxel-slice technique for simulation of samples with strongly varying X-ray optical properties. The tools were evaluated and validated for SBI and PBI settings in a cone-beam geometry with polychromatic X-ray spectra. The demonstrated examples are the proof-of-concept of the capability of the proposed approaches in accurate and efficient simulation of laboratory experimental setups.
The proposed MC framework can be linked with CT modules designed in Gate platform for a comprehensive simulation of various phantoms. Additionally in this platform, one can design a specific readout system for the imaging detector, which would be useful for experimental tests. The program can incorporate digital voxelized phantoms and simulate a complete tomography system. However, applying the refraction process to a voxelized phantom in Geant4/Gate comes with extra challenges to modify the source codes because the particle transport method in a voxelized phantom mode is different from a normal designed volume and should be addressed in future studies. The computation time in the MC framework in some cases is longer than for the numerical approach but through parallelization on a GPU, it can be considerably shortened. For the MC method, an advantage is that the run time is not dependent on the thickness of the sample, the number of image pixels, or the number of energy points in the spectrum. Furthermore, we have shown that the MC approach captures e.g. incoherent scatter contribution in the signal and, therefore, has a higher accuracy to reproduce actual experiments.
The numerical framework based on TIE-MS method has revealed a good potential for virtualization of X-ray SBI systems both in parallel and cone-beam implementations. While it poses no limitation on the data sampling rate within the paraxial approximation, it can efficiently simulate thick/inhomogeneous objects. The computation time for generating a projection image for polychromatic radiation in this method is short (in the order of several seconds to several minutes for a 500–1000 length of the 3D image matrix). In the case of a tomography simulation, with very thick samples or larger image volumes at the order of 1010 voxels, the computation time could be a day or more on a parallelized GPU. One way to decrease the simulation time is to make thicker slices in the object for the interpolation and propagation or reduce the number of sub-voxel interpolation points that don't add much accuracy to the calculations. However, the TIE-MS approach is computationally very efficient for reasonably sized images, sufficient to be used for virtual parametric or tomographic studies in optimizing SBI systems. This is a topic of future works.
In summary, having the two approaches in our framework enables a great flexibility in choosing and purposefully applying them for a broad range of different SBI systems and applications. In this way one can utilize the different benefits of the two approaches regrading e.g. image volume size, computational efficiency, accuracy, dosimetry, etc., according to a given simulation goal. MC approach has the advantage of providing a more realistic and accurate simulation of photon interactions within the imaging system, and therefore, is very practical in the design of diffusers and dosimetry studies of SBI techniques for biomedical applications. However, the TIE-MS is very time-efficient for parametric and tomography studies. Both approaches can be employed to evaluate and optimize SBI systems.
Funding
H2020 Marie Skłodowska-Curie Actions (754364); Eidgenössische Materialprüfungs- und Forschungsanstalt (EMPAPOSTDOCS-II programme).
Acknowledgments
We acknowledge the support from the institutions mentioned in the funding block. We appreciate the useful comments of the reviewers that contributed to improving the manuscript.
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.
Supplemental document
See Supplement 1 for supporting content.
References
1. A. Bravin, P. Coan, and P. Suortti, “X-ray phase-contrast imaging: From pre-clinical applications towards clinics,” Phys. Med. Biol. 58(1), R1–R35 (2013). [CrossRef]
2. M. Stampanoni, Z. Wang, T. Thüring, C. David, E. Rössl, U. Van Stevendaal, T. Köhler, M. Trippel, G. Singer, R. A. Kubik-Huch, M. K. Hohl, and N. Hauser, “Toward clinical differential phase contrast mammography: Preliminary evaluations and image processing schemes,” J. Instrum. 8(05), C05009 (2013). [CrossRef]
3. S. C. Irvine, K. S. Morgan, Y. Suzuki, K. Uesugi, A. Takeuchi, D. M. Paganin, and K. K. W. Siu, “Assessment of the use of a diffuser in propagation-based x-ray phase contrast imaging,” Opt. Express 18(13), 13478 (2010). [CrossRef]
4. H. Wang, Y. Kashyap, and K. Sawhney, “Hard-X-ray directional dark-field imaging using the speckle scanning technique,” Phys. Rev. Lett. 114(10), 103901 (2015). [CrossRef]
5. H. Wang and K. Sawhney, “Hard X-ray omnidirectional differential phase and dark-field imaging,” Proc. Natl. Acad. Sci. U. S. A. 118(9), e2022319118 (2021). [CrossRef]
6. H. Wang, S. Berujon, J. Herzen, R. Atwood, D. Laundy, A. Hipp, and K. Sawhney, “X-ray phase contrast tomography by tracking near field speckle,” Sci. Rep. 5, 8762 (2015). [CrossRef]
7. M. C. Zdora, “State of the Art of X-ray Speckle-Based Phase-Contrast and Dark-Field Imaging,” J. Imaging 4(5), 60 (2018). [CrossRef]
8. T. Zhou, F. Yang, R. Kaufmann, and H. Wang, “Applications of laboratory-based phase-contrast imaging using speckle tracking technique towards high energy X-rays,” J. Imaging 4(5), 69 (2018). [CrossRef]
9. M.-C. Zdora, P. Thibault, W. Kuo, V. Fernandez, H. Deyhle, J. Vila-Comamala, M. P. Olbinado, A. Rack, P. M. Lackie, O. L. Katsamenis, M. J. Lawson, V. Kurtcuoglu, C. Rau, F. Pfeiffer, and I. Zanette, “X-ray phase tomography with near-field speckles for three-dimensional virtual histology,” Optica 7(9), 1221 (2020). [CrossRef]
10. S. Berujon, H. Wang, and K. Sawhney, “X-ray multimodal imaging using a random-phase object,” Phys. Rev. A 86(6), 063813 (2012). [CrossRef]
11. D. M. Paganin, H. Labriet, E. Brun, and S. Berujon, “Single-image geometric-flow x-ray speckle tracking,” Phys. Rev. A 98(5), 053813 (2018). [CrossRef]
12. K. M. Pavlov, H. Li, D. M. Paganin, S. Berujon, H. Rougé-Labriet, and E. Brun, “Single-Shot X-Ray Speckle-Based Imaging of a Single-Material Object,” Phys. Rev. Appl. 13(5), 054023 (2020). [CrossRef]
13. S. C. Irvine, D. M. Paganin, R. A. Jamison, S. Dubsky, and A. Fouras, “Vector tomographic X-ray phase contrast velocimetry utilizing dynamic blood speckle,” Opt. Express 18(3), 2368 (2010). [CrossRef]
14. A. Olivo and R. Speller, “A coded-aperture technique allowing x-ray phase contrast imaging with conventional sources,” Appl. Phys. Lett. 91(7), 074106 (2007). [CrossRef]
15. F. Pfeiffer, M. Bech, O. Bunk, T. Donath, B. Henrich, P. Kraft, and C. David, “X-ray dark-field and phase-contrast imaging using a grating interferometer,” J. Appl. Phys. 105(10), 102006 (2009). [CrossRef]
16. M.-C. Zdora, I. Zanette, T. Walker, N. W. Phillips, R. Smith, H. Deyhle, S. Ahmed, and P. Thibault, “X-ray phase imaging with the unified modulated pattern analysis of near-field speckles at a laboratory source,” Appl. Opt. 59(8), 2270 (2020). [CrossRef]
17. S. Saghamanesh, S. M. Aghamiri, A. Kamali-Asl, and W. Yashiro, “Photon detection efficiency of laboratory-based x-ray phase contrast imaging techniques for mammography: A Monte Carlo study,” Phys. Med. Biol. 62(18), 7394–7406 (2017). [CrossRef]
18. M. C. Zdora, P. Thibault, F. Pfeiffer, and I. Zanette, “Simulations of x-ray speckle-based dark-field and phase-contrast imaging with a polychromatic beam,” J. Appl. Phys. 118(11), 113105 (2015). [CrossRef]
19. T. U. Z. Hou, M. A. H. Z. Dora, I. R. Z. Anette, J. E. R. Omell, and H. A. N. S. M. H. Ertz, “Noise analysis of speckle-based x-ray phase-contrast imaging,” Opt. Lett. 41(23), 5490–5493 (2016). [CrossRef]
20. I. A. Aloisio, D. M. Paganin, C. A. Wright, and K. S. Morgan, “Exploring experimental parameter choice for rapid speckle-tracking phase-contrast X-ray imaging with a paper analyzer,” J. Synchrotron Radiat. 22(5), 1279–1288 (2015). [CrossRef]
21. S. Meyer, S. Z. Shi, N. Shapira, A. D. A. Maidment, and P. B. Noël, “Quantitative analysis of speckle-based X-ray dark-field imaging using numerical wave-optics simulations,” Sci. Rep. 11(1), 1–9 (2021). [CrossRef]
22. T. Faragó, P. Mikulík, A. Ershov, M. Vogelgesang, D. Hänschke, and T. Baumbach, “Syris: A flexible and efficient framework for X-ray imaging experiments simulation,” J. Synchrotron Radiat. 24(6), 1283–1295 (2017). [CrossRef]
23. R. Zboray, “Optimizing and applying high-resolution, in-line laboratory phase-contrast X-ray imaging for low-density material samples,” J. Microsc. 282(2), 123–135 (2021). [CrossRef]
24. H. Rouge-Labriet, L. Quenot, S. Bohic, B. Fayard, D. M. Paganin, E. Brun, and S. Berujon, “Comparison of X-ray speckle-based imaging deflection retrieval algorithms for the optimization of radiation dose,” Phys. Med. Biol. 66(6), 065005 (2021). [CrossRef]
25. E. Abadi, W. P. Segars, B. M. W. Tsui, P. E. Kinahan, N. Bottenus, A. F. Frangi, A. Maidment, J. Lo, and E. Samei, “Virtual clinical trials in medical imaging: a review,” J. Med. Imaging 7(4), 042805 (2020). [CrossRef]
26. J. Fu, X. Hu, A. Velroyen, M. Bech, M. Jiang, and F. Pfeiffer, “3D algebraic iterative reconstruction for cone-beam X-Ray differential phase-contrast computed tomography,” PLoS One 10(3), 1–13 (2015). [CrossRef]
27. O. Chubar and P. Elleaume, “Accurate and efficient computation of Synchrotron Radiation in the near field region,” EPAC Proc., 1177–1179 (1998).
28. M. Sanchez Del Rio, N. Canestrari, F. Jiang, and F. Cerrina, “SHADOW3: A new version of the synchrotron X-ray optics modelling package,” J. Synchrotron Radiat. 18(5), 708–716 (2011). [CrossRef]
29. J. Bahrdt, U. Flechsig, S. Gerhardt, and I. Schneider, “PHASE: a universal software package for the propagation of time-dependent coherent light pulses along grazing incidence optics,” Proc. SPIE 8141, 81410E (2011). [CrossRef]
30. E. Bergbäck Knudsen, A. Prodi, J. Baltser, M. Thomsen, P. Kjær Willendrup, M. Sanchez Del Rio, C. Ferrero, E. Farhi, K. Haldrup, A. Vickery, R. Feidenhans, K. Mortensen, M. Meedom Nielsen, H. Friis Poulsen, S. Schmidt, and K. Lefmann, “McXtrace: A Monte Carlo software package for simulating X-ray optics, beamlines and experiments,” J. Appl. Crystallogr. 46(3), 679–696 (2013). [CrossRef]
31. M. S. Del Rio and L. Rebuffi, “OASYS: A software for beamline simulations and synchrotron virtual experiments,” AIP Conf. Proc. 2054(January), (2019).
32. L. Quénot, E. Brun, J. M. Létang, and M. Langer, “Evaluation of simulators for x-ray speckle-based phase contrast imaging,” Phys. Med. Biol. 66(17), 175027 (2021). [CrossRef]
33. H. Itoh, K. Nagai, G. Sato, K. Yamaguchi, T. Nakamura, T. Kondoh, C. Ouchi, T. Teshima, Y. Setomoto, and T. Den, “Two-dimensional grating-based X-ray phase-contrast imaging using Fourier transform phase retrieval,” Opt. Express 19(4), 3339 (2011). [CrossRef]
34. E. R. Shanblatt, Y. Sung, R. Gupta, B. J. Nelson, S. Leng, W. S. Graves, and C. H. McCollough, “Forward model for propagation-based x-ray phase contrast imaging in parallel- and cone-beam geometry,” Opt. Express 27(4), 4504 (2019). [CrossRef]
35. D. G. Voelz and M. C. Roggemann, “Digital simulation of scalar optical diffraction: Revisiting chirp function sampling criteria and consequences,” Appl. Opt. 48(32), 6132–6142 (2009). [CrossRef]
36. C. Zuo, J. Li, J. Sun, Y. Fan, J. Zhang, L. Lu, R. Zhang, B. Wang, L. Huang, and Q. Chen, “Transport of intensity equation: a tutorial,” Opt. Lasers Eng. 135, 106187 (2020). [CrossRef]
37. A. V. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” J. Opt. Soc. Am. A 19(3), 472 (2002). [CrossRef]
38. S. Agostinelli and E. Al, ““GEANT4 - A simulation toolkit,” Nucl. Instruments Methods Phys. Res. Sect. A Accel. Spectrometers,” Nucl. Instrum. Methods Phys. Res., Sect. A 506(3), 250–303 (2003). [CrossRef]
39. S. Jan, G. Santin, D. Strul, S. Staelens, K. Assié, D. Autret, S. Avner, R. Barbier, M. Bardiès, P. M. Bloomfield, D. Brasse, V. Breton, P. Bruyndonckx, I. Buvat, A. F. Chatziioannou, Y. Choi, Y. H. Chung, C. Comtat, L. Simon, T. Y. Song, J.-M. Vieira, D. Visvikis, R. Van De Walle, E. Wieërs, and C. Morel, “GATE -Geant4 Application for Tomographic Emission: a simulation toolkit for PET and SPECT,” Phys Med Biol,” Phys. Med. Biol. 49(19), 4543–4561 (2004). [CrossRef]
40. A. Deresch, “Modellierung von Röntgenspektren für technische Anwendungen,” Universität Potsdam (2015).
41. M. Langer, Z. Cen, S. Rit, and J. M. Létang, “Towards Monte Carlo simulation of X-ray phase contrast using GATE,” Opt. Express 28(10), 14522 (2020). [CrossRef]
42. S. Cipiccia, F. A. Vittoria, M. Weikum, A. Olivo, and D. A. Jaroszynski, “Inclusion of coherence in Monte Carlo models for simulation of x-ray phase contrast imaging,” Opt. Express 22(19), 23480 (2014). [CrossRef]
43. I. Zanette, T. Zhou, A. Burvall, U. Lundström, D. H. Larsson, M. Zdora, P. Thibault, F. Pfeiffer, and H. M. Hertz, “Speckle-based x-ray phase-contrast and dark-field imaging with a laboratory source,” Phys. Rev. Lett. 112(25), 253903 (2014). [CrossRef]
44. M. R. Teague, “Deterministic Phase Retrieval: a Green’S Function Solution,” J. Opt. Soc. Am. 73(11), 1434–1441 (1983). [CrossRef]
45. K. A. Nugent, T. E. Gureyev, D. F. Cookson, D. Paganin, and Z. Barnea, “Quantitative phase imaging using hard X rays,” Phys. Rev. Lett. 77(14), 2961–2964 (1996). [CrossRef]
46. D. M. Paganin and D. Pelliccia, “Tutorials on X-ray Phase Contrast Imaging: Some Fundamentals and Some Conjectures on Future Developments,” arXiv (2019).
47. K. Li, M. Wojcik, and C. Jacobsen, “Multislice does it all—calculating the performance of nanofocusing X-ray optics,” Opt. Express 25(3), 1831 (2017). [CrossRef]
48. D. Paganin, S. C. Mayo, T. E. Gureyev, P. R. Miller, and S. W. Wilkins, “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” J. Microsc. 206(1), 33–40 (2002). [CrossRef]
49. T. E. Gureyev, A. Pogany, D. M. Paganin, and S. W. Wilkins, “Linear algorithms for phase retrieval in the Fresnel region,” Opt. Commun. 231(1-6), 53–70 (2004). [CrossRef]
50. T. Schoonjans, A. Brunetti, B. Golosio, M. Sanchez Del Rio, V. A. Solé, C. Ferrero, and L. Vincze, “The xraylib library for X-ray-matter interactions. Recent developments,” Spectrochim. Acta, Part B 66(11-12), 776–784 (2011). [CrossRef]
