A metal artifact reduction algorithm in CT using multiple prior images by recursive active contour segmentation

We propose a novel metal artifact reduction (MAR) algorithm for CT images that completes a corrupted sinogram along the metal trace region. When metal implants are located inside a field of view, they create a barrier to the transmitted X-ray beam due to the high attenuation of metals, which significantly degrades the image quality. To fill in the metal trace region efficiently, the proposed algorithm uses multiple prior images with residual error compensation in sinogram space. Multiple prior images are generated by applying a recursive active contour (RAC) segmentation algorithm to the pre-corrected image acquired by MAR with linear interpolation, where the number of prior image is controlled by RAC depending on the object complexity. A sinogram basis is then acquired by forward projection of the prior images. The metal trace region of the original sinogram is replaced by the linearly combined sinogram of the prior images. Then, the additional correction in the metal trace region is performed to compensate the residual errors occurred by non-ideal data acquisition condition. The performance of the proposed MAR algorithm is compared with MAR with linear interpolation and the normalized MAR algorithm using simulated and experimental data. The results show that the proposed algorithm outperforms other MAR algorithms, especially when the object is complex with multiple bone objects.


Introduction
Metal artifacts are one of the most common problems in computed tomography (CT). Various metal implants in the human body such as dental fillings [1], orthopedic implants [2], hip prostheses [3], and implanted marker bins [4], generate dark and bright streaks owing to the high attenuation of metals, which degrade the image quality and diagnostic value of CT images. To reduce metal artifacts, several metal artifact reduction (MAR) algorithms have been developed [5]. These are classified into two groups: projection completion methods and iterative methods. Projection completion methods treat projections through the metal trace region as missing data and estimate them using neighboring measured data by linear [6], higher-order [7][8][9], wavelet [10,11], prior knowledge [12,13], and diffusion [14,15] approximations. Interpolation-based MAR algorithms are simple, but produce additional artifacts when large metal implants are present. To prevent this, MAR algorithms utilizing a prior image have been suggested [16][17][18][19][20]. In these algorithms, the metal trace region is replaced by the forward projection of the prior image obtained by applying thresholding on the pre-corrected CT images [16][17][18][19] with the assumption that pixel values do not vary substantially within materials. Instead of using simple thresholding, Zhang et al. [20] used a TV minimization regularized algebraic reconstruction technique (ART) to find the prior image. In general, MAR algorithms with prior images demonstrate better performance in artifact reduction, but acquiring a high quality prior image is essential.
Iterative methods attempt to iteratively reduce the mismatch between measured data and forward projection of the object [21][22][23][24][25][26][27]. By excluding data in the metal trace region of the sinogram, MAR with iterative methods becomes an ill-posed problem that is stabilized by using a priori knowledge of the image as a regularization term [21][22][23][24]. Wang et al. [21] proposed two iterative deblurring algorithms that used an expectation maximization (EM) formula and ART adapted for MAR. Williamson et al. [22] proposed statistical iterative reconstruction using prior information about the geometry of metal objects and the detector response model. Choi et al. [23] and Zhang et al. [24] proposed iterative MAR based on constrained optimization with total variation (TV) and a quadratic smoothness function as a regularizer, respectively. Verburg and Seco [25] used a regularized iterative method only applied on a high-Z implants traced sinogram. Lemmens et al. [26] used maximum a posteriori (MAP) reconstruction with thresholding-based multimodal prior, and Nasirudin et al. [27] used material decomposed images of spectral CT as priors, and reduced the metal artifacts by regularized maximum likelihood approach. Compared to the projection completion methods, MAR with iterative methods are robust to noise because they are controlled by a regularization term. However, the high computational cost is a practical obstacle.
The proposed MAR algorithm is categorized in the projection completion group. To reduce estimation errors in the metal trace region, we generate multiple prior images by applying a recursive active contour (RAC) segmentation technique to the MAR image with linear interpolation [6], where the number of prior image is controlled by the RAC depending on the object complexity. A sinogram basis is acquired by conducting forward projection of prior images. The metal trace region of the original sinogram is then replaced by the linearly combined sinogram of prior images. When the object is heterogeneous and the data acquisition condition is not ideal, residual errors occur between the original sinogram and the linearly combined sinogram. Additional correction is performed to compensate. We iteratively continue the process until the norm of the residual errors is minimized. As a result, the proposed algorithm finds good prior images and effectively reduces metal artifacts. To evaluate the performance of the proposed MAR algorithm, both numerical and physical experiments are conducted. Quantitative evaluations are performed and compared with the MAR algorithm with linear interpolation (LIN) [6] and normalized MAR (NMAR) algorithm [18].

Materials and methods
In this section, we briefly review the RAC segmentation technique and then describe the proposed MAR algorithm.

