A new metal artifact reduction algorithm based on a deteriorated CT image

Abstract

BACKGROUND: Computed tomography (CT) is an established imaging technology primarily used as a non-invasive diagnostic tool that reconstructs axial images. However, significant problems with metal artifacts remain. A metal artifact is a strong radial noise in an image, which makes it difficult to diagnose patients and inspect products containing metal implants. Historically, studies related to metal artifact reduction used projection data, though the data is not typically saved after processing.

OBJECTIVE: This study proposes a new metal artifact reduction algorithm that does not require projection data, for new applications and for accurate diagnostic techniques.

METHODS: The algorithm utilizes reconstructed images and is based on iterative reconstruction. By reproducing an accurate forward projection on simulation and combining it with an iterative calculation, discrepancies causing metal artifacts are eliminated.

RESULTS: Validation was completed with numerical phantom models. Our results indicate that the proposed algorithm effectively reduces metal artifacts, even if numerous complex-shaped metal pieces were embedded in the cross-section.

CONCLUSIONS: We developed a novel reconstruction algorithm for metal artifact reduction in CT imaging that only requires reconstructed images and projection conditions. Any historical CT data containing metal artifacts could be improved with this method.

Keywords

Computer tomography (CT)metal artifact image reconstruction image processing iterative reconstruction

1 Introduction

Computed tomography (CT) is a method for generating cross-sectional images of an object. Generally, two steps are involved in this process. First, an X-ray source rotates around an object and diametrically positioned sensors record the projection data. The second step involves computer-aided reconstruction. By superimposing the detection data along the projected beam to an empty image on a computer, a cross-sectional image can be obtained. The most common reconstruction algorithm is filtered-back projection (FBP).

CT is commonly used as a diagnostic tool in medicine such as in the early detection of cancer, as well as for non-destructive inspection in industrial engineering applications. Despite its widespread use, technical challenges remain. A pervasive problem is an imaging artifact caused by radial noise resultant of objects containing high-Z materials; typically metals (Fig. 1). Metal artifacts make it difficult to diagnose patients and inspect products containing metal implants, and this continues to be problematic in CT imaging; numerous studies have been dedicated to this subject.

Alvarez et al. noted that the absorption of X-rays by materials depends on the incident energy spectrum, and proposed a reconstruction algorithm based on projection data collected at two different energies [1]. Methods using this routine are described as Dual Energy (DE) methods, and many studies related to DE methods have been performed [2 –7]. DE methods have been applied to various types of artifact reduction including metal artifacts. However, these methods often require hardware modification. Furthermore, because DE methods require two scans, patient dose and robustness of the imaging routine are of concern.

E. Meyer et al. proposed normalize metal artifact reduction (NMAR) [8], which is divided into a segmentation operation and a normalizing operation. This method partially removed metal artifacts, but generated other streaking artifacts. In addition, it was assumed that metal regions could be extracted accurately by thresholding.

Other methods of metal artifact reduction include Zhao et al.’s method of using a wavelet-based weighting function in the field of view [9]. Zhang et al. modified metal regions on CT images through projection data from two directions [10]. Yazdi et al. interpolated missing projections in raw CT data using unaffected projections [11]. Abdoli et al. corrected projection data using weighted virtual sinograms [12]. Iterative reconstruction of back projection data is a widely used CT imaging technique and several studies have attempted to exploit this aspect of the process to reduce metal artifacts [13, 14]. Although other modifications of projection data exist, they are limited to application on simple shapes. Additionally, the metal artifact reduction may prove insufficient for widespread use. Furthermore, almost all of the proposed methods use projection data. In general, CT imaging systems save reconstructed image data, although they do not save projection data as the file size is prohibitively large. Therefore, even if metal artifact reduction algorithms are established, they cannot be applied to historical CT data if based on projection data.

There are also some methods that are based on reconstructed images. M. Bal et al. conducted adaptive filtering algorithm [15], and segmented reconstructed images into different material classes. Y. Chen et al. enhanced metal artifacts with a large-scale non local mean filter [16], and segmented images into metal parts and artifact parts. However, these segmentation methods require to adjust parameters by hand, and the result is highly dependent on the amount of metals and the material structure of the cross-section.

In this paper, we propose a new reconstruction algorithm for metal artifact reduction using reconstructed images that does not require projection data and applicable in any situation.

