1 Introduction

Computed tomography (CT) is an imaging technique that has been widely applied in medical and industrial applications such as for medical diagnosis and non-destructive testing. However, CT usually suffers from undesirable artifacts such as scattering, and beam hardening (BH) due to the physical nature of the X-rays [1]. One of the major artifacts that influence image quality is the BH artifact, which results from the polychromatic nature of the X-ray source. When photons of different energies penetrate through the object, low-energy photons are more easily absorbed than high-energy photons. This shifts the mean of the spectrum to a higher value and results in the hardening of the beam. Since the attenuation value of the material is energy-dependent, it decreases as the beam is hardened. This results in a non-linearity between the attenuation of the beam and the length of the propagation path. Using a linear reconstruction method such as the filtered back projection (FBP), which assumes that the object attenuates the X-ray linearly with the path length, results in cupping and streak artifacts which degrade the quality of the reconstructed image.

Several methods have been implemented to correct the BH artifact, such of these methods are filtration, dual-energy, linearization, post-reconstruction, statistical polychromatic reconstruction, and iterative reconstruction methods [2, 3]. Although the filtration method is simple and widely applied, it is only able to reduce the BH artifacts and does not perform well in multi-material applications [4, 5]. The dual-energy method requires two scans, this is in general expensive. In addition, it is limited to certain applications [2]. The linearization and the post-reconstruction methods require the material to be known which is not always the case. Although some linearization methods do not require material knowledge [6,7,8], they have some limitations. These methods are either limited to a certain class of objects or a specific geometry [2]. Moreover, the main burden of the statistical and iterative reconstruction methods is the long computation time required, which makes them not suitable for industrial applications [9].

In this work, a multi-material BH correction algorithm is proposed. This method does not require prior knowledge of the materials, which is denoted blind, and can correct first-order (cupping) and high-order (streak) BH artifacts by accurately mapping the polychromatic projection to an equivalent monochromatic one through a correction term that is added to the projections from the scanner. The simulation of the projections, which are needed for the calculation of the correction term, is performed using two implemented fast multi-GPUaccelerated forward projection models. As a consequence, the BH correction can be performed in near real-time using the proposed method on four GPUs. Even a real-time BH correction can be achieved by the use of more GPUs, i.e, a shorter processing time is required in comparison to the acquisition time of the complete number of projections from the real-world scanner [13]. The first forward projection model used is a fast and sufficiently accurate GPU-based Monte Carlo (MC) photon transport model (FPM) which can simulate the scatter projection accurately in addition to the primary projection, this model is introduced in [10]. However, the forward projections that are involved in the proposed BH correction method are only primary projections without scatter. The second forward projection model (FPLS) is based on the line-simple method [11]. The FPLS model is only used to provide a reference for the proposed material estimation method. The experimental projections are first corrected from the scatter artifacts prior to the correction of the BH artifacts. The scatter estimation and correction were performed using the GPU-based MC model (FPM) and an iterative scatter correction algorithm in near real-time following the work in [10]. The scatter correction results from this method match the experimental reference near scatter-free results from the collimator. After the correction of the scatter artifacts, the scatter-corrected projections are imported to the BH correction method and used by the implemented materials estimation method, the polychromatic and monochromatic projections estimation methods.

The motivation of this work is to introduce a fast BH correction algorithm that can correct higher-order BH artifacts in near real-time and at the same time does not require the materials’ knowledge. This algorithm suits the industrial applications as in industry the materials are not always known, especially for compound materials, and fast artifact correction is necessary. Existing BH correction methods which do not require the materials’ knowledge either work iteratively, which requires a long execution time, or their applications are limited to a specific geometry or objects. In the proposed algorithm, novel methods to estimate the materials and the BH-corrupted projection from the scanner are implemented. The materials are estimated using a look-up table, while the polychromatic projection is derived using least square estimation (LSE) fitting from eight simulated polychromatic projections using spectra of different filtration. The introduction of these methods facilitates the derivation of an accurate non-linearity correction term that can accurately correct higher-order BH artifacts.

Our approach represents a time-efficient alternative in near real-time to the iterative BH correction algorithms in industrial computed tomography. Moreover, the proposed work is more efficient than the linearization methods, as these methods require the materials’ knowledge and they follow an iterative approach in the case of multi-material objects. In comparison to the post-reconstruction methods, the proposed BH correction method, supported by the LSE fitting approach, has the potential to derive an accurate BH correction term. As the proposed fitting approach overcomes the possible mismatch between the simulated polychromatic projection using any forward projection simulator and the BH-corrupted projection from the scanner. Such a mismatch is possible due to the inaccuracy of the attenuation tables and the spectrum used in such simulators. Finally, the statistical approaches are very slow and also require additional system information. Using the proposed BH correction method, the BH artifacts were corrected in near real-time and the correction results have been evaluated by several experimental examples from real-world data. It is shown that this algorithm can accurately correct the BH artifacts without prior knowledge of the materials.

