Imaging-system modeling

Before starting the MC-based physical-interaction simulations, the FPXS, imaging geometry, and detector in the FPXS imaging system were modeled, and a unified coordinate system was established. In addition, the detailed geometric parameters and relative positions of each component were defined and determined.

To obtain the PHSP and analyze the beam characteristics of the X-rays generated by the FPXS, the BEAMnrc package (version 1.25) [35] was used to perform a structure modeling of the source. Since the cathode and anode of the FPXS include millions of patterns, directly simulating the structure of the entire FPXS to generate a PHSP is difficult and time-consuming. Under the premise that the fabrication process is stable and that all Al-Mo-Al patterns have exactly the same structure and output the same PHSP data, we propose an approach to represent the PHSP of the entire source by duplicating the pattern PHSP and sampling the pattern position. Under a strong electric field between the cathode and anode plates, the direction of motion of the electrons produced by the cathode is almost perpendicular to the anode target; therefore, a parallel electron beam is used to represent the electron beam emitted by the ZnO nanowires. Considering that the FPXS anode pattern had a diameter of 30 μm and spacing of 80 μm, numerous electrons were emitted from the ZnO nanowires bombarding the ITO in the gap of the Al-Mo-Al patterns, and the resulting X-rays significantly affected the FPXS spectrum. A square electron beam with a side length of 80 μm was used to completely cover each Al-Mo-Al pattern and the surrounding ITO. The kinetic energy of the simulated electron beam was 35 kV. Figure 3(a) shows the structure of the FPXS modeling. The FLATFILT module in BEAMnrc was used to simulate the Al-Mo-Al pattern of the anode with a diameter of 30 μm, and the outer region of the module was vacuum. Three square SLABS modules with side lengths of 4.80 cm were used to simulate the ITO film, glass, and virtual detector. A total of 20 billion electrons were simulated, and the 280 million photons generated were recorded by the virtual detector and saved as PHSP data. The states of each photon included energy, motion direction, and relative distance to the center of the pattern.

Fig. 3Full-screen View

(Color online) (a) Schematic of FPXS modeling; (b) schematic of FPXS imaging system

Imaging-geometry modeling and detector modeling were also performed, and a unified coordinate system for the FPXS imaging system was established. To model the imaging geometry, parameters such as the matrix size, resolution, and offsets of the geometric phantom were defined in the imaging system. The required parameters for the detector including pixel size, resolution, and offsets were defined, and the detector response function was defined to represent the probability of a photon being attenuated by the realistic detector material and the efficiency of converting the optical signal into an electrical signal. The input parameters of the detector response were the photon energy and position of the photon arriving at the detector. Figure 3(b) shows the unified coordinate frame of the entire imaging system. The geometric coordinate frame (X_g, Y_g, Z_g) was established by setting the center of the geometry as the origin G. For the FPXS, a coordinate frame (X_s, Y_s, Z_s) was defined, with its origin S at the center of the FPXS. For the detector, a coordinate frame (X_d, Y_d, Z_d) was established with its origin D at its center. Notably, the centers of the FPXS, geometry, and detector were all on the same line; the X_s axis was parallel to the X_g and X_d axes, and the Y_s axis was parallel to the Y_g and Y_d axes. Through the transformation of different coordinate frames, the exact coordinates at which the photons were located in the process of physical interaction in the geometry can be determined. Furthermore, the relative positions of the FPXS, imaging geometry, and detector within the imaging system were determined by defining the source-to-axis distance (SAD) and source-to-imager distance (SID) values over three coordinate frames.

Photon initialization in gPSMC

Photon initialization determines the energy, position, and direction of motion of the initial photons in the FPXS imaging system. Photon initialization in gPSMC was performed by obtaining the initial photon state in the PHSP, including pattern sampling and PHSP sampling techniques in the FPXS coordinate frame (X_s, Y_s, Z_s), as indicated in Fig. 4. A pattern sampling technique was proposed to generate the pattern position to which the emitted photon belongs in the entire patterned array and calculate its center coordinates. Specifically, random sampling was used to determine the row and column indices of the sampled pattern in the patterned array, and the row and column indices were used to calculate the distance between the center of the sampled pattern ($x_0\mathrm{,}{y}_0\mathrm{,}{z}_0$) and that of the FPXS in the (X_s, Y_s, Z_s) coordinate system, as shown in Fig. 4(a). The PHSP sampling technique was proposed to randomly sample photons from the simulated PHSP data on the lower surface of the pattern to further determine the state (${\alpha }_0\mathrm{,}{\beta }_0\mathrm{,}{\gamma }_0\mathrm{,}{v}_0\mathrm{,}{e}_0$) of the sampled initial photons, as depicted in Fig. 4(b), where (${\alpha }_0\mathrm{,}{\beta }_0\mathrm{,}{\gamma }_0$) denote the coordinates of the sampled photons in the PHSP, and v₀ and e₀ denote the direction and energy of the sampled photons, respectively. The energy and motion directions obtained by the PHSP sampling were utilized as the energy and direction of the current photon, whereas the final photon position (x, y, z) was calculated from the sum of the sampled photon coordinates (${\alpha }_0\mathrm{,}{\beta }_0\mathrm{,}{\gamma }_0$) and sampled pattern coordinates ($x_0\mathrm{,}{y}_0\mathrm{,}{z}_0$). Thus, the latest state ($x\mathrm{,}y\mathrm{,}z\mathrm{,}{v}_0\mathrm{,}{e}_0$) of the initial photon was updated.

