Beam hardening correction based on image noise statistics

Beam hardening artefacts in x-ray computed tomography (CT) result from energy-dependent attenuation of x-rays in matter and cause inaccuracies in reconstructed 3D volume data. Effects due to beam hardening can easily be corrected if the detected spectrum of x-rays is known after having passed through an object. Conventional scintillator x-ray detectors, however, are incapable of measuring spectral information directly. The innovative idea of this paper is to extract information about the detected spectrum from image noise statistics and to estimate the spectrum using a new semi-empirical model function for a partially absorbed x-ray spectrum depending on a single unknown variable. The beam hardening correction factor is thus determined for each image pixel prior to the 3D reconstruction, does not require knowledge of the material of the CT-scanned object and is determined by modelling physical effects directly, without relying on an iterative approach or elaborate image processing.


Introduction
X-ray computed tomography (CT) imaging has emerged as a powerful tool in metrology for mechanical components.This non-destructive imaging technique allows for detailed 3D visualisation and measurement of internal structures and features within mechanical parts.High-resolution x-ray CT scans can reach sub-micron accuracy and provide precise data on dimensions [1], detect defects [2] and material properties, aiding in quality control, failure analysis, and design optimisation.X-ray CT is an invaluable asset, especially for modern additive manufacturing, in ensuring the reliability and performance of mechanical components in various industries.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Nevertheless, the process from capturing two-dimensional (2D) x-ray images (projections) of an object to reconstructing the grey value (gv) of each voxel within the scanned volume and to finally determine the position of the surface of the object to render an accurate virtual 3D volume of a mechanical part is complex and does neither take into account all physical interactions between x-ray radiation and matter nor the exact geometry of the system.Reconstruction algorithms like the widely used filtered back projection algorithm of Feldkamp et al [3] assume that the images result from cone-beam (CB) projections taken from a circular orbit of a single point x-ray source and that the x-rays travelling through matter experience an exponential attenuation with material thickness.In addition, the algorithm is strictly exact only within a plane containing the source and perpendicular to the rotation axis.Thus, well-identified artefacts appear in the reconstructed volume like artefacts due to the CB geometry outside the plane of exact reconstruction, artefacts due to the wave nature of x-rays known as phase contrast, or artefacts due to the beam hardening which is the subject of this paper.A more comprehensive summary of artefacts can be found in [4].

Beam hardening
Beam hardening artefacts in x-ray CT imaging arise from the poly-energetic nature of x-rays emitted from x-ray tubes which are commonly applied as sources for CT systems: since the radiation is issued from colliding electrons onto a target, a broad continuous x-ray spectrum is generated.The attenuation coefficient of x-rays in matter, however, varies with the xray energy and material composition and is generally strong at low energies and weaker at higher energies.Therefore, the attenuation of a continuous spectrum in matter is not strictly exponential with the material thickness as required by the reconstruction algorithm, causing inaccuracies in reconstructed 3D volume data.Addressing beam hardening is thus crucial for enhancing the accuracy and reliability of x-ray CT imaging.
Several advanced techniques have been developed to mitigate the effects of beam hardening.One common approach called linearisation consist in correcting the gv of the x-ray images before reconstructing the 3D volume.This correction, which usually requires prior knowledge about the part, can be performed on the basis of simulations or calibration measurements on a thickness standard (according to VDI 2 / VDE 3 guideline 2630-1.3section 6) made out of the same material as the mechanical part being analysed [5,6].While this is an easy approach for single material parts, a more complex material decomposition technique is being increasingly utilised for multi-material parts [7,8]: by identifying and isolating individual materials according to the CAD (computeraided design) model of a part, material decomposition helps in accurately estimating the material-specific attenuation coefficients and in correcting the gv of each x-ray projection.
Additionally, iterative reconstruction algorithms with or without prior knowledge of the part being analysed have gained popularity in recent years for beam hardening correction [9].These algorithms iterate between the measured data and a forward projection model, incorporating knowledge about the x-ray source beam spectrum and the materials being imaged.By iteratively updating attenuation parameters, these algorithms optimise the image contrast, thereby compensating for beam hardening effects and producing more accurate reconstructions.Their main drawback is the vast amount of computing time required and the lack of traceability of the method.
Advanced hardware solutions can also be an interesting but more expensive solution in mitigating beam hardening.Spectral x-ray detectors capable of measuring the energy spectrum of the x-ray beam can provide valuable information for beam hardening correction.Dual-energy or multi-energy CT systems utilise such detectors to acquire data at multiple energy levels, enabling more accurate material decomposition and attenuation coefficient estimation [10][11][12].
In this paper, we propose a new linearisation method for correcting beam hardening artefacts based on the analysis of the image noise.This method has the advantage of being 2 Association of German Engineers (Verein Deutscher Ingenieure) 3 Verband der Elektrotechnik Elektronik Informationstechnik e. V. physically traceable and does not need prior calibration with a step thickness standard.