This paper is organized as follows. The rest of Sect. 1 presents the contribution of this work and the related works. Section 1.1 describes the BH problem, the proposed methods to estimate the materials, the BH-corrupted projection from the scanner, and the monochromatic projection. Section 3 shows the experimental results on real-world data. In the same section, discussions about the robustness of the proposed methods to estimate the material and the BH-corrupted projection are also provided. Finally, the discussion and the conclusion are presented in Sects. 4 and 5, respectively.

1.1 Contributions

1.1.1 Fast and robust material estimation based on look-up table

To overcome the limitation of many other BH correction methods, which require the knowledge of the materials used, we have created a look-up table that contains experimentally measured polychromatic linear attenuation values of a wide range of well-known materials versus the polychromatic energy. Based on the look-up table and the measured linear attenuation value of each material in the uncorrected volume, the proposed algorithm can iteratively find the best-matched material by means of minimum mean square error (MMSE). Experimental results show that the proposed material estimation is robust to strong BH effects and a small amount of deviation of the estimated material has a negligible impact on the BH-corrected results.

1.1.2 Accurate regression of the measured projections by aggregation of weighted MC simulations

To correct the non-linearity of the acquired projections, a correction term is required. One determinant part of the correction term is the estimated polychromatic projection. Estimation of such a projection using forward projection simulators leads to an inaccurate result since the linear attenuation tables of the materials used in such simulators are not very precise, the same is also applied to the polychromatic spectrum used. To enable an accurate fitting to the measured data, the estimated projection is constructed by a linear combination of eight simulated polychromatic projections of different spectra using diverse filtration. It is shown that the constructed projection very well matches the raw data.

1.2 Related work

In the literature, the methods of BH correction can be classified into six categories which are based on hardware filtration, dual-energy, linearization, post-reconstruction, statistical polychromatic reconstruction, and iterative reconstruction [2, 3]. Some other methods combined two classes to produce an effective correction method.

The first approach implies the use of filters [12, 14] in which a thin plate of metal, aluminum, copper, or other material is placed between the source and the object to filter out the low-energy photons in the spectrum before their penetration through the object. This method has been widely used in the CT field due to its simplicity and the limited requirement of this method. However, it is not completely able to eliminate this artifact, especially in the case of multi-material objects with different high attenuation characteristics [4, 5]. Besides, the use of these kinds of filters results in a lower detected signal-to-noise ratio and reduces the contrast between the materials [15].

In the dual-energy methods [4, 16], two images from the low and high energies are acquired. Then, a material decomposition is performed by representing the linear attenuation coefficient by two bases functions, i.e., Compton scattering and photoelectric absorption [6, 16, 17]. The contributions from these two bases functions are calculated which enables the reconstruction of a monochromatic image. Similar to the approach in [16], the authors in [18] introduced a method to decompose the linear attenuation value of a material \(\mu\) into density and atomic number basis. Then the determination of both at every point in the sample is achieved by the use of two energies. The coefficients are first estimated by the use of a known materials phantom, then this method can be used with a phantom of unknown materials. The main flaw of this method is the need to acquire two separate distinct measurements which is in general expensive.

In the linearization methods [4, 6, 19, 20], the non-linear relationship of the experimental polychromatic attenuation values and the object thickness is fitted with a polynomial. In which they are mapped to a linear line to eliminate the BH artifacts. For small BH artifacts, the fit of a second-order polynomial is enough [6, 20]. However, for large BH artifacts, the fit of a higher-order polynomial is required [6]. For multi-material objects, the first reconstructed data is used to provide prior information about the intersection of the materials in the object with each ray path. The linearization can be then performed using the iterative post reconstruction (IPR) approach. The main limitation of the linearization methods is that the materials should be known, only a few methods have been proposed that do not require the materials’ knowledge [6,7,8]. The method proposed in [6] assumes that the tissues have energy dependence similar to that of the water, such an assumption performs badly in the presence of high-Z materials [34]. On the other hand, the authors in [7] proposed a method that is limited to objects made from cylinders with a high-Z material outside and a low-Z material inside. In [3], the authors proposed a BH correction method that is iterative in nature and similar to the IPR methods. The proposed method does not require the materials’ knowledge nor the spectrum and uses a hypersurface and a hyperplane to represent the polychromatic and the monochromatic attenuation values plotted with respect to the ray path through each material. The difference between the hypersurface and the hyperplane is calculated which is then added to the original projection from the scanner to correct the BH artifact. This method works well in case a proper segmentation can be performed on the original uncorrected volume [2].