Recursive active contour segmentation
RAC segmentation is based on classical active contour segmentation proposed by Chan-Vese [28]. The active contour model uses a level set function to divide an image into two separate regions.
Let I be a given image on O & R 3 such that I : O ! R. The active contour model minimizes the following equation [28]: where H is the Heaviside function and D is the distributional derivative. A function φ is the levelset function which separates an image into two subregions, and c 1 and c 2 are mean values of corresponding subregions; μ is a regularization parameter to control the length of the active contour and segmentation sensitivity. Image segmentation into multiple subregions can be performed by recursively applying the active contour model. Suppose image I 0 has two objects with different intensities as shown in Fig 1(a). By applying the active contour model to I 0 , two objects are separated together from the background by the levelset function. To further segment subregions, the background values in I 0 are replaced by the minimum value of the foreground regions in I 0 . Then a new image I 1 is acquired as shown in Fig 1(b). By continuing the segmentation process until the entire image is filled with a constant as shown in Fig 1(c), multiple subregions are segmented. Note that red lines in Fig 1 indicate a region where the values of the levelset function are zero. The levelset function has positive values in the foreground regions.
The main advantage of the RAC technique is that the number of subregions is determined adaptively during the segmentation process. Although we describe image segmentation using a piecewise linear object, the RAC technique is still applicable to heterogeneous objects by controlling the regularization parameter μ in Eq (1) [29]. Since the regularization parameter controls the segmentation sensitivity [29], small μ is needed for the segmentation of low contrast intensity object, which increases the number of segmented subregions. Fig 2 illustrates the segmentation results depending on the regularization parameter, where the background intensity increases gradually from zero to one toward the right direction.

Proposed method
The proposed MAR algorithm consists of three steps as described in Fig 3. In the first step, the initial image is reconstructed by the Feldkamp, David and Kress (FDK) algorithm using the original sinogram p 0 , and metal regions are segmented by the active contour [28]. Note that the active contour technique is more effective in metal segmentation than the simple thresholding technique when multiple metals are present. Then, the metal trace region, T M , is identified by the forward projection of the segmented metals within the sinogram domain O. The initial corrected image, f k for k = 1, is acquired by applying linear interpolation on the metal trace region T M [6].
In the second step, f k is divided into N subregions, I 1 , . . ., I N by RAC. In our implementation, we chose the regularization parameter in Eq (1) as 0.1, and the number of subregions N is determined accordingly. The segmented subregions are used as prior images, which have the intensity as one; thus, finding the proper coefficient d j for each subregion I j is necessary to satisfy the following equation: where R is the Radon transform operator and b j ¼ RðI j Þ.
The optimal coefficients d j are estimated by solving the following minimization problem: and @T M is the boundary of T M . The first term in Eq (3) minimizes the mismatch between the original sinogram and the linearly combined sinogram on T c M , and the second term regulates the smoothness along the boundary with the regularization parameter α > 0. In this work, the minimization problem is solved using Nelder-Mead simplex direct search [30] with α = 0.1. The mean value of each segmented subregion is used as an initial guess for solving Eq (3). With computed coefficients d j , a new linearly combined sinogramp ¼ P j d j b j is generated. If the segmentation is performed without errors and the coefficients d j are estimated correctly,p is equivalent to p 0 on T c M with ideal data acquisition condition. However, real projectin data are polychromatic and contain some amount of scatter, and thus residual errors (i.e., d ¼ p 0 Àp) still exist on T c M ; thus, using data in the metal trace region ofp is suboptimal for metal artifact reduction. In step three, the residual errors in the metal trace region are compensated by adding a linearly interpolated sinogramd of δ on the metal trace region andp.
The new sinogramp þd coincides with the original sinogram on T c M sinced is equivalent to δ on T c M . Then, k is updated by k + 1, and steps two and three are repeated using the reconstructed image f k ofp þd. This process continues until the norm of the residual errors k d k T c M is less than tolerance τ.
In NMAR, interpolation errors in the metal trace region are minimized by using a normalized sinogram. If the object is piecewise constant and the X-ray spectrum is monochromatic, NMAR and the proposed method produce the same MAR results. When the object is heterogeneous and the X-ray spectrum is polychromatic, the denormalization step of NMAR may amplify the interpolation errors. This effect will be more significant when the metals exist on the surface of the object. In contrast, our approach is more stable since the interpolation errors are minimized by error subtraction. This effect will be described more clearly in our result section.