Noise statistics
The noise in an image recorded with a scintillator x-ray detector will be dominated by different numbers of photons hitting different detector pixels, causing varying pixel brightness.The thermal background noise originating from electronic circuits is usually very small and can be subtracted with a noise offset calculated from dark images if needed.
For an idealised system consisting of an x-ray detector with uniform pixel size irradiated with monochromatic x-ray radiation of uniform intensity, the number X of incident xray photons per detector pixel is described by the Poisson distribution with the property that expectation value E (X) and variance Var (X) are equal and given by where the Poisson parameter λ corresponds to the product of irradiation intensity and exposure time.The equality of expectation value and variance of a Poisson distributed signal leads to noise behaviour commonly referred to as shot noise.
If scattering processes are ignored, photons of different energies ε are statistically independent and the above expressions remain valid for a polychromatic x-ray beam with a continuous photon spectrum λ (ε) instead of a single parameter λ.
For a scintillator x-ray detector, the statistical modelling of the signal intensity created by an incident photon spectrum λ (ε) needs to take into account two more quantities, both of which are energy-dependent: the interaction probability p sc (ε) of the photons with the scintillator layer and the detector response φ (ε) caused by a detected photon.
Using the notation λ ε := λ (ε) p sc (ε), the probability of detecting k photons of energy ε is given by Taking into account the detector response, the expectation value of the signal intensity of photons of energy ε is and similarly, the variance is obtained as The crucial conclusion of this expression is that the detector response φ (ε) breaks the equality of expectation value and variance which is the key to extracting useful information from image noise statistics.
Ignoring scattering, the expectation value and variance of the total signal intensity of the complete spectrum are given by

Example noise statistics measurements
To illustrate the influence of the spectrum on the noise statistics, a measurement of the image noise of projections of a titanium step wedge was performed using an x-ray tube voltage of 150 kV and a 0.03 mm aluminium filter.Twentyfive images of the same projection were recorded, and for every pixel, the variance of the gv across all images was calculated to obtain a variance image as shown in figure 1. Gaussian (normally distributed) background noise of the detector was calculated from a set of dark images (x-ray source turned off) and subtracted as a constant variance offset.The influence of the spatial intensity profile of the CB on the variance image was corrected by calculating the variance of a set of flat images (no absorbing object).
The relative projection gv and the variance to gv ratio are illustrated in a 2D histogram in figure 2 (with a Gaussian image filter applied to reduce the data spread).It can clearly be seen that with decreasing transmission, the variance to gv ratio increases.This can be explained by the hardening of the transmitted spectrum with increasing thickness of the Ti object.Monochromatic radiation would exhibit a constant variance to gv ratio.
At very low gvs, the variance to gv ratio appears to decrease, although the spectral hardening is expected to increase further.This discrepancy may be due to a background irradiation of scattered low-energy photons which are not related to the absorption in the Ti object.
To show that the observed change of variance is materialdependent, the same analysis has been performed for projections containing step wedges of three different materials  (titanium, aluminium and polyoxymethylene) as depicted in figure 3. The different behaviour of these materials can easily be seen qualitatively.