In the post-reconstruction methods, monochromatic and polychromatic projections are simulated relaying on the reconstructed image, the spectrum used, and material characteristics such as the attenuation coefficient, mass density, and atomic number. The difference between the two projections is added to the original corrupted projection to eliminate the BH artifacts [21, 22]. In [23], the authors have proposed a method to correct the CT artifacts including the scatter and the beam hardening artifacts. The BH is corrected by deriving a non-linearity correction term with the assumption that the materials are known. The scatter is estimated using a hybrid approach, i.e., using MC simulation and an analytic scatter convolution model. The proposed BH correction method in this work belongs to this category of BH correction methods, however, in contrast to these methods, our approach does not require material knowledge, and the scatter artifact is estimated using full MC simulation utilizing the fast multi-GPU model in [10]. The scatter artifact is then corrected using an iterative scatter correction method [10]. Moreover, we propose the use of the LSE method to derive an estimation of the polychromatic projection which has the potential to provide an accurate estimate of the experimental polychromatic projection.

Statistical correction methods [24,25,26,27,28,29,30, 33] exploit multiple information to correct the BH artifacts in the X-ray images. Such information includes the polychromatic nature of the source, the detector model, the noise distribution, the measurement non-linearity, and the scatter effect are all used in the statistical methods in which they are all incorporated into the maximum likelihood (ML) algorithm [2]. In [24], an objective function is formulated from the Poisson likelihood, then a quadratic approximation to this likelihood is derived, this quadratic approximation leads to a penalized weighted least square (PWLS) estimate, and as a consequence, a simpler objective function is derived. By minimizing this simplified objective function, an updating step for this algorithm is obtained which can estimate the unknown densities in each voxel with reduced beam hardening artifacts. The same authors proposed an algorithm in [25] that follows the same approach of [24]. In addition to the formulation of the penalized-likelihood function, as in the former method, they developed an ordered-subsets iterative algorithm for the estimation of the unknown densities in each voxel. The methods in [24, 25] assume that the object contains non-overlapping materials and the segmentation is possible. These works are then extended further in [26] to allow for overlapping materials in the same voxel. A statistical method to reduce the BH artifact has been introduced in [34]. This method uses the same physical model proposed in [33] and they use a reconstruction algorithm that ensures a fast monotonic convergence using the non-linear conjugate gradient method. The major limitation of the statistical approaches is the long computational time required.

On the other hand, iterative algorithms improve the images iteratively by minimizing the cost function formed from the difference between the measured polychromatic images from the scanner and the simulated one. The spectrum, which is already involved in the minimization method along with the materials attenuation values, is approximated by a small number of energy bins [2]. In [31] the authors proposed a new reconstruction algorithm that corrects the BH artifacts by involving the polychromatic characteristic of the X-ray beam. The algorithm in [32] corrects the BH effect by considering a polychromatic acquisition model and taking into account the energy dependency of the materials. In this method, the linear attenuation value of the material is decomposed into photoelectric and a Compton scatter components which are weighted based on a prior assumption. A model-based BH correction algorithm has been introduced in [35]. This algorithm neither requires the knowledge of the X-ray spectrum nor the attenuation values of the materials. A polynomial function of both low- and high-density materials is first formulated which represents a poly-energetic X-ray forward projection model. Then, an alternating minimization algorithm is then developed which is used to estimate the reconstructed image, the material segmentation, and the two coefficients of the polynomial function which models the BH function.

Recently, the energy-sensitive detector, also known as a photon-counting detector, is introduced. These kinds of detectors provide multiple information about each detected photon including the energy of this photon. In [36], it is shown that simple energy thresholding applied in such a detector allows the reduction of the BH effect and enhances the quality of the CT images.

2 Materials and methods

2.1 Description of the beam hardening problem

Multiple material specifications influence the attenuation of the source intensity \(I_0\). Such specifications include the material atomic number, the density of the material, the energy of the photon, and the source–detector propagation path. According to the Beer–Lambert law, the intensity on the detector from a monochromatic source is given by Eq. 1 [22],

$$\begin{aligned} I_{mono}=N\eta (E) exp\Big (-\int ^L_0\mu (E,{\hat{x}} ,\rho )d{{\hat{x}}}\Big ),\ \end{aligned}$$
(1)

N is the number of incident photons, \(\eta (E)\) is the energy response of the detector at energy E, \(\mu (E,{\hat{x}},\rho )\), where \({\hat{x}}\) is spatial position, and \(\rho\) is the density of the material, is the energy-dependent linear attenuation coefficient, and L is the propagation path length. For a monochromatic source, the intensity of the rays has the same amount of attenuation for every slice of thickness dx as in Eq. 2,

$$\begin{aligned} dI/I=-\mu dx. \end{aligned}$$
(2)