2 Metal artifacts

Metal artifacts are largely due to improper consideration of the X-ray energy spectrum during image reconstruction. The Beer–Lambert law describes the attenuation of X-rays passing through an object, $I_{out} = I_{in} exp (- μ t),$ (1) where I_in and I_out are the incident and transmitted X-ray intensities, respectively, μ denotes the material-specific X-ray attenuation coefficient, and t is the path length. From Equation (1) we define the projection data p as: $p = - ln \frac{I_{out}}{I_{in}} = μ t .$ (2)

We see that p is linear with respect to an object thickness t when the object is homogeneous, which is only true for monochromatic X-rays. Modern CT imaging systems use polychromatic X-ray fields. Because attenuation depends on energy, the transmission intensity of polychromatic X-rays, I_out, is given by: $I_{out} = \int_{E} I_{in} (E) exp (- μ (E) t) dE .$

The projection data p becomes: $p = - ln \frac{I_{out}}{I_{in}} = - ln \frac{\int_{E} I_{in} (E) exp (μ (E) t) dE}{\int_{E} I_{in} (E) dE} .$

Thus for polychromatic X-rays, p is a non-linear function of object thickness t. However, in CT the polychromatic X-rays are assumed to be monochromatic, using Equation (2). This assumption makes differences for scanning targets, and makes beam hardening artifacts. Some of these artifacts can be reduced by simple operations when the targets are non-metals, such as human tissue, because the energy dependence of the absorption coefficient is not high. However, X-ray absorption coefficients of high-Z materials are extremely energy dependent and produce erroneous projection data, which is the source of metal artifacts.

2.1 Development of CT simulation software

CT simulation software was developed to generate metal artifacts for evaluation of the efficacy of our proposed algorithm. A numerical phantom was used as a virtual cross-sectional image in X-ray projection simulation. In this study, numerical phantoms were represented by bitmap images of 16-bit gray scale, where brightness corresponded to actual CT values. To reproduce metal artifact on simulations, forward projection calculations using continuous X-rays assuming actual scan environments are necessary. However, it is difficult to obtain X-ray intensities for each energy contained in a continuous X-ray spectrum. Thus, assuming a continuous source of X-rays, we defined eight discrete energy bins (10, 20, 30, 40, 50, 60, 80, and 100 keV) (Table 1). This discrete energy spectrum can be described by the Birch–Marshall equation [17]. X-ray attenuation coefficients for each element were calculated from mass attenuation coefficients found in the National Institutes of Standards and Technology (NIST) database [18]. Note that the most important factor in generating discrepancies that cause metal artifacts is the simultaneous irradiation of an object with both low- and high-energy X-rays, and it is generally accepted that a rough discretization of irradiation energy is sufficient.

2.2 Simulation experiments

We irradiated the Shepp–Logan phantom (Fig. 3(a)) with discrete X-rays from 800 directions using the simulation software developed in this study. The FBP result is shown in Fig. 3(b). Although some beam hardening artifacts were observed, the reconstructed image is almost identical to the original image. Figure 3(c) shows the Shepp–Logan phantom with the locations of six simulated iron pieces indicated, and Fig. 3(d) is the FBP result of Fig. 3(c). The results of FBP on an image containing metal clearly show the artifact and image quality degradation. More specifically, these results imply that metal artifacts are energy-dependent as described by the attenuation coefficients, validating our simulation methods.

3 Method

To reduce metal artifacts, we focused on iterative reconstruction algorithms. Commonly, iterative reconstruction includes metal artifact effects during forward projection calculations that do not consider the energy spectrum. To solve this calculation discrepancy, we propose the following iterative method, accounting for energy information:

Prepare an original reconstructed image μ₀ (x, y) and an initial image for iterative reconstruction μ′ (x, y).

Perform forward projection that includes energy dependence on μ′ (x, y), generating projection data p′ (X, θ).

p' (X, θ) = - \ln {\frac{\int_{E} I_{in} \exp (- \int_{- \infty}^{\infty} μ' (x, y) d Y) E d E}{\int_{E} I_{in} E d E}} .

(3)

Perform filtered back projection to p′ (X, θ), and generate the reconstructed image μ″ (x, y).

