Reduce beam hardening artifacts of polychromatic X-ray computed tomography by an iterative approximation approach

1 Introduction

Today, many clinical diagnoses and treatments depend on CT images. The metallic implants or other high density artificial prostheses usually degrade CT quality in form of metal artifacts [1–4]. The beam hardening artifact is one of most important modalities of metal artifacts, which is very difficult to be removed by general approaches for improving CT quality such as higher dose scans [5–9]. When the effect of beam hardening has not been suppressed efficiently, for example, FBP algorithm is employed directly, the beam hardening artifacts will become very serious in the form of dark and bright streaking artifacts in reconstructed images, which impair the image quality and thereby hamper the correct diagnosis.

Since more and more artificial prostheses and metallic implants are used along with an aging society is drawing near, many researchers have made great efforts to reduce metal artifacts [10–14]. The major challenge for metal artifact reduction (MAR) roots in the highly-nonlinear variation of attenuation coefficients of metallic objects along with beam energy. The proposed approaches include dual energy correction [15], projection completion and statistical iterative correction schemes (iterative image reconstruction algorithms). The latter two do not need extra dose of radiation and become more common choices. In the projection completion approach, the missing projections are synthesized from neighboring projections using linear or polynomial interpolations [16–20], Poisson inpainting [21], wavelet interpolation [22, 23], adaptive filtering, tissue-class model techniques [24–27]. The performance of these algorithms depends on the redundancy of available projection data and/or a priori information about missing projections. Recently, the iterative model-based image reconstruction methods are proposed to estimate missing or incomplete data using statistical methods [28–31]. However, numerous study results show these approaches cannot solve the problem very efficiently. For example, these algorithms cannot entirely compensate for severe projection corruptions, or might blur the image details, or induce new artifacts sometimes, and some need huge calculation power, or need other extra techniques [32].

In this paper, the idea behind the proposed approach is to work out the synthetical geometry projections from the (original) projections. The synthetical geometry projections are the nearly linear weighted sums of element geometry projections with respect to mean of each attenuation coefficient, where the element geometry projections are the theoretical projections that are produced only by density distribution of certain element(component) in underlying images. According to Lambert-Beer law, the projections can be expressed as monotonic nonlinear functions of synthetical geometry projections. Newton-Raphson algorithm is employed to work out the solutions, the synthetical geometry projections, from the projections. In the process, Newton-Raphson formula is modified to make Newton-Raphson function satisfy the convergence conditions of FPI so that the solutions are fine approximations of the synthetical geometry projections. In theory, since the element geometry projections or their weighted sums, the synthetical ones, are free or almost free from beam-hardening effect, the beam-hardening artifacts would be suppressed efficiently.

2 Method: Obtaining the synthetical geometry projection from original projection

According to Lambert-Beer law, for X-ray source at energy level E, I ₀ (E), the X-ray beam intensity passing through a two-dimensional object along a line L_φ,s : = { x  = (x₁, x₂) : x₁ cos φ + x₂ sin φ = s} can be described by following formula I ( E , φ , s ) = I 0 ( E ) exp { - R f ( E , φ , s ) }(1) where R f ( E , φ , s ) = ∫ L φ , s f ( E , x ) d l and f (E, x ) denotes the attenuation coefficient at position x and energy level E. If the attenuation coefficient is constant, for example, the object consists of only one component(element) with attenuation coefficient at energy level E, the integration can be simplified as R f ( E , φ , s ) = α E R f ( φ , s ) = α E ∫ L φ , s f ( x ) d l, where f ( x ) is the grayscale image. For the polychromatic X-ray source and multicomponent object, the normalized projection becomes I ( φ , s ) = ∫ I ( E , φ , s ) dE ∫ I 0 ( E ) dE = ∫ η ( E ) exp { - R f ( φ , s ) } dE(2) where η (E) =  I ₀ (E)/∫ I ₀ (E′) dE′ represents the fractional energy of photon energy E in the spectrum of X-ray source, or the normalized distribution of energy spectrum of X-ray beam.

