Scatter Correction for Industrial Cone-Beam Computed Tomography ( CBCT ) Using 3 D VSHARP , a fast GPU-Based Linear Boltzmann Transport Equation Solver

The use of cone-beam computed tomography (CBCT) to inspect objects containing aluminum, Inconel and other metals is a challenge due to high amounts of scattered radiation in the projection data. Following in the footsteps of VSHARP, our kernelbased two-dimensional (2D) software scatter correction algorithm that operates on x-ray projection data, we introduce threedimensional (3D) VSHARP designed to offer improved accuracy and flexibility. The core of 3D VSHARP is fast finite element Linear Boltzmann Transport Equation (LBTE) solver that models photon transport within the object-of-interest and computes the scatter component of the projection data with Monte Carlo-like accuracy but greatly reduced computation times. 3D VSHARP utilizes an object model consisting of 3D material density maps based on a first pass reconstruction, and an imaging chain model describing the system geometry, the x-ray source components (polychromatic spectrum, beam filtration, and collimation hardware), and the x-ray detector components (flat-panel imager in this case). The 3D material density maps making up the object model may be further informed by Computer Aided Design (CAD) files. Validation against Monte Carlo (MC) simulations of overlapping metal blocks showed that 3D VSHARPand MC-computed scatter fractions agreed to within 3% RMS. The 3D VSHARP execution time was less than 1 sec per projection using a single NVidia TitanXp GPU (NVidia Inc., Santa Clara CA) while the MC simulations took 1566 CPU-hours using two Intel-Silver-Xenon, 2.2 GHz, 10x core CPU (Intel Inc., Santa Clara CA), in a Dell 7920 workstation (Dell Inc., Austin TX). Experimental CBCT scans were taken of an aluminum motorcycle cylinder head and Inconel turbine blade at 450kVp. Results showed that high quality reconstructions and volume renderings were obtained which may open up the possibility of in-line CBCT applications.


Introduction
Cone-beam CT (CBCT) has proven to be an invaluable tool for industrial imaging.The acquisition of extensive volumetric data in a single rotation removes the burden of object translation and allows for higher throughput.However, the wide-area beam generates a large amount of scatter in the projections that, if not corrected for, can produce images that suffer from cupping, shading, streaks, false inhomogeneities, and quantification inaccuracies [1].Kratz et.al. [2] evaluated the performance of three different hardware-based scatter correction techniques and a software-only method on several phantoms.The hardware-based methods included a beam-stop array [3], a primary modulator [4,5] and a linecollimator.The software-only method was Varex's 2D VSHARP algorithm utilizing adaptive kernels to deconvolve scatter from the projection data [6].The potential advantages of the software-only approach are that no extra hardware is required and throughput may be higher since pre-scans are not needed.The results showed that for low or medium absorbing objects, 2D VSHARP image quality was comparable to that provided by the hardware-based approaches, but that image quality was limited for the highly absorbing objects.The aim of this study was to improve the VSHARP algorithm to offer high quality results for all objects.
A weakness of the 2D VSHARP algorithm is that, since it operates on projection data and lacks prior knowledge of exact geometric shapes and material compositions, it may be less effective for heterogeneous, irregularly shaped objects of the type seen in industrial applications.A more accurate approach is to perform a first pass reconstruction and then estimate the scatter component of the projection data by performing x-ray transport calculations through a model of the object derived from the first pass reconstruction.The computed scatter signal is then subtracted from the original projection data before a second pass reconstruction is performed.In the past, traditional Monte Carlo (MC) methods that stochastically solve the Linear Boltzmann Transport Equation (LBTE) have been employed for this purpose and have produced excellent results.However, execution times have been prohibitively slow thus limiting the approach to academic and laboratory systems.Here we introduce 3D VSHARP, which utilizes finite element tools that rapidly solve the LBTE deterministically [7,8] with accuracies equal to that provided by traditional MC methods.We first describe the 3D VSHARP method, present validation results against traditional Monte Carlo techniques, and then show experimental CBCT reconstructions and volume renderings of an aluminum motorcycle cylinder head and Inconel turbine blade.