Experiments
Both numerical and physical experiments were conducted to evaluate the performance of the proposed MAR algorithm.
Numerical simulations. The numerical simulations were performed with a fan beam CT geometry using Shepp-Logan, jaw and abdomen phantoms. A Shepp-Logan phantom (shown in Fig 4(a)) contains two rectangular gold inserts, a jaw phantom (shown in Fig 4(b)) contains three gold implants, and an abdomen phantom (Fig 4(c)) has two gold implants. For each phantom, polychromatic projection data were acquired with a 120 kVp tube voltage using the Siemens X-ray spectrum [31]. Attenuation coefficients of muscle, bone, and gold were obtained from the NIST X-ray Attenuation Database (ICRU-44) [32]. Each detector cell received 10 6 photons without the object, and Poisson noise was simulated. Each phantom consists of 512 × 512 pixels with a 0.2 × 0.2 mm 2 pixel size. To obtain the reference images, noiseless projection data of tissues and bones were reconstructed using a polychromatic X-ray spectrum. The simulation parameters are summarized in Table 1.  Physical experiments. For the physical experiments, a tooth and polymethylmethacrylate (PMMA) phantoms were used. As shown in Fig 5(a), the tooth phantom consists of an orthodontic bracket with one silver clip at the fore tooth and two gold rings at the cheek tooth. A PMMA phantom consists of a 10 cm-diameter cylinder with five holes and a 16 cm diameter extension annulus with four holes (Fig 5(b)). A metal bar was placed at the center of the PMMA phantom, and two screws were placed at the edge of the PMMA phantom. Projection data were acquired from the bench-top cone beam CT (CBCT) system. The bench-top system includes a generator (Indico 100, CPI Communication & Medical Products Division, Georgetown, ON, Canada), a tungsten target X-ray source (Varian G-1592, Varian X-ray Product, Salt Lake City, UT, USA) with a 0.6 × 0.6 mm 2 focal spot, and a 400 × 300 mm 2 flat-panel detector (PaxScan 4030CB, Varian Medical Systems, Salt Lake City, UT, USA) with a tube voltage of 120 kVp and current of 5 mA. The experimental parameters are summarized in Table 2.