Fig. 4Full-screen View

(Color online) Schematic of the photon-initialization process in gPSMC. (a) Pattern sampling and center-coordinates calculation; (b) PHSP sampling; (c) PHSP optimization used in gPSMC-o

Owing to the large range of the initial directions of the photons in the PHSP, preliminary experiments revealed that most of the sampled initial photons (up to 90%) obtained through the above steps could not reach the detector when the motion direction was constant without scattering; thus, the simulation of these photons was impractical. To improve the utilization of the initial photons and accelerate the subsequent simulation process, a PHSP optimization technique was proposed to screen the sampled initial photons, as depicted in Fig. 4(c). A new simulator gPSMC-o was proposed by incorporating the PHSP optimization technique into gPSMC. Considering that the direction of scattering changes when a photon interacts with the geometry, the screening process is based on whether the photon intersects the detector or the surface of the geometry when moving in the initial direction. In the unified coordinate system, the initial direction and coordinates of the initial photon were used to determine whether the photon could enter the geometry or reach the detector surface after an outward expansion of 2 cm. Only photons that can reach the geometry surface or detector region were used for subsequent physical-interaction simulations, and unreachable photons were discarded.

Photon initialization in gFSMC

For gFSMC, photon initialization was implemented based on the energy spectrum and fluence map, rather than on the PHSP. Because of the identical structure of all the patterns in the FPXS, the energy spectrum of a single pattern can represent the energy spectrum of the entire FPXS. The energy spectrum curve was calculated by counting the number of photons at different energy levels for all the photons in the PHSP generated as described in Sect. 2.2.1, shown in Fig. 5(a). In this study, the fluence map refers to the energy distribution of the photons reaching the detector without imaging geometry, and the values of different pixels in the distribution represent the probability density of the photons arriving at different locations on the detector. The fluence map can be obtained by acquiring or simulating an FPXS air-scan image, where the air-scan image corresponds to the energy-deposition distribution of all photons emitted by the FPXS reaching the detector in the absence of geometry. The projection of the air phantom simulated by gPSMC was used as the fluence map to ensure the consistency of the simulation conditions between the two MC simulators.

Fig. 5Full-screen View

(Color online) Schematic of the photon-initialization process in gFSMC. (a) Energy-spectrum sampling; (b) fluence-map sampling and initial-direction calculation

The proposed photon-initialization approach in gFSMC includes energy-spectrum, pattern, and fluence-map sampling techniques. Figure 5 presents a schematic of the photon-initialization process. The energy e₀ of the initial photons was determined by performing Metropolis sampling [36] in 235 energy bins of the energy spectrum curve, e.g., the sampled energy corresponding to the red line in Fig. 5(a). To determine the direction of motion of the initial photon, a pattern-sampling technique was used to generate pattern positions from the entire patterned array. Notably, each pattern in the FPXS arrays was approximated as an infinitesimal spot source, and the center coordinates ($x_0\mathrm{,}{y}_0\mathrm{,}{z}_0$) of the sampled pattern were used to represent the coordinates (x, y, z) of the initial photon, as depicted in Fig. 5(b). A pixel in the detector was randomly selected by the Metropolis sampling of the fluence map, which represented the location (x_fmp, y_fmp), where the photon would hit the detector if it moved in the initial direction. The initial direction v₀ of the photon was determined by connecting this detector pixel to the center of the sampled pattern in the unified coordinates of the imaging system, as depicted in Fig. 5(b). Thus far, all parameters of the initial photon have been calculated, and photon initialization in gFSMC has been completed.

Physical-interaction simulation in phantom

