Efficient CT Metal Artifact Reduction Based on Fractional-Order Curvature Diffusion

Abstract

We propose a novel metal artifact reduction method based on a fractional-order curvature driven diffusion model for X-ray computed tomography. Our method treats projection data with metal regions as a damaged image and uses the fractional-order curvature-driven diffusion model to recover the lost information caused by the metal region. The numerical scheme for our method is also analyzed. We use the peak signal-to-noise ratio as a reference measure. The simulation results demonstrate that our method achieves better performance than existing projection interpolation methods, including linear interpolation and total variation.

1. Introduction

Metal artifact reduction (MAR) is still one of the major challenges in X-ray computed tomography (CT) imaging [1–6]. For high density objects, such as a metal object, the severe attenuation of X-rays allows only a limited number of photos to reach the receiving CT array of sensors. As a result, streak artifacts appear in the reconstructed image after filtered back projection (FBP). The artifacts spread through the whole image, thereby contaminating the imaging quality.

Since the projection data of the metal regions are much larger than ordinary tissue projection data sets, we can assume that the segments of the sinogram with the metal projection data sets are dominated by the metal component only. Based on this assumption, we can deal with projection data of metal objects as lost information, using one of the two main categories of methods: projection interpolation methods [7–12] or iterative reconstruction methods [13–17]. After theoretical analysis, MAR methods based on iterative reconstruction have better reconstruction performance than projection interpolation methods, but they often incur high computational costs and are difficult to implement in current CT imaging systems. In this paper, we concentrate on projection interpolation methods for MAR. In 1978, Lewitt and Bates used a Chebyshev polynomial to implement interpolation [7]. Kalender et al. employed linear interpolation (LI) [8], while Crawford added various assistant processes based on it [9]. Also based on linear interpolation, Gu et al. presented a more accurate metal region segmentation method using the differences in neighboring pixels [10]. Zhao et al. proposed interpolating the wavelet coefficients of the projection data [11]. To obtain a better visual effect, inpainting based on partial differential equations (PDEs) was introduced in [12, 13]. The inpainting method proposed by Gu et al. is based on Euler's elastica and curvature [12]. Duan et al. introduced a classical image inpainting method based on total variation (TV) for MAR [13]. However, the existing PDE image inpainting method cannot connect a wide inpainting region smoothly, and therefore, if wide regions exist, it does not achieve satisfactory results.

In this paper, we propose a fractional-order curvature-driven diffusion (FCDD) MAR model based on our generalized image regularization framework [18]. First, we introduce the FCDD inpainting model. Then, we give the numerical scheme for our method. After presenting our simulation, we conclude this paper.

2. Method

The main steps in our algorithm are similar to those in conventional projection interpolation algorithms as shown in Figure 1.

Figure 1 

Flowchart of projection interpolation algorithms.

The main difference between our method and the conventional projection interpolation methods lies in the step “correction for metal projection data.” When the metal damages the original projection data in the form of a gap, conventional methods apply interpolations [7, 8] or PDE inpainting [12, 13] algorithms to restore the data gap. But these conventional methods have some drawbacks. First, the gap boundaries after inpainting lack smoothness. Second, if the gap is wide, the inpainting results do not achieve satisfactory visual effects. The proposed CT MAR method is based on FCDD that can fix these flaws. This method achieves smoother results and can also cope with the artifacts caused by wide data gaps. In this section, we first introduce the well-known TV inpainting model, and then present our FCDD model. Finally, the numerical algorithm is given for our model.

2.1. Review of Classic TV Inpainting Model

Assume a standard image model defined as𝑢(𝑥)=𝑢₀(𝑥)+𝑛(𝑥),(1) where 𝑢₀ is the original image, 𝑛 is additive noise, and 𝑢 is the contaminated image with noise.