Performance comparison
We compared the performance of the proposed algorithm with LIN [6] and NMAR [18]. The performance of NMAR was evaluated using three different prior images: prior image by simple thresholding (NMAR th ), prior image by RAC segmentation of LIN (NMAR seg s ), and the same prior image with the proposed algorithm (NMAR seg e ). In NMAR th , the thresholding value was chosen by the attenuation coefficient of the material (i.e., air = 0, soft tissue = 0.28, and  For the quantitative evaluation, the mean square error (MSE) and standard deviation were used. For the simulation data, the MSE between the reference image and each MAR image was computed using the entire image except the metal region. For the experimental phantoms, the region of interest (ROI) was selected within homogeneous regions, and the standard deviation was calculated for each MAR image. Fig 6(a) shows the MAR results of the Shepp-Logan phantom. Due to significant beam hardening artifacts in the uncorrected image, high quality prior images could not be acquired by simple thresholding (shown in Fig 6(b)). As a result, the NMAR th image contained severe streak artifacts. Using good prior images improved the image quality as shown in the NMAR seg s image, and further improvement was observed in the NMAR seg e image. It is also shown that the proposed MAR image shows similar performance to the NMAR seg e image since the error amplification in the denormalization step of NMAR is not significant because of the negligible beam hardening effect of the soft tissues in the original sinogram. Fig 7(a) shows the MAR results of the jaw phantom and Fig 7(b) compares the prior images acquired by a simple thresholding and RAC technique. Unlike the Shepp-Logan results, the proposed MAR method shows better performance than the NMAR seg e image. Fig 8(a)-8(f) show the sinogram of the reference image, prior images by simple thresholding, the NMAR th image, prior images by the RAC technique, the NMAR seg e image, and the proposed MAR image, respectively. Profiles of the sinogram images (indicated by the magenta line in Fig 8(a)) are also compared. As shown in Fig 8(g) and 8(h), the NMAR th produces significant errors within the metal trace region. Using good prior images reduces the estimation errors in the metal trace region as shown in Fig 8(i). However, residual errors still exist in the profile of NMAR seg e as shown in Fig 8(j). Note that the original sinogram contains beam hardening effects generated from multiple bones; thus, the denormalization step of NMAR amplifies the interpolation errors on the metal trace region. The proposed MAR algorithm reduces the residual errors by the additional residual error correction step (i.e., Step 3 in Fig 3), and thus, further improvement in MAR can be achieved. Fig 9(a) shows the MAR results of the abdomen phantom. Compared to the previous numerical phantoms, the abdomen phantom is more heterogeneous, and thus more accurate prior images would be required to reduce metal artifacts effectively. However, due to the severe noise and metal artifacts in the uncorrected image, prior images generated by simple thresholding miss details of inner structures in the abdomen phantom (shown in Fig 9(b)). As a result, the NMAR th image contains severe streak artifacts. In contrast, prior images generated by the RAC technique contain more details, which improve the performance of the NMAR. While all other MAR algorithms contain residual streaks artifacts between two metals, the proposed algorithm reduces them effectively.

Numerical simulations
In simulations, the number of prior images was chosen by RAC automatically, and the tolerance τ was set to 10% of the first norm of the residual errors. Table 3 summarizes the number of prior image and coefficients for numerical simulations, and Table 4 summarizes the MSE values for different MAR images. For all numerical phantoms, the proposed algorithm shows the best performance in metal artifact reduction.

Physical experiment
For the tooth phantom, three slices from the reconstructed three-dimensional volume were chosen. Fig 10 shows MAR results and prior images of slice 17. Notice that the NMAR th image contains significant streak artifacts owing to the poor quality of the prior images. Although the streak artifacts are reduced in the NMAR seg s and NMAR seg e images, residual streak artifacts are still observed. Since the tooth phantom contains the metals on the surface of the object, small interpolation errors in the metal trace region are significantly amplified during the For the quantitative evaluation, two ROIs were chosen on the homogeneous region, which is indicated by the red lines in Figs 10(a)-12(a). The standard deviation of each ROI was calculated and summarized in Table 5. For the tooth phantom, the proposed MAR algorithm has the smallest standard deviation for all three slices. Fig 13 shows MAR results and prior images of the PMMA phantom with three metal inserts, where the central slice from the reconstructed three-dimensional volume is selected. Notice that the prior images by simple thresholding and the RAC technique are very similar and thus, all NMAR algorithms show similar MAR results with the proposed MAR algorithm.
For the quantitative evaluation, four ROIs were selected on the homogeneous region, which are indicated by red lines in Fig 13(a). The standard deviation of each ROI was calculated and summarized in Table 5. Notice that the three NMAR images and the proposed MAR image show similar standard deviation owing to the simple structure of the PMMA phantom. Table 6 summarizes the number of iterations for each phantom. More iterations are required when the image quality of LIN is poor. Table 7 summarizes the number of prior image and coefficients for each experimental phantom.

Discussion and conclusion
We presented a new MAR algorithm using multiple prior images with subsequential residual error correction. Traditional interpolation-based MAR algorithms replace the metal trace region with neighboring data depending on the distance, which introduces streak artifacts owing to interpolation errors in the metal trace region. The proposed method minimizes interpolation errors by utilizing multiple prior images generated by the RAC technique and additional residual error correction method in sinogram space. The performance of the proposed method was compared with LIN and NMAR with different prior images. For the homogeneous object, NMAR and the proposed method show similar performance in metal artifact reduction. However, when the object is complex with multiple bone objects, the proposed method shows better MAR results. For a complex object, we also showed that the performance of the NMAR algorithm can be improved significantly using prior images generated by the RAC technique.
In NMAR, interpolation errors in the metal trace region were reduced by performing the interpolation on the normalized sinogram. However, using poor quality prior images can degrade the performance of NMAR since the interpolation errors are magnified by the denormalization step of the NMAR algorithm. This effect is more significant when the metals lie on the surface of the objects. Using good prior images improves the performance of the NMAR algorithm significantly. As shown in the results of NMAR seg s , RAC segmentation improves NMAR compared to NMAR th , but the residual errors caused by the denormalization step still exist, especially for complex objects containing multiple bone objects. In contrast, the proposed algorithm does not amplify interpolation errors since the metal trace region is filled in by the summation of linearly combined sinograms of prior images and subsequential residual error correction.     The proposed method iteratively reduces interpolation errors. Each iteration consists of prior image generation, sinogram basis generation, and residual error correction. Each iteration takes approximately 15 s for a two-dimensional image in MATLAB (MathWorks, Natick, MA) on a work station with an Intel XEON E5-2630V2 processor at 2.6GHz. The number of iterations depends on the object's complexity and the quality of the LIN image. Since the proposed algorithm can be applied independently for each slice, the computation time for threedimensional volume data would be similar to the two-dimensional case by parallel computing using a graphics processing unit (GPU). The computation time can be more accelerated by C/ C++ implementation, which is essential for real clinical applications.
In this work, we treated the projection data in the metal trace region as missing or completely unreliable. However, when the density of metals is not high enough, the metal trace region may contain available information of the object, and thus, further improvement of the MAR can be achieved by frequency split MAR (FSMAR) [19]. Comparison of the image quality improvement with FSMAR with NMAR images and proposed MAR images is a subject of future research. In conclusion, we presented a new MAR algorithm using multiple prior images with additional residual error correction. The results show that the proposed algorithm outperforms other MAR algorithms, especially when the object is complex with multiple bone objects.