A transmission measurement-based spectrum estimation method incorporating X-ray tube voltage fluctuation

Zengxiang Pan; Xu Li; Guoqin Lin; Peishan Qin; Zun Piao; Shuang Huang; Wangjiang Wu; Mengke Qi; Linghong Zhou; Yuan Xu

doi:10.21037/qims-22-1055

Original Article

A transmission measurement-based spectrum estimation method incorporating X-ray tube voltage fluctuation

Zengxiang Pan^#, Xu Li^#, Guoqin Lin, Peishan Qin, Zun Piao, Shuang Huang, Wangjiang Wu, Mengke Qi, Linghong Zhou, Yuan Xu

School of Biomedical Engineering, Southern Medical University, Guangzhou, China

Contributions: (I) Conception and design: Y Xu, Z Pan; (II) Administrative support: Y Xu, L Zhou; (III) Provision of study materials or patients: Y Xu, L Zhou; (IV) Collection and assembly of data: Z Pan, G Lin, X Li; (V) Data analysis and interpretation: Z Pan, G Lin, X Li; (VI) Manuscript writing: All authors; (VII) Final approval of manuscript: All authors.

^#These authors contributed equally to this work.

Correspondence to: Yuan Xu; Linghong Zhou. School of Biomedical Engineering, Southern Medical University, Guangzhou Avenue North 1838, Guangzhou 510515, China. Email: yuanxu@smu.edu.cn; smart@smu.edu.cn.

Background: The energy spectrum is the property of the X-ray tube that describes the energy fluence per unit interval of photon energy. The existing indirect methods for estimating the spectrum ignore the influence caused by the voltage fluctuation of the X-ray tube.

Methods: In this work, we propose a method for estimating the X-ray energy spectrum more accurately by including the voltage fluctuation of the X-ray tube. It expresses the spectrum as the weighted summation of a set of model spectra within a certain voltage fluctuation range. The difference between the raw projection and the estimated projection is considered as the objective function for obtaining the corresponding weight of each model spectrum. The equilibrium optimizer (EO) algorithm is used to find the weight combination that minimizes the objective function. Finally, the estimated spectrum is obtained. We refer to the proposed method as the poly-voltage method. The method is mainly aimed at the cone-beam computed tomography (CBCT) system.

Results: The model spectra mixture evaluation and projection evaluation showed that the reference spectrum can be combined by multiple model spectra. They also showed that it is appropriate to choose about 10% of the preset voltage as the voltage range of the model spectra, which can match the reference spectrum and projection quite well. The phantom evaluation showed that the beam-hardening artifact can be corrected using the estimated spectrum via the poly-voltage method, and the poly-voltage method provides not only the accurate reprojection but also an accurate spectrum. The normalized root mean square error (NRMSE) index between the spectrum generated via the poly-voltage method and the reference spectrum could be kept within 3% according to above evaluations. There existed a 1.77% percentage error between the estimated scatter of polymethyl methacrylate (PMMA) phantom using the two spectra generated via the poly-voltage method and the single-voltage method, and it could be considered for scatter simulation.

Conclusions: Our proposed poly-voltage method could estimate the spectrum more accurately for both ideal and more realistic voltage spectra, and it is robust to the different modes of voltage pulse.

Keywords: Spectrum estimation; voltage fluctuation; transmission measurement; ray-tracing algorithm; cone-beam computed tomography (CBCT)

Submitted Sep 30, 2022. Accepted for publication Mar 24, 2023. Published online Apr 07, 2023.

doi: 10.21037/qims-22-1055

Introduction

The X-ray energy spectrum, which displays the relative fraction of photons at different energies during an exposure, is a fundamental distinctive feature of the X-ray tube. It has a direct or indirect impact on correction for beam-hardening artifacts, Monte Carlo (MC) simulation, and others. For example, MC simulation (1,2), which is regarded as the gold standard in the field of particle transport, utilizes the source energy spectrum and its fluence map to initialize the photons of the source. However, the gold standard will no longer exist if the physical data loaded at the beginning are inaccurate.