Let Ω be the inpainting (open) domain with its boundary 𝜕Ω, and 𝐸 an extended domain surrounding the 𝜕Ω such that 𝜕Ω lies within 𝐸∪Ω. The image inpainting model based on TV proposed by Chan and Shen [19] is followed here:min𝐽_𝜆[𝑢]=_𝐸∪Ω||∇𝑢||𝑑𝑥𝑑𝑦+𝜆^Ω2_𝐸||𝑢−𝑢₀||²𝑑𝑥𝑑𝑦.(2) The first term is the regularizing term for inpainting damaged domains, while the second term in the energy function is the data fidelity term that can keep important features and sharp edges when noise exists. 𝜆_Ω is a scale function to tune the weights of the two terms. According to variational theory, the Euler-Lagrange equation corresponding to (2) is−∇⋅||∇𝑢||⁻¹∇𝑢+𝜆_Ω𝑢−𝑢₀=0(3) with the Neumann boundary condition, 𝜕𝑢/𝜕𝑛=0 on 𝜕Ω, where 𝜆_Ω=𝜆⋅𝟏_𝐸(𝑥,𝑦)=𝜆0,,(𝑥ot,h𝑦e)r∈wi𝐸se.(4) This model is inspired by the classic TV denoising model [20].

2.2. Proposed FCDD Inpainting Model

In the classical TV inpainting model, the conductivity coefficient of the diffusion strength, which only depends on the numerical value of the isophotes, is𝐹=||∇𝑢||⁻¹.(5) The geometric information of the isophotes is not considered, which is why a wide gap cannot be restored perfectly. Motivated by the methods in [21, 22], to recover from this situation, we use a new fractional-order conductivity coefficient instead of the old one:𝐹=𝑓||𝜅^𝛼||||∇^𝛼𝑢||⁻¹,(6) where 𝑓(⋅) is a function with the following property:𝑓(𝑠)=⎧⎪⎪⎨⎪⎪⎩0∞b,e,tween0and∞,0𝑠𝑠==<0∞𝑠<∞.(7) If the isophote has a large fractional-order curvature, the diffusion strength is also large. In this paper, we choose 𝑓(𝑠)=𝑠. Thus, the Euler-Lagrange equation for the FCDD model is𝜕𝑢𝜕𝑡=∇⋅𝑓||𝜅^𝛼𝑢=𝑢||₀||∇^𝛼𝑢||⁻¹∇𝑢,in𝐷,,in𝐷^𝑐.(8) Here, the inpainting domain 𝐷 is a mathematical open set, 𝐷^𝑐 denotes the outer area of 𝐷, and 𝑢₀ is the available part of the image. The fractional-order curvature 𝜅^𝛼 is defined as𝜅^𝛼=∇^𝛼⋅∇^𝛼𝑢||∇^𝛼𝑢||⁽⁻¹⁾.(9)

2.3. Numerical Scheme

In this section, we apply a time marching scheme to our model. Assuming a time step size of Δ𝑡 and a space grid size of ℎ, we let 𝑥𝑡_𝑖=𝑖ℎ,𝑦_𝑗=𝑗ℎ,𝑖,𝑗=0,1,…,𝑁,with𝑁ℎ=1,_𝑛=𝑛Δ𝑡,𝑛=0,1,…(10)

The explicit scheme iterates as𝑢^𝑛+1=𝑢^𝑛+Δ𝑡Δ𝑢||𝜅^𝛼||||∇^𝛼𝑢||⁻¹.(11)

First, we describe the discretization of the fractional-order gradient operator ∇^𝛼 and fractional-order Laplacian operator Δ^𝛼 using the Grümwald-Letnikov definition [23] in fractional calculus as the mask we proposed in [24, 25].

The 𝛼-order Grümwald-Letnikov definition-based fractional differential can be expressed as𝐷^𝛼_G-L𝑠(𝑥)=𝑑^𝛼[𝑑(𝑥−𝑎)]^𝛼𝑠(𝑥)||||_G-L=lim_𝑁→∞((𝑥−Γ(𝑎−)/𝛼𝑁))^−𝛼^𝑁−1_𝑘=0Γ(𝑘−𝛼)×𝑠𝑥−𝑘𝑥−𝑎Γ(𝑘+1)𝑁,(12) where the duration of signal 𝑠(𝑥) is [𝑎,𝑥], 𝛼 is any real number, and 𝑠(𝑥−𝑘((𝑥−𝑎)/𝑁)) is the discrete sampling.