Physical files including differential cross-section (DSC) data, mass attenuation coefficient, and density-mapping data of various simulated materials need to be initialized to ascertain how the photons interact with multiple materials and determine the photon-path length within the phantom. For the initial photons that have completed initialization on the central processing unit (CPU), we propose a GPU-based MC approach for physical-interaction simulation in the phantom; the gPSMC and gFSMC simulators have the same simulation process.

The parallel-computing advantage of the GPU was exploited to accelerate the physical-interaction simulation, as shown in Table 1. The total histories were divided into multiple batches based on the energy of the initial photons sampled from the photon initialization, with each batch having similar energies. Batch division is used only for categorizing photons of different energies to speed up the simulation process rather than for sampling the energies of the initial photons. Ten batches were used in this study. For photons within the same batch, multiple GPU kernels were launched to simulate the transport of all the photons simultaneously, with one kernel tracking the history of multiple photons sequentially. The actual energy deposited on each pixel of the detector was calculated using the detector response function. By accumulating the energy deposition of the same batch of photons arriving at the detector, the energy deposition of different batches and total energy deposition were obtained and transferred to the CPU to calculate the final projection.

A detailed physical-interaction simulation workflow for a single photon is shown in Fig. 6. The Woodcock tracking [37] method was utilized in this study to avoid calculating the intersection of the photon paths with each voxel boundary to speed up the simulation. During the simulation, the forward distance of the initial photon was sampled and a new interaction point was calculated. Next, random sampling was used for the Woodcock method to determine the interaction type at the current position, including the photoelectric effect, Compton scattering, Rayleigh scattering, and virtual interactions. If the photoelectric effect occurred, the simulation of this photon was terminated and the generated electrons were deposited in situ. When Compton or Rayleigh scattering occurred, the scattering direction and energy of the scattered photons were determined by sampling the DSC data. The photons remained unchanged in the virtual interaction, which was utilized to speed up the simulation process. After photon simulation at the current position was completed, the forward distance was sampled again, and the above steps were repeated until the photon was absorbed or moved out of the phantom. We also determined if the photons moving outside the phantom could reach the detector. The position at which a photon reached the detector was calculated and transformed to the detector coordinate system (X_d, Y_d, Z_d) to determine the detector pixel that was hit, and the deposited energy of the photon was recorded. Another stopping condition was set to terminate the simulation of the photons with energies below 1 keV. Through the simulation of physical interactions in the phantom, the final projection consisted of two parts: the primary signal, which underwent attenuation during the propagation of X-rays, and the scatter signal, which underwent a change in the direction of the photons.

Fig. 6Full-screen View

(Color online) Workflow of physical-interaction simulation in phantom

Water phantom experiment

A simulated water phantom was created, and the projections (primary signal) of the phantom under the same simulation conditions were calculated using gPSMC, gPSMC-o, gFSMC, and the ray-tracing approach. Fig. 7(a) shows a schematic of the water-phantom imaging system. The ray-tracing approach superimposed the projections calculated using different patterns to obtain the final projection, which was suitable for comparison with the MC simulator proposed in this study. The water phantom was a three-dimensional matrix with a Housfield unit (HU) value of 0, the size of the phantom was 6.00 cm × 6.00 cm × 10.00 cm, and the dimensions were 120 × 120 × 200. In the imaging system, two small patterned arrays of 4 × 4 and 6 × 6 were utilized for projection calculation and comparison, with a distribution area of 4.80 cm × 4.80 cm and spacings of 1.60 cm and 0.96 cm, respectively. The size of the detector was 12.00 cm × 12.00 cm and the resolution was 1548 × 1548. The digital water phantom was placed close to the detector, the SAD was set to 35 cm, and the SID was set to 40 cm. The inputs of gFSMC and the ray-tracing approach were the same energy spectrum and fluence map of FPXS, and the fluence was calculated based on gPSMC using an air phantom with a size of 12.00 cm × 12.00 cm × 40.00 cm and dimension of 60×60×200. To evaluate the projections calculated using the three methods quantitatively, the mean absolute percentage error (MAPE), contrast-to-noise ratio (CNR), and full width at half maximum (FWHM) were used to calculate the differences between the projections, which are defined as follows: $\text{MAPE}=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|\frac{I_{\text{rt}}^i-I_{\text{MC}}^i}{I_{\text{rt}}^i}\right|}\mathrm{,}$ (1) $\text{CNR}=\frac{|{\mu }_{\text{r}}-{\mu }_{\text{b}}|}{\sigma }\mathrm{,}$ (2) $\text{FWHM}=V_{50\%}^{\text{max}}-V_{50\%}^{\text{min}}\mathrm{,}$ (3) where $I_{\text{MC}}^i$ and $I_{\text{rt}}^i$ are the values at voxel i of the projections calculated by MC simulators and the ray-tracing approach, and n is the number of voxels in the projection. μ_r and $\sigma$ are the mean value and standard deviation of the selected region of interest (ROI) in the border region of the projection, respectively, and μb is the mean value of the selected central region of the projection. $V_{50\%}^{\text{max}}$ and $V_{50\%}^{\text{min}}$ are the maximum and minimum voxel coordinates at 50% of the height of the wave crest, respectively.

