Semi-analytic off-axis X-ray source model

Spectral computed tomography (CT) systems are employed with energy-resolving photon counting detectors. Incorporation of a spectrally accurate x-ray beam model in image reconstruction helps to improve material identification and quantification by these systems. Using an inaccurate x-ray model in spectral reconstruction can lead to severe image artifacts, one of the extreme cases of this is the well-known beam-hardening artifacts. An often overlooked spectral feature of x-ray beams in spectral reconstruction models is the angular dependence of the spectrum with reference to the central beam axis. To address these factors, we have developed a parameterized semi-analytical x-ray source model in the diagnostic imaging range (30-120 kVp) by applying regression techniques to data obtained from Monte Carlo simulations (EGSnrc). This x-ray beam model is generalized to describe the off-axis spectral information within ±17o along θ (vertical direction), ±5o along ϕ (horizontal direction) of the central axis, and can be parameterized for specific x-ray tube models. Comparisons of our model with those generated by SpekCalc, TOPAS, and IPEM78 at central axis show good agreement (within 2 %). We have evaluated the model with experimental data collected with a small animal spectral scanner.


Introduction
Spectral CT is an emerging x-ray imaging technique that has the potential to revolutionize x-ray imaging techniques as the energy of photons in an x-ray beam can be resolved [1,2]. This energy information can then be used for identifying and quantifying materials based on their spectral signatures [3]. This spectral information is typically analyzed through spectral reconstruction techniques which uses some model of the energy distribution of the x-ray beam. For such techniques to be successful, they require this model to be an adequate representation of the true x-ray beam. If an incorrect energy distribution of the beam is assumed, then the optimization techniques used in the spectral reconstruction may produce biased solutions or, in the worst case, may not converge [4].
-1 -others used general purpose Monte Carlo code such as EGS4, MCNP, ITS (Integrated TIGER Series) [25,26,28,29]. Ay et al. presented a Monte Carlo model for diagnostic radiology and mammography using MCNP4C. It provides information along the cathode-anode direction but not along the vertical axis (perpendicular to cathode-anode direction). They used anode angle ranging from 6 • to 18 • . Bhat et al. used EGS4 code to estimate the off-axis spectra along the 6 • cathode-anode direction [26]. Bontempi et al. used the Monte Carlo code to estimate the spectra for the tungsten anode for different anode angles ( 10 • , 12 • , 15 • , 17 • , and 20 • ) for 50-140 kVp [30]. Although these simulations give off-axis spectral information, the spectral information is not commercially available and is not parameterized for efficient use. This limitation can be covered either starting with an existing established model and then to extend it to include off-angle energy dependence or developing a new x-ray model. The future use of this off-axis spectral information is to incorporate it into the spectral image reconstruction which favours the development of a new parameterised x-ray model.
In this paper we present a parameterized semi-analytical x-ray source model that we have developed to describe the output from a Source Ray SB-120-350 x-ray tube. These x-ray tubes are installed in several of our small animal Medipix All Resolution Systems (MARS). We developed this model by applying regression methods to data from Monte Carlo simulations that we have run for this particular x-ray tube as configured within the MARS scanner. We compare the spectrum of our model with other x-ray models such as SpekCalc, IPEM78, and TOPAS. We also present beam profile comparisons between this model and data collected from the MARS scanner using a CdTe Medipix3RX energy discriminating photon counting detector. This x-ray model is referred as MARS x-ray model and can be used for all those x-ray tubes having 20 • anode angle and 50 µm effective focal spot. The spectral information from the MARS x-ray model can be made available on request.

Spectral system
The MARS small animal scanner shown in figure 1 is a spectral system employed with hybrid Medipix3RX detector. High-Z sensor layers such as CZT and CdTe are being used in the MARS spectral system. The Medipix3RX detector has a pixel pitch of 110 µm in charge summing mode (CSM) [31]. The detector inside the MARS scanner can be moved up/down perpendicular to the rotation axis (θ) to enlarge the field of view. In the current MARS setup, maximum vertical (along y axis) and horizontal (along z axis which is anode-cathode direction) coverage of the detector at a source to detector distance (SDD) of 11 cm ranges from ±17 o along vertical direction (θ), ±5 o along horizontal direction (φ) respectively.