This is because there is no change in the energy in such a source and the attenuation value \(\mu (E,{\hat{x}},\rho )\), which is energy dependent, does not change. However, the source used by the real-world scanners is a polychromatic source in which the energy of the photons emitted by such a source varies over several energies. When the rays from such a source pass through an object, the ratio of the high to the low energy of the photons increases. This is because the low-energy photons easily got absorbed during their penetration and the mean of the beam is then shifted. The attenuation of the materials inside the object, as its energy-dependent, will have then a lower value due to the increase in the mean energy. Thus, the expression of the intensity from Eq. 1 is no longer valid. The intensity scored on the detector from such sources can be expressed by Eq. 3 [22],

$$\begin{aligned} I_{poly}=N\int ^{E_m}_0\eta (E)S(E)dE\,exp\Big (-\int ^L_0\mu (E,{\hat{x}},\rho )d{\hat{x}}\Big ). \end{aligned}$$
(3)

Here, \(E_m\) and S(E) are the maximum energy and the spectrum of the source, respectively. In the absence of the object, the monochromatic and the polychromatic intensities are given by Eqs. 4 and 5, respectively,

$$\begin{aligned} I_{mono,0}&=N\eta (E). \end{aligned}$$
(4)
$$\begin{aligned} I_{poly,0}&=N\int ^{Em}_0\eta (E)S(E)dE. \end{aligned}$$
(5)

The non-linear relation between the intensity and the distance described in Eq. 3, causes the BH artifacts. Such a non-linear relation could be compensated by finding a proper correction term that can map this non-linear relation to a linear one. A certain class of correction methods derives a correction term by calculating the difference between estimated monochromatic and polychromatic projections. This term is then added to the original projections to correct the BH artifacts. If Eqs. 6 and 7 represent the logarithm of the monochromatic and the polychromatic projections, i.e., \(\chi _{mono}\) and \(\chi _{poly}\), respectively,

$$\begin{aligned} \chi _{mono}&=ln\Big (\frac{I_{mono,0}}{I_{mono}}\Big ), \end{aligned}$$
(6)
$$\begin{aligned} \chi _{poly}&=ln\Big (\frac{I_{poly,0}}{I_{poly}}\Big ). \end{aligned}$$
(7)

The correction term is then calculated as in Eq. 8.

$$\begin{aligned} \Phi _{CT}=\chi _{mono}-\chi _{poly}, \end{aligned}$$
(8)

where \(\Phi _{CT}\) is the correction term. The correction of the BH-corrupted projections is then performed according to Eq. 9.

$$\begin{aligned} \chi _c=\chi _u+\Phi _{CT}, \end{aligned}$$
(9)

where \(\chi _c\) is the corrected attenuation projection, and \(\chi _u\) is the BH-corrupted attenuation projection.

Fig. 1
figure 1

The flowchart of the proposed BH correction algorithm

Full size image

2.2 The proposed BH correction method

The flowchart of the proposed BH correction algorithm is shown in Fig. 1. This algorithm is divided into four parts. The first part is to estimate the materials of the BH-corrupted volume. The estimated materials are then used in the following two parts of the algorithm to estimate the polychromatic projection described in Sect. 2.2.2 and the monochromatic projection in Sect. 2.2.3, based on which we can yield the BH-corrected volume as formulated in Eqs. 8 and 9.

2.2.1 Materials estimation

Most of the available BH correction methods require prior knowledge of the materials. However, such information in many practical cases is not available [2]. The existing methods which do not require material information either follow an iterative correction approach, which is computationally expensive [9], or perform a linearization with some restrictions [2], which limits the performance of BH correction. To circumvent these issues, the proposed multi-material blind BH correction method first performs a simple but effective estimation of the materials and based on which, we can obtain the BH correction term using the GPU-accelerated forward projection model FPM. Before we explain the method in detail, we first present a look-up table as shown in Table 1, containing polychromatic linear attenuation values of a wide range of well-known materials such as copper, and aluminum are measured experimentally for different energies. These linear attenuation values are derived from cylindrical objects with fixed diameters and different materials. The look-up table is measured using our in-house Nikon XT H 225 CT scanner equipped with the PaxScan 4030E detector from Varian. A collimator was used to suppress the scatter artifact in all the experiments. The experimental setups used to acquire the attenuation values of the table were with a current of 150 \(\mu\)A, a tungsten source target, and 1 mm Cu filtration. Since CT systems with different sources, materials of source target, beam spectra, and filtration would give a different response, a new look-up table should be constructed for individual CT systems.

Table 1 A Look-up table of the materials attenuation values \((m^{-1})\) versus the energy in keV
Full size table

The material estimation starts from the scatter-corrected raw data as shown in the upper part of Fig. 1. A volume is reconstructed from these scatter-free projections using the FBP reconstruction method. Different materials in the reconstructed volume can be separated by the segmentation method Otsu [37]. These materials are assigned individually into volumes of the same size as the original input, i.e., each volume contains a single material. Since the reconstructed volumes are almost scatter-free, they are generally of adequate quality to perform the segmentation by Otsu. For each volume, a single attenuation value is initialized by comparing with the linear attenuation values of all the collected materials in the look-up table by Eq. 10.