Fig. 7Full-screen View

(Color online) Experiments conducted on the water phantom and PMMA phantom. (a) Schematic of the water-phantom imaging system; (b) schematic of the PMMA-phantom imaging system; (c) photograph of the projection-acquisition devices for PMMA phantom

PMMA phantom experiment

To further verify the accuracy of the two MC simulators in the projection (primary signal + scatter signal) simulation of the entire FPXS (with patterned arrays of 600×600 and a center spacing of 80 μm), a PMMA phantom experiment was designed, and the projections calculated by gPSMC, gPSMC-o, and gFSMC for the two PMMA phantoms were compared with the realistic acquired projections. Figure 7(b) shows a schematic of the PMMA-phantom imaging system. The two simulated PMMA phantoms had the same volume and density as the realistic PMMA slabs, where the 4-cm-thick PMMA phantom comprised four PMMA slabs with a size of 10.00 cm × 15.00 cm× 1.00 cm, and the 6-cm-thick phantom comprised one PMMA slab with a size of 10.00 cm× 10.00 cm× 2.00 cm and four PMMA slabs with a size of 10.00 cm × 15.00 cm × 1.00 cm. The sizes of the two simulated PMMA phantoms were 10.00 cm × 15.00 cm × 4.00 cm and 10.00 cm × 15.00 cm × 6.00 cm with dimensions of 200 × 300 × 80 and 200 × 300× 120, respectively. A photograph of the PMMA projection-acquisition devices is shown in Fig. 7(c). Given the size of the PMMA slab, a flat-panel CMOS X-ray detector (Xineos-2329, Teledyne DALSA, USA) with an active area of 22.80 cm × 29.20 cm and resolution of 4608 × 5888 was used for projection acquisition. The detector for the PMMA projection simulation had the same dimensions and parameters as those of a realistic detector. During the PMMA projection simulation and acquisition, the two PMMA phantoms were placed close to the detector with the SID set to 40 cm and SAD set to 38 cm and 37 cm.

Line-pair phantom experiment

A simulated line-pair phantom was generated with reference to a real line-pair testing card (Type18D, IBA, Germany), and the MC-calculated projections (primary signal + scatter signal) were compared with the actual acquired projections. A schematic of the line-pair phantom imaging system is shown in Fig. 8(a). The line-pair testing card used in this study had a resolution range of 0.50 lp/mm to 20.00 lp/mm, and the simulated line-pair phantom comprised 21 types of line-pairs with different resolutions (0.50 lp/mm to 5.00 lp/mm), with a size of 4.00 cm× 4.00 cm× 0.68 cm and dimension of 4000× 4000× 8. Figure 8(b) shows the simulated line-pair phantom. A line-pair phantom was utilized to investigate the variation in imaging resolution with the distance between the FPXS and detector, where the axis-to-imager distance (AID) remained unchanged (0.34 cm) and the SID varied from 16 to 90 cm, with a total of 13 groups. In the process of collecting realistic projections, a line-pair testing card and flat-panel CMOS X-ray detector (Xineos-1515, Teledyne DALSA, USA) were used for the projection acquisition of the FPXS. Figure 8(c) shows a photograph of the projection-acquisition devices, including the detector, lifting platform, and metal frame. The detector had an active area of 15.30 cm × 15.30 cm, a pixel size of 99.00 μm × 99.00 μm, and a resolution of 1548× 1548. The line-pair card was fixed at the center of the detector surface, and realistic projections at different SIDs from 16 to 90 cm were acquired and compared with the projections calculated using the three MC simulators.

Fig. 8Full-screen View

(Color online) Experiment on the line-pair phantom. (a) Schematic of the line-pair phantom imaging system; (b) simulated line-pair phantom with resolutions of 0.50 lp/mm to 5.00 lp/mm; (c) photograph of the projection-acquisition devices