Correction strategy outline
Beam hardening can be corrected based on knowledge of the detected spectrum as a function of the x-ray absorption which is reflected in the detected intensity (gv).As can be seen from equations ( 1) and ( 2), the ratio Var (X) /E (X) is related to the spectrum through φ (ε) and may therefore contain some information about the spectrum: It is interesting to note that the incident spectrum λ (ε) on the detector is only revealed by the spectral sensitivity of the detector φ (ε) as a weighting function: in the case of an ideal photon counting detector with φ (ε) = const, no information about the spectrum can be extracted since the integral over the photon energy is cancelled out.Thus, measuring the ratio of variance to gv ν of one pixel can provide one single quantity of information about the incident spectrum on the detector λ det (ε) at this pixel.If this spectrum is known and using the 'air spectrum' λ air (ε) from pixels around the attenuating measuring object as a reference spectrum not subjected to beam hardening, one can compute for every pixel a correcting factor f corr for its measured gv that compensates the hardening of the detected spectrum.This pixel-specific correction factor f corr is illustrated in figure 4 and obtained from: where the numerator is a constant scaling factor and the denominator varies according to the magnitude of beam hardening.This correction factor must be applied to each pixel of all projections prior to reconstruction in order to correct for beam hardening effects.

Auxiliary quantities
Solving the single equation ( 3) can provide one quantity only as information regarding the incident spectrum on the detector.Obtaining an estimate of the spectrum from noise statistics alone therefore requires a spectrum model function with only one unknown parameter.Due to the lack of a generally applicable formula for a partially absorbed x-ray spectrum, a semi-empirical expression will be introduced in paragraph 4.1.
Additionally, solving equation ( 3) requires prior knowledge of both the interaction probability p sc (ε) of photons with the scintillator layer and the detector response φ (ε).The latter requires a careful characterisation which will be detailed in paragraphs 4.2 and 4.3.

Spectrum model function
The following semi-empirical spectrum model function is based on Kramers' law of the bremsstrahlung spectrum (1/ε behaviour) and multiplied with an exponential absorption term for low energies and a logistic tail function for modelling the behaviour close to the maximum photon energy (as illustrated in figure 5): Illustration of the correction factor as the ratio of the air spectrum integral (blue area) to the hardened spectrum integral (hatched area).Note that the derivative of both curves must be equal at the maximum photon energy in order to account for pure hardening and to avoid a faulty attenuation "correction".where b is the maximum photon energy given by the tube voltage and a defines a lower energy limit of the model spectrum and is the only unknown parameter.As an intensity model function, it must be strictly positive, thus it is defined in the interval [a, b] and 0 elsewhere.For simplicity, a linear scaling factor corresponding to the beam power has been omitted and all parameters (ε, a, b) are treated as dimensionless quantities in units of keV.Using ε/ ( ε 2 + 1 ) as the first term instead of simply 1/ε solely serves the purpose of avoiding division by zero in numerical applications and is irrelevant for the energy range of interest.
The functions f a (ε, a) and f b (ε, b) are used for tuning the low and high energy ends of the model function and may depend on the x-ray source.For the authors' system consisting of a 2 µm tungsten transmission target on a 300 µm diamond substrate and a maximum electron acceleration voltage of 190 kV, the following functions were empirically found to be convenient: This empirical match has been achieved by comparing our model to simulations of spectra attenuated by various materials and thicknesses.A qualitative illustration of the flexibility of the model function is shown in figure 6 by fitting it to simulated spectra using aRTist [13] for the source spectrum and xraylib [14] for the filter attenuation.Note that characteristic tungsten emission peaks are not included in the model and only filter materials up to the first transition metal period are shown.For calculating hardening correction factors using this model function, care must be taken at low x-ray tube voltage or very strong hardening because the derivatives of the source spectrum and the hardened spectrum may deviate at the maximum photon energy (see figure 4): a modification of the correction factor is required if the absorption term ( 1 − e −fa(ε, a) ) significantly deviates from 1 at the maximum energy.

Detection probability
The detection probability of x-ray photons in the detector scintillator is assumed to correspond to the total interaction probability of photons with the scintillator layer as shown in figure 7. The interaction cross section was calculated using the Python library xraylib [13] and is strongly dominated by the photoelectric effect in the voltage range of the authors' system.The scintillator of the authors' system is made of CsI and has a thickness of approximately 200 µm.