$μ^{″} (x, y) = \frac{1}{2 π} \int_{0}^{2 π} p^{'} (X, θ) d θ .$ (4)

Compare the reconstructed image μ″ (x, y) with the original image μ₀ (x, y), and calculate the differences dμ.

$d μ = μ_{o} (x, y) - μ^{″} (x, y) .$ (5)

Divide the dμ by a coefficient α, and add it to the initial image μ′ (x, y).

$μ^{'} (x, y) = μ^{'} (x, y) + \frac{d μ}{α} .$ (6)

Repeat steps 2–5 until the feedback value is sufficiently small.

The equation for the iterative calculation is:

$μ^{k + 1} = μ^{k} + \frac{1}{α} [μ_{o} + \frac{1}{2 π} \int_{0}^{2 π} In {\frac{\int_{E} I_{in} exp (- \int_{- \infty}^{\infty} μ^{k} dY) EdE}{\int_{E} I_{in} EdE}} d θ] .$ (7)

Following this process, there is no need to prepare projection data; our method can be applied as long as there is a reconstructed image. The primary goal of this method is to reproduce accurate CT images with metal artifacts. If the proposed method generates accurate metal artifacts equivalent to those of an actual CT image on the simulation, the numerical phantom can be considered equal to the actual cross-sectional image.

3.1 Assumption of X-ray attenuation coefficients

For the proposed algorithm to be fully utilized, CT values (and thus X-ray attenuation coefficients) of all the materials in an image must be determined. In this study, we prepared a distribution of X-ray attenuation coefficients for 5000 materials in the range of – 1000 to 4000 HU. The distribution was calculated from X-ray attenuation coefficients of six materials shown in Table 1, using linear interpolation. If a CT value was below – 1000 or exceeded 4000, the CT value would be modified to – 1000 or 4000, respectively.

4 Application to simulation data

In this section, the simulation experiments performed to confirm the efficacy of the proposed algorithm are described.

4.1 Simulation conditions

Shepp–Logan phantoms embedded with metal targets were prepared. Here, we supposed that the kind of metals on the phantoms was only iron. The image size was 512×512, and the field of view was 300 mm. In all phantom images and reconstructed images, the window width was 1500, and the window level was 300. The FBP images, used as the original images of our algorithm, were rounded down to 14 bit assume actual data format of CT images. The algorithm was applied until the feedback coefficient α was 10.

As an indicator of quantitative performance evaluation, we use the root mean squared error (RMSE), which is expressed as: $\sqrt{\frac{1}{n} \sum_{i = 1}^{n} {({\hat{μ}}_{i} - μ_{i})}^{2}},$ (8) where ${\hat{μ}}_{i}$ is a predicted value, and μ_i is an observed value. In this study, two types of evaluation values were used; RMSE and RMSE_m. RMSE indicates the error between a phantom and a reconstructed image, while RMSE_m indicates the error between an original reconstructed image and a reconstructed image in the iterative calculation process.

4.2 Shepp–Logan phantom with six metal pieces

The result of applying the proposed algorithm to the Shepp–Logan phantom with six metal objects is shown in Fig. 4. The artifacts generated from metal objects on the conventional FBP CT image made the cross-sectional shape unclear. Upon applying the proposed algorithm, it is clear that the metal artifacts are reduced. Owing to the nature of the proposed algorithm, regions with intense differences between μ_o and μ″ are preferentially reflected. In the early stages of application, coefficients were extracted from the approximate locations of metal regions and other regions that had high X-ray attenuation. Then, other regions that had low X-ray attenuation coefficients were gradually modified to generate the same image as the original one. As seen in Fig. 4, metal artifacts were still observed at iteration 30 but were significantly reduced at iteration 1000. Corresponding changes in the RMSE are shown in Fig. 5. It was confirmed that the RMSE decreases in accordance with a fractional function. From this result, it was expected that the RMSE value would continue to decrease as the iteration number increases, but the reduction rate would also decrease. The number of iterations can be adjusted depending on the required image quality.

The results of varying the feedback coefficient α are shown in Fig. 6. At α = 2, the results are of generally good appearance; however, some streaks of metal artifacts remain, and the image is noisy overall. Changes in RMSE at α = 2 are shown in Fig. 7. The reduction ratio of the RMSE is high; however, the value does not continually decrease. In contrast, the reduction ratio of the RMSE increased as the feedback coefficient decreases; however, the value did not continue to decrease. This is thought to occur because it is difficult to reflect small differences of CT values when the feedback value is large. At α = 100, the risk of diverging values decreases, but metal artifacts again remain at iteration 1000. This image and the RMSE value are quite similar to the case at α = 10, and iteration 100. This implies that 10,000 times the iterations would be needed to generate an image is similar to the image at α = 10, resulting in 1000 iterations. To efficiently produce reconstructed images, it is necessary to select an appropriate coefficient α.

4.3 Shepp–Logan phantom with a lot of metal pieces

To verify the performance of our algorithm, another Shepp–Logan phantom including numerous metal pieces including complex-shaped pieces was modeled. Figure 8(b) shows the very large metal artifacts observed in the FBP image. By applying the proposed algorithm, metal artifacts were gradually reduced, similar to the first phantom case. The image at iteration 3000 shows the metal and organ shapes clearly, although some metal artifacts remain especially in regions with sharp-angled pieces.

5 Discussion

5.1 Image quality

For the phantom cases studied, the proposed algorithm effectively reduced metal artifacts. Although the efficacy was demonstrated both qualitatively and quantitatively, it was also confirmed that the performance of the algorithm was dependent on the feedback coefficient α. When α was too small, the reconstructed image tended to diverge. When α was too large, the convergence speed increased dramatically. This study confirmed that the resulting images were essentially the same in the range of 5 ≤ α ≤ 50, indicating that there is almost no need to select a specific value for α for cases within this range. However, in cases including a lot of metals or complex shaped metals, information of attenuation coefficients in some regions disappears and it becomes hard to converge the calculation and remove remained metal artifacts, even if a large α value is used.

5.2 Comparison with another metal artifact reduction

The comparisons of the NMAR (normalized metal artifact reduction) and our proposed method are shown in Fig. 9 and Fig. 10. Note that the NMAR requires raw projection data. When the number of metal pieces was small (Fig. 9), both methods reduced metal artifacts effectively. Other streaking artifacts were generated in the NMAR result, while metal artifacts remained slightly in our method. When the number of metal pieces was large (Fig. 10), the NMAR generated a lot of streaking artifacts and could not reduce metal artifacts sufficiently. The reason seems due to the difficulty to extract metal regions by thresholding. On the other hand, although metal artifacts remained in our proposed method, the result is qualitatively superior to the NMAR. Furthermore, you can see that the RMSE values of the proposed method also superior to those of the NMAR. However, the kinds of remained artifacts are different to each other. By combining these approaches, the effectiveness could be further improved.

5.3 Calculation time

The time required for reconstructions is an important consideration for clinical and industrial applications. All measurements were performed with the same PC, and the average time was calculated from 1000 trials. The conventional FBP time in this study took an average of 1.60 seconds. In contrast, the calculation time for a single iteration of the proposed algorithm was an average of 20.6 seconds. This means that 1000 iterations would require approximately 5.72 hours. This length of time is obviously prohibitive for widespread use of the proposed algorithm. Two methods to reduce calculation time have been considered; general-purpose computing on graphics processing units (GPGPU) and dynamically changing the feedback coefficient α. Because the compatibility between projection calculations and parallel processing is good, it is expected that the calculation time will be dramatically improved with GPGPU. In contrast with this study, which used a fixed α for each experiment, being able to actively select an optimal coefficient during the iterative calculation will reduce the number of iterations needed to achieve an acceptable result.

5.4 Application to actual data

All experiments in this study were based on simulated geometries and experimental parameters. Specifically, the incident X-ray intensity distribution was assumed to be known. To experimentally verify the proposed method, the actual CT X-ray source spectrum and material properties need to be known. Future work could utilize standard materials with known X-ray properties and measurement of the X-ray spectrum. In addition, the proposed algorithm could be affected by noise in the system and will need to be evaluated.

6 Conclusion

We have developed a novel CT image reconstruction algorithm that dramatically reduces the metal artifacts characteristic of FBP reconstruction algorithms on simulation. Discrepancies that cause metal artifacts were eliminated by using iterative reconstruction based on energy information. Because the iterative calculation proposed in this study only requires reconstructed images and projection conditions, it has great potential for widespread use. Once CT projection conditions are digitized accurately, the proposed algorithm can be directly applied to preexisting CT data.

Footnotes

Acknowledgments

We wish to acknowledge valuable discussions with Mr. Kazuo Kikuchi and Mr. Takashi Tomitsuka of Comscantecno Co., Ltd, and Dr. Tomohiro Iwashita of White Rabbit Co.

References

Alvarez

and Macovski

, Energy-selective reconstructions in X-ray computerized tomography, Physics in Medicine and Biology 21(5) (1976), 733–744.

Coleman

and Sinclair

, A beam-hardening correction using dual-energy computed tomography, Physics in Medicine and Biology 30(11) (1985), 1251–1256.

Kalender

W.A.

, Perman

, Vetter

and Klotz

, Evaluation of a prototype dual-energy computed tomographic apparatus. I. Phantom studies, Medical Physics 13(3) (1986), 334–339.

Stenner

, Berkus

and Kachelriess

, Empirical dual energy calibration (EDEC) for cone-beam computed tomography, Medical Physics 34(9) (2007), 3630–3641.

Bamberg

, Dierks

, Nikolaou

, Reiser

M.F.

, Becker

C.R.

and Johnson

T.R.

, Metal artifact reduction by dual energy computed tomography using monoenergetic extrapolation, European Radiology 21(7) (2011), 1424–1429.

, Leng

and McCollough

C.H.

, Dual-Energy CT–Based Monochromatic Imaging, American Journal of Roentgenology 199(5) (2012), 9–15.

and Zhao

, A practical material decomposition method for x-ray dual spectral computed tomography, Journal of X-ray Science and Technology 24(3) (2016), 407–425.

Meyer

, Raupach

, Lell

, Schmidt

and Kachelrieß

, ß, Normalized metal artifact reduction (NMAR) in computed tomography, Medical Physics 37(10) (2010), 5482–5493.

Zhao

, Bae

K.T.

, Whiting

and Wang

, A wavelet method for metal artifact reduction with multiple metallic objects in the field of view, Journal of X-ray Science and Technology 10(1-2) (2002), 67–76.

10.

Zhang

, Zhang

, Zhu

X.R.

, Lee

A.K.

, Chambers

and Dong

, Reducing metal artifacts in cone-beam CT images by preprocessing projection data, International Journal of Radiation Oncology Biology Physics 67(3) (2007), 924–932.

11.

Yazdi

and Beaulieu

, A novel approach for reducing metal artifacts due to metallic dental implants, IEEE Nuclear Science Symposium Conference Record 4 (2006), 2260–2263.

12.

Abdoli

, Ay

M.R.

, Ahmadian

, Dierckx

R.A.

and Zaidi

, Reduction of dental filling metallic artifacts in CT-based attenuation correction of PET data using weighted virtual sinograms optimized by a genetic algorithm, Medical Physics 37(12) (2010), 6166–6177.

13.

De Man

, Nuyts

, Dupont

, Marchal

and Suetens

, An iterative maximum-likelihood polychromatic algorithm for CT, IEEE Transactions on Medical Imaging 20(10) (2001), 999–1008.

14.

Wang

, Frei

and Vannier

M.W.

, Fast iterative algorithm for metal artifact reduction in X-ray CT, Acad Radiol 7 (2000), 607–614.

15.

Bal

and Spies

, Metal artifact reduction in CT using tissue-class modeling and adaptive prefiltering, Medical Physics 33(8) (2006), 2852–2859.

16.

Chen

, Li

, Guo

, Hu

, Luo

, Yin

, Gu

and Toumoulin

, CT metal artifact reduction method based on improved image segmentation and sinogram in-painting, Mathematical Problems in Engineering 2012 (2012), 786281.

17.

Birch

and Marshall

, Computation of Bremsstrahlung x-ray spectra and comparison with spectra measured with a Ge-Li detector, Physics in Medicine and Biology 24(3) (1979), 505–517.

18.

Hubbell

J.H.

and Seltzer

S.M.

, Tables of X-Ray Mass Attenuation Coefficients and Mass Energy-Absorption Coefficients from 1keV to 20MeV for Elements Z=1 to 92 and 48 Additional Substances of Dosimetric Interest, NIST, http://www.nist.gov/pml/data/ (2015).