To quantitatively evaluate the difference between the acquired and MC-calculated projections as well as the difference in the exact line-pair values, the contrast transfer-function (CTF) curve and maximum line-pair values corresponding to the 10% CTF curve were calculated, where the CTF curve reflects the square-wave response of the imaging system at different frequencies [38]. To calculate the CTF, the modulation value M₀ at zero frequency was first calculated using high attenuation and background regions of the same size. Then, M₀ was used to normalize M(ω) to obtain the discrete CTF(ω). ω is the spatial frequency (lp/mm) in one direction. The value corresponding to 10% CTF is the desired maximum line-pair value. CTF(ω) was calculated as follows: $CTF(\omega )=\frac{M(\omega )}{M_0}=\frac{\frac{I_{\text{max}}-I_{\text{min}}}{I_{\text{max}}+I_{\text{min}}}}{\frac{Ave{I}_{\text{N}}-Ave{I}_{\text{H}}}{Ave{I}_{\text{N}}+Ave{I}_{\text{H}}}}\mathrm{,}$ (4) where AveIN and AveIH are the average values of the areas with no attenuation and high attenuation in the projection, respectively, and I_max and I_min are the maximum and minimum values of the line-pair at the corresponding frequency, respectively.

Qualitative comparison

The gPSMC, gPSMC-o, and gFSMC simulators were used to calculate the projections of the simulated line-pair phantom with simulated histories of 1×10¹², and the projections were compared with the actual acquired projections of the line-pair testing card. In the line-pair phantom experiment, the simulated line-pair phantom and actual testing card were fixed close to the detector, and the projections were calculated and acquired under an SID from 16 to 90 cm. Figure 11(a)-(d) show the acquired and calculated projections using gPSMC, gPSMC-o, and gFSMC at SID = 35 cm. The maximum line pairs that can be recognized in the projections are marked in red boxes, and the line-pair values for all four projections reached up to 1.25 lp/mm. The profiles of the four projections represented by the blue line are shown in Fig. 11(e)-(h), from which five distinct peaks (red boxes) at the maximum line pair can be clearly recognized. Figure 11(i)-(l) present the projection comparison at SID = 50 cm, with maximum line-pair values of 1.60, 1.61, 1.61, and 1.61 lp/mm, respectively. The profiles in Fig. 11(m)-(p) exhibit five distinct peaks. The maximum line-pair values of the calculated and acquired projections for the different SIDs were the same. The projection and profile comparisons also indicated that the gPSMC-o- and gFSMC-calculated projections were less noisy than that of gPSMC, and the edges of the line pairs in the projection were sharper.

Fig. 11Full-screen View

(Color online) Comparison of projections and profiles at different SIDs. (a)-(d) Acquired projections and projections calculated using gPSMC, gPSMC-o, and gFSMC simulators at SID=35 cm; (e)-(h) profiles of (a)-(d) on the blue line; (i)-(l) acquired projections and projections calculated using gPSMC, gPSMC-o, and gFSMC simulators at SID=50 cm; (m)-(p) profiles of (i)-(l) on the blue line

Quantitative evaluation

For each calculated and acquired projection, the maximum line-pair values were quantitatively evaluated. Specifically, the CTF values were calculated for all line pairs at the corresponding frequencies in each projection, and the CTF curves were plotted. Figure 12(a) shows the CTF curves for the acquired and simulated projections at an SID of 50 cm. The line-pair values of the acquired and calculated projections start from 0.50 lp/mm, and the CTF values decrease gradually with an increase in line-pair values. For the acquired projections, M₀ was calculated based on the specific no-lead and background regions in the projections, whereas for the simulated projection, the CTF of the 0.50 lp/mm line pair was taken as M₀. For line pairs with low frequencies, the CTF values of the gPSMC-, gPSMC-o-, and gFSMC-calculated projections were close to 100% and decreased gradually as the line-pair values increased. Moreover, the CTF curves of gPSMC, gPSMC-o, and gFSMC were consistent. The inset (green box) in Fig. 12(a) presents the 10% CTF values for the acquired and calculated projections at an SID of 50 cm, which were 1.61, 1.65, 1.63, and 1.62 lp/mm. Figure 12(b) summarizes the maximum line-pair values (10% CTF values) of the acquired and calculated projections for different SIDs. The line-pair values of the acquired and calculated projections were consistent, with errors of 4.13%, 3.86%, and 4.15%. Moreover, all maximum line-pair values increased with increasing SID, which indicated that the imaging resolution increased with increasing imaging distance while keeping the AID unchanged.

Fig. 12Full-screen View

(Color online) Comparison of CTF curves and maximum line-pair values of acquired and calculated projections using gPSMC, gPSMC-o, and gFSMC simulators at SID=50 cm. (a) Comparison of CTF curves of the CTF values varying with line-pair values; (b) comparison of maximum line-pair values varying with SIDs