$$\begin{aligned} mse=\mathop {{\mathbb {E}}}\Big [\parallel \mu -{\hat{\mu }}\parallel _2^2\Big ], \end{aligned}$$
(10)

where \(\mu\) represents the linear attenuation value of the material in the individual volume, and \({\hat{\mu }}\) represents the linear attenuation values of all the materials in the table for the selected energy. The best fit according to the MMSE is then selected as an initial guess of the material for each volume. To iteratively refine the estimation, an L2-norm minimization between the forward projection of the volume with assigned material using FPM [10] and the forward projection of the reconstructed volume by FPLS [11] is carried out. The first projection \((\chi )\) is derived by the Monte Carlo based FPM which is simulated using the primary photons only and based on the initial guess of the material. The linear attenuation values of the estimated material are converted to the density since FPM supports only density-based simulation according to Eqs. 3, 5, and 7 [38]. The second projection \(({\hat{\chi }})\) is derived using the volume which has the separated material. This projection is acquired using the forward projection model FPLS by simulating rays from the source to each pixel of the detector. The pixel value of the detector is approximated by summing the voxels’ linear attenuation values at equidistant points along the ray as formulated in Eq. 11.

$$\begin{aligned} {{\hat{\chi }}}_{j_{l}}=\sum _{v=1} \mu _{v} d_v, \end{aligned}$$
(11)

where \({{\hat{\chi }}}_{j_{l}}\) is the value recorded in the \(l_{th}\) pixel of the detector from the volume which contains the separated material j, \(\mu _{v}\) is the linear attenuation value of the voxel v, and \(d_v\) is the distance that the ray travels in voxel v. The simulation of \(({\hat{\chi }})\) is performed without converting the linear attenuation values into densities. The projection \(({\hat{\chi }})\) is closer to reality as it is unlike the projection \((\chi )\), it neither relies on the tables of the attenuation values nor on the spectrum. Thus, it overcomes the possible mismatch between the simulated projection \((\chi )\) and the experimental projection \((\chi _u)\).

The method used to calculate the projection \(({\hat{\chi }})\) can provide an exact match with the BH-corrupted projection of the scanner (the scatter-corrected projection) if the volume used in the simulation is the reconstructed volume that contains all the materials. Therefore, the projection of FPLS is used as the reference to determine the accuracy of the estimated material. Although FPLS has the potential to generate an exact match with the BH-corrupted projection of the scanner, it cannot be used to generate the polychromatic projection for the correction term in Eq. 8. This is to ensure that the algorithm does not end up simply replacing the measured projections with the monochromatic projections, which are highly dependent on the accuracy of the segmentation. Thus, we propose the material estimation and the polychromatic regression to get rid of the influence of inaccurate segmentation. To determine the optimal estimation of the material, the simulation and minimization are iteratively performed until a best fit is found. The same process is repeated to estimate the rest of the materials. Section 3.2.2 discusses the robustness of the material estimation and the effect of the incorrect material estimation on the result of the BH correction.

Fig. 2
figure 2

Different spectra simulated using different filtration

Full size image

2.2.2 Simulation of the polychromatic projections

As the materials in the BH-corrupted volume are determined using the proposed materials estimation method, the linear attenuation values of the segmented volume are replaced by the density value of the materials. To derive a correction term for the BH-corrupted projections, an accurate estimation of the simulated polychromatic projection is required. The use of forward projection simulators can provide inaccurate estimations which do not match the BH-corrupted projection. This is mainly because the linear attenuation tables of the materials and the polychromatic spectrum used in such simulators are not very precise. In addition, some effects in the actual scanners are ignored in these simulators. Thus, there is always a noticeable mismatch between these simulated polychromatic projections and the experimental ones.

To construct an accurate match with the raw projection, we propose to perform a regression from eight polychromatic projections simulated by FPM as presented in the middle part of Fig. 1. These eight projections are simulated using primary photons only and with different filtration under spectra of the same energy. Figure 2 shows the spectra used to acquire the different projections. These spectra are all simulated using SpekCalc software [39]. Then, the LSE algorithm is used to minimize the difference between the BH-corrupted projection and the weighted sum of the eight simulated projections. This is achieved according to Eqs. 12 and 13.

$$\begin{aligned} c=arg\; min_{c} \mathop \parallel \chi _u-\chi _{p}\parallel _2^2, \end{aligned}$$
(12)

where c are the coefficients to be estimated, \(\chi _u\) is the BH-corrupted attenuation projection, and \(\chi _{p}\) is the weighted sum given by Eq. 13,

$$\begin{aligned} \chi _{p}=\sum _{i=1}^{k} c_{i}\chi _{p_{i}}, \end{aligned}$$
(13)

