Suppression of the CT Beam Hardening Streak Artifact Using Predictive Correction on Detector Data

The purpose of the research was to develop an automated program incorporating a predictive artifact correction technique (PACT) to correct for the signal deviations from metal beam hardening artifacts in Computed Tomography (CT) detector raw data. Thin-slice sequential CT scans were performed on a dosimetry head phantom using a Somatom Sensation 16 scanner to establish a ground truth image. Metal pins were then affixed to either side of the phantom at the three and nine o’clock positions to cause streak artifact in detector raw data and a subsequent streak image. The program automatically detected the extent of the overlap peaks in the detector raw data causing the artifact. It profiled a correction using adjacent projections so that the peak error could be corrected rather than simply being removed or smoothed by interpolation. The PACT algorithm modified raw data was then reconstructed on a SYNGO CT reconstruction workstation. This image was then compared against ground truth and that produced by commercially available metal artifact reduction projection completion and also a research based iterative technique. Qualitative results illustrate superior suppression of streak artifact in images using PACT when compared directly to tested projection completion methods but inferior to iterative reconstruction. Recovery of voxel data underlying the streak is also demonstrated to be quantitatively superior with PACT when referenced to the original ground truth image. Limitations were however detected with the threshold technique for initial localisation of the streak sources. The work still demonstrates the feasibility of this predictive artifact correction technique in correcting beam hardening affected voxel data without recourse to expensive additional options such as iterative reconstruction or dual energy that are not so commonly available in the clinical setting.