In order to simplify analysis calculation, the equation (2) has to be discretized. Suppose a polychromatic X-ray source consists of K energy levels (marked by i = 1, ⋯ , K), and an object consists of M components (metals and other organic tissues, marked by j = 1, ⋯ , M), which occupy certain regions in the underlying image. Let f _j denote an image that only consists of j-th component, i.e, if a pixel (n, m) is of j-th component, f _j  (n, m) =1, f _j  (n, m) =0 otherwise. For example, let f₁ denote an image that only consists of copper, f₂ an image that only consists of calcium. In this paper, R f j ( φ , s ) is named as the element geometry projection of component j, which is the theoretical projection and produced by the pixel intensities (component density distribution) only. Let α_i,j denotes the attenuation coefficient of j-th component at i-th energy level, the normalized X-ray beam intensity passing through a two-dimensional object can be approximately expressed as I ( φ , s ) = ∑ i = 1 K η i exp { - ∑ j = 1 M α i , j R f j ( φ , s ) }(3) where η _i  = η (E _i ).

Generally, the metallic objects consists of few elements, such as only titanium or calcium or copper, which are all known. Except the large difference of the attenuation capacities between organic tissues and metals (or other high attenuation components), the variation pattern of attenuation properties are similar: the higher the X-ray beam energy, the smaller the attenuation coefficient. The attenuation properties of metals in common use, such as titanium, calcium, iron and copper are rather similar to each other, and are much bigger than that of organic components such as water [33], which are depicted in Fig. 1. Although the attenuation coefficients per unit volume among the metals make a huge difference, the variation rule are rather same, which can be demonstrated by Fig. 2. In Fig. 2, when the attenuation coefficients are multiplied by proper numbers, they vary similarly.

OPEN IN VIEWER

Fig.1

The attenuation coefficients of titanium, calcium, iron, copper and water.

OPEN IN VIEWER

Fig.2

The attenuation coefficients of titanium, calcium, iron, copper, where are multiplied by some numbers.

The mechanisms of beam hardening artifact for metals and organic tissues are similar except the huge difference in magnitude. For organic tissues, the variations of attenuation coefficients are slight and can be ignored, however, for metallic objects, the variations cannot be ignored, or serious beam hardening artifacts will arise. For the convenience of further analysis, an averaged or representative attenuation coefficient takes the place of several attenuation coefficients, i.e., the equation (3) is simplified as the followings I ( φ , s ) = ∑ i = 1 K η i exp { - α i R f ( φ , s ) }(4) where R f ( φ , s ) is named as the synthetical geometry projection relative to element geometry projection R f j ( φ , s ); α _i is the averaged or the representative attenuation coefficient of metals.

It is reasonable to think that the beam hardening effect in R f ( φ , s ) will be much smaller than that in projections P (φ, s) (i.e. - ln  I  (φ, s)). Therefore, an effective approach for reducing beam hardening effect is to work out R f ( φ , s ) from I _d (the practical version of I  (φ, s)) first, and then reconstruct the underlying image by it. In this paper, Newton-Raphson algorithm is employed to work out R f ( φ , s ) from I _d . The algorithm is an iterative operation for finding successively better approximations to the root of a real-valued function from initial approximation, which is rewritten as following [34].

Newton-Raphson Theorem. Assume that f ∈ C² [a, b] (the function set whose second derivative is continuous in [a, b]) and there exists a number p ∈ [a, b] so that f (p) =0. If f′ (p) ≠0, then there exists a δ > 0 such as the sequence { p k } k = 0 ∞ defined by the iteration p k = p k - 1 - f ( p k - 1 ) f ′ ( p k - 1 ) , for k = 1 , ⋯(5) will converge to p for any initial approximation p₀ ∈ [p - δ, p + δ].

Obviously, Newton-Raphson Theorem can be extended to high dimension problems when variables are independent of each other. Since the elements of R f ( φ , s ) in the equation (4) are independent each other, Newton-Raphson algorithm is available for these questions.

For convenience, let x = R f ( φ , s ), the equation (4) becomes I ( x ) = ∑ i = 1 K η i exp { - α i x }(6)