If 𝑁 is big enough with 𝑎=0, the limit symbol can be dispensed with and (12) is rewritten as𝑑^𝛼𝑑𝑥^𝛼𝑠(𝑥)|||_𝐺-𝐿≅𝑥^−𝛼𝑁^𝛼Γ(−𝛼)^𝑁−1_𝑘=0Γ(𝑘−𝛼)Γ(𝑘−1)𝑠×𝑥+2𝛼𝑁𝑥−𝑘𝑁𝑥.(13) To obtain the value of 𝑠(𝑥+𝛼𝑥/2𝑁−𝑘𝑥/𝑁), we use Lagrange 3-point interpolation with 𝑠(𝑥+𝑥/𝑁−𝑘𝑥/𝑁), 𝑠(𝑥−𝑘𝑥/𝑁), and 𝑠(𝑥−𝑥/𝑁−𝑘𝑥/𝑁). Then, we obtain 𝑑^𝛼𝑑𝑥^𝛼𝑠≅(𝑥𝑥)_−𝛼𝑁^𝛼Γ(−𝛼)^𝑁−1_𝑘=0Γ(𝑘−𝛼)×𝑠Γ(𝑘+1)_𝑘+𝛼4𝑠^𝑘−1−𝑠_𝑘+1+𝛼²8𝑠^𝑘−1−2𝑠_𝑘+𝑠_𝑘+1.(14) If 𝑘=𝑛≤𝑁−1, from (13), the anterior 𝑛+2 approximate backward differences of the fractional partial differentials on the negative 𝑥- and 𝑦-axes, respectively, are expressed as 𝜕^𝛼𝑠(𝑥,𝑦)𝜕𝑥^𝛼≅𝛼+4+1𝛼−8²𝛼₂𝑠(𝑥+1,𝑦)×2−𝛼8³×𝑠(𝑥,𝑦)+Γ(1−𝛼)_𝑛−2_𝑘=1⋅Γ1((𝑘𝑘−−+𝛼𝛼1₂)+!1)⋅𝛼4+𝛼8²+Γ(𝑘𝑘−!𝛼)×𝑠(𝑥−𝑘,𝑦)+4+Γ(Γ𝑛(−(𝑘𝑘𝛼−−−𝛼11)−)!1)⋅−𝛼4+𝛼8²+Γ(𝑛(𝑛−−1)𝛼!−Γ(2−)𝛼)⋅1−𝛼4²×𝑠(𝑥−𝑛+1,𝑦()𝑛+−2Γ)(!𝑛Γ−(−𝛼𝛼−)⋅1)−𝛼4+𝛼8²⋅−𝛼(𝑛−1)!Γ(−𝛼)𝜕4+𝛼8²×𝑠(𝑥−𝑛,𝑦),_𝛼𝑠(𝑥,𝑦)𝜕𝑦^𝛼≅𝛼+41+−𝛼8𝛼²₂𝑠(𝑥,𝑦+1)×2−𝛼8³×𝑠(𝑥,𝑦)+Γ(1−𝛼)_𝑛−2_𝑘=1⋅Γ(1(𝑘𝑘−−+𝛼𝛼1₂)+!1)⋅𝛼4+𝛼8²+Γ(𝑘𝑘−!𝛼)×𝑠(𝑥,𝑦−𝑘)+4Γ+(𝑛Γ−((𝑘𝑘𝛼−−−𝛼11)−)!1)⋅−𝛼4+𝛼8²+Γ(𝑛(𝑛−−1)𝛼!−Γ(2−)𝛼)⋅1−𝛼4²×𝑠(𝑥,𝑦−𝑛+1()𝑛+−2Γ)(!𝑛Γ−(−𝛼𝛼−)⋅1)−𝛼4+𝛼8²⋅−𝛼(𝑛−1)!Γ(−𝛼)4+𝛼8²𝑠(𝑥,𝑦−𝑛).(15)

For simplicity, we only use fractional order masks in four directions for calculation, including positive 𝑥- and 𝑦-coordinates and negative 𝑥- and 𝑦-axes. Let 𝐷^𝛼𝑥+, 𝐷^𝛼𝑥−, 𝐷^𝛼𝑦+, and 𝐷^𝛼𝑦− denote the calculations in the four directions, as illustrated in Figure 2.

(a)

(b)

(c)

(d)

Figure 2 

Masks of the four directions: (a) 𝐷𝛼𝑥+, (b) 𝐷𝛼𝑥−, (c) 𝐷𝛼𝑦+, and (d) 𝐷𝛼𝑦−.