For obtaining the energy spectrum of the X-ray imaging system, various methods have been proposed. Due to the limited count rate of the spectrum detector, it is difficult to measure the spectrum directly. Most of the methods are indirect, which are approximately classified into four categories: Compton-scattering measurement, MC simulation, spectra modeling, and transmission measurements. For compton-scattering measurement (3,4), the spectrum is calculated using the scattered spectrum and the corresponding scattering angle. The scattered spectrum can be generated via measuring the Compton-scattering photons directly. System calibration is complicated due to the absorption in the scatterer and the limited absorption efficiency of the detector. Most of the MC packages (5,6) provide comprehensive physics modeling and are well validated. However, in addition to the energy spectrum, there are other factors that influence the simulation results. Spectra modeling methods (7,8) usually adopt an empirical or semi-empirical approach that fits a parameterized model based on physical processes. Transmission measurements (9,10) estimate the spectrum using projections of a calibration phantom which requires dedicated hardware or workflow. In order to avoid such drawbacks, an indirect transmission measurement-based spectrum estimation method (11) was proposed, which expresses the estimated X-ray spectrum as the weighted summation of a set of model spectra. In this method (11), each model spectrum is the spectrum that the X-rays penetrate with a filter of a certain thickness and voltage. Then, an estimated projection can be computed via the ray-tracing method using the estimated spectrum. The final spectrum is calculated with these calibrated weights which minimize the difference between the reprojection and the raw projection. The method does not need any special phantom or high-precision measuring tools and can reduce the complexity of solving.

However, the method introduced above has a drawback in that it does not take the voltage fluctuation (12,13) into consideration. It only utilizes the model spectra at a preset voltage, and we therefore refer to it as the single-voltage method. In reality, it is difficult to generate an ideal voltage. The X-ray tube will not stabilize at a voltage during the entire exposure and will instead fluctuate around the preset voltage (14). There will be a process of rising, oscillating, and falling (15). The amplitude of the fluctuation is related to the quality of the machine and the X-ray tube. Different from the single-voltage method, we take the voltage fluctuation into account to estimate the spectrum as accurately as possible. Specifically, the model spectra that cover a voltage range are utilized to estimate the real spectrum of the X-ray tube, which expands the voltage range of the model spectra. For the sake of simplicity, we refer to the method proposed in this study the poly-voltage method.

Methods

Figure 1 shows the proposed workflow of the poly-voltage method. The first step is to acquire the raw projection data exposed by the X-ray tube. A modified Feldkamp-Davis-Kress (FDK) algorithm (16), one of the classical analytic reconstruction methods, integrated with geometric calibration results (17), is used to reconstruct a three-dimensional (3D) image. The second step is to load the reconstructed image and the initial estimated spectrum into the ray-tracing algorithm for calculating the estimated projection, a process referred to as reprojection. In order to minimize the difference between the raw projection and the reprojection, we set the difference as the objective function, and it is optimized for iteratively updating the weights of the model spectra. Finally, the best estimated spectrum is calculated with the model spectra and the corresponding weight.

Figure 1 The workflow of the poly-voltage method. 3D, three-dimensional.

In the poly-voltage method, we take the voltage fluctuation into account. The model spectra are produced from various voltages, and the spectrum can be represented as the weighted sum of the model spectra as follows:

$S_{r e a l} = \sum_{i = 1}^{M} \sum_{j = 1}^{N} w_{i j} S_{i j}$ [1]

Here, S_real is the real spectrum produced by the X-ray tube during exposure; M is the number of voltages, which represents the voltage range of the model spectra; N is the number of model spectra at each voltage, which represents different thicknesses of the filter; S_ij is the model spectrum at i^th voltage with j^th thickness of the filter; and w_ij is the weight of S_ij. At the beginning, the initial estimated spectrum can be calculated with the weights set randomly, and the weights are updated with an iterative optimization algorithm.

Voltage fluctuation

From the technical standard (18,19), we know that the deviation of the tube voltage should not exceed ±10% or ±10 kV. In addition, according to the relevant literature (20,21) and the tube voltage pulse measured in our in-house cone-beam computed tomography (CBCT) system, we can know that the phenomenon of voltage fluctuation is real. In an ideal situation, the exposure voltage will remain stable at the preset value during the entire exposure, and the voltage will increase and decrease very quickly at the beginning and end of the exposure, respectively. Focusing on the tube voltage pulse measured in our in-house CBCT system, we can observe that the tube voltage will rise to a value that is higher than the preset value at first and then gradually decrease to the preset value. The difference between the ideal voltage pulse and the real voltage pulse can be clearly seen in Figure 2. Of the pulses shown in Figure 2, the first pulse is obtained by fitting the actually measured pulse. For convenience, the second value is copied from the first. Based on the above phenomenon, a single real exposure can be interpreted as many short exposures at different voltages. In other words, the actual spectrum of the X-ray tube can be expressed as the weighted summation of numbers of model spectra from different voltages.

Figure 2 The ideal voltage pulse and the real voltage pulse.