Detector response
The energy dependence of the response φ (ε) of a detector to photons of a given energy ε cannot be measured directly with a polychromatic x-ray source and has proven to be challenging, because any measurement corresponds to the integral over the incident spectrum.Since for a monochromatic x-ray beam, the detector response could directly be obtained from the variance to gv ratio, the strategy of choice was to apply filters with varying thicknesses to reduce the width of the spectrum for a range of source voltages and measure the variance to gv ratio.
An initial approach of extrapolating measurements obtained with increasing filter thickness to 0 gv corresponding to infinite filtering (i.e.limiting the spectrum to the source voltage) did not yield satisfactory results, since the variance to gv ratio changes rapidly towards 0 gv (shown in figure 8), rendering a numerical fit unstable, and attempts to find a suitable linearisation were unsuccessful.
Instead, a combined measurement and simulation approach was chosen: for a set of filter thicknesses, the variance to gv ratio was measured (figure 9) and the spectrum was simulated in order to obtain the mean energy of detected photons (figure 10).Both the measured variance to gv ratio and the simulated mean photon energy were fitted with the empirical model function which was chosen for its flexibility and interpolated to a set of reference relative gvs to allow consistency across different source voltages and plausibility checking.The interpolated  variance to gv ratio is then taken as an approximate measurement of the detector response φ (ε) at the interpolated mean photon energy.The relation between the x-ray tube acceleration voltage and the simulated mean energy of detected photons is shown in figure 11.

Detector response model function.
For the detector response model, two different functions were combined for a lower and higher photon energy range with an interpolated range in-between and fitted to the measured detector response  as shown in figure 12.While this decision was empirical, it can be motivated by the jump in the detection probability due to interaction with lower shell electrons in the CsI scintillator (figure 7).The model chosen for the detector response is with an interpolating transition function where a, b, c, d and the model transition energies ε 1 , ε 2 are fitted to the interpolated measurement data.The parameters obtained for our detector are (treating ε as a dimensionless quantity in units of keV again): 4.25 0.020 18.9 1.55 31.6 49.6

Consistency test.
To check the consistency of the determined models, the variance to gv ratio was simulated for a wide range of source voltages (20-180 kV) and filter thicknesses (0-4 mm stainless steel with 0.03 mm Al pre-filter) and compared with measurements.While agreement between simulation and measurement was found to be reasonable for strong filtering, the simulation strongly overestimated the variance to gv ratio for weak filtering with respect to the x-ray source voltage, as shown in figure 13, indicating an overestimation of the contribution of high energy photons in broad and intense spectra.One possible explanation for this effect might be saturation of the scintillator due to high x-ray intensity.However, changing the power of the x-ray source for otherwise constant parameters does not influence the measured variance to gv ratio as expected for a Poisson process (except for a slight decrease at very weak signal) as shown in figure 14.Saturation due to high x-ray intensity can therefore not explain the overestimation of the variance to gv ratio at weak filtering.Instead, our interpretation of this observation is that the presence of low energy photons causes a 'spectral saturation'  of the scintillator and alters its response behaviour to high energy photons, leading to a smoothing of the detected signal strength.
To model this effect, a 'saturated spectrum' was introduced as ) with an empirical saturation factor f sat and weighting function Assuming this spectral saturation only affects the variance but not the gv, the modified simulated variance to gv ratio ν sat is then given by Reasonable agreement between simulation and measurement was found for the authors' system for f sat = 0.2 and ε sat = 140 keV, yielding the results shown in figure 15.