where k denotes the amount of the basis polychromatic projections for the regression and \(\chi _{p_{i}}\) is the individual polychromatic projection simulated by FPM. From Eqs. 12 and 13 we get:

$$\begin{aligned}\mathop\parallel \chi_u-\chi_{p}\parallel_2^2&=\mathop\parallel\chi_p\parallel_2^2-2(\chi_p,\chi_u)+\mathop\parallel\chi_u\parallel_2^2,\\ &=\sum_{i=1}^{k} \sum_{j=1}^{k} c_{i}(\chi_{p_{i}},\chi_{p_{j}})c_{j}-2\sum_{i=1}^{k}(\chi_{p_{i}},\chi_u)c_i \\ & \qquad +\mathop\parallel\chi_u\parallel_2^2,\\ &=c^TG_{\chi_{p}}c-2r_{\chi_{p}\chi_{u}}^Tc+\mathop\parallel\chi_u\parallel_2^2,\end{aligned}$$
(14)

where the cross-correlation vector \(r_{\chi _{p}\chi _{u}}\) and the matrix of inner product \(G_{\chi _{p}}\) are given by;

$$\begin{aligned}{} & {} r_{\chi _{p}\chi _{u}}= \begin{vmatrix} (\chi _u,\chi _{p_{1}})\\ (\chi _u,\chi _{p_{2}})\\ \vdots \\ (\chi _u,\chi _{p_{k}}) \end{vmatrix},\\{} & {} G_{\chi _{p}}= \begin{vmatrix} (\chi _{p_{1}},\chi _{p_{1}})&(\chi _{p_{1}},\chi _{p_{2}})&\cdots&(\chi _{p_{1}},\chi _{p_{k}})\\ (\chi _{p_{2}},\chi _{p_{1}})&(\chi _{p_{2}},\chi _{p_{2}})&\cdots&(\chi _{p_{2}},\chi _{p_{k}})\\ \vdots&\vdots&\vdots&\vdots \\ (\chi _{p_{k}},\chi _{p_{1}})&(\chi _{p_{k}},\chi _{p_{2}})&\cdots&(\chi _{p_{k}},\chi _{p_{k}}) \end{vmatrix}.\end{aligned}$$

To obtain the coefficients c, the gradient of the distance function with respect to c in Eq. 14 is calculated and we get the solution in the following equation:

$$\begin{aligned} c=G_{\chi _{p}}^{-1}r_{\chi _{p}\chi _{u}}. \end{aligned}$$
(15)

The optimized coefficients are used to weight the eight simulated projections and construct the best estimation of the acquired projections. Further discussion about the robustness of the polychromatic projection estimation can be found in Sect. 3.2.3.

2.2.3 Simulation of the monochromatic projections

To complete the correction term, an estimation of the monochromatic projection is required as formulated in Eq. 8. The energy of the monochromatic projections is determined by minimizing the difference between the simulated projections using FPM following Eqs. 1, 4, and 6 and the BH-corrupted projection from the scanner in terms of MSE. The energy that results in the lowest MSE is selected. The simulated monochromatic and polychromatic projections are then used to correct the BH artifacts by Eq. 9. The monochromatic estimation method and the correction part can be seen in the lower part of Fig. 1. It should be noted that the aforementioned correction procedure can be repeated iteratively to achieve superior performance.

3 Experiments and results

In this section, comparisons between the results of the BH correction from the proposed method and the state-of-the-art empirical BH correction method [9, 40] are first presented using real-world datasets. These datasets were acquired using an in-house Nikon XT H 225 scanner for different test objects. Moreover, analyses of the BH correction results were also performed.

Fig. 3
figure 3

Experimental results of the proposed BH correction method. The first column shows the slices from the scatter and BH-corrupted volumes. The second column shows the slices from the scatter-corrected volumes. The third and fourth columns show the slices from the BH-corrected volumes using the proposed method and the empirical method, respectively. The last column shows the profile lines of the slices in columns one–four. The profiles are marked by white dashed lines in these sub-figures. The profile lines were averaged over multiple lines to suppress the noise. The red squares represent the ROI and the background used in the calculation of the contrast-to-noise ratio (CNR)

Full size image

3.1 Evaluation of the proposed BH correction method