The model spectra used in the single-voltage method are generated at the preset voltage with different thicknesses of the filter. The high-energy photon that exceeds the preset voltage will not be simulated, and the estimated spectrum does not match the real spectrum well. Therefore, we propose the poly-voltage method to cover a wide voltage range of the model spectra.

Model spectra

There are a variety of ways to generate model spectra, such as the MC simulation (5,22,23) and SPEKTR 3.0 software (24). In our work, we generate model spectra using SPEKTR3.0 software, which was developed to calculate X-ray spectra using the TASMICS (tungsten anode spectral model using interpolating cubic splines) algorithm (25) based on the TASMIP (tungsten anode spectral model using interpolating polynomials) spectral model (26). As for the process of generating model spectra in the spectrum estimated method, the tube voltage, material type, and thickness of the filter need to be defined. The tube voltage is included in a certain range in the “Model spectra mixture evaluation” section. The model spectra are filtered with 3, 4, 6, 8, 10, 12, and 18 mm of aluminum at each voltage, respectively. If one can use a broad range of the thicknesses of the filter, and sample more intensively, the accuracy of the estimated spectrum will improve; however, the calibration time may also increase. The requirement of the model spectra is that the softest model spectrum should be softer than the true spectrum while the hardest model spectrum should be harder than the true spectrum. For simplicity, we used the filter thickness (3, 4, 6, 8, 10, 12, 18 mm) from the reference which can obtain an acceptable spectrum.

Ray-tracing algorithm

After phantom image reconstruction, the ray-tracing algorithm is used to obtain the reprojection. The ray-tracing algorithm is based on Beer-Lambert law (27). It is expressed in Eq. [2]:

$I = I_{0} e^{- μ L}$ [2]

where I and I₀ are the incident and transmitted X-ray intensities, L is the length of the X-ray through the phantom, and µ is the linear attenuation coefficient of the material, which is a function of the incident X-ray photon energy. In Eq. [2], the X-ray beam is single energy and the phantom is composed of a single homogeneous material, which is inconsistent with reality. In reality, X-rays are multienergy, and a phantom is composed of a variety of materials. We therefore apply Beer-Lambert law as in Eq. [3]:

$I = \sum_{i} I {}_{i}e^{- \sum_{j} μ_{i j} L_{j}}$ [3]

where I_i is the normalized flux of the i^th energy segment, µ_ij is the linear attenuation coefficient of the photon in the i^th energy segment when passing through the j^th voxel, and L_j is the length that the X-ray penetrates the j^th voxel. For implementing Eq. [3], Siddon’s algorithm (28,29) is firstly used to calculate the length of each X-ray passing through each voxel. The attenuation coefficients of the different energy rays in each voxel are from the table of the National Institute of Standards and Technology (NIST) database (30). Finally, the reprojection is calculated via Eq. [3].

Updating weights

To make the reprojection match the raw projection well, we optimize the weights of the model spectra constantly until a preset number of iterations has been completed. Prior to this study, we compared particle swarm optimization (PSO) (31,32), genetic algorithm (GA) (33,34), and the equilibrium optimizer (EO) (35-37) algorithm. We found that EO yielded the best results based on the stability of the solution and iteration speed. Therefore, we used EO in this study. EO is inspired by the control volume mass balance model for estimating both dynamic and equilibrium states. In EO, the concentration of each particle is a solution, and the change in concentration of each particle acts as a search agent. The search agents update the concentration randomly from the best-so-far solutions, namely equilibrium candidates, and reach the equilibrium state as the optimal result. In order to keep the solution physically meaningful, a nonnegative constraint is applied to the weights. In this study, we execute EO in MATLAB (MathWorks, Natick, MA, USA). The pseudocode is presented below and can be briefly summarized as follows.

First, we initialize the number of independent runs R_n=1; the number of random number P_n=25, the maximum number of iterations M_i=360, and the number of weight W_d is equal to the number of model spectra. Then, line 2 of the pseudocode initializes the weight combination set C, and each weight value is in the interval [lb, ub], lb=0 and ub=1. From lines 3 to 7, we obtain the optimal four groups (C_eq1, C_eq2, C_eq3, C_eq4) from the weight combination set (C) according to fitness. The fitness fit(C(i, :)) at line 5 represents the difference of the middle row profiles between the estimated projection (p_est) and the raw projection (p_raw) for the i^th weight combination. The specific expression is shown in line 5 and is equal to the sum of the absolute values of the difference. We build the equilibrium pool (C_{eq_pool}) with C_eq1, C_eq2, C_eq3, C_eq4, and C_{eq_ave}at line 9. C_{eq_ave} is the arithmetic mean of the mentioned four particles. Finally, line 12 generates the new weight combination set (C) according to the equilibrium pool (C_{eq_pool}). See Faramarzi et al. (35) for the specific generating process.