Boltzmann Transport Equation Solver
The flow of photons during a CBCT projection can be organized into three steps as shown in Fig. 1.In each step, the steadystate LBTE governs how photons behave: where is the angular fluence and is a source of photons that are inserted into position �⃗ with energy E traveling along direction � .The linear directional scatter coefficient � �⃗, ′ → , � ′ → � � describes the fraction of photons with energy ′ and direction � ′ that scatter into a new direction � with a new (lower) energy E. encompasses all scattering events and is the total linear attenuation coefficient accounting for all scattering and absorption events.The LBTE states that for a given position �⃗ and direction � , the streaming rate of photons (first term) is equal to the angular fluence generated by external sources (second term) plus that which scatters into �⃗ with direction � (third term) minus that which is removed by collisions (last term).The 3D VSHARP algorithm discretizes the problem in the spatial, energy and streaming direction domains maintaining the accuracy of a continuous solution by solving the LBTE through finite element means.The highly parallel implementation was in CUDA (NVidia, Santa Clara, CA) for execution on a GPU.

VSHARP Pre-Calculated Data
The 3D VSHARP scatter correction process requires three sets of pre-calculated information: 1) a discretized incident x-Ray spectrum, 2) total attenuation and scattering cross-sections, and 3) a detector response model.

X-Ray Spectra and Energy Grouping
To solve the Boltzmann transport equation using finite element means, the incident x-ray beam spectrum needs to be characterized and discretized into finite energy groups which form the basis for the subsequent transport calculations.Our experiments at 450kVp used 13 energy groups for both the primary and scatter calculations, and the 125 kVp experiments used 6 groups for the primary and 7 energy groups for scatter calculations.In both cases, the energy groupings were non uniformdenser at lower energies near the location of the source's characteristic x-rays and further apart over the higher-energy, smoother bremsstrahlung tail region of the spectrum.Detector

Cross-Sections
For each energy group E', an average total absorption cross section (′) and scattering cross-section � �⃗, ′ → , � ′ → � � is pre-computed to be used by the transport solver.The scattering cross section includes both interaction and directional ( � ′ → � ) probabilities of a Compton event.Since the transport computations for the first photon scattering event ("first collided") and subsequent scattering events are handled differently, different scattering cross-sections are required for each of those two cases [7].All possible energy group-to-energy group ( ′ → ) downscattering combinational possibilities are included in the scattering cross sections.

Detector Response
To obtain an accurate object scatter estimate, 3D VSHARP incorporates the flat-panel's detector response, which depends on the incident energy and angle of an x-ray photon, into its transport calculation.The relative energy deposited was determined using Geant4.9.6 simulations of photon pencil beams incident onto a scintillator at set angles and energies [9].In Geant4, the Livermore EM physics package was used and optical transport was not modeled.Shown in Figure 3 is the detector response to the average energy of each group at an incident angle of 0 degrees (perpendicular to the panel surface) for a 600um Cesium Iodide (CsI) scintillator.Up to about 60keV, the signal increases monotonically with incident photon energy.Afterwards, the signal decreases because the increased incident energy does not compensate for the reduced scatter and absorption of x-rays.

The 3D VSHARP Scatter Correction Process
Figure 4 diagrams a typical Varex Cone Beam Software Tools (CST) Pipeline layout of processing step "Plugins" used to run the 3D VSHARP physics-based object scatter correction.A first-pass CBCT reconstruction, depicted in Column 1, is performed to estimate the spatially dependent material properties and density of the object, which are then passed to the 3D VSHARP Boltzmann equation transport solver to estimate the primary and scatter components of each projection.For the first-pass reconstruction, typical steps include lag correction, 2D VSHARP scatter correction, normalization, beam hardening correction, filtering, and back-projection.Depicted in the second column are the 3D VSHARP scatter correction steps that include voxel material segmentation, volume down-sampling, physics-based forward projection, and projection-by-projection scatter subtraction.Due to the low-frequency nature of the scatter signal, only every 5-10 projections need to be calculated, with the scatter signal interpolated between those projections.Shown in the last column is the second-pass reconstruction that is run after scatter is subtracted from each original projection.These steps may include normalization, filtering, beam hardening correction, back-projection, Hounsfield unit scaling, and ring correction.

