Kun Yang | Ziteng Yang | Weina Yan | Jiankang Zhao | Yu Du | Shuang Liu | Kun Liu*
© 2019 IIETA. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).
OPEN ACCESS
This study aims to investigate a post-construction method that can correct hardening artifacts in computed tomography (CT) images without measuring the X-ray spectrum. Reconstructed images may include several hardening artifacts, which weaken the image quality and affect diagnosis. Therefore, herein, we adopt a reconstruction algorithm based on the equivalent tissue length approach. Firstly, the image pixels were divided into different equivalent tissues based on their CT values. Then, the equivalent tissues were projected to obtain their lengths in different ray directions. Considering the equivalent tissue length as the independent variable, a projection model was constructed and the correction coefficients were computed. Next, erroneous projection was identified and removed from the image. Finally, an artifact-free image was reconstructed from the corrected projection data. The results of a phantom model experiment and four patient data experiments show that the proposed method can effectively remove beam hardening artifacts and improve image accuracy. We believe the findings will be significant in clinical applications.
beam hardening, computed tomography (CT), equivalent tissue length, proportional guidance
Conventional reconstruction algorithms for computed tomography (CT) often assume X-rays to be monochromatic. However, in practice, X-rays are polychromatic and contain photons with various energies. When an X-ray beam passes through an object, low-energy photons are more likely to be absorbed than high-energy photons. Consequently, beam attenuation is not related to object thickness linearly; ignoring this nonlinear relation during image reconstruction produces a resulting image containing several hardening artifacts, thereby reducing image quality and affecting the diagnosis. In the last few decades, many methods have been proposed for reducing such hardening artifacts [1, 2].
Existing artifact removal methods fall into two broad categories: hardware and software correction methods. The hardware correction method filters low-energy photons by placing a metal filter between the X-ray source and the object. Although this approach is easy to implement and effectively reduces artifacts, the filtering process also reduces the signal-to-noise ratio (SNR) [3]. Typical software correction strategies include dual-energy correction, single-energy correction, iterative correction, and polynomial fitting. We briefly describe each of these approaches below.
As proposed by Alvarez et al., theoretically, dual-energy correction is the most effective method to remove hardening artifacts [4, 5]. This method considers linear attenuation to be the combined result of the photoelectric and Compton scattering effects. However, it requires two radiation sources and two detectors, and these additional instruments lower the scanning efficiency and increase the patient’s radiation dose, making the method less practical.
Single-energy correction is based on a hypothesis by Nalcioglu and Lou [6]. However, their original method does not produce good results if the scanning area contains more than two materials. This method was later improved by Joseph et al. [7] to handle cases where the scanning area contains three materials; however, the computation becomes increasingly complex with the increasing number of materials.
Iterative correction is a relatively new method that iteratively recurses new values with old ones. For example, Van Gompel et al. [8] introduced an iterative correction method based on a cost function. However, their nonconvex cost function yields multiple solutions; thus, the solution quality depends on the rational selection of initial parameter values. Zhao et al. [9] presents another iterative correction strategy that can be applied to various hybrid materials, but it can’t remove all hardening artifacts. In summary, iterative correction currently has considerable practical limitations and thus cannot be widely used.
Polynomial fitting essentially converts the measured polychromatic projection data into monochromatic data, and then uses these data to reconstruct the object image. This is currently the most popular approach for artifact correction owing to its simplicity, ease of use, and good correction ability [10-12]. For instance, Kachelrieß et al. [11]; Kyriakou et al. [12] construct polynomial equations using the segmented tissue image as the independent variable, and then compute the correction coefficients by minimizing the total variation. Although this produces good correction, the calculations are complicated. Another approach is to correct the hardening artifacts post reconstruction using a ladder-like experimental database of equivalent tissues [13-16]. Although this method is effective and convenient for simple objects, it is not sufficiently flexible to handle objects with complex components, where the tissue database is unfeasibly large.
In our previous work [17], we presented an algorithm based on the mathematical perspective of constructing a projection correction model using the projection data as the independent variable. This method is suitable for hardening artifact correction with a single homogeneous material or a combination of two such materials. It can be used for water pre-correction in conventional CT, as well as hardening artifact correction for both single parts in industrial CT images and two-component splicing materials.
Polynomial fitting is the most popular approach for correcting beam hardening artifacts in CT images. However, existing correction methods are either based on image space or projection space; i.e., they consider the base image or base projection as the independent variable. In addition, such need to mathematically solve for the correction coefficients to complete the correction. For example, Kyriakou et al. [12] solves for the coefficients by minimizing the total variation, but the solution process is rather complicated. Yang et al. [15] obtains the coefficients by minimizing the entropy. Although such methods can remove the hardening artifacts, they require complex computation. In practice, we cannot accurately represent the imaging law by considering the base image or base projection as the independent variable.
To resolve this issue, we develop a simple and universal beam hardening correction (BHC) method that corrects beam hardening artifacts under the proportional guidance of the equivalent tissue lengths. Specifically, we construct an innovative projection model that uses the equivalent tissue length as the independent variable.
In the proposed model, correction coefficients reflect the correlations between the attenuation coefficients, with low-order terms describing the actual attenuation projection and high-order terms representing the hardening process. (Each term in the projection model carries a specific physical meaning.) Thus, the entire projection model reflects the actual polychromatic projection. Next, we employ an established model to approximate the collected polychromatic projection and obtain the correction coefficients via the least-squares method. Compared with previous polynomial fitting methods [18, 19], our approach is simple and fast, and every term has a clear physical meaning, making the approach easy to understand and utilize.
During image segmentation, the equivalent tissues are often divided inaccurately because of the improper operation or lack of the experience. This inaccuracy may cause errors in the calculated equivalent tissue lengths, which in turn cause discrepancies in the calculated attenuation coefficients (i.e., differences in the attenuation coefficient ratio of water to bone obtained from the equivalent single-energy process). In such cases, the hardening artifacts cannot be completely corrected. This led us to propose the proportional guidance method for obtaining the correction coefficients. Our analysis results show that, under proportional guidance, our correction method can overcome the common practical issue of incomplete segmentation, demonstrating the universal applicability of our method.
Considering the equivalent tissue length as the independent variable, we propose a BHC algorithm that removes hardening artifacts by eliminating the erroneous component of the projection. Figure 1 shows a detailed flowchart of the proposed method. Based on the physical laws of ray projection, we constructed the projection model through the following five steps.
First, we acquired the original image f0 (with hardening artifacts) from the CT image. Next, we decomposed each image pixel according to the proportions of air, water, and bone using the proportional division function. Following that, we separately projected the water and bone tissue regions to determine the total length of each tissue type for each ray. Then, we removed the nonlinear error term from the original projection, yielding the corrected projection. Finally, we reconstructed an artifact-free image based on the corrected projection. While in existing methods [14], erroneous components must be removed from the original artifact-containing image, our algorithm can correct all artifacts in a single reconstruction and output quantifiable CT values.
Figure 1. Flowchart of the proposed method
2.1 Equivalent tissues
We roughly determined the CT distribution of each tissue type from a histogram of the CT value distribution. Each pixel in the artifact-containing image was decomposed according to the proportions of the three equivalent tissue types (air, water, and bone) using the proportional division function (Figure 2). Here, the fractions of the three tissue types always sum to one.
Figure 2. Proportional division function
We selected −1,000 Hounsfield units (HU), denoted by T1, as the upper limit for air, and adopted 1,300 HU, denoted by T4, as the lower limit for bone tissue and 0–100 HU, denoted by T2–T3, as the range for water. Thus, a pixel with a CT value below T1, a CT value above T4, or a CT value between T2 and T3, is considered pure air, bone, or water, respectively. Furthermore, a pixel with a CT value between T3 and T4 is deemed a mixture of water and bone. Herein, we denote the equivalent tissue proportions of water, bone, and air as Ww, Wb, and Wa, respectively, and the pixel CT value is denoted as Z.
The proportions of the three equivalent tissues in each pixel are determined using the proportional division function (i.e., formulas (1) and (2)), which were established by considering that X-rays are attenuated more significantly by passing through bone than by passing through the same thickness of water. The proportional relationships for water and bone can be expressed as follows:
$W_{w}=\left\{\begin{array}{cc}{\frac{z-T_{1}}{T_{2}-T_{1}}} & {T_{1} \leq z<T_{2}} \\ {1} & {T_{2} \leq z<T_{3}} \\ {\cos ^{2}\left(\frac{\pi}{2} * \frac{z-T_{3}}{T_{4}-T_{3}}\right)} & {T_{3} \leq z<T_{4}}\end{array}\right.$ (1)
$W_{b}=\left\{\begin{array}{cc}{0} & {z<T_{3}} \\ {\sin ^{2}\left(\frac{\pi}{2} * \frac{z-T_{3}}{T_{4}-T_{3}}\right)} & {T_{3} \leq z<T_{4}} \\ {1} & {z \geq T_{4}}\end{array}\right.$ (2)
For the original image f0 in Figure 3(a), we can use these division functions to obtain the equivalent tissue images for water and bone as Figure 3(b) and 3(c), respectively.
Figure 3. Original image (a), and equivalent tissue images for water (b) and bone (c)
2.2 Projection model
Let E0 be the energy of the monochromatic X-ray source. When an object is irradiated by this source, the projection q is linearly correlated with the thickness through which the rays penetrate [15]. In this case, the projection can be described as:
$q=\int_{0}^{D} \mu\left(E_{0}, x\right) d x$ (3)
where, x is the total thickness penetrated by the X-rays and μ(E0,x) is the cross-sectional distribution of the object’s linear attenuation coefficient for an energy of E0.
In the previous subsection, we obtained the equivalent tissues for air, water, and bone using the proportional division function. Hardening is mainly due to the attenuation of low-energy X-rays in water and bone because radiation is only slightly attenuated in air. Hence, herein, we primarily consider the hardening effects of water and bone, and ignore the effect of air. Then, we perform separate forward projections for the equivalent tissues of water and bone Ww and Wb to obtain the equivalent tissue lengths Lw and Lb for each ray at each angle. Based on these projections, we construct a projection model via polynomial fitting to correct the considered image.
During the scan, each ray simultaneously passes through water, bone, and air, leaving a projection P on the detector. If the beam was purely monochromatic, then P would be linearly correlated with Lw and Lb. However, because actual beams are polychromatic, P’s relationship with Lw and Lb is actually nonlinear. This implies that the individual or combined effects of the high-order terms for bone and water create the hardening artifacts. Specifically, high-order terms for bone produce bone hardening artifacts, whereas those for water produce water hardening artifacts. Here, we neglect higher-order terms, such as third- and fourth-order terms, and construct the following theoretical model for polychromatic projection:
$P=c_{1} L_{w}+c_{2} L_{b}+c_{3} L_{b}^{2}+c_{4} L_{w b}+c_{5} L_{w}^{2}$ (4)
Thus far, we determined Ww and Wb from the proportional division function and acquired Lw and Lb via forward projection. Now, using (4), we calculated the Lb×Lb, Lw×Lb, and Lw×Lw terms, denoted as $L_{b}^{2}$, Lwb, and $L_{w}^{2}$ , respectively. Zhang et al. [17] obtained artifact-free images by removing the $L_{b}^{2}$ error terms. Accordingly, we modified the polychromatic projection model to:
$P=c_{1} L_{w}+c_{2} L_{b}+c_{3} L_{b}^{2}$ (5)
2.3 Attenuation coefficient determination
Here, we discuss (5) further as the main aim of the linear fitting correction process is to convert polychromatic projection data into monochromatic data for the equivalent single-energy source. We divided the projection terms into linear (c1Lw+c2Lb) and nonlinear ( $c_{3} L_{b}^{2}$ ) terms. Then, we obtained monochromatic projection (c1Lw+c2Lb) for the equivalent single-energy source by eliminating the nonlinear terms from the projection data. According to the Beer–Lambert Law, the projection data is given by the product of the linear attenuation coefficient and the thickness penetrated by the ray in the monochromatic projection, as shown in (3). Coefficients c1 and c2 in the linear term (c1Lw+c2Lb) in (5), which was modified from (3), should be proportional to the linear attenuation coefficients of water and bone. Therefore, to obtain c1 and c2, we used the proportional guidance method and thus ensured that the linear attenuation coefficients of water and bone matched the actual law in the single-energy case.
Referring to the table of X-ray mass attenuation coefficients provided by NIST [20], we obtained the mass attenuation values of water and bone at different energies. The relationship between mass and linear attenuation is as follows:
$\mu_{L}=\mu_{m} \rho$ (6)
where, μL and μm are the linear and mass attenuation coefficients, respectively, and ρ is the density. we obtained the ratios of the linear attenuation coefficients of bone to water at different energies, as shown by the fitted curve. In Figure 4, we can observe that although the ratio is different at different energies, it has a fixed value for any given energy, such as that at point A.
Figure 4. Decay ratios of bone to water at different energies
The linear correction method involves correcting the polychromatic projection data to match that for the equivalent single-energy case. In that case, the linear attenuation coefficient ratio has a certain value. Our results indicate that the proportional guidance method can achieve good correction for energies of 40-60. For example, in Experiment 1 (Section 3.1), when we do not apply the proportional guidance method, the ratio may substantially differ from the true value, resulting in artifacts such as those observed in Figure 5.
Figure 5. Original and corrected images for Experiment 1
The proportional value, denoted by t, indirectly depends on the equivalent energy:
$\frac{c_{1}}{\mu_{\text {water}}}=\frac{c_{2}}{\mu_{\text {bone}}}=t$ (7)
Equation (7) can be rewritten as:
$\frac{c_{2}}{c_{1}}=\frac{\mu_{\text {bone }}}{\mu_{\text {water }}}=t$ (8)
where, μwater and μbone are the linear attenuation coefficients of bone and water, respectively. Eq. (8) shows that the ratio of c2 to c1 can be approximated by the linear attenuation coefficient ratio of bone to water. We therefore introduced the proportional guidance method to obtain more accurate values for the coefficients c1, c2, and c3. First, we determined the ratio of c2 to c1 based on the ratio t for the monochromatic case:
$t=\frac{c_{2}}{c_{1}}$ (9)
For example, the linear attenuation coefficient ratio of polyvinyl chloride (PVC) to water is 2.8250, i.e., t=2.8250, at an equivalent energy of 50 kVp. Here, we can obtain the value of t from the image’s CT value as follows:
$C T=\frac{\mu-\mu_{\text {water }}}{\mu_{\text {water }}} \times 1000$ (10)
Eq. (10) shows the relationship between the linear attenuation coefficient ratio and the CT value. Combining (8) and (10), we obtain:
$t=\frac{c_{2}}{c_{1}}=\frac{\mu_{\text {bone }}}{\mu_{\text {water }}}=\frac{C T_{\text {bone }}}{1000}+1$ (11)
We determined the value of t via self-adaptive proportional guidance as follows. First, we extracted N points with high CT values (i.e., bone pixels) from the original image f0 without water hardening artifacts, and calculated their mean CT value. Herein, N was determined based on the image’s size and actual bone content. Then, we substituted the mean CT value for bone into (11) to obtain a more representative t value that is more in line with physical values.
Substituting this t value into (5), we can simplify the polynomial as follows:
$P=c_{1} L_{w}+t c_{1} L_{b}+c_{3} L_{b}^{2}=c_{1}\left(L_{w}+t L_{b}\right)+c_{3} L_{b}^{2}$ (12)
Thus, only two fitting coefficients were left to be found.
2.4 Coefficient calculation
First, we calculated the projection data p0 by applying a Radon transform to the original image f0. Next, we approximated p0 by fitting the projection data P using (12). Following that, we determined the unknown coefficients c1 and c3 via the least-squares method:
$E=\min \iint\left(P-p_{0}\right)^{2} d x d y$ (13)
where, E represents the fitting error. Then, we determined the unknown coefficients in (13) using:
$\nabla_{\mathrm{c}} E=\nabla_{\mathrm{c}} \iint\left(P-p_{0}\right)^{2} d x d y=0$ (14)
After obtaining the coefficients c1 and c3 by solving (14), we calculated the value of c2 from c1 and t, and the erroneous projection Perror was calculated as:
$P_{\text {error}}=c_{3} L_{b}^{2}$ (15)
To approximate the single-energy projection, we removed Perror from the original projection $p_{0}$ to obtain the final corrected monochromatic projection Pfinal:
$P_{\text {final}}=p_{0}-P_{\text {error}}$ (16)
Finally, we reconstructed the artifact-free image ffinal from Pfinal.
2.5 Evaluation
We evaluated the proposed beam hardening method’s performance on both phantom data and patient data obtained using a cone-beam spiral CT scanner (Beijing Arrays Medical Imaging Corporation, China). The bulb voltage, bulb current, medium thickness, source–detector distance, and source–patient distance was 120 kVp, 280 mA, 1.25 mm, 988 mm, and 560 mm, respectively. In addition, the reconstructed convolution kernel was adopted as the standard convolution kernel.
2.5.1 Phantom model
To simulate the water tissue in a human head, we filled a cylindrical bucket with water. The bucket’s diameter was 20 cm, which is close to that of a human head. In addition, we used two 30-mm-diameter pure PVC rods (SIMONA, Germany), referred to as bone rods hereafter, to simulate human bones as their attenuation effect on X-ray is similar to that of human bones. We placed the two bone rods symmetrically around the water-filled bucket, then placed the bucket on the scanning table and conducted a spiral CT scan.
The main difference between the two experimental materials is that, under irradiation by X-rays of the same energy, the X-ray attenuation caused by water is less than that caused by bone; i.e., compared with water, bone is a strong attenuator. We built this phantom model to simulate the hardening effect of bone. Specifically, because the cylinder was similar in size to a human head, we used it to simulate the hardening effect of a human head.
2.5.2 Patient data
Four volunteers participated in Experiment 2. During the experiment, they were asked to lie flat on the scanning table, following which they were subjected to spiral CT scans of their heads.
3.1 Experiment 1
Figure 5 compares the original and corrected images. The uncorrected image (Figure 5(a)) shows clear black stripe-shaped artifacts (high-attenuation objects) between the models. However, post correction using the proposed method, we obtained a corrected image (Figure 5(b)) in which the artifacts were qualitatively removed. This image is in good agreement with previous correction results [12, 16].
To further investigate the correction effect, Figure 6 compares the CT values of pixels from horizontal and vertical stripes in the middle of the original and corrected images. Figure 6 shows that the original image contains a cupping artifact for the bone (Figure 6(a)), which is largely absent from the corrected image. In addition, Figure 6(b) and the magnified view in Figure 6(c) show that the CT values in the water tissue region have been fully restored, corresponding to the elimination of the black stripe artifact in Figure 5(b).
To evaluate the benefits of the proportional guidance method, we also corrected the results of Experiment 1 by directly applying the traditional polynomial correction method in (5), without proportional guidance. Figure 7 shows the results of applying both methods to correct the same CT image. The original image (Figure 7(a)) contains an obvious black stripe artifact, which was reduced but not eliminated via polynomial fitting (Figure 7(b)). In contrast, the black stripe artifact has clearly been removed by our proposed method, as observed in Figure 7(c).
Figure 6. Comparison of CT values for pixels from vertical (a) and horizontal (b) stripes in the middle of the original and corrected images. Additionally, (c) shows a magnified view of (b)
Next, we compared the pixels from horizontal and vertical stripes in the middle of the original and corrected images to further evaluate the CT value variations produced by correction.
Figure 7. Comparison of the original image (a) with the results obtained via standard polynomal correction (b) and the proposed method (c)
Figure 8 shows the CT values from vertical and horizontal stripes in the middle of the original image and the image corrected by the standard polynomial method. Figure 8(a) and 8(b) compare the pixel contrast values for the vertical and horizontal stripes, respectively.
Next, we selected regions of interest (ROIs) from the images. Table 1 shows the mean CT values in a water tissue region (ROI1) for the original and corrected images, with and without proportional guidance.
(a)
(b)
Figure 8. Comparison of CT values for pixels from vertical (a) and horizontal (b) stripes in the middle of the original image and the image corrected via the standard polynomial method
Table 1. Comparison of mean CT values in ROI1
|
Image |
Mean CT value (of ROI1) |
|
Original image |
-22 |
|
Image corrected by our method |
-3 |
|
Image corrected by the polynomial method |
-19 |
As shown in Table 1, the mean CT value in ROI1 was -22 HU before correction, well outside the normal range of 0 ± 4 HU. Without correction, this image could produce serious errors in diagnosis. In contrast, the mean CT values become -19 and -3 after being corrected via the polynomial and proposed methods, respectively. Both values are higher than the mean CT value for the original image, indicating that both methods can improve image quality. However, only our method resulted in a value within the normal mean CT value range, indicating that proportional guidance is an important artifact removal tool.
To make these comparison results more intuitive, Figure 9 shows the CT values of pixels in the middle of Line 1, comparing the original image with those corrected via the proposed and polynomial methods. Here, unlike the polynomial method, the proposed method managed to restore the bone rods’ CT values and suppress the cupping artifacts.
Figure 9. CT values for pixels in the middle of Line 1, for the original image and the images corrected via the proposed and polynomial methods
In addition, the polynomial method yielded a c2/c1 ratio of 2.7479, which is smaller than that yielded by our method (2.8586); hence, the result differed from the linear attenuation coefficient ratio of bone to water under monochromatic conditions. Although both correction methods could suppress the image artifacts, our method’s CT value is closer to the original true value. Therefore, our correction method, which uses proportional guidance, can achieve good correction results.
3.2 Experiment 2
Figure 10 shows the test results for one of the volunteers. Here, there is a significant stripe artifact between the bone rods in the original image (Figure 10(a)). However, we cannot observe this artifact in the images corrected by or method or the polynomial method (Figure 10(b)), demonstrating the good performance of both methods.
To further compare the effects of the two methods, we selected two ROIs in the image, denoted as ROI2 and ROI3 and representing water tissue and reference regions, respectively. Table 2 shows the mean CT values in both ROIs, both before and after correction.
Figure 10. Original and corrected images for Experiment 2
Table 2. Comparison of mean CT values in ROI2 and ROI3
|
Image |
Mean CT value (ROI2) |
Mean CT value (ROI3) |
|
Original |
20 |
41 |
|
Corrected (our method) |
39 |
41 |
|
Corrected (polynomial method) |
37 |
41 |
As shown in Table 2, correction by our method increased the mean CT value in ROI2 from 20 HU to 39 HU, close to that of water tissue in the other parts of the brain. Meanwhile, correction by the polynomial method increased the value from 20 HU to 37 HU. Although both methods improved the CT value, our method came closer to the mean reference value (ROI3), suggesting that it is more reliable and precise for image correction.
In addition, the c2/c1 ratios were 2.4575 and 2.4376 for our method and the polynomial method, respectively. Introducing proportional guidance thus reduced the ratio. As mentioned above, correction by our method yielded a higher mean CT value in ROI2, suggesting the resulting data is more accurate and reliable and demonstrating that it can achieve good correction results.
Next, we applied our method to all four volunteers’ scans to verify its validity and universality; Figure 11 shows the results. For the first patient, the original image (Figure 11(a1)) shows a hardening artifact, due to the presence of bones. The dark region (indicated by an arrow) has a low grayscale value, while the bright region (circled) has a high value. The CT values are not uniform in the soft tissue region of the brain. By contrast, after correction by our method (Figure 11(a2)), the distribution is now uniform in the soft tissue region, and the hardening artefact has been removed.
For the second patient, the original image (Figure 11(b1)) shows three obvious stripe artifacts (arrows), which can no longer be seen in the corrected image (Figure 11(b2)). For the third patient, the original image (Figure 11(c1)) shows an obvious stripe artifact (arrows) and a relatively dim region (circled). In the corrected image (Figure 11(c2)), the stripe artifact has disappeared, the values have returned to normal in the circled area, and the CT values are more uniform across the soft tissue region. For the fourth patient, the original image (Figure 11(d1)) shows a black hardening artifact (arrow) between the bones and a low-brightness area (circled) in the soft tissue region. In the corrected image (Figure 11(d2)), the stripe artifact has disappeared, the brightness is normal in the circled region, and the values in the soft tissue region are highly uniform. These four examples demonstrate our algorithm’s wide applicability to artifact removal and CT value restoration.
Finally, we investigated the applicability of our method further by correcting CT images of human heads. Specifically, we used the proposed method to correct 104 CT images obtained from different layers of the human head in Experiment 2. Due to space constraints, we cannot present all these images in this paper. Instead, Figure 12 compares how the CT values varied with the layer, with and without correction.
Figure 11. Original and corrected images for all four patients in Experiment 2
Figure 12. Contrast maps of a human head before and after CT correction
As shown in Figure 12, correction does not change the average CT values of the target images. In other words, the correction process only corrects the CT values in hardening artefacts, leaving the CT values in other tissues unchanged and avoiding any increase in the CT values across the image. Thus, our algorithm is both widely applicable and robust.
In order to more fully demonstrate the broad applicability of the proposed correction method, we constructed a model of human length. Specifically, we constructed a model of the human body (30 cm in diameter) consisting of a cylindrical bucket full of water with two PVC bone rods distributed evenly around the bucket. Figure 13 compares the proposed correction method with the scanner’s native correction method. Figure 13(a) includes obvious black stripe artifacts between the high-attenuation regions. The scanner’s native correction method (Figure 13(b)) reduced but does not remove the stripe artifacts. By contrast, our method (Figure 13(c)) entirely eliminated all the stripe artifacts. Thus, our method achieved clearly better correction than the scanner’s native method.
Figure 13. Original image (a), compared with the results of the scanner’s native correction method (b) and the proposed algorithm (c)
The scanner’s native correction method is based on the projection. This is essentially a fitting strategy with the projection data as the independent variable, and is similar to the approaches in [14, 16]. By contrast, our approach is based on the equivalent tissue length, and uses proportional guidance to calculate the correction coefficients, which is why it achieved better results.
The differences in the two methods’ results come mainly from the following facts. In the scanner’s method, even though some bones are in direct contact with the air in the model, the segmentation function neglects the bone–air interface and the resulting segmentation errors lead to discrepancies in the projection calculations. Since the projection acts as an independent variable, this greatly affects fitting results, meaning that some artefacts remain in the corrected image (Figure 13(b)). In our method, we instead take the equivalent tissue length as the independent variable for fitting, and we use proportional guidance to ensure our model obeys the true projection law. This enabled our approach to effectively eliminate the artifacts, as shown in Figure 13(c). In summary, proportional guidance can solve the correction errors induced by inaccurate segmentation.
Figure 14. Original image (a), compared with the results of the scanner’s native correction method (b) and the proposed algorithm (c)
In order to illustrate the proposed algorithm’s applicability to different CT machines, we also scanned a volunteer using another model of spiral CT scanner, and again compared the proposed correction method with the machine’s native correction approach. The results are shown in Figure 14.
The original image (Figure 14(a)) includes obvious black stripe artifacts between the high-attenuation regions. After correction by the scanner’s native method (Figure 14(b)), the stripe artifacts have been reduced but not removed. By contrast, our method (Figure 14(c)) was able to completely eliminate all the stripe artifacts. Thus, our method again achieved better correction results than the native method.
The results of the phantom model and patient data experiments show that our correction method, i.e., beam hardening correction with proportional guidance based on equivalent tissue length, can effectively reduce hardening artifacts and enhance image precision.
Without proportional guidance, the polynomial method could eliminate the artifacts, improve image quality, and approximate the actual pixel values in the patient data experiments; however, it could not fully eliminate the hardening artifacts in the phantom model experiment because the proportional division function ignores the air–bone interface, even though they in direct contact with each other in the phantom model. Thus, the equivalent tissues were not accurately divided, and hence the equivalent tissue lengths could not reflect the actual contributions of the water and bone lengths to the projection. Importing these inaccurate data into the correction model introduced errors to the correction results. In contrast, the bone was not in contact with the air in the patient data experiments. Thus, the proportional division function could accurately divide the equivalent tissues. Consequently, the equivalent tissue lengths, linear correction coefficients, corrected projection, and final image were all accurately obtained.
In summary, the proposed method successfully uses proportional guidance to prevent inaccurate division between equivalent tissues. Since there are frequent errors in this division process in real-world scenarios, this shows that the proposed correction method is widely applicable.
This research was supported by the Natural Science Foundation of Hebei Province (H2019201378), the key R&D Program of the Ministry of Science and Technology (2016YFC0104203), Foundation for High-level Talents in Higher Education of Hebei Province (B20190030010).
[1] Bracewell, R.N., Riddle, A.C. (1967). Inversion of fan-beam scans in radio astronomy. The Astrophysical Journal, 150(2P1): 427-434. http://dx.doi.org/10.1086/149346
[2] Kak, A.C., Slaney, M., Wang, G. (2002). Principles of computerized tomographic imaging. Medical Physics, 29(1): 107-107. http://dx.doi.org/10.1118/1.1455742
[3] Gao, H.W., Zhang, L., Chen, Z.Q., Xing, Y.X., Li, S.L. (2006). Beam hardening correction for middle-energy industrial computerized tomography. IEEE Transactions on Nuclear Science, 53(5): 2796-2807. http://dx.doi.org/10.1109/TNS.2006.879825
[4] Marshall, W.H., Alvarez, R.E., Macovski, A. (1981). Initial results with prereonstruction dual-energy computed tomography. Radiology, 140(2): 421-430. http://dx.doi.org/10.1148/radiology.140.2.7255718
[5] Maass, C., Meyer, E., Kachelriess, M. (2011). Exact dual energy material decomposition from inconsistent rays (MDIR). Medical Physics, 38(2): 691-700. http://dx.doi.org/10.1118/1.3533686
[6] Nalcioglu, O., Lou, R.Y. (1979). Post-renconstruction method for beam hardening in computed tomography. Physics in Medicine and Biology, 24(2): 330-340. http://dx.doi.org/10.1088/0031-9155/24/2/009
[7] Joseph, P.M., Ruth, C. (1997). A method for simultaneous correction of spectrum hardening artifacts in CT images containing both bone and iodine. Medical Physics, 24: 1629-1634.
[8] Van Gompel, G., Van Slambrouck, K., Defrise, M., Batenburg, K.J., de Mey, J., Sijbers, J., Nuyts, J. (2011). Iterative correction of beam hardening artifacts in CT. Medical physics, 38(S1): S36-S49. http://dx.doi.org/10.1118/1.3577758
[9] Zhao, Y.S., Li, M.F., Zhang, P. (2015). Iterative beam hardening correction for multi-material objects. PLoS One, 10(12): e0144607. http://dx.doi.org/10.1371%2Fjournal.pone.0144607
[10] Kijewski, P.K., Bjärngard, B.E. (1978). Correction for beam hardening in computed tomography. Physics in Medicine and Biology, 5(3): 209-214. http://dx.doi.org/10.1118/1.594429
[11] Kachelrieß, M., Sourbelle, K., Kalender, W.A. (2006). Empirical cupping correction: A first-order raw data precorrection for cone-beam computed tomography. Medical Physics, 33(5): 1269-1274. http://dx.doi.org/10.1118/1.2188076
[12] Kyriakou, Y., Meyer, E., Prell, D., Kachelriess, M. (2010). Empirical beam hardening correction (EBHC) for CT. Medical Physics, 37(7): 5179-5187. http://dx.doi.org/10.1118/1.3477088
[13] Chen, C.Y., Chuang, K.S., Wu, J., Lin, H.R., Li, M.J. (2001). Beam hardening correction for computed tomography images using a post reconstruction method and equivalent tissue concept. Journal of Digital Imaging, 14(2): 54-61. http://dx.doi.org/10.1007/s10278-001-0003-2
[14] Shanghai Lianying Medical Technology co, Ltd. (2013). A correction method of bone hardening artifacts in CT image reconstruction: in Chinese, CN 103186883 A.
[15] Yang, K., Xue, L.Y., Yin, K., Liu, S., Meng, J. (2017). Microbubble generation and trapping induced by femtosecond laser and acoustic signal analysis. Traitement du Signal, 34(1-2): 33-44. http://dx.doi.org/10.3166/TS.34.33-44
[16] Schüller, S., Sawall, S., Stannigel, K., Hülsbusch, M., Ulrici, J., Hell, E., Kachelrieß, M. (2015). Segmentation-free empirical beam hardening correction for CT. Medical Physics, 42(2): 794-803. http://dx.doi.org/10.1118/1.4903281
[17] Zhang, X.K., Xue, L.Y., Yin, K., Wang, T.Y., Yang, K. (2019). Computed tomography beam hardening correction based on non-linear segmentation. International Journal Bioautomation, 23(1): 87-96.
[18] Elbeltagy, A.E.H.M., Youssef, A.M., Bayoumy, A.M., Elhalwagy, Y.Z. (2018). Fixed ground-target tracking control of satellites using a nonlinear model predictive control. Mathematical Modelling of Engineering Problems, 5(1): 11-20. http://dx.doi.org/10.18280/mmep.050102
[19] Xie, Q., Zhang, C.Z. (2006). CT imaging technology – principle, design, artifacts and progress. Science press.
[20] Hubbell, J.H., Seltzer, S.M. (2004). X-ray mass attenuation coefficients. Physical Measurement Laboratory. https://dx.doi.org/10.18434/T4D01F