Introduction
CT since its inception has now become one of the mainstays of Diagnostic Imaging departments [1] replacing many imaging procedures that were previously performed radiographically and changing the practice of medicine by substantially reducing the need for exploratory surgery [2]. The significant growth of population doses from the 1980s to this decade is evidence if increase is the use of CT [3,4]. However this technology suffers significantly from the physical nature of the x-ray beam production, namely its polychromaticity. The problem with polychromaticity is that the CT number assigned to a voxel is not purely dependent on the tissue occupying that voxel but also on the location of the voxel in the slice [5]. This can result in deviations in measurements along a beam path depending on the materials the beam traverses, which can manifest as artifacts in the image. This artifact is well known in anatomical areas such as the posterior fossa and in fact the classical presentation of the dark band in this region has been titled the Hounsfield bar [6].  Figure 1 is an example of a beam hardening streak artifact caused by metal pins on either side of a head anthropomorphic phantom. These image defects can often be a result of more than one cause as this artifact can have both a beam-hardening component due to the aforementioned polychromaticity related signal error but can also contain a partial volume artifact due to the fan nature of the x-ray beam in the z direction [7]. They can also be further complicated by the use of outer slices in multi-slice detectors that causes smearing due to the nutating or wobbling slice acquisition [8]. CT manufacturers do have inbuilt correction to minimize beam hardening effects [9] however this is commonly defeated by the alignment of high-density objects within the scan field. Earlier efforts to address this directly in the detector data profiles by linear interpolation Metal Artifact Reduction (LI-MAR) [10] to remove the offending structures demonstrated early promise. Later additions of adaptive filtering and multi-dimensional adaptive filtering [11] in 2001 were combined with LI (LI-MAF) in 2004 [12] to good effect [13] but yet again the use of these metallic focused corrections still do not see common use. There have been many other approaches at implementing beam-hardening and metallic artifact suppression including Iterative [14,15,16,17] and Dual Energy (DE) reconstruction [18,19,20,21] techniques.
Iterative reconstruction (IR) is currently a hot topic research and development focus area in CT [22] but in fact CT iterative reconstruction is as old as CT itself. The original CT scanner produced by EMI used an iterative reconstruction method known as Algebraic Reconstruction Technique (ART) [23] but this was replaced by the computationally faster Filtered Back Projection (FBP) method in later CTs. Iterative reconstruction is computationally intensive and while it is common in nuclear medicine applications such as single photon emission computed tomography (SPECT) [1] it is only now seeing use in clinical CT applications with the introduction of fast modern computers. An example of IR uses forward projection of an initial or blank image to create an estimate raw data set that is then compared to the actual raw data. An error is computed and then this error term is back projected onto the image and the cycle repeats [24].
Variations include systems that go beyond simple noise or error iterations to systems that model the system geometry and noise. Different products operate in image domain, raw data domain or in the most modern commercial implementations, a combination of both. Commercial IR dedicated iterative computing engines primarily focus on producing low noise images enabling lower dose scanning [25] however there has been significant research into different iterative reconstruction techniques that have been specifically applied to beam hardening.
These methods do require segmentation of the initial reconstruction to facilitate processing without the metal for example, in order that a constrained image can be constructed to control the effect of adding back in the metal [26]. Iterative techniques can also use segmentation techniques to classify different materials [27] to isolate the respective polychromatic errors and iteratively reconstruct the data to minimise the errors or converge to an expected value. The number of iterations defines the level of convergence and the results clearly outperform filtered back projection and demonstrate the ability to suppress the artifact to varying degrees. Using the Selective Algebraic Reconstruction Technique (SART) is another iterative technique whereby non-metal voxels are selectively reconstructed using projection data that does not pass through metal [28].
The metal deletion technique (MDT) is another example of an iterative technique that uses thresholding of an initial image to localise the metal voxels and then a linear interpolation process removes the corrupted metal from the sinogram data to create a prior image. This prior image is then put through four successive iterations of forward projection followed by back projection, each cycle improving the estimation of the linear interpolated data to produce an image without the metal object. The final image produced is then free of metal artifact and the metal voxels from the initial reconstruction are added to the metal free image to complete the image restoration [28]. While there is no doubt about the efficacy of the approach in suppressing the beam hardening artifact it is computationally intensive in that it requires 9 instances of projection followed by an image voxel fusion step. The authors do note the limitation that by discarding the metal projection data results in a loss of spatial resolution near the metal objects. In simulations, the MDT technique scored best in qualitative comparisons with Linear Interpolation (LI or LI-MAR), Selective Algebraic Reconstruction Technique (SART) and Filtered Back Projection (FBP). MDT also scored best when quantitative average absolute errors were calculated with MDT, LI (LI-MAR), SART and FBP incorporating not only streak artifact suppression but also scores for noise reduction. In an interesting observation from Wang et al., the prior images at the start of these iterative cycles to produce streak free images can be a source of defects to the iterative corrective process even with thresholding or linear interpolation MAR (LI-MAR) techniques used to suppress the metal. LI-MAR images although free of metal still contain streak artifact elements in the image data and iterative processing of this data is therefore seeded with incorrect information [29].
The research in this work sees a return to direct manipulation of the detector raw data profiles while avoiding the use of DE and IR methods, as these require expensive additional hardware and software. These experiments also were designed to incorporate automated thresholding techniques to remove operator interaction or image based segmentation steps [13] prior to final image reconstruction. Rather than using interpolation to replace the offending data, a predictive approach using the information from adjacent non-overlapping rotational sample points was used calculate the expected value without using iterative processing. The objective being to recover the lost signal represented by the streak artifact rather than interpolating or smoothing out the streak, as this would provide more clinical information about the area under the streak. The overall goal was to investigate the possibility of a computationally simple but effective step that could be implemented into the current CT data processing chain with minimal impact and minimal modification of the CT reconstruction engine. This algorithm will be referred to as the Predictive Artifact Correction Technique (PACT) in this paper and will be compared to non-iterative LI-MAR & LI-MAF corrections and also the well-documented successful MDT iterative reconstruction technique.

Materials and Methods
In order to accurately test the efficacy of the correction algorithm, not only were images required demonstrating beam hardening artifact but also equivalent images without artifact to determine the effects of the correction compared to a ground truth.

Experiment Design
The phantom used was the Computerized Imaging Reference Systems, Incorporated (CIRS) Model 701 (sections 1-10) Atom Dosimetry Verification Phantom ( fig. 2). This anthropomorphic phantom was chosen as the skull detail and the apertures in the phantom for dose monitoring devices gave sufficient detail to demonstrate the effect of the processing on the underlying voxel data.
Hardened steel rods were placed on opposite lateral sides of the adult head phantom in order to cause a measurable beam hardening streak artifact. The rods chosen were 5.5mm diameter steel (Fe3C) rods, as these would cause clearly visible artifact while being practical to handle. It was determined by beam spectrum attenuation simulation, utilising the SpekCalc program [30] created by Poludniowski et al, that the 5.5mm Fe3C rods equate approximately to 29mm of Aluminium (Al) at the 120keV levels produced from Tungsten tube targets as typically used for such CT scans.