To evaluate the proposed BH correction method, a comparison has been conducted between the BH correction results of the proposed method and the BH correction results using the state-of-the-art empirical BH correction method. The empirical method is very robust and one of the widely used methods as it does not require prior knowledge of the materials nor the spectrum [22]. The parameters of the empirical method were derived according to [9, 41]. The BH correction results of six experimental objects using the two methods and the scatter correction results are shown in Fig. 3. The measurement parameters used in these examples are shown in Table 2. All the experimental datasets were acquired using a detector resolution of \(2304\times 3200\) pixels with a pixel size of 0.127 mm. Before the application of the BH correction method, these datasets were downsampled to \(576\times 800\) pixels resolution. The first object in this figure is a water pipe distributor of 7 cm in length. The material of this object is unknown. The second object is an electrical wire connector strip of 11.5 cm\(\times\) 2 cm \(\times\) 1.5 cm. The exact materials of this object are also unknown. The third object is a motorcycle cylinder head made from aluminum and iron. The dimension of this object is 12 cm\(\times\)11 cm\(\times\) 5.5 cm. In the fourth row, the object is a cement cylinder with an 11 cm outer diameter and a 5 cm inner diameter with eight iron rods inserted in this cylinder. The object in the fifth row is a 7 cm diameter cement cylinder with eight iron rods. The last object represents three aluminum cylinders positioned away from each other in space. The biggest cylinder has a 2 cm diameter, while the two small cylinders have a 1 cm diameter each. The volumes from the scanner, shown in the first column of this figure, were reconstructed using 1500 projections for the object in the first row, 2000 projections for the object in the second row, and 3000 projections for the rest of the objects.

Table 2 Measurements parameters for the objects shown in Fig. 3, the exposure time is 1 s for all the examples
Full size table
Table 3 Time required for the scatter and the BH correction methods using four GPUs and the acquisition time required by the scanner
Full size table

The volumes from the scanner were first corrected from the scatter artifacts before the application of the proposed BH correction and the empirical BH correction algorithms. Although the scatter estimation is performed using the MC method which requires the knowledge of the materials, the results of the scatter correction of the first two rows in Fig. 3 were corrected from the scatter artifacts without the exact knowledge of the materials. The same method to estimate the materials for the BH correction method was used to estimate the materials for the scatter simulation. Thus, our assumption of the blindness of the proposed BH correction method does not contradict the approach we propose for the correction of the scatter artifacts. The scatter correction results and the results of the proposed BH correction method are shown in the second and third columns of this figure, respectively. While the BH correction results using the empirical method are shown in the fourth column. The correction of the scatter artifacts results in improved images in which the cupping and streak artifacts in the corrected images were highly suppressed. However, these images still suffer from noticeable cupping and streak artifacts. The application of the proposed BH correction method results in much-improved image quality in which the BH-corrected images are now better visualized and they are almost artifacts-free images. From the same figure, it is also shown that the proposed BH correction method produces better BH correction results in comparison to the results derived using the empirical method. As the corrected images using the empirical method still suffer from cupping and streak artifacts. The results in Fig. 3 were further analyzed in Sect. 3.2.1 in which the images are quantitatively evaluated in terms of CNR.

Table 3 shows the execution time required for the scatter and the proposed BH correction methods. With respect to the acquisition time from the scanner, i.e., 0.8 h for 3000 projections, the time required for the BH correction is near real-time using four GPUs. Since the number of projections, that are needed to be simulated, are distributed among multi-GPU. This scale the simulation time almost linearly with the number of GPUs. As a result, the simulation time for the full set of projections is divided by the number of GPUs. Real-time BH correction can be achieved with a higher number of GPUs used. Scatter correction was not performed for the last example, as it has only a minor contribution. Alternatively, we attached the near scatter-free image from the collimator for this example which is shown in the second column.

Fig. 4
figure 4

Results of the BH correction using the proposed method for different cases of materials estimation. The first column shows the slices from the BH-corrupted volumes. The second column shows the BH-corrected result when the materials are known. The third and fourth columns show the results when the materials are estimated with \(20\%\) lower and higher attenuation values than the values on the look-up table, respectively. The last column shows the comparisons between the profile lines of the slices in columns one to four. The profile lines were averaged over multiple lines to suppress the noise

Full size image

3.2 Analysis of the proposed BH correction method

3.2.1 Quantitative evaluation of the results

To quantitatively evaluate the results of the proposed BH correction and the empirical methods, the CNR is used, which is calculated by Eq. 16. The CNR was measured in a region-of-interest (ROI) of \(16\times 16\) pixels within the target images. The ROIs of the six examples are labeled with red squares in Fig. 3.

$$\begin{aligned} CNR=\frac{|\mu _{ROI}-\mu _{Background}|}{\sigma _{Background}},\ \end{aligned}$$
(16)

where \(\mu _{ROI}\) and \(\mu _{Background}\) are the mean values of the ROI and the background, respectively, while \(\sigma _{Background}\) is the standard deviation of the background.

Table 4 Quantitative evaluation of the scatter and the BH correction results shown in Fig. 3 using the CNR
Full size table

Table 4 shows the CNR results of the volumes shown in Fig. 3. For all the experimental examples, it is shown that the CNR values are improved after the application of the proposed BH correction method in comparison to the CNR values of the volumes in the case without any artifact correction or the case with only scatter artifact correction. Moreover, it can be noticed that the CNR results of most of the volumes corrected using the proposed BH correction method are higher than the CNR values derived from the volumes corrected using the empirical BH correction method. Therefore, better image visualization and quality have been achieved using the proposed BH correction method in comparison to the empirical method.

