In kV cone beam CT, the 2D projection on the planar detector is contaminated by significantly increased scatter due to the large irradiated object volume compared to conventional CT modalities with fan-beam geometry. The increased scatter can be comparable to the primary X-ray transmission in magnitude and thus severely degrades image quality by introducing quantification errors and image artifacts   .
The X-ray scatter consists of any secondary photons other than the incident primary photons and is mainly produced by direct photon interactions with medium, including Compton (incoherent) and Rayleigh (coherent) scattering. In Compton scattering, an incident photon transfers part of its energy to an outer shell atomic electron and is deflected from its incident direction. In Rayleigh scattering, the incident photon only changes its direction with no energy transferred. Analytic formulas have been discovered to calculate the distribution of the first-order Compton scatter using the Klein-Nishina formula and Rayleigh scatter using the Thomson formula . Besides, the scattered photon may experience multiple successive interactions with medium before arriving at the detector, which is referred as multiple scatter. In kV-CBCT, the distribution of multiple scatter is generally smooth and can be approximately calculated as a constant background     . As Compton scattering dominates photon interactions in kV-CBCT, in general the spatial distribution of scatter is mainly determined by the Compton single scatter (CSS) component  .
Many methods have been proposed and developed for determining the spatial distribution of scatter in CBCT including: measurements using beam-stopper-array  , moving blocker  or primary modulator  ; Monte Carlo simulations with optimized code package   , or GPU-based high computational efficiency platform   ; and analytical calculations with projection-based pencil-beam-scatter-kernel (PBSK) models   , or image-based single scatter calculation combined with multiple scatter estimation     . For the mathematical model based approaches for scatter correction (Monte Carlo simulations and analytical calculations), computational cost and quantification accuracy are the two primary competing factors under consideration for performance evaluation. In general, PBSK based superposition or convolution methods have the highest computational efficiency but lowest accuracy for scatter correction, while Monte Carlo simulations and analytical single scatter calculations are on the other end of the spectrum.
In this study, we propose an analytical model to calculate the spatial distribution dominant CSS component of the PBSK aiming to improve the computational accuracy of the PBSK based superposition method for scatter correction. The CSS component of the PBSK can be analytically formulated as line integration over the pencil-beam range. In the proposed method, the attenuation term is separated from the line integration and approximated by using an average attenuation path length and an effective attenuation coefficient. For the line integration of the unattenuated CSS, we derived a compact formula for easy implementation (Section 2.1). We also introduced the effective scattering center on the PBSK beam line as the point source of the unattenuated CSS so that the average attenuation could be efficiently calculated (Section 2.2). Exact calculations of the CSS of the PBSK with varying parameter settings and the integrated CSS from point scattering targets in a slab phantom were benchmarked to evaluate the performance of the proposed model (Section 2.3).
2. Methods and Materials
Figure 1 shows the schematic of calculating the scatter kernel from a pencil-beam. The distance from the X-ray source to the isocenter of the imaging system is s and the height from isocenter to the detector plane is h. The range of the pencil-beam is from x− to x+ with the origin defined at the isocenter. The pencil-beam has a finite size of cross section area Apb when projected at the detector plane. Denoting as the differential fluence of the CSS produced from a beam segment at x to point P with off-axis distance r on the detector plane, the spatial distribution of the CSS can be written as a line integration as Equation (1).
The attenuation path length and attenuation coefficient of the scattered photons when traversing the object are dependent of the external contour and heterogeneity of the object, which makes the scatter kernel object-dependent and spatially variant.
We assume that for a particular PBSK, the variations of the attenuation path length and attenuation coefficient of the CSS are smooth with respect to the beam segment x and the differential attenuation term can be equivalently replaced by an average value that is independent of the integration variable x. Mathematically, the above line integration is approximately calculated using Equation (2).
Hereinafter, we denote the first part in the bracket, on the right side of Equation (2), as the “unattenuated CSS component of PBSK” and the second term as the “post-scattering attenuation of CSS”. In the following, we will derive analytical formulas as approximate solutions for the two parts.
Figure 1. The schematic diagram for calculating the Compton single scatter (CSS) component of the pencil-beam scatter kernel (PBSK). The scatter kernel represents the scatter distribution on the detector plane from the line integration of scatter production along the beam line (range (x−,x+) as shown in the diagram).
2.1. The Unattenuated CSS Component of PBSK
Considering the beam segment as a point target, the number of Compton photons produced from segment B to pixel P, as shown in Figure 1, can be given by:
The parameter is the fluence of incident primary photons at the segment x, is the Compton interaction cross section (given by the
Klein-Nishina formula  ), is the linear electron density and is the solid angle subtended by the pixel area (Apxl) with respect to segment B. They are given by the following equations.
In above, is the projection fluence of the primary beam at the detector plane, µ is the linear attenuation coefficient of the primary photons, ρe is the volumetric electron density, r0 is the classic electron radius, hν and hv' are the energies of incident and scattered photons, respectively. The relation between hν and hv' is  :
where is the electron rest energy (~511 keV).
We define short denotations for the relative photon energy and the cosine of scattering angle as following:
Substituting Equations (4)-(9) into Equation (3), the differential fluence of CSS can be given by Equation (10).
The term g(p) is a function of p and relative photon energy E, given by Equation (11).
We derived that under the conditions of and , the integration of as in Equation (10) over can be approximated to be a compact form (Equation (12)) as a quartic function of off-axis distance r (see the details in the Appendix).
The coefficients Ck (k = 0, 2, 4) are independent of off-axis distance r and given by the following:
Upon given incident photon energy and imaging system geometry, these coefficients are specified by the pencil-beam range and can be pre-calculated as a lookup table with varying pencil-beam lengths.
In the coordinate system of the detector plane, we define as the calculating point of a detector pixel and as the projection position of a pencil-beam. The general formula to calculate the unattenuated CSS at contributed by pencil-beam at is:
2.2. The Post-Scattering Attenuation of CSS
We assume the average post-scattering attenuation term (as in Equation (2)) can be equivalently calculated using the attenuation path length from a specific point on the beam line to the calculating point r and the corresponding attenuation coefficient determined by the scattering angle formed between them. We define this specific point as the “effective scattering center” of the integrated unattenuated CSS. The position of the effective scattering center on the beam line can be solved from Equation (15) as an approximate solution (see the details in the Appendix):
Similar to the quartic formula coefficients Ck (k = 0, 2, 4) (Equation (13)), the position of effective scattering center as a function of r can be pre-calculated as a lookup table with varying pencil-beam lengths for given incident photon energy and imaging system geometry.
2.3. Model Evaluation
To evaluate the performance of the derived analytical solutions as described in the preceding sections, we used the specifications of the Varian On-Board Imager (OBI) system (Varian Medical Systems, Palo Alto, CA)  in the following calculations. For the OBI kV-CBCT system, the X-ray tube is capable of generating spectra in the range of 40 - 130 kVp and the most often used voltages in clinical protocols are 100 kVp and 125 kVp with the effective mean energies 51 keV and 58 keV, respectively. The isocenter-to-detector height (h) is in the range of 40 - 70 cm with the standard h being 50 cm. In addition, the source-to-isocenter distance (s) is 100 cm and the active detector area is 40 × 30 cm2.
2.3.1. The Quartic Formula
The accuracy of using the quartic formula (Equation (12)) to approximate the exact line integration of the unattenuated CSS is dependent of three parameters including the primary photon energy (hν), the isocenter-to-detector height (h) and the size of the object (x−, x+). The profiles of the unattenuated CSS were calculated and compared between the quartic formula solutions and the exact line integrations with varying hν (40 keV, 60 keV, 100 keV and 130 keV), h (40 cm, 50 cm, 60 cm and 70 cm) and (x−, x+) ((−5, 5) cm, (−10, 10) cm, (−15, 15) cm and (−20, 20) cm).
2.3.2. The Effective Scattering Center
As the object external contour plays a major role in determining the post-scattering attenuation of the CSS, we select two representative types of external contour―round and flat―for evaluating the performance of the effective scattering center method for calculating the average post-scattering attenuation. The profiles of CSS as a function of off-axis distance r along the radial and axial cross sections of a cylinder (corresponding to circular and flat contours, respectively) with three different sizes (10, 20 and 40 cm in diameter) are calculated using the exact differential attenuation term and the effective scattering center method, with comparisons made.
2.3.3. Phantom Experiment
A preliminary evaluation of applying the proposed model as an analytical solution for the CSS distribution from a volumetric object was performed on a slab water phantom (30 cm in thickness) with an incident cone beam X-rays (26 × 20 cm2 at the isocenter plane and ~40 × 30 cm2 at the detector plane). The primary photon energy is 58keV, which is the mean energy of the 125 kVp X-ray tube voltage used in the Varian OBI system. The SAD is 100 cm and the isocenter-to-detector height is 50 cm.
The benchmarked CSS distribution was calculated by conventional analytical method using point scattering target, in which the irradiation volume is discretized as voxels and the CSS contribution from each voxel to each detector pixel is calculated analytically with exact implementation of the Klein-Nishina formula followed by exact calculation for the attenuation term .
The relative root-mean-square-error (RMSE) in the CSS distribution is used as the metric to assess the accuracy of the proposed model. The RMSE is calculated using Equation (16), with the CSS distribution on the detector (pixel dimension is ) calculated by the proposed model and the benchmarked method, respectively.
The computations were performed in MATLAB (MathWorks, Inc., Natick, MA) with a single processor (Intel® i5, 2.4 GHz, 4 GB RAM).
3.1. Evaluation of the Quartic Formula
3.1.1. Dependence of Photon Energy
The profiles of CSS calculated by the quartic formula (Equation (12)) and the exact line integration for photon energies 40 keV, 60 keV, 100 keV and 130 keV are shown in Figure 2. The other parameter settings are s = 100 cm, h = 50 cm and (x−, x+) = (−10, 10) cm. The deviation of the quartic formula compared to the exact line integration increases with the off-axis radius and slightly increases with photon energy. The maximum differences at r = 20 cm are +1.3%, +1.6%, +1.9% and +2.2% for hν =40 keV, 60 keV, 100 keV and 130 keV, respectively.
3.1.2. Dependence of Isocenter-to-Detector Height
Figure 3 shows the CSS profiles calculated by the quartic formula and the exact line integration for four different h values (40 cm, 50 cm, 60 cm and 70 cm), with the same primary photon energy 58 keV and beam range (−10, 10) cm. The deviation between the quartic formula and the exact line integration rapidly decreases with h. The maximum differences at r = 20 cm are +6.1%, +1.5%, +0.5% and +0.2% for h = 40 cm, 50 cm, 60 cm and 70 cm, respectively.
3.1.3. Dependence of the Length of Pencil-Beam
Figure 4 shows the CSS profiles calculated by the quartic formula and the exact line integration for four different pencil-beam ranges (−5, 5) cm, (−10, 10) cm,
Figure 2. The profiles of the unattenuated CSS component of PBSK calculated by exact line integration (solid) and the quartic formula (dashed) for incident primary photon energy at 40 keV, 60 keV, 100 keV and 130 keV, with h = 50 cm and beam length 20 cm.
Figure 3. The profiles of the unattenuated CSS component of PBSK calculated by exact line integration (solid) and the quartic formula (dashed) for h = 40 cm, 50 cm, 60 cm and 70 cm, with incident primary photon energy 58 keV and beam length 20 cm.
Figure 4. The profiles of the unattenuated CSS component of PBSK calculated by exact line integration (solid) and the quartic formula (dashed) for pencil-beam ranges of (−5, 5) cm, (−10, 10) cm, (−15, 15) cm and (−20, 20) cm, with incident primary photon energy 58 keV and h = 50 cm.
(−15, 15) cm and (−20, 20) cm with the same hν = 58 keV and h = 50 cm. The ranges of pencil-beam correspond to nominal object diameters 10 cm, 20 cm, 30 cm and 40 cm, respectively. With larger pencil-beam range, the quartic formula increasingly matches the exact line integration. The maximum differences at r = 20 cm are +2.1%, +1.5%, +1.0% and +0.7% for beam length = 10 cm, 20 cm, 30 cm and 40 cm, respectively.
3.2. Evaluation of the Model of Effective Scattering Center
The effective scattering center position on the pencil-beam as a function of off-axis distance was determined using Equation (15) for three different pencil-beam ranges (−5, 5) cm, (−10, 10) cm and (−20, 20) cm, with the corresponding cylinder radius (R) of 5 cm, 10 cm and 20 cm, respectively. Figure 5 shows the attenuated CSS profiles calculated using the exact attenuation term and the average attenuation based on the effective scattering center model. The primary photon energy is 58 keV and h = 50 cm. The differences are in the ranges of [−0.5%, 0.5%], [+0.9%, +1.6%] and [+1.5%, +9.3%] along the radial cross section and [−1.4%, +0.5%], [−5.0%, +0.9%], [−11.3%, +1.5%] along the axial cross section for cylinder radius of 5 cm, 10 cm and 20 cm, respectively. The deviation of using the effective scattering center for average attenuation increases with the cylinder size.
Figure 5. The profiles of attenuated CSS component of PBSK along the radial and axial directions of long cylinders with different radius of 5 cm, 10 cm and 20 cm, calculated by exact attenuation term (solid) and using the effective scattering center model (dashed).
3.3. Total CSS Distribution of Incident Cone Beam on the Slab Phantom
In the benchmarked calculations, the phantom volume within the X-ray field consists of [26, 20, 30] voxels with voxel size 1 × 1 × 1 cm3 and in the model calculations, the irradiated phantom volume consists of [26, 20] pencil-beams with beam size 1 × 1 cm2. The CSS distribution on the detector plane consists of [80, 60] pixels with pixel size 0.5 × 0.5 cm2. Figure 6 shows the CSS distributions by the benchmarked calculation and the proposed model, respectively. The RMSE is 3.8%, indicating high agreement achieved. The computational time by the model was about 30 times less than the benchmarked method.
Figure 6. The CSS distribution on the projection plane (40 × 30 cm2, pixel size 0.5 × 0.5 cm2) with cone beam field (26 × 20 cm2 at SAD = 100 cm) incident on the slab phantom (thickness 30 cm) calculated by (a) the point target method as the benchmarked result and (b) the proposed analytical framework. The pixel value represents the scatter-to-primary ratio. And the RMSE is 3.8%.
The general form of PBSK applied by superposition or convolution methods for scatter calculation is usually obtained by measurements or Monte Carlo simulations using large slabs or disks and has symmetric, spatially invariant formulation  . However, the realistic objects in question for scatter correction usually deviate from slabs or disks in contours, size and medium homogeneity, causing asymmetric deformation of the PBSK. The method proposed in this study explores an analytical approach for correcting the spatial variation of the PBSK, yet only focusing on the CSS component which is generally dominant in spatial distribution over the other components of the PBSK.
The quartic formula as an approximation of the line integration of unattenuated CSS was derived under the conditions of and . In
kV-CBCT, the maximum photon energy is 130keV and the maximum value of E is 130/511 = 0.25. As shown in Figure 2, the accuracy of the quartic formula varies slightly with different photon energies. Because the high order terms of
are truncated in the approximation, the error of the quartic formula increases with the off-axis distance r. Assuming h = 50 cm and the object size less than 40 cm in diameter, the condition of may be violated at off-axis
distance r beyond 30 cm for beam segment (x) at the lower end of the pencil-beam and increased error can be caused by using the quartic formula. However, as the detector active area has a maximum size of 40 × 30 cm2, such increased error only happens at the periphery of the detection area. Using larger isocenter-to-detector height (h) can increase the accuracy of the quartic formula, as shown in Figure 3. Also the quartic formula matches the exact line integration better for larger beam length, as shown in Figure 4. With increasing beam
length, the value of is reduced for beam segments at the upper end
which have more contribution to scatter due to larger incident primary photon fluence, and thus the truncation error is reduced as well.
We introduced the concept of effective scattering center as the equivalent scatter point source of the PBSK for calculating the average post-scattering attenuation. An analytical equation (Equation (15)) was obtained to determine the position of the effective scattering center with a first-order approximation made. Figure 5 shows that the deviation of such method from the exact results increases with the size of the object. The subtle inaccuracy in determining the location of the effective scattering center introduces increased error in larger object because the attenuation path length is increased and plays an increasing role in modulating the magnitude of the scatter. By comparing the two different types of object contours (round and flat), the performance of the effective scattering center model has insignificant variations, as shown in Figure 5. It indicates that this model may be applied as a general solution with negligible dependence on the object contour.
The scope of the current study has not considered medium heterogeneity in analytical calculation of the CSS component of PBSK. However, as the object-specific variance of the PBSK is determined by the post-scattering attenuation term, the impact of medium heterogeneity may be accounted for by ray tracing the attenuation path length from the effective scattering center. In addition, the electron binding effect in Compton interaction was considered as a secondary effect to the distribution of CSS and a correction factor may be included in future work.
The method proposed in this study shows highly increased computational efficiency compared to the conventional analytical calculation method based on point scattering model. It is also potentially useful for correcting the spatial variant PBSK in adaptive superposition calculation for the purpose of scatter correction in kV-CBCT.
Conflicts of Interest Disclosure
This work awarded US patent US 9,615,807 B2, April 2017 .
1.1. Derivation of the Quartic Formula of φS(r) (Equation (12))
From Equation (11), g(p) can be rewritten as:
The coefficients fj (j = 1, 2, 3, 4, 5) are
Applying Taylor series expansion for the last three terms in Equation (A1):
By substituting Equations (A2) and (A3) into Equation (A1), g(p) then becomes:
The coefficients are given as following:
Because of , the series in Equation (A4) is convergent when and the high order terms (n > 2) are truncated hereinafter.
The Taylor series expansion for p (the cosine of the scattering angle, as in Equation (9)) is:
Then g(p) can be written in the form of as following:
Under the condition of and neglecting the high order terms (≥6) in Equation (A9), the differential form of Compton single scatter as in Equation (10) is given as a quartic function of r as following:
The coefficients ck(x) with k = 0, 2, 4 are given as following:
1.2. Derivation of the Position of the Effective Scattering Center (Equation (15))
First consider the scatter reaching the detector point on the pencil-beam’s axis (i.e. r = 0). The scatter produced from segment x is
. The attenuation path length is given by and the scattering angle is 0 such that the attenuation coefficient is µ. The effective scattering center position can be accurately determined by:
For r ≠ 0, the general form of is given by Equation (10) and the attenuation path length is determined by intersection of the ray with the object contour. It becomes complicated to accurately calculate the effective scattering center position using Equation (2). However, we can approximate the effective scattering center position for r ≠ 0 by introducing a correction factor based on the above equation for r = 0. Assume the attenuation path length is linear with x and change in the attenuation coefficient is negligible (for small scattering angle, the scatter photon energy by Equation (8) is very close to the primary photon energy). Meanwhile, approximate the term g(p) to the first order. Then the effective scattering center position can be approximately determined by Equation (15).
 Rührnschopf, E. and Klingenbeck, K. (2011) A General Framework and Review of Scatter Correction Methods in X-ray Cone-Beam Computerized Tomography. Part 1: Scatter Compensation Approaches. Medical Physics, 38, 4296-4311.
 Rührnschopf, E. and Klingenbeck, K. (2011) A General Framework and Review of Scatter Correction Methods in Cone Beam CT. Part 2: Scatter Estimation Approaches. Medical Physics, 38, 5186-5199.
 Siewerdsen, J.H. and Jaffray, D.A. (2001) Cone-Beam Computed Tomography with a Flat-Panel Imager: Magnitude and Effects of X-ray Scatter. Medical Physics, 28, 220-231.
 Kyriakou, Y., Riedel, T. and Kalender, W.A. (2006) Combining Deterministic and Monte Carlo Calculations for Fast Estimation of Scatter Intensities in CT. Physics in Medicine & Biology, 51, 4567-4586.
 Ingleby, H.R., Elbakri, I.A., Rickey, D.W. and Pistorius, S. (2009) Analytical Scatter Estimation for Cone-Beam Computed Tomography. Proceeding of SPIE 7258, Physics of Medical Imaging, Physics of Medical Imaging, 725839.
 Gong, H., Yan, H., Jia, X., Li, B., Wang, G. and Cao, G. (2017) X-ray Scatter Correction for Multi-Source Interior Computed Tomography. Medical Physics, 44, 71-83.
 Chen, Y., Song, Y., Ma, J. and Zhao, J. (2016) Optimization-Based Scatter Estimation Using Primary Modulation for Computed Tomography. Medical Physics, 43, 4753-4767.
 Thing, R.S. and Mainegra-Hing, E. (2014) Optimizing Cone Beam CT Scatter Estimation in egs_cbct for a Clinical and Virtual Chest Phantom. Medical Physics, 41, Article ID: 071902.
 Watson, P.G.F., Mainegra-Hing, E., Tomic, N. and Seuntjens, J. (2015) Implementation of an Efficient Monte Carlo Calculation for CBCT Scatter Correction: Phantom Study. Journal of Applied Clinical Medical Physics, 16, 216-227.
 Jia, X., Yan, H., Cervino, L., Folkerts, M. and Jiang, S.B. (2012) A GPU Tool for Efficient, Accurate, and Realistic Simulation of Cone Beam CT Projections. Medical Physics, 39, 7368-7378.
 Xu, Y., Bai, T., Yan, H., Ouyang, L., Pompos, A. and Wang, J. (2015) A Practical Cone-Beam CT Scatter Correction Method with Optimized Monte Carlo Simulations for Image-Guided Radiation Therapy. Physics in Medicine & Biology, 60, 3567-3587.
 Sun, M. and Star-Lack, J.M. (2010) Improved Scatter Correction Using Adaptive Scatter Kernel Superposition. Physics in Medicine & Biology, 55, 6695-6720.
 Spies, L., Evans, P.M., Partridge, M., Hansen, V.N. and Bortfeld, T. (2000) Direct Measurement and Analytical Modeling of Scatter in Portal Imaging. Medical Physics, 27, 462-471.
 Yao, W. and Leszczynski, K.W. (2009b) An Analytical Approach to Estimating the First Order Scatter in Heterogeneous Medium. II. A Practical Application. Medical Physics, 36, 3157-3167.
 Ding, G.X. and Coffey, C.W. (2010) Beam Characteristics and Radiation Output of a Kilovoltage Cone-Beam CT. Physics in Medicine & Biology, 55, 5231-5248.