Application
Calculating the pixel-resolved variance of an x-ray image requires recording the same image multiple times with the resulting variance image being noisy, unless a large number of images is recorded.Repeating this procedure for every projection of a CT scan would be possible but extremely timeconsuming and is therefore not suitable for applications.For CT scans of objects consisting of only one material or a statistically homogeneous material mix (e.g.powder mixture) however, the relation between variance and gv is unique and can be calculated from a set of images of a single projection and fitted with a model function.An empirical model function which proved to be suitable for this task is where a, b, c, d, e are fit parameters and x is the relative gv.In order to improve the stability of the fit against bad data quality at very low gv and flat field correction artefacts at very high gv, the data points for fitting were limited to a range of typically 15%-95% relative gv, as indicated in figure 16.
It should be noted that the orientation of the object used for measuring and fitting the variance to gv ratio model does not need to correspond to an actual projection of the CT scan.However, it is important for the fit to be reliable that a sample orientation is chosen whose projection includes a broad range of thicknesses (i.e.gvs) and especially thin parts, since the behaviour of the model close to a relative gv of 1 is crucial for the result of the correction.This is due to the fact that beam hardening effect is strongest in the outermost layer of material as the x-ray beam penetrates the object.Also, detector afterglow can have a significant influence on the image variance (especially in dark areas) and should be avoided by letting the detector stabilise for ideally a few minutes before recording the fit model projections.
In summary, this model fit method using multiple projections of a single orientation can be applied to any object of a statistically homogeneous material composition with a continuous thickness distribution in at least one orientation.Apart from this requirement, the usability of the method is not restricted to specific applications but can reduce surface determination errors and improve contrast in internal structures for any type of geometrical analysis.
The computation time required for calculating a model as shown in figure 16 is less than a minute on a normal computer.The correction of an individual projection takes several seconds on one CPU (implemented in Python) but can be run in parallel and performed during the measurement.The total correction computation time for a scan of 3000 images with 4096 × 4096 pixels was around 30 min on a workstation with 20 CPU cores and limited by simultaneous read and write operations on the same solid state drive, which could easily be improved by using different drives for reading and writing.

Example results
As an example to demonstrate the correction method, a machined part of an unknown brass alloy with an outer diameter of 9 mm was scanned with an x-ray source voltage of 160 kV and a 0.1 mm Cu filter, representing a CT scan of unfavourable parameter choices.From a set of 16 images of a single sample orientation, the variance to gv ratio was determined and fitted as shown in figure 16.Equation ( 4) was then numerically solved for the spectrum model parameter a for a range of several thousand gvs, allowing the calculation of a lookup table for the correction factor f corr shown in figure 17 (normalised to 1 at relative gv 1).This correction factor was then applied to all projections of the CT scan prior to reconstruction as shown in figure 18.A qualitative comparison of the reconstructed object with and without correction is presented in figure 19.
Note that even though the model assumes the object to be composed of a single material, a structure resembling grains can be seen in the sections and line profiles which may arise from inclusions of heavy elements such as Sn or Pb which are sometimes added to brass alloys for improved machinability.Knowledge of such alloy elements is not required for applying the correction method.

Limitations
The method described in this paper could be applied to a general case of multi-material objects but would require measuring the noise statistics for every pixel at every projection angle which would require excessive measurement time.Nevertheless, smart denoising algorithms might reduce the number of images required for obtaining usable results.
Reliable application of the method to single-material objects requires object geometries which can be positioned in a way that a projection image with a wide range of gvs is obtained, resulting in a suitable distribution of data points for fitting the material-specific variance to gv ratio (figure 16).Geometries such as spheres, cylinders or cubes with large relative wall thickness are therefore not well-suited.
The method does not take into account scattering of x-ray photons which enhances the number of photons at low energies detected as an additional background illumination.This effect has not been investigated but is not expected to have a relevant influence on the results.
In addition, the method only works for suitable combinations of x-ray source voltage and sample material where the absorption cross section is continuous across the x-ray spectrum, leading to continuous hardening of the spectrum towards higher energies.While the method has worked well for elements up to the first transition metal period, its application to objects commonly composed of heavy elements such as jewellery or tungsten carbide tools is expected to be less accurate and has not been investigated.