According to Newton-Raphson Theorem, the convergence depends on the initial approximation p₀ and radius of convergence δ. The former should be selected carefully so that the initial approximation is near the convergence center p. The latter depends on the properties of function f (x). According to our experience on using Newton-Raphson Theorem, it is quite difficult to select the initial approximations so that the proper solutions are obtained finally. In order to extend the convergence region, i.e, enlarge δ to ensure the proper solution p is obtained finally, a skill is necessary. According to FPI theorem [34], for Newton-Raphson iteration function g (x) = x - f (x)/f′ (x), if |g′ (x) |<1, i.e. satisfy the convergence conditions of FPI, δ becomes sufficiently large. Consider (6), Newton-Raphson iteration function can be rewritten as | g ′ ( x ) | = | ( ( x − I ( x ) − I d I ′ ( x ) ) ′ | = | ( ∑ i = 1 K η i exp { − α i x } − I d ) ( ∑ i = 1 K a i 2 η i exp { − α i x } ) | ( ∑ i = 1 K a i η i exp { − α i x } ) 2(7) where the division in the right side is the element-by-element division. Generally, |g′ (x) | < 1 cannot be satisfied for all the coefficients α _i , η _i and x . In order to enlarge the convergence region, let x _c  = c x with a large positive constant c > max {α _i }. The equation (6) becomes (8) I ( x c ) = ∑ i = 1 K η i exp { - α i x c c } = ∑ i = 1 K η i exp { - b i x c } = I d(8) where b _i  = α _i /c < 1. Let b _m  = max {b _i }, then | g ′ ( x c ) | = | ∑ i = 1 K η i exp { - b i x c } - I d | · ( ∑ i = 1 K b i 2 η i exp { - b i x c } ) ( ∑ i = 1 K b i η i exp { - b i x c } ) 2 < b m | ∑ i = 1 K η i exp { - b i x c } - I d |(9)

For the initial x _c , suppose max | ∑ i = 1 K η i exp { - b i x c } - I d | = ɛ, there always exists c (or b _m ) so that for any initial x _c it results |g′ ( x _c ) | < 1. It means the solution of (8), R f ( φ , s ), can be finally obtained by an iterative approach from projection I _d .

Consider the similarity of attenuation coefficients among the metals and the large difference of attenuation coefficients in magnitude between organic and metals, we think the error caused by the simplification in (4) will be limited. For example, there are four different metallic objects and organic tissue in the second example of Example 1 in the next section, the simulation results show the normalized error between two I  (φ, s) calculated using five different attenuation coefficients (four metals and one organic tissue) and using the mean of four metal attenuation coefficients is smaller than 0.01%. The normalized error between the synthetical geometry projection R f ( φ , s ) and linear weighted sum of four element geometry projections R f j ( φ , s ) with respect to attenuation coefficient means is smaller than 1.7%, where the means are the average values of attenuation coefficient with respect to X-ray beam energy distribution. So, the synthetical geometry projection can be thought to be nearly linear weighted sum of element geometry projections with respect to attenuation coefficient means.

When it is reconstructed using the synthetical geometry projection R f ( φ , s ), the underlying image f should be adjusted to obtain proper CT numbers. Suppose the grayscale of an image is in [0, 1], the CT number of pixel with intensity p _i will be transformed as C T n u m b e r = p i a m e t a l a H 2 O · 1000 − 1000(10) where a_metal and a_H2O are the attenuation coefficient means of the metal selected (or an average or representative one) and water. Newton-Raphson algorithm forms an iterative approximation approach for MAR in the sinogram domain, which can be summarized by the flowchart in Fig. 3, and is named as Newton-Raphson Iteration Technique (NRIT) in this paper.

OPEN IN VIEWER

Fig.3

The flowchart for calculating the synthetical geometry projection from polychromatic projection.

The performance of the proposed approach is evaluated by two numerical simulation and one practical experiment examples. All examples are achieved on a Windows machine with Intel Xeon E5-2620 2.1GHz CPU (2 processors) and 32GB RAM.

Example 1. In the first part of this example, the original image is the head phantom with 512 × 512 pixels, whose size is suppose to be 128 × 128mm and component is supposed to be brain (grey/white matter, ICRU-44). The image is added four circles whose component is supposed to be copper with mass density 8.9gram/cm³. The image is shown in Fig. 4 (a). In the simulation of projection, the two components are performed parallel projection transform to produce two versions of element geometry projections separately, where the angle interval is selected as 0 . 1 ^o . The original projection is calculated by equation (3), where η _i is selected as the simulation X-ray spectra distribution of Siemens tube with the peak voltage 80kev [35], α _i is selected as the mass attenuation coefficient supported by National Institute of Standards and Technology (NIST) [33].

OPEN IN VIEWER

Fig.4

The original and reconstructed images. (a) The original image, which is the head phantom image together with four circles(copper), (b) and (c) are the images reconstructed by FBP and the proposed algorithm, respectively(WW/WL=200HU/100HU).

Since the mean of mass attenuation coefficient of copper in 8 ∼ 80kev (the valid spectra range of Siemens tube) is 32.78 times of H₂O and the mass density of copper is 8.9 times of that of H₂O, the intensity of copper is set as 291.7 times of H₂O, so does the CT number. The synthetical geometry projection is calculated by the flowchart in Fig. 3.

When the image is reconstructed by FBP directly, the reconstructed image is shown in Fig. 4 (b). When the proposed approach in Fig. 3 is employed, the reconstructed images are shown in Fig. 4 (c), where the stop criterion is selected as ε = 1.0E - 8.

The normalized root mean square difference (NRMSD) and mean absolute deviation(MAD) are employed as criteria in accessing the correction performance [4, 22]. The NRMSD in certain region between the reconstructed and reference image (the original image) is computed according to following formula NRMSD = 1 x max - x min 1 N ∑ i = 1 N ( x r , i - x c , i ) 2 × 100 % where N is the pixel number in the region selected, x _max and x _min are the maximum and minimum intensities in the region of reference image, respectively; x_r,i and x_c,i are the intensities of the corresponding pixels in the reconstructed and reference images, respectively. Unlike NRMSD, MAD is defined as the mean absolute difference between the reconstructed and reference image in the term of CT number. In this example, ROI3 is selected as water in calculated the CT numbers. The results are shown in Table 1.

In the second part of this example, the elements of four circles are supposed to be copper, iron(7.8 gram/cm³), calcium(1.5 gram/cm³) and titanium(4.5 gram/cm³), respectively. The image is shown in Fig. 5 (a). At the first, the five components (one tissue and four metals) are performed parallel projection transform to produce five versions of element geometry projection separately. Then the simulation projections are calculated by equation (3). When the image is reconstructed by FBP directly, the reconstructed image is shown in Fig. 5 (b), while the reconstructed image using the proposed approach is shown in Fig 5 (c). In the latter process, the mean of attenuation coefficients of four metals is employed in working out the synthetical geometry projection. The criteria results are also shown in Table 1.

OPEN IN VIEWER

Fig.5

The original and reconstructed images. (a) The original image, which is the head phantom image added with four circles of different metals. The circles at left-top, right-top, left-bottom and right-bottom corners are supposed to be of copper, iron, titanium and calcium, respectively. (b) and (c) are the images reconstructed by FBP and the proposed algorithm, respectively. (WW/WL=200HU/100HU).

Example 2. In this simulation example, the original image shown in Fig. 6 (a) is a jaw phantom of 512 ∗ 512 pixels (128 ∗ 128mm) with six disk-shaped gold objects [36]. In the first step of simulation, the three components (soft tissue, bone and gold) are performed fan-beam projection transform to produce three versions of element geometry projection separately. For the projection, the rotation angular increment and the sensor spacing are all selected as 0 . 1 ^o , the distance from the fan-beam vertex to the rotation center is selected as 1000 pixels. The simulation projections are calculated by equation (3). The reconstructed image is obtained by the proposed approach is shown in Fig. 6 (f). When the image is reconstructed by FBP directly, the reconstructed image is shown in Fig. 6 (b). Several MAR algorithms, LI, NMAR and the approach in [4], are also employed to reduce beam-hardening artifacts, the results are shown in Fig. 6 (c), (d) and (e), respectively.

OPEN IN VIEWER

Fig.6

The reconstructed images using different approaches. (a) the original image, (b), (c), (d), (e) and (f) are the images reconstructed by FBP, LI, NMAR, the approach in [4] and NRIT, respectively. (WW/WL=600HU/200HU).

NRMSD and MAD are also employed in accessing the correction performance. The regions of interest (ROIs) for calculation is shown in Fig. 6 (a) marked as ROI₁ and ROI₂ (excludes the metal regions by threshold segmentation). All results of accession criteria are shown in Table 2. The results demonstrate the proposed approach has satisfactory performance. For example, it shows NRIT has reduce NRMSD 3.49(%) compared to the algorithm in [4], or NRMSD has been reduced by 17.52% compared to [4].

Example 3. In this experiment, the objects are four power cables that are cut and bound together, which is shown in Fig. 7 (a), whose cores are made of copper. The CT image is obtained by Siemens Inveon Micro-CT/PET machine, which is shown in Fig. 7 (b). Owing to the reasons well known, only the reconstructed image is available, and the projections have to be obtained by reprojection. In the process, the synthetical geometry projection is calculated using the reprojection data, and then reconstructed using the same scan settings used by the machine (that can be obtained in the log file). The results are shown in Fig. 7 (c). The CT numbers are calculated according to the mean of attenuation coefficient of copper roughly (since there is not water in this example). The mean of CT numbers at the four regions in Fig. 7 (b) is 1137 HU, while it becomes 1392 HU after MAR in (c). The latter is more near the CT number of PVC (the power cable sheath), 1717HU [37].

OPEN IN VIEWER

Fig.7

The object, the initial and postprocessing images. (a) The object, four power cables bound together, (b) the initial (WW/WL=2000HU/1000HU), (c) the image after NRIT postprocessing. (WW/WL=2000HU/1000HU).

The results show the image quality has been improved effectively, for example, the contour of PVC material become more clear and full. However, it is still not enough satisfactory and need to be improved furthermore. We think the reason is twofold: (1) there might are some extra processing procedures that are hidden by the machine, which make the reprojection data is different from the true original projections; (2) the scatter (might comes with Compton effect) and Poisson noise are also the causes of metal artifacts, which should be represented by some different models other than (8), so do the MAR algorithms. In the future study, we will focus on these questions.

In this study, an iterative approach is proposed for reducing beam hardening artifacts in the sinogram domain. In the approach, the idea of element geometry projections and synthetical geometry projections are introduced. An element geometry projection is the theoretical projection of certain component, the Radon transform of the image that only consists of the component, and the synthetical one is the linear weighted sum of element geometry projections with respect to means of attenuation coefficients in certain energy range. As a result, the projections can be expressed as monotonic nonlinear functions of synthetical geometry projections. With help of prior knowledge on spectrum parameters and attenuation coefficients, the functions have explicit expressions. Newton-Raphson algorithm is employed to work out the synthetical geometry projections. By some techniques, Newton-Raphson functions can satisfy the convergence conditions of FPI, which ensures the solutions will be obtained steadily.

Since the similarity of attenuation coefficients among metals and the large difference of attenuation coefficients in magnitude between organic and metals, the solutions (the approximations of synthetical geometry projections) are nearly linear weighted sums of the element geometry projections with respect to means of attenuation coefficients, which has been demonstrated by the numerical simulation results. The underlying images are reconstructed using the solutions by general algorithms, and then are adjusted by linear magnification to obtain proper CT numbers. Some examples demonstrate the proposed approach is efficient in reducing beam hardening artifacts and has satisfactory performance. When the original projection is absent, the approach can also alleviate metal artifacts somewhat. The essential conditions for the approach consist of the spectrum parameters and energy-dependent attenuation coefficients, which are usually available in the practical applications. Therefore, the approach has the potential applications in practical CT systems, or can be used in postprocessing CT images.