The coefficients of the masks in Figure 2 are given below:𝐶_𝑠−1=𝛼𝐶_𝑠04+𝛼8²=1−𝛼²𝐶_𝑠12−𝛼8³=−5𝛼𝐶⋮4+51𝛼6³+1𝛼6⁴_𝑠𝑘=1Γ(−𝛼)+ΓΓ((𝑘(𝑘𝑘−+−𝛼1𝛼)+)!1)⋅𝛼4+𝛼8²+Γ(𝑘𝑘−!𝛼−⋅1)1−𝛼4²𝐶⋮(𝑘−1)!⋅−𝛼4+𝛼8²_𝑠𝑛−2=1Γ(−𝛼)+ΓΓ((𝑛(𝑛𝑛−−−𝛼1𝛼)−−!12))⋅𝛼4+𝛼8²+Γ((𝑛𝑛−−𝛼2)−!3)⋅1−𝛼4²𝐶(𝑛−3)!⋅−𝛼4+𝛼8²_𝑠𝑛−1=Γ(𝑛−𝛼−1)+(𝑛−Γ(1𝑛)−!Γ𝛼(−−𝛼2))⋅1−𝛼4²𝐶_𝑠𝑛(𝑛−2)!Γ(−𝛼)⋅−𝛼4+𝛼8²=Γ(𝑛−𝛼−1)(𝑛−1)!Γ(−𝛼)⋅−𝛼4+𝛼8².(16) Thus, we get a discrete representation of each itemΔ𝑢𝑢_𝑥=𝑢_𝑥+𝑢_𝑦,𝑢=minmod𝑢^𝑐𝑥,minmod2𝑢^𝑏𝑥,2𝑢^𝑓𝑥_𝑦,||∇=minmod𝑢^𝑐𝑦,minmod2𝑢^𝑏𝑦,2𝑢^𝑓𝑦_𝛼,𝑢||⁻¹=1||𝐷^𝛼𝑥+𝜅𝑢||²+||𝐷^𝛼𝑦+𝑢||²+𝜀_𝛼,=∇^𝛼∇^𝛼𝑢=𝐷_𝛼𝑥−||∇^𝛼𝑢||⎛⎜⎜⎜⎜⎝𝐷_𝛼𝑥+𝑢||𝐷^𝛼𝑥++𝐷_𝛼𝑦−𝑢||²+||𝐷^𝛼𝑦+𝑢||²+𝜀⎞⎟⎟⎟⎟⎠⎛⎜⎜⎜⎜⎝𝐷_𝛼𝑦+𝑢||𝐷^𝛼𝑥+𝑢||²+||𝐷^𝛼𝑦+𝑢||²+𝜀⎞⎟⎟⎟⎟⎠(17) where 𝜀 is a small positive number to prevent dividing by zero, the superscript indexes 𝑐, 𝑏, and 𝑓 denote central, backward, and forward differences, respectively, and the minmod function satisfies minmod(𝑥,𝑦)=sign(𝑎)⋅max(0,min(|𝑎|,𝑏⋅sign(𝑎))).(18)

3. Results

In this section, we present the experimental results of our FCDD model compared with the LI and TV models.

As there is no quantitative method to measure the performance of CT MAR [11], we apply the peak-signal to-noise ratio (PSNR), which is commonly used in image inpainting [26], as the available criterionPSNR=10×log₁₀‖‖𝑢2−55𝑢²₀‖‖²².(19)

Here, 𝑢 is the image after inpainting and 𝑢₀ is the original image. The greater the value of the PSNR is, the better is the performance. PSNR is usually used to measure the similarity between an inpainted image and a real image, and as such it is a suitable reference.

Several concrete steps are followed in our method. First, because in this paper we focus on the inpainting algorithm, for simplicity, we only use the threshold method to extract the metal region. It is easy to say that a more accurate segmentation algorithm would enhance the performance of MAR. Second, we locate the corresponding metal regions in the projection data set. Third, we employ the FCDD algorithm to inpaint the metal regions. Finally, we reconstruct the image from the inpainted sinogram and insert the metal regions.

In the first experiment, while dealing with multiple metals, we compared the performance of FCDD against other methods. Five metal regions with much higher attenuation were added into the Shepp-Logan (S-L) phantom (256 × 256) to simulate the metal artifacts. The parameters are given in Table 1.