Monte Carlo (MC) Validation
Comparisons were made between 3D VSHARP scatter estimates and traditional Monte Carlo (MC) estimates.Transport was computed through a heterogeneous phantom (Figure 2) containing multiple shapes and materials.The 2 cm-thick aluminumwalled box with two solid metal rods inside (one aluminum and the other either Inconel or stainless steel) was chosen as a difficult case as there are high amounts of scatter along with steep changes in scatter-fraction across the phantom.The phantom's center was placed 100 cm from a point x-ray spectrum source, and 125kVp and 450kVp beams were tested.The 40x30cm 2 detector, with a 0.4x0.4cm 2 pixel size, was placed 50cm from the center of the object.Comparisons were made between 3D VSHARP output and that of the Monte Carlo packages, Geant4.10.4 and MCNP6.1 (Los Alamos National Laboratory, New Mexico, USA).For each test, total photon flux (scatter + primary) and primary-only flux at the detector were computed.The scatter-fraction (SF) was then calculated as the SF = (Total Flux -Primary Flux) / Total Flux.
In MCNP, the total photon flux per-pixel was calculated using an F2 tally placed at the center line of the U-axis of the detector surface.The F2 tally includes both scattered and primary photons and is a photon-surface flux counter that is normalized by the area of the pixel and the starting number of photons.An FIR5 tally was used to calculate the primary, uncollided photon flux at each of the pixel locations of the F2 tally.The scatter-only photon flux was then the result of the (F2 -FIR5primary) tally at each pixel along the center-line U-axis of the detector.

Industrial CBCT Measurements
To experimentally verify the new 3D VSHARP scatter-correction tool, industrial CBCT scans of an aluminum motorcycle cylinder head and a small Inconel turbine blade where performed at 450kV (Figure 5).The scans were taken on a Varex Scanning Services 450 kV system using a 1620-AN3 panel which has an active area of 41cm x41cm and a 200µm pixel pitch.For all scans, the source-to-imager distance (SID) was 1275 mm and 720 projections were taken.For the aluminum motorcycle cylinderhead, the source-to-object distance (SOD) was 850 mm while, for the turbine blade, the SOD was 1060 mm.The aluminum cylinder-head was reconstructed into 0.4 mm 3 voxels while the Inconel turbine blade was reconstructed into 0.1mm 3 voxels.In addition to object scatter correction, a "scanning system scatter correction" was implemented which accounted for scatter within the detector housing [6] and the effects of focal spot blurring.

Monte Carlo Validation Results
Computational times for the metal blocks phantom were around 10 6 time shorter for 3D VSHARP compared to traditional Monte Carlo (MC): A typical single projection MC simulation took 1566 CPU-hours using 2 Intel-Silver-Xenon, 2.2 GHz, 10x core CPU in a Dell 7920 workstation (Dell Inc., Austin TX); while a 3D VSHARP projection took about 1 second to finish.

Case 1: Aluminum and Stainless Steel at 125 kV
The MCNP-based and 3D VSHARP-based scatter fractions (SFs) are plotted in Figure 6 for the case where one internal metal block is stainless steel.It is seen that the peak SF value of 0.95 occurs behind the 14cm long aluminum box edges (pixel locations = +/-17cm), with a local maximum of 0.75 occurring at location -5 cm where the aluminum rod and Inconel-rod overlap.Excellent agreement is achieved between the two simulation methods with a root mean squared (RMS) relative difference between the VSHARP and MCNP estimates of 3.3%.

Case 2: Aluminum and Inconel at 450 kV
The MCNP-based and 3D VSHARP-based SFs are plotted in Figure 7 for the case where one internal metal block is Inconel.It is seen that the peak SF value of 0.7 occurs behind the 14cm long aluminum box edges (pixel location = +/-17cm), with a local maximum of 0.5 occurring at location -5 cm where the aluminum rod and Inconel-rod overlap.Excellent agreement is achieved between the two simulation methods with a root mean squared (RMS) relative difference between the VSHARP and MCNP estimates of 3.8%.