Development of MARS x-ray model
The development of the MARS x-ray model comprises of three steps include (A) running the Monte Carlo simulations (B) extracting the data from phase space files, and (C) applying the regression on the Monte Carlo extracted data and mapping dependency of coefficients of regression using an independent variable order polynomial to get the final MARS x-ray model.

Monte Carlo simulation
In this work, we chose to use EGSnrc (a Monte Carlo simulation code) as the modular design of this code allows us to use any combination of various types of the sources with different thickness of material as a filter or collimators [32]. Furthermore, using the advanced features like variance reduction capabilities provides a good accuracy with processing the less value of history number. A thickness of 1 mm tungsten was used for anode surface which was mounted on a copper holder plunged in a vacuum container. The anode angle was set to 20 • . Since the exit window of actual x-ray tube is made of glass material which is equivalent to 1.8 mmAl filter, I chose this value for the intrinsic filter as this is the best reliable value provided by the x-ray tube's manufacturer. There are manufacturing variations (non-uniform thickness) in the glass window which are not accounted for in these simulations. This may affect the simulations results in the low energy range. To define electron beam, we selected monochromatic "parallel circular beam inside from side" with the kinetic energy of 0.12 MeV and the radius of 0.0073 cm. We used 0.521icru.pegs4dat (cross section library). AE (inelastic collision threshold) and AP (photon production threshold) were set to 0.521 and 0.01 MeV, respectively. The parameter of ECUT (electron transport cutoff) and PCUT(photon transport cutoff) was set to 0.521 and 0.01 MeV, respectively. We used NIST's XCOM photon cross sections for scattering and photoelectric absorption. Bound Compton scattering, atomic relaxation, electron impact ionization, and rayleigh scattering were switched on. We used directional beam splitting (DBS) as a variance reduction technique. The beam splitting number was chosen to 1000 and radius of a splitting field were set to 5 cm. A scoring plane of 7 × 7 cm 2 was defined at the source to surface distance of 11 cm. The number of particle histories was 2 × 10 9 . We ran ten simulations to cover the tube voltages ranging from 30 to 120 kVp with the step size of 10 kVp. An example of a simulated spectrum by BEAMnrc with statistical uncertainties is shown in figure 2. The maximum statistical uncertainty of the photon fluence at each energy bin (1 keV), is less than 1% in the head and tail of the spectrum, and in the mid energies, it is even less than 0.3%.

Monte Carlo data extraction and unit transformation
An in-house custom application was used to extract a subset of the information in the form of text file for further processing (regression) in MATLAB. This custom application extracts the number of photons (for a particular energy band) deposited in each pixel. The energy resolution and pixel size can be changed as desired. For the work presented here, photons were divided by energy into 1 keV bins, and the scoring plane was divided into pixels with a size of 1mm 2 . The pixel size of 1mm 2 was used for better data resolution and less statistical noise for the said number of histories. We transformed the "counts/(pixel×incident particle×keV)" generated from the Monte Carlo simulations to "counts/(µsr×µA× s)". Where "sr" is for steradian (SI unit of solid angle), "A" is for Ampere (tube current), and "s" is for second (exposure time).

Source model regression
A surface fitting approach was applied to the Monte Carlo data using linear-quadratic function regression in MATLAB. This regression provides an equation describing the photon distribution along the angular dimensions. Figure 3 demonstrates a fitted surface to simulated data for the integrated energy spectrum of 120 kVp in four views. As it can be observed in this figure, the fitted surface formed a symmetrical shape along the θ direction (figure 3b), while in the φ direction, high counts concentrate more in the right side of the curve (cathode side), consequently fitted curve tend to that side (figure 3c). Residual values were obtained in both directions by subtracting fitted values from the simulated values at the same points and are shown in figures 4a, and 4b. The equation (3.1) extracted from this fitted curve for broad bin includes six terms.