Conclusion
We present a new method for beam hardening correction in xray CT by using image noise statistics as a measurable quantity to obtain an estimate of the incident x-ray photon spectrum on a scintillator detector.The key element for the calculation of the correction is a versatile spectrum model function with only one unknown parameter which can numerically be solved for by simulating the corresponding noise statistics.The method does not require knowledge about the material composition or dimensions of the object being scanned since the beam hardening is measured from the image noise statistics of projections of the object directly.
While the measurement and calculation of the beam hardening correction for a specific object can be performed quickly, the method requires a carful characterisation of the scintillator detector of the CT system which is complex and relies on several measured quantities and assumptions.Although a sophisticated setup with a tuneable x-ray spectrum could enable a more precise characterisation of the detector properties, the approach presented in this paper has provided satisfactory results without the requirement for additional hardware except for a stack of sheet metal as a variable thickness filter.
The obtained beam hardening correction gives qualitatively very good results attested by the flattening of the background noise level and cupping effect.The edge contrast of the part is greatly enhanced avoiding fundamentally erroneous or impossible surface determination.Contrast improvement is most pronounced for internal structures such as highlighted in figure 19, which is additionally expected to improve the sensitivity for porosity or inclusion detection.The magnitude of beam hardening correction is traceable to the physical property of noise.Nevertheless, the quantitative influence of the correction method on dimensional measurements of reconstructed objects still needs to be systematically assessed.
A strategy for the application of the method to x-ray CT scans of objects of general combinations of multiple different materials has not been developed.For 2D x-ray imaging however, analysing the variance to gv ratio from multiple images may be beneficial for material segmentation by providing an extra parameter in addition to the gv alone.
Aside from being used for beam hardening correction, detailed understanding and modelling of attenuationdependent noise statistics may also improve CT simulations by making image noise more realistic, especially in highly attenuating areas.

Outlook
Next steps will involve comparing the accuracy of the correction method both with a calibrated sample and with other approaches, in particular with linearisation using calibration measurements on a thickness standard.Regarding multimaterial applications, tests are intended to find out under what circumstances the model fit method might still be applicable, for instance for combinations of a strongly and a weakly absorbing material like a part composed of one metal and synthetics.The general case of application to any multi-material combination using multiple recorded images of every projection of a CT scan is not a priority, since the required scan time for useful results is expected to be unsuitable for practical applications.

Figure 1 .
Figure 1.One projection of a Ti step wedge (left) and a detail of the variance image calculated from 25 projections (right).

Figure 2 .
Figure 2. Relation between the variance to grey value (gv) ratio and the gv in projections of a Ti step wedge.

Figure 3 .
Figure 3. Relation between the variance to gv ratio and the gv in projections containing step wedges of three different materials.

Figure 5 .
Figure 5. Composition of the spectrum model function.

Figure 6 .
Figure 6.Comparison of simulated and model spectra for 90 kV x-ray tube voltage without filtering (left), 120 kV with 1 mm Al filter (centre) and 150 kV with 4 mm Cu filter (right).

Figure 8 .
Figure 8. Simulated variance to gv ratio for 150 kV source voltage and different materials.

Figure 9 .
Figure 9. Measured variance to gv ratio with a fitted interpolation function for 150 kV source voltage, 0.03 mm Al pre-filter and stainless steel filters of variable thicknesses.

Figure 10 .
Figure 10.Simulated mean energy of detected photons for 150 kV source voltage, 0.03 mm Al pre-filter and stainless steel filters of variable thicknesses.

Figure 11 .
Figure 11.Relation between the x-ray tube acceleration voltage and the simulated mean energy of detected photons filtered with stainless steel.

Figure 12 .
Figure 12.Data points obtained as measurements of the detector response and fit of an empirical model function.

Figure 13 .
Figure 13.Comparison of measurements (squares) and simulation (dotted lines) of the variance to gv ratio for different voltages and stainless steel filter thicknesses, showing large discrepancies at low filtering and high voltages.

Figure 14 .
Figure 14.Measurement of the variance to gv ratio with 120 kV source voltage, 0.03 mm Al filter and 3 s exposure for varying source power.

Figure 15 .
Figure 15.Simulation (dotted lines) of the variance to gv ratio taking into account a model for spectral scintillator saturation, showing better agreement with measurements (squares) than figure 13.

Figure 16 .
Figure 16.2D histogram and model fit of the variance to gv ratio vs. relative gv of a projection of a brass object.

Figure 17 .
Figure 17.Relation between the measured and corrected gv obtained from the fit model of the variance to gv ratio.For relative grey values ⩾1 (noise in the air background), the correction factor was extrapolated with constant derivative.

Figure 18 .
Figure 18.Projection of the example object as recorded (left) and with the gv correction applied (right).

Figure 19 .
Figure 19.Comparison of uncorrected (top) and corrected (bottom) reconstruction volume: 3D view with surface determination (left), horizontal section (centre) and line profiles (right, plotted with arbitrary offsets).Surface determination was performed using a gv threshold as starting point with gradient-based local optimisation.