Scanner Protocol
In order to simplify the application of the algorithm to a single slice of raw data, sequential rather than spiral/helical scanning was preferred. Scanning was carried out on a Siemens Sensation 16 slice scanner. Scanning parameters used from a base of brain protocol were 12 slices, 120kVp, 1 second, 0.75mm slice thickness, 320mAs with a H31s reconstruction kernel. The only modification to the clinical head protocol was the use of the narrowest available slice thickness to minimize any partial volume effect in the sampled slice. The central slice was used to minimize as much as possible any nutation effects in the data acquisition from peripheral detector rows in the multi-slice detector array.  The sinogram plot of the detector data versus the rotational view is depicted in fig. 4. The plot demonstrates the position of the pins shifting across the detector as the system rotates and the largest peaks represent the overlap points. These maximal peaks are the point at which the strongest streak artifact is created. The scan started at a lateral overlap position thus the first and last peaks seen in the plot are in fact the same contiguous rotational overlap peak.

Correction Algorithm
The data was subsequently analysed and processed in MATLAB (The MathWorks Inc., Natick, MA, 2010) using the PACT predictive algorithm. This algorithm was designed to be a "zero click" pre-processing step that adds additional beam hardening correction to the raw data projections in advance of image reconstruction in a single non-iterative step. The DFD in figure 5 shows the processing steps of the predictive (beam hardening) artifact correction technique. In order to analyse and modify the raw data, thresholding of the sinogram data is required but attempting to use an absolute threshold on the projection attenuation profile has drawbacks especially when trying to localize peak structures not in the centre of the scan field as they become superimposed on the underlying oval head or body structures typically in the scan field. A median filtration technique is utilized in the PACT algorithm to remove the underlying body structure so that the peaks can be operated on in isolation. The algorithm uses an empirical threshold for normal adult attenuation that was determined from the study of the mean of the raw data from 120 adult head slices that were determined to be free of streaks by expert viewers. The algorithm can then determine the number of overlap or coincidence peak events in the raw data, which equates to twice the number of identified streaks in the image. However it must be noted that this method did exhibit limitations that will be addressed in the discussion.
The coincidence event as the high density structures come in and out of alignment does not affect just one projection and therefore the algorithm determines the range of projections in the sinogram that are affected. The sinogram in figure 4 demonstrated that the peak locations from the pins are not static as the system rotates. Due to the rotational sampling the pin alignment overlap range spanned almost 2 degrees (i.e. 12 views) during rotation and the close up representation of the sinogram in figure 6 demonstrates that this resulted in the peak travelling across three central detectors in this range.
The algorithm also uses the view-spread information in order to locate an adjacent projection where the high-density structures have just become divergent. The predicted correction factor was determined by comparing the maximum overlap peak view of an alignment of structures as shown in figure 7a and then using the neighbouring divergent view shown in figure 7b in order to calculate the correction factor to address the beam hardening deficit in the coincidence attenuation measurement.
A further challenge that required addressing was the changing non-symmetrical pattern of the overlap data peak.  The corrected raw data sets were processed on a Siemens Syngo workstation equipped with the VAMP Syngo Explorer CT reconstruction software [31]. The metal artifact reduction (MAR) algorithm and also MAR in conjunction with adaptive filtration (AF) were used for comparison purposes in evaluating the efficacy of streak suppression by the PACT algorithm. No additional corrections were enabled in the VAMP package when producing the PACT images in order to evaluate the new algorithm in isolation.

Data Analysis
The images were evaluated by analysing the voxel data at a location where the artifact would manifest as a horizontal streak.   figure 9(b). This particular slice and voxel profile location were selected due to the wide range of detail available at this scan position. The centre of the profile captures the wide range of air to bone CT values whereas the light grey circles to the left and right of the bone represent more subtle density changes against the tissue equivalent background. The phantom was not moved during the addition of the metal rods to ensure that the following scans and voxel value analyses would be consistent. Therefore Figure 9(a) and the voxel measurements of 9(b) represent the gold standard or ground truth before the artifact and any ensuing corrections affect subsequent image voxel data.