Aluminum Motorcycle Cylinder Head
Figure 8 shows a central axial slice of both a 3D VSHARP-corrected and a non-scatter-corrected volume.As can be seen, the 3D VSHARP reconstruction is crisper and has improved material homogeneity.The red box outlines a zoomed-in region where a thin aluminum wall is only visible following scatter correction.When the CBCT volumes are rendered in VG Studio (Volume Graphics, Heidelberg, Germany) [10], we see improved surface definition following 3D VSHARP scatter correction reflecting the improvement in the material uniformity (Figure 9).Artifactual holes appear in the side walls of the uncorrected volume when none are, in fact, present.The improved quality of the 3D VSHARP reconstruction is also illustrated in Figure 10, which are histograms of the CT units from the slices shown in Figure 8.As shown, following scatter correction, there is enhanced separation between aluminum and air voxels with the aluminum voxels also displaying a more symmetric profile.

Inconel Turbine Blade
Another challenging set of objects for CBCT industrial imaging are Inconel turbine blades.Representative axial slices are shown in Figure 11 for a 3D VSHARP-corrected and an uncorrected reconstruction respectively.The most prominent artifact in the uncorrected image is the bleeding of signal into the air region adjacent to the concave wall.Figure 12 plots the voxel intensity values along the line displayed in Figure 11 that extends from the concave wall through a central cooling channel and through the convex side.It is seen that in the corrected volume the concave wall is better defined with a steeper fall-off.Moreover, in the line profile from the uncorrected image, there is a 24% difference in CT numbers on each side of the cooling channel compared to only a 3% difference in the corrected image line profile.When the turbine blade is surface-rendered in Volume Graphics, the 3D VSHARP volume shows improved surface definition.As seen in Figure 13, the concave wall is well defined in the corrected volume, whereas in the uncorrected volume the cooling vents blend with the definition of the concave wall.

Discussion
The results from this preliminary study suggest that 3D VSHARP may produce sufficiently high quality CBCT reconstructions to perform inspection and other demanding industrial tasks including detecting porosity, defining surfaces and measuring interior features.The software-only approach has the inherent advantages of not requiring extra hardware components such as a modulator or beam stop array, nor needing pre-scans, and may open up the possibility of in-line CBCT applications.
The deterministic LBTE solver was shown to offer equivalent accuracy to traditional stochastic Monte Carlo algorithms with run times that are many orders-of-magnitude faster.3D VSHARP and Monte Carlo-based scatter fraction estimates agreed to  within 3.5% RMS for both the 125 and 450kVp cases.Experimental studies at 450 kVp of an aluminum motorcycle cylinder head and Inconel turbine blade showed that excellent image quality was achieved.VSHARP scatter correction led to reconstructions with improved homogeneity with more accurate CT numbers, and better separation between object voxel intensities and air voxel intensities.Subsequent surface renderings of scatter-corrected volumes showed significantly fewer artifacts, such as false holes, than the renderings of uncorrected images.Particularly encouraging was the reduced artifact adjacent to the concave wall of the Inconel turbine blade which is a region previously identified by Kratz et.al. [2] as needing improvement.
Further testing will be conducted to assess 3D VSHARP capabilities with comparison to 2D VSHARP performance.Although preliminary results (not shown) back up the previous study by Wang et.al. [8] using medical CBCT data showing that the 3D VSHARP approach leads to improved accuracy of reconstruction compared to a kernel-based method, more investigations are still required.This includes demonstrating the ability to achieve good reconstructions of heterogeneous objects that have components with different densities and atomic numbers.Aside from offering improved accuracy, another advantage that 3D VSHARP may hold over 2D VSHARP is that it may be more flexible since it does not require custom convolution kernels to be pre-computed for each use case.The parameters that do have to be pre-computed, including incident spectral shapes, energy groupings, scattering and absorption cross sections, detector responses, are fairly ubiquitous and should encompass many use cases The main optimization steps that need to be completed involve finding the best balance between run times and accuracy.Increased discretization in terms of energy groupings, object voxelization and image dexelization, generally increases accuracy but also increases run times.As demonstrated by Wang. et. al. [8] for medical applications, Pareto optimization coupled with a genetic optimization algorithm can be used to find a set of operating points that achieve both high image quality and adequate throughput.We intend to explore a similar approach cognizant of the fact that the higher diversity of objects encountered in industrial applications may produce a more complicated Pareto front.