3.2.2 Effect of inaccurate estimation of materials on BH correction

The proposed materials estimation method depends on the polychromatic linear attenuation values of the different materials in the uncorrected volume and the polychromatic linear attenuation values of the relevant materials in the look-up table. If the volume, that is needed to be corrected, is derived from the same object and geometry setup which were used to fill the look-up table, the estimation of the material will be accurate. Otherwise, the linear attenuation values of the different materials in the uncorrected volume will deviate from the ones assumed in the table. This deviation could be a consequence of a stronger BH and scatter artifacts. Thus, the proposed materials estimation method would wrongly estimate these materials. Due to the proposed estimation scheme of the polychromatic projection, even when the materials are wrongly estimated, the proposed approach can estimate the polychromatic projection in a way that will be close to the BH-corrupted polychromatic projection. However, this is only possible when the materials estimated are not far away from their original types.

Figure 4 shows the effect of the wrong estimation of the materials on the BH correction results for two different objects. In this figure, three cases of materials estimation were tested. In the first case, the materials were assumed to be known. The results of the BH correction, in this case, are shown in the second column of this figure. In the second case, the materials from the uncorrected volume were assumed to have lower linear attenuation values than the values of the same materials on the look-up table. The deviation is assumed to be \(20\%\) between the two values. For the first object, the material was estimated as titanium instead of the assumed material, i.e., iron. For the second object, the materials were estimated as cement and titanium instead of aluminum and iron, respectively. The BH correction results of this case are shown in the third column of Fig. 4. In the last case, the materials were assumed to have higher attenuation values, i.e., \(20\%\) deviation was assumed again. For the first object, the material was selected as copper, while for the second object the materials were estimated as titanium and copper. The results of the BH correction, in this case, are shown in the fourth column of this figure. It is shown that the results of the BH correction for all the cases are almost the same and the wrong estimation of the materials has almost no effect on the quality of the correction using the proposed BH correction method. Table 5 shows the actual estimation of the materials by the proposed materials estimation method for the six experimental objects shown in Fig. 3.

Table 5 Materials estimation results
Full size table
Fig. 5
figure 5

Comparison of the BH-corrupted projections from the scanner and the projections estimated using the proposed polychromatic projection estimation method. The first column shows the BH-corrupted projections from the scanner. The second column shows the estimated projections. The last column shows the profile lines of the projections in the first two columns. The profiles are marked by white dashed lines in these sub-figures

Full size image

3.2.3 Robustness of the polychromatic projection estimation

As stated before, to derive an accurate correction term that can correct the BH artifact of the experimental scatter-corrected projections, the estimation of the polychromatic projections should be as close as possible to these projections. Figure 5 shows a comparison between the experimental scatter-corrected projections and the estimated polychromatic projections using the LSE method for the objects in rows 2–5 of Fig. 3. It is shown that the proposed LSE method can estimate a polychromatic projection which is very close to the experimental scatter-corrected projection in these examples.

4 Discussion

The application of the proposed BH correction algorithm improves the CNR ratio of the corrected images. The detectability of low-contrast materials in the corrected images is also improved which facilitates an accurate geometry extraction. The BH-corrected results, using the method in this work, were derived in near real-time based on limited available information from each scan. Based on the aforementioned, the method proposed provides a time-efficient alternative in near real-time to the classic iterative BH correction algorithms in industrial computed tomography.

5 Conclusion

In this work, a blind and robust BH correction method with unknown materials is proposed. The materials within the uncorrected volume are first estimated based on a look-up table of the measured polychromatic linear attenuation values. A BH artifacts correction term is then obtained by simulating the polychromatic and the monochromatic projections using an implemented GPU-based forward projection model. The polychromatic projection is estimated accurately utilizing the LSE method, while the monochromatic projection is estimated by selecting the energy bin which produces the lowest MSE with the acquired projection. Unlike most of the available BH correction methods, the proposed one corrects the BH artifacts in near real-time and does not require prior knowledge of the materials. Moreover, as the accuracy of the BH correction highly depends on the estimation accuracy of the polychromatic projection, it is shown that the proposed estimation method can predict a precise polychromatic projection that is very close to the used scatter-corrected projection from the scanner. Thus, a successful non-linearity adjustment of the BH-corrupted projections has been achieved which allows significant suppression of BH artifacts. We have validated the proposed BH correction method and compared with the representative empirical BH correction method. Experimental results show that the proposed method outperforms the empirical method with a significantly suppressed streak and cupping artifacts. This work shows that the proposed method is very effective and robust for the removal of the BH artifacts without prior knowledge of the materials which makes it very promising and practical for industrial applications.