Results
The qualitative or subjective image recovery will be discussed now in the context of the actual quantitative corrected values when compared with the known truth-values. Figure 10 presents the artifact and four individual methods correction including PACT, and one combination method of correction including PACT in order to provide subjective assessment of the corrections.   Figure 10(b) illustrates the reconstruction of the streak raw data using the LI-MAR metal artifact reduction algorithm. Figure 10(c) illustrates the reconstruction using LI metal artifact reduction followed by MAF algorithms combined. Figure 10(d) illustrates this reconstruction using only the PACT algorithm. Figure 10 While linear interpolation clearly suppresses the impact of the primary streak artifact on the image, it causes the addition of streak effects in other peripheral areas of the image, particularly the frontal sinus area of the skull phantom, for both the LI-MAR 10(b) and LI-MAF 1o(c) images. Linear interpolation not only affected the data from the primary structures treated but also from other structures that aligned with these structures throughout the complete rotation which caused this additional image degradation. The PACT algorithm produced better visual streak correction as evidenced in 10(d) in the area originally affected by the primary streak. The target structures are visibly distinguishable and comparable to the original truth image unlike the LI-MAR and LI-MAF results. The localization of the correction also gave better overall performance on the image as its effects were limited to the primary streak affected area only resulting in no changes or deterioration elsewhere in the image. The MDT images using the original streak image or the PACT image as the initial seed for the iterative correction both look very similar. While some additional minor streaks are present extending from the central bone structure, the MDT corrected images demonstrate the best visual improvement in terms of overall artifact suppression. Note that while all the algorithms detected the primary streak between the metal rods, the PACT algorithm failed to detect the weaker artifacts extending from the metal rods to the long bone paths extending from the Parietal bones back to the Occipital bone.

Hounsfield number result
Each of the images was then quantitatively assessed using the data analysis method of figure 10. This visualised and quantified the effect of not only the artifact but also each of the subsequent correction methods.  Figure 11 illustrates the corresponding voxel value profile across the reference line in the streak-affected area of the image. The red dashed line in each profile represents the baseline truth determined from the initial scan in advance of application of the steel rods to the phantom as in Figure 9(b).
In figure 11(a) the solid black line represents the negative shift in Hounsfield numbers due to the streak and it should be noted that the pattern of the new profile maintains reasonable shape correlation to the original despite the offset. The other point of interest here is the amplitude of the offset is non-linear along the black line as it is material dependant, which is to be expected from the beam hardening phenomenon itself.
In figures 11(b) and 11(c), the black lines represent the effects of correction from the application of the LI-MAR and LI-MAF algorithms respectively. While these two results seem quite similar, closer scrutiny did reveal that there was a subtle numerical difference in the amount of modification of image voxel data with those two techniques. Linear interpolation clearly corrects for the negative offset of the streak artifact; however, plots 11(b) and 11(c) demonstrate that there is also a loss of information fidelity when compared to the original data along the profiles. The first loss of information evident is the insertion of false voxel values in the region where the artifact suppression takes place, which is particularly evident at the minimum and maximum data points in the truth profile. This is in addition to the voxel value distortion elsewhere in the image from the secondary streaks as evidenced by images 10(b) and 10(c). However, plots 11(b) and 11(c) also demonstrate that object displacement took place in the image with both these techniques when examining the location of the peaks and troughs along the detector axis raising further concerns over the validity and diagnostic efficacy of these techniques. Figure 11(d) demonstrates that the PACT algorithm not only corrects for the negative offset but also preserves the shape in terms of correct voxel HU values and also spatial position of the corrections of the voxel data along the profile. Air values are underestimated by approximately 10%.
The MDT image Hounsfield numbers using the original streak image for seeding are illustrated in figure 11(e) and not only is a 10% underestimation of air values noted but also a 100% overestimation of the bone Hounsfield numbers. Soft tissue values demonstrate reasonable correlation and only exhibit errors at transition boundaries to other materials. This signal ringing is a known issue of iterative reconstruction and is a result a mismatch between the original continuous data and the discretization of the data in the voxels through the use of iterative forward and back projection [32]. This can be mitigated by the use of finer grids than the voxel matrix for the iterative process [33] and while this can add further time required for reconstruction, it is a clear focus of research and development as indicated by the growth of model based commercial iterative reconstruction products.
In figure 11(f) the MDT image generated using the PACT image as an initial seed demonstrates similar behaviour to the MDT image using the streak image as a seed except for air values that now no longer demonstrate the 10% underestimation.
PACT demonstrates the closest visual correlation between the profile of the recovered voxel data and the original truth voxel data. In the next Figure, 13, the errors between the sampled profile voxel values in the reconstructed images when referenced to the original truth values are plotted and the impact of the offset and ringing errors are evident.
Each of the images was then quantitatively assessed using the data analysis method of figure 10. This visualised and quantified the effect of not only the artifact but also each of the subsequent correction methods.