Conclusion
Presented here is an experimental demonstration of 3D VSHARP, a CBCT scatter correction algorithm that calculates the Boltzmann transport through a first-pass CBCT reconstructed volume.3D VSHARP was demonstrated to accurately subtract the scattered X-ray component improving material uniformity throughout the CBCT reconstructed volume and increasing the quality of surface renderings.The algorithm is currently planned to be implemented into Varex's Cone Beam Software Toolset (CST).

Source 1 )
Trace photons from source to all voxels in the object 2) Photons are absorbed and/or scattered within object 3) Trace scatter flux from all voxels to each pixel Voxelized Object

Figure 3 :
Figure 3: Flat-Panel x-ray relative energy-dependent response for a 600um CsI scintillator.The incident angle was perpendicular to the surface of the detector.

Figure 4 :
Figure 4: CST Pipeline flow-chart showing the processing plugins used for the physics-transport based scatter correction 3D VSHARP.The blocks colored in orange will be key components of an upcoming Varex software release, Cone Beam Software Tools (CST2.1)

Figure 5 :
Figure 5: Objects scanned for the CBCT reconstructions at 450kV.Inconel turbine blade (left) and aluminum motorcycle cylinder head (right).

Figure 6 :
Figure 6: Scatter fraction (SF) comparison between 3D VSHARP (blue) and MCNP (red) for the phantom in Figure 2 for a 125 kVp x-ray spectrum source.

Figure 7 :
Figure 7: Scatter fraction (SF) comparison between 3D VSHARP (blue) and MCNP (red x) for the phantom in Figure 2 for a 450 kVp x-ray spectrum source.

Figure 8 :
Figure 8: Aluminum motorcycle cylinder head central axial-slice from the 3D VSHARP scatter-correction (left) and with no scatter-correction (right).Zoomed-in region on the top right showing restoration of a thin aluminum wall following scatter correction.For both slices: reconstruction voxel sizes are 0.4 mm 3 .CT units are calibrated so that aluminum intensities are 3500 with air set to 0 (display window = 6451, level = 3226).

Figure 9 :
Figure 9: Surface renderings (VG Studio Max 3.1) of the aluminum motorcycle cylinder head with 3D VSHARP correction (left) and without correction (right).The red arrows point to missing aluminum walls in the non-scattercorrected volume.

Figure 10 :
Figure 10: Histograms of the intensities of all voxels in 3D VSHARP slice 275 (orange) and the uncorrected (blue) reconstructions of the aluminum cylinder head.With 3D VSHARP, there is improved separation between aluminum and air voxels with the aluminum voxels also displaying a more symmetric distribution, and the air voxels displaying a tighter distribution.

Figure 11 :
Figure 11: Left: 3D VSHARP-corrected Inconel turbine blade axial slice.Right: Uncorrected slice where bleeding of the signal into air region adjacent to the concave wall is evident.The voxel size 0.1 mm 3 .The Inconel voxels were calibrated to 10000 CT units (display window = 1184, and level = 5970).The red arrows indicate the line profile plotted in Figure 12.

Figure 12 :
Figure 12: Line profiles across the corrected and uncorrected Inconel turbine blade reconstructions.The concave wall on the left part of the profile is better defined in the 3D VSHARP reconstruction and CT numbers are more uniform.

Figure 13 :
Figure 13: Inconel turbine blade surface renderings (VG Studio Max 3.1) following 3D VSHARP correction (left) and without correction (right).The concave wall is well defined in the corrected volume, whereas in the uncorrected volume the cooling vents blend with the definition of the concave wall.