Line 1: Initialize R_n, P_n, M_i, W_d;

Line 2: Calculate C = rand(P_n, W_d) × (ub − lb) + lb;

Line 3: While iteration < M_i;

Line 4: For i=1:P_n;

Line 5: Fit(C(i, :)) = sum(abs(p_est - p_raw));

Line 6: Update C_eq1, C_eq2, C_eq3, C_eq4;

Line 7: End;

Line 8: C_{eq_ave} = (C_eq1 + C_eq2 + C_eq3 + C_eq4)/4;

Line 9: C_{eq_pool} = {C_eq1, C_eq2, C_eq3, C_eq4, C_{eq_ave}};

Line 10: For i=1:P_n;

Line 11: C_eq = C_{eq_pool} {rand(size(C_{eq_pool}, 1), 1)};

Line 12: Generate the new weight combination set C;

Line 13: End;

Line 14: Iteration = iteration + 1;

Line 15: End while.

Experiments and evaluation

In this section, we described how the model spectra mixture evaluation and projection evaluation verified the correctness of the poly-voltage method and identified the voltage range of the model spectrum. The aim of using phantom evaluation was to not only prove that the poly-voltage method not only provides accurate reprojection, but can also the accurately estimate the spectrum. The aim of the scatter evaluation was to assess the difference between the estimated scatter signals using the estimated spectrum via the two methods and to determine whether it is necessary to load the spectra via the poly-voltage method in MC simulation for scatter simulation.

Model spectra mixture evaluation

In this evaluation, three modes of X-ray tube voltage pulse including the ideal pulse mode, long pulse mode, and short pulse mode were used to highlight the difference among the spectra estimated by the different methods. As shown in Figure 3, the ideal pulse mode was set to 90 kV during the whole exposure, and the other two were the same as the real voltage pulse. Both the preset voltage of the long and the short pulse mode were 90 kV. The former slightly oscillated around 93 kV, whereas the latter slightly oscillated around 90 kV. In this evaluation, the reference spectrum was composed of the model spectra filtered with 10 mm of aluminum.

Figure 3 The voltage pulse of the three X-ray tube modes.

Voltage range refers to the voltage range covered by the model spectra. For instance, we set 90 kV as the exposure voltage and 10% as the voltage range. This means that the voltage of the model spectra covered from 81 to 99 kV. Figure 4 shows a sample of the model spectra in the poly-voltage method and single-voltage method. One line in the figure indicates as one model spectrum, and the model spectra which have the same end point are generated at the same voltage. For the model spectra, there were 6 model spectra at 1 voltage, which were filtered with 3, 4, 6, 8, 12, and 18 mm of aluminum, respectively. The number of model spectra used in the poly-voltage method was larger than that in the single-voltage method. Moreover, the voltage range was sufficiently wide to cover the peak of the voltage pulse.

Figure 4 The model spectra in this study. A single line represents a single model spectrum; 6 model spectra at 1 voltage were filtered with 3, 4, 6, 8, 12, and 18 mm of aluminum, respectively. PV spe. (blue), the model spectra used in the poly-voltage method; SV spe. (black), the model spectra used in the single-voltage method. PV, poly-voltage method; SV, single-voltage method.

Projection evaluation

To verify the difference between two reprojections obtained via the two methods, a 3D phantom image of a pelvis was used in the ray-tracing algorithm. The phantom matrix and the size of the voxel were 512×512×70 and 0.1172×0.1172×0.3 cm³, respectively. The source to axis of rotation distance (SAD) and axis of rotation to image distance (AID) were set to 100 and 50 cm, respectively. Then, we generated reprojections via the ray-tracing algorithm and compared the difference between the generated projections.

Phantom evaluation