Discussion
The changing rotational nature of the raw data requires that the algorithm continually realign itself to this data. The automation of the correction includes enhancements to analyse each individual view profile so that the location, peak amplitude, peak width and overall peak shape be accurately determined for correct application of the coincidence peak signal amplification to compensate for the beam hardening signal deficit.
This research focused purely on the streak artifact and the goal was to provide signal recovery so that any information underneath the streak artifact could be correctly visualized. The PACT algorithm employs a relatively simple automated technique that does not involve a large computational burden thus negating the need for additional complex and expensive iterative computing hardware and software. It also remains completely in the polychromatic domain and therefore has no requirement for additional expensive dual energy hardware and software. In our experiments, the predictive algorithm approach demonstrably outperforms other linear interpolation projection completion techniques in both terms of quality of the visual image output and also accuracy of the voxel signal recovery leading to improved detail in the image.
The results of the qualitative tests confirm that the PACT algorithm approach produces more improvement to the image than the other projection completion LI-MAR and LI-MAF corrective measures tested but the MDT iterative corrections clearly produce the greatest amount of visual artifact suppression. Quantitatively, the PACT algorithm produced better fidelity in the recovery of the voxel data HU values compared to all the other techniques and also better spatial accuracy of the detail recovered in the image particularly when compared to the LI-MAR and LI-MAF corrective measures tested while maintaining equivalent performance to MDT in this regard as the ringing errors noted were of an amplitude nature rather than a spatial one. Clearly the PACT method of not only assessing attenuation signal amplitude loss but also of profiling the shape of the amplitude correction to the coincidence event raw data profiles yields optimal results in terms of voxel HU value recovery.
The work focussed on adult head imaging and thresholds would need to be calibrated for Paediatric and Body imaging so that the detection system could scale to these body areas. It must also be noted however that MAR, MAR-AF and MDT detection did affect other streaks as metal aligned with long bone edges, the metal rod to Parietal to Occipital Bone streaks evident in figure 10(a) for example. The empirical threshold detection method employed failed to detect these additional weaker beam hardening streaks. This empirical detection method is clearly insufficient to use in a final commercial product and while the goal was to create a computationally efficient (i.e. time wise) algorithm that did not need to use iterative techniques, the author does acknowledge the use of at least one iteration to provide improved metal detection and then subsequent targeting of PACT error correction on the sinogram data using this initial detection would be optimal for raw data based corrective techniques as evidenced in other corrective work used in the majority of other corrective approaches [28,34,31].

Conclusions
The efficacy of image modification using automated PACT has been presented in this proof of concept work. The empirical threshold metal detection method initially employed failed to detect weaker streaks that fell below the threshold detection level and the literature repeatedly demonstrates that the use of at least one iteration to detect the origin of the streaks would be a far superior approach.
The corrective element of the PACT algorithm demonstrated superior recovery of the voxel information to all the other methods however the non-iterative metal detection was a significant weakness. Therefore the next logical step, whilst remaining in the raw data domain, is to take the corrective technique and apply it using a stronger detection method such as a single stage iterative process. A blended approach of this mathematical correction technique together with an Iterative Reconstruction approach could also yield the higher image quality gains of the IR artifact and noise reduction while leveraging the improved voxel value restoration of the PACT approach. However the initial attempts to use PACT images as the side for the iterative reconstruction MDT process did not yield significant improvement over the use of the conventional seed and this approach would need more investigation. Finally the experience with visualising the impact of the artifact on the voxel data creates another clear opportunity to use these profiling and correction techniques completely within the image domain resulting in an even higher ease of use and speed gain.
While the results of this work are admittedly less than perfect, the research does take a new approach to the problem that does not require the addition of expensive hardware that shows promise and the corrective element of the algorithm demonstrates a novel method of calculating superior beam hardening correction values in CT detector raw data.