JINST 12 P10013
Where S θφ E is the spherical count at given angle θ and φ and S 00 E is the spherical counts at the central beam axis. At the extreme values for θ and φ, the last two terms of this equation varied ±6 × 10 −5 and ±2.72 × 10 −4 , respectively. Since these two terms do not contribute significantly (less than 0.01%), they are negligible, and the equation (3.1) can be rewritten as: A quadratic surface fitting was applied to the Monte Carlo data to derive an equation for the photon distribution in the energy range of 10 to 120 keV. R-square value for this surface fitting was 0.986. Residual values were under 0.5% of total counts in both θ and φ direction. In this model, characteristic component was fitted separately from the bremsstrahlung. The general form of the source function for bremsstrahlung part with 1 keV energy bins for range of tube voltage is shown in the equation (3.3).
S 00 EV is the spectrum at the central axis of the beam and it varies with the energy and tube voltage. A EV , B EV , and C EV are the regression coefficients depends on the energy and tube voltage. ξ θ and ξ φ are the geometrical offset corrections along θ and φ direction respectively which are determined experimentally.
The sub-equations (3.4) are tungsten characteristics in each solid angle, where the values of α and β vary by various tube voltages which are represented in tables 7 to 10 (appendix).
Each dependency of the S 00 EV and A EV , B EV , C EV , α and β coefficients was mapped using an independent variable order polynomial regression to get the final x-ray model. The S 00 EV , A EV , B EV , C EV , α and β coefficients are given in the appendix.

Comparison with other models
The accuracy of the spectral information obtained from the MARS x-ray model is evaluated by comparison with three other spectral models. The comparisons made include shape of spectrum, beam quality, and off-axis photon distribution (with TOPAS only).

Visual assessment of spectrum shape
The spectrum shape is the best parameter for qualitative visual assessment of differences between spectra of different models. It is a common assessment technique for the evaluation of new models. To compare spectra we used the central axis of the MARS x-ray model as it is the only data provided -7 -by IPEM78 and SpekCalc. Spectra at 117, 60, and 90 kVp tube voltage with 1.8 mm Al filter were used for this assessment. These spectra were normalized to the total number of photons. TOPAS model data [33] is only available for 117 kVp therefore TOPAS is not used in the 60 and 90 kVp comparisons.
Off-axis spectra at three different position between MARS x-ray model and TOPAS for 117 kVp are compared. These three positions are selected where the maximum difference with respect to the central beam axis is expected. These position are (a) Figure 5a shows 117 kVp spectra from the MARS x-ray model and other spectral models. The MARS x-ray model has higher intensity in low energies while it has equal intensity in the high energy range (greater than 35 keV) in comparison with the IPEM78. The MARS x-ray model matches well with SpekCalc in the entire energy range. TOPAS and the MARS x-ray model produce the same intensity at lower energies while TOPAS has lower intensity at higher energies (greater than 70 keV). Characteristic peaks of IPEM78 spectra are larger than the others.     shows that the MARS x-ray model, SpekCalc, and IPEM78 are in line with each other at 60 and 90 KVp except the characteristic x-ray at 90 kVp. The difference between tungsten characteristic lines is the main reason for the overall discrepancy between the source model presented in this work and others at higher tube voltages. This difference is likely to be due to using different fundamentals to create the characteristic lines in each model. For example, the characteristic lines of IPEM 78 were manually added to the spectra on the basis of the published values [34]. In the MARS source model, they were fully simulated in the Monte Carlo simulation phase by enabling "electron impact ionization". The MARS x-ray model was found closer to the TOPAS off-axis spectra with minor difference near the bremsstrahlung peak as shown in 6a and 6b.

Beam quality
Beam quality is the penetrating power of the x-ray beam. Commonly used metric of beam quality is HVL. It is used to characterize the shape of the spectrum [35]. The central axis photon spectra are used here for the HVLs computations. The HVLs are computed by iterating the material thickness until the intensity of the beam becomes half.
-9 -HVLs are calculated with filtration of 3.8 mm Al and 5.8 mm Al for tube voltage of 50, 60, 70, 80, 90, 100, 110, and 120 kVp. Only 117 kVp spectrum data is available from TOPAS for these calculations. Two filtration are used to determine the effect of filtration on the HVL. Mass attenuation coefficients from NIST XCOM are used [36]. Tables 1 and 2 show the computed HVLs and their relative difference with respect to MARS x-ray model data for given tube voltages with the filtration of 3.8 mm Al and 5.8 mm Al, respectively. The HVL comparison of the MARS x-ray model against IPEM78 showed that there is an increase in difference with the increase in tube voltage with 3.8 mm Al filtered spectra. On the other hand, this trend is opposite for 5.8 mm Al filtered spectra. The maximum relative difference was found to be 4 % at 120 kVp for 3.8 mm Al filtered spectra and 1.9 % for 5.8 mm Al filtered spectra. The comparison against SpekCalc showed a good agreement with a maximum difference of -3.4% for 50 kVp and 3.8 mm Al filtered spectra. The MARS x-ray model was found closer to the other spectral models with high filtration. This is due to a slight difference present between them in the low energy range.