Table 1 

Parameters of the modified S-L phantom with five metal regions: 𝜏 is the attenuation coefficient; 𝐴, 𝐵 are the lengths of the semiaxes of each ellipse; 𝑥, 𝑦 are the central coordinates of each ellipse; 𝜃 is the rotation angle (in degrees).

Figure 3 shows the results of the different algorithms, where Figure 3(a) gives the original phantom with a metal region while Figure 3(b) shows the phantom reconstructed from the projection data and containing severe metal artifacts. Figures 3(c) to 3(e) illustrate the results of LI, TV, and FCDD, respectively. In these figures, we can see that metal artifacts are suppressed to different degrees. Compared with LI and TV, FCDD achieves a better visual effect. The structure information (edges and shapes) in Figure 3(e) (𝛼=1.8) is the clearest. In particular, the shapes of the virtual organs are almost maintained, except that the TV method causes an artificial effect. Near the metal regions, the staircase effect is marked. The results of FCDD show the best performance with the highest PSNR. To explain the reason for this, we also give the sinograms after inpainting with these three methods in Figure 4.

(a)

(b)

(c)

(d)

(e)

Figure 3 

Comparison of different sinogram inpainting methods with multiple metal regions: (a) phantom with five metal regions, (b) filtered back projection with metal artifacts, (c) filtered back projection after LI (PSNR = 26.4205), (d) filtered back projection after TV (PSNR = 25.9497), and (e) filtered back projection after FCDD (PSNR = 26.4399).

(a)

(b)

(c)

Figure 4 

Sinogram inpainting results for different methods with multiple metal regions: (a) sinogram after LI, (b) sinogram after TV, and (c) sinogram after FCDD.

When there are multiple metal regions, the gaps to be inpainted are much wider than a single metal region and there is much more missing information that needs to be interpolated. The classic interpolation methods only use the information in the same column for interpolation. The useful information is too scant to obtain an accurate result, as shown in Figure 4(a). As the TV and FCDD methods are 2D inpainting methods, they make use of not only the information in the columns, but also that in the rows and thus they obtain better results. However, TV has a flaw in dealing with wide regions, and the staircase effect occurs. The reason for this is that the order of TV is 2 and to obtain an accurate result, the order must be greater than 2 [27] but less than 4 [28]. In this experiment, we set 𝛼=1.8 in FCDD, which means that the order approximates to 3 and thus the visual effects of inpainting sinograms of FCDD are superior.

Figure 5 gives the variation of the PSNR value with the order of the algorithm. If 𝛼=1, the model is the same as in [22]. As the order increases, the PSNR value also increases. If 𝛼≈1.8, the PSNR peaks and subsequently declines.

Figure 5 

Variations in PSNR with fractional order.

Figure 6 shows the clinical case of a patient with metal implants in both femurs. In this experiment, we do not know the original image, so we cannot choose 𝛼 according to the PSNR. Based on the analysis in the previous experiment, we approximately set 𝛼=1.8. Figure 6(a) shows the original image with dark, board streaks radiating from the metals. Figures 6(b) to 6(d) illustrate the reduction in the most severe artifacts using LI, TV, and FCDD, respectively. The dark streaks are not adequately suppressed in Figures 6(b) and 6(c), while fictitious artifacts caused by TV are visible in Figure 6(c). These disturbing artifacts are reduced to a large degree in Figure 6(d). Structural information previously invisible because of artifacts or the incomplete correction becomes visible.

(a)

(b)

(c)

(d)

Figure 6 

Comparison of different methods dealing with a clinical case: (a) original image with metal artifacts, (b) FBP after LI, (c) FBP after TV, and (d) FBP after FCDD.

4. Conclusion

A new MAR algorithm has been proposed based on the classic TV inpainting model. We replaced the conditional conductivity coefficient for TV with a new fractional-order curvature, and in comparison with linear interpolation and TV, our method obtained better quantitative results and visual effects.

To achieve the best performance for different images, the fractional order also needs to be different. Thus, future work will focus on the adaptive selection of the order of the algorithm.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China (60972131). The authors thank the anonymous reviewers for their comments, which have substantially improved the paper.