A digital phantom simulation and a polymethyl methacrylate (PMMA) phantom experiment were performed. A digital cylindrical water phantom was used in the numerical simulation. The phantom matrix was 780×780×620, and the size of voxel was 0.2×0.2×0.2 cm³. The reference projection is the reprojection of the phantom via the ray-tracing algorithm using the reference spectrum, which is combined with model spectra with consideration to the measured voltage fluctuation. We executed the PMMA phantom experiment with a 6-degrees of freedom (DoF) robot-based in-house CBCT system as shown in Figure 5A. In this configuration, there were two robotic manipulators upon which the X-ray tube and the detector were respectively mounted. The optional voltage range was from 40 to 140 kV, and the optional tube current range was from 10 to 120 mA. The detector resolution was 1,920×1,536 with a physical size of 24.4×19.5 cm². The PMMA phantom consisted of an upper cone and a lower cylinder, as shown in Figure 5B. The height of the cone and the cylinder were both 90 mm, and the diameter of the bottom of the cone was 140 mm, equal to the diameter of the cylinder. The raw projections were obtained from the 6-DoF robot-based in-house CBCT system. The SAD and AID were 100 and 50 cm, respectively. A total of 360 projections were equally acquired in 360 degrees with the tube voltage of 85 kV. In order to minimize the impact of scattering, the scanner was equipped with a narrow collimator. Therefore, the contribution of the scatter signal to the raw projections could be neglected. The 3D image was reconstructed via the FDK algorithm. Then, we calculated the two spectra estimated via the poly-voltage method and single-voltage method according to the workflow presented in Figure 1. The resulting spectra were utilized for beam-hardening artifact correction (38-41). The artifacts could be corrected using the water correction method, which transforms the total attenuation values from polychromatic to monochromatic using polynomial-based nonlinear mapping.

Figure 5 CBCT System and PMMA phantom used in the phantom evaluation. (A) Six-DoF robot-based in-house CBCT system established in the laboratory. (B) The PMMA phantom consists of an upper cone and a lower cylinder. CBCT, cone-beam computed tomography; PMMA, polymethyl methacrylate; DoF, degrees of freedom.

Scatter evaluation

In this evaluation, the MC simulations of CBCT projections were developed using the GPU-based MC package (gMCDRR) by Jia et al. (42). In gMCDRR, the detector can collect and distinguish the primary signal and the scatter signal. Multiple GPU threads are launched to transport photons simultaneously. Each photon is transported from the X-ray source until it exits the phantom or is absorbed in the phantom. Within each thread, the state of the X-ray source photon is initialized by the energy spectrum and fluence. The source photon energy must be generated according to the spectrum. As for the fluence, it is used to specify the probability density of a photon traveling toward the detector after emission from the source. If a photon exits from the phantom, its energy is recorded at a corresponding detector pixel in either the primary counter or the scatter counter depending on whether a scattering event has occurred. In this study, we used a digital water phantom and the gMCDRR to simulate the scatter signals with two estimated spectra calculated in the PMMA phantom experiment. The phantom matrix was 128×128×96, and the voxel size was 0.1×0.15×0.15 cm³. The detector resolution was 512×384 with size of 40×30 cm². The SAD and AID were set as 100 and 50 cm, respectively. The total number of photons we simulated in gMCDRR was 2×10¹².

Performance comparison

In the model spectra mixture evaluation, for quantitative comparisons, absolute error and percentage error of each energy bin were calculated to show the difference between the reference spectrum and the estimated spectrum. Furthermore, the mean energy difference (MED) and normalized root mean square error (NRMSE) were used to evaluate the effect of the voltage range on the accuracy of the estimated spectrum. MED and NRMSE can be defined as follows:

$M E D = \sum_{n = 1}^{N} [E_{n} \cdot (S_{e} (n) - S_{r e f} (n))]$ [4]

$N R M S E = \sqrt{\frac{\sum_{n = 1}^{N} {(S_{e} (n) - S_{r e f} (n))}^{2}}{\sum_{n = 1}^{N} {(S_{r e f} (n))}^{2}}}$ [5]

where S_e(n) is the normalized flux of bin n in the estimated spectrum, S_ref(n) is the normalized flux of bin n in the reference spectrum, E_n is the photon energy of bin n, and N is the total bin number of the spectrum. In the projection evaluation, profiles of the projections and the percentage error graph were compared in a quantitative analysis. In the phantom evaluation, beam-hardening artifacts appeared as cupping artifacts on reconstructed images of the PMMA and water phantom. The image appeared dark in the middle and light in the periphery. With the profile of the middle row of the image, a concave curve appeared with the two sides high and the middle low. The deeper the concave curve was, the more serious the artifact. We determined the severity of artifacts qualitatively by observing the degree of sag of the profile (i.e., flatness). The spectra generated by the two methods were used to perform beam-hardening correction on the projection, and the severity of cupping artifacts in the two reconstructed images were compared qualitatively to verify which method is more accurate. We also calculated the NRMSE of the middle row profile between the estimated projection and the reference projection, as well as the NRMSE indices between the estimated spectrum and the reference spectrum. The accuracy of the estimated spectrum from the two methods and the deviation of the reprojection were then compared using the NRMSE index. In scatter evaluation, we compared the profiles of the scatter signals generated using two spectra directly, and the percentage error was used as the quantification index of the difference.