Off-axis photon distribution
We compared the off-axis photon distribution from the MARS x-ray model against another Monte Carlo simulation tool, TOPAS. We used the off-axis photon distributions from both models at 11 cm with angular dimensions of   -11 -values are evidence of a good match between the energy spectrum as a function of angle of photon distribution of the MARS source model and TOPAS. A large part of this residual is due to the Monte Carlo noise in the TOPAS simulation. The contribution of low energy photons towards the cathode side (+φ) is more as compared to the anode side (-φ) and can be seen in low energy bins. On the other hand, high energy photon contributions are more towards the anode side as compared to the cathode side and can be found in the high energy window (71-117 keV). This effect is seen in both the MARS x-ray model and TOPAS off-axis photon distributions. These results indicate a good match between the MARS x-ray model and TOPAS generated off-axis photon distributions.

Comparison with experimental data
The MARS x-ray model is also evaluated against the experimental data collected in a MARS scanner by comparing with the photon counts and off-axis photon distribution.

Photon counts
We compared the photon counts estimated by the MARS x-ray model with experimental data collected in the MARS scanner. Tube voltage and tube current of x-ray tube were set to 80 kVp and 30 µA respectively. Measurements were taken at a source to detector distance of 18.68 cm. Exposure time was set to 120 ms for each open-beam image. An additional filter of 3.1 mm Al was used for these measurements. 720 open beam images were taken at each of five camera positions.
Pixel masking was applied to avoid the contributions from the noisy and malfunctioning pixels [7,37]. The mean counts of 720 open-beam images at each camera position were calculated. We divided the mean image into 32 groups having 4 rows in each group along the θ direction. We calculated the mean counts of each group. The midpoint of these groups at φ= 0 was considered for θ angle calculations. Similarly, with five camera positions, 160 data points along the θ direction were obtained. The experimental data was corrected for co-incident photon pulse pile-up [38]. The MARS x-ray model data was corrected for detector absorption efficiency based on 2 mm CdTe.
The experimental data points were least square fitted by three different degrees of polynomials linear, quadratic, and cubic. These polynomial fittings were assessed using sum of squares of error (SSE), R-square value, adjusted R-square, and root mean square error. The best fit of these polynomials was used for applying the fitting curve on the experimental data sets. These fitted curves were then compared with the data obtained from the MARS x-ray model. The deviation between the MARS x-ray model and fitted experimental data was quantified by the RMSD. Figure 8a displays the experimental mean counts (red dots) measured with the CdTe-Medipix3RX detector along the vertical direction (θ). Figure 8b shows the comparison of the MARS x-ray model counts along the θ direction with three experimental data sets. The RMSD values calculated for these experimental data sets are 0.0119, 0.0122, and 0.0116. An average RMSD value of 0.012 (approximately within 1 %) shows that our model has good agreement with the experimental data.

Off-axis photon distribution
Off-axis photon distribution was compared against the experimental measurement. As there is no experimental off-axis photon distribution data available in the literature (with desired specifications), we collected the experimental data with the Medipix3RX detector inside the MARS scanner.  Experimental measurements were taken at 45 camera position along the vertical direction (θ) with the CdTe-Medipix3RX detector. The tube voltage was 80 kVp and tube current was 9 µA. The exposure time was 275 ms. Source to detector distance (SDD) was set to 160 mm. 360 open-beam images were acquired at each camera position. The pixel masking technique is applied to remove the noisy and malfunctioning pixels from the measurement data. The best pixel was selected in each column of the chip having minimum absolute count difference with the mean of that column. Measurement at each camera position represents the one data point for the best-selected pixel. Quadratic polynomial regression applied to this data created a fitting surface. The fitting surface is normalized at the central beam axis to compare with the MARS x-ray model. The percentage difference relative to MARS x-ray model was calculated to quantify the difference between the MARS x-ray model and the fitting surface. Figure 9 shows the comparison between the MARS x-ray model and experimental results. There is an offset of 1.4 • along the vertical directions in the experimental results as shown in figure 9a. This variation could be due to off-set in the anode angle, tube assembly or the mechanical offset of the MARS scanner setup. Figure 9 (b) shows the experimental result with offset adjustment. The photon density is higher towards the cathode side, and it reduces from cathode +φ to anode -φ. There is also an off-set along the φ direction but it is hard to quantify because of measurement limitation along this direction. There is a maximum difference of 3 % at the extreme edges (θ) between model and experimental data. This difference could be due to the collimator shadow effect (absorption of photons at the collimator edge of the x-ray tube) or different amount of scattering inside the scanner or different off-focal radiations or error in the SDD. The small angle of φ is used for the comparison because (1) technical restrictions in the experimental setup (cannot move the chip along the φ direction) and (2) almost half of pixels in the chip were poorly bonded and not used.

Discussion
The accurate estimation of x-ray spectrum emitted from an x-ray source is an important study problem in medical physics. Especially, it is extensively critical in computed tomography applications such as correction of beam hardening artifacts. A parameterized semi-analytical x-ray source model presented here simulates the x-ray spectrum emitted from the x-ray tubes (Source Ray: SB-120-350) used in MARS scanners. The accuracy of this model in application depends on the accuracy of the geometrical parameters used in the Monte Carlo simulation, appropriate use of physics principles, computational data capacity, and interaction cross-sections databases. Equally important is the application of polynomial interpolation functions to describe the output as a function of tube voltage, energy, and spectral dimensions.
There is a chance of variation in geometrical parameters such as non-uniformity of the glass window, anode angle, but this assumption is not incorporated in this model. Added filtration of any type can be included easily using the Lambert-Beers law and materials attenuation coefficients. As the x-ray tube output is integrated over the image exposure time during the image acquisition in MARS scanner, therefore, voltage ripple is ignored in this study.
The generated spectra using the MARS x-ray model was assessed visually against IPEM78, SpekCalc, and TOPAS for 60, 90, and 117 kVp. They all matched with the other models with a minor difference of characteristics x-rays. The MARS x-ray model was observed closer to the TOPAS off-axis spectra with minor difference near the bremsstrahlung peak. The evaluation of the MARS x-ray model beam quality was conducted with the spectra ranging from 50 kVp to 120 kVp with two filtrations, 3.8 mm Al and 5.8 mm Al. There is a maximum difference of 4.0 % in HVL -14 -between the model and IPEM78 at 120 kVp with 3.8 mm Al filtration. However, this difference reduced to 1.2% for 5.8 mm Al filtration.
Off-axis photon distribution comparison between the MARS x-ray model and TOPAS showed a good agreement with a maximum RMSD value of 0.0075 (less than 1 %). The reason for this small difference could be due to different Monte Carlo code, underlying physics, and variance reduction technique.
The experimental evaluation of the MARS x-ray model showed a good agreement with the experimental measurements in terms of photon counts and off-axis photon distribution. Off-axis distribution comparison between the MARS x-ray model and experiment revealed an off-set of 1.4 • in the experimental data. This offset could be due to the mechanical offset of the MARS setup or there may be a tilt in the anode angle of the x-ray tube. A maximum difference of 3% was observed at the extreme vertical positions between model and experimental data after adjustment. This variation may be due to the collimator shadow effect. The absorption or the scattering from the edge of the collimator or penumbra affect or SDD error. The measured distribution includes the detector effects which need to be exclude for one-to-one comparison. The MARS x-ray model helps to compare theoretical knowledge with experimental results. As this model is not implemented yet, it is unlikely to determine the effects of this variation in material identification.
In conclusion, these comparisons showed that the MARS x-ray model has a good agreement with other spectral models and experimental measurements. The model can give a satisfactory estimate of the realistic beam in the diagnostic imaging energy range not only along the central axis but also off-axis.
This source model can be extended to different anode angles if required. In the planned future goals, the MARS project will develop a breast scanner and this model is applicable provided the same x-ray tube is used. With conventional x-ray tubes where Molybdenum/Rhodium targets are used this model can be reparameterized following the same methodology which will require further Monte Carlo simulations. The future study will focus on finding the ways to incorporate this model into the spectral reconstruction to achieve the ultimate benefit of better material identification and quantification.

Conclusion
A parameterized semi-analytical x-ray source model has been presented. This model has a good agreement with other models and measurement taken on MARS spectral scanner. We have applied this model to identify the geometrical offset in the beam profile taken by MARS scanner. In the future, this model will be used to improve the quantitative accuracy of spectral reconstruction techniques.
A Bremsstrahlung S 00 EV tables