CT Metal Artifact Reduction Method Based on Improved Image Segmentation and Sinogram In-Painting

The streak artifacts caused by metal implants degrade the image quality and limit the applications of CT imaging. The standard method used to reduce these metallic artifacts often consists of interpolating the missing projection data but the result is often a loss of image quality with additional artifacts in the whole image. This paper proposes a new strategy based on a three-stage process: (cid:2) 1 (cid:3) the application of a large-scale non local means ﬁlter (cid:2) LS-NLM (cid:3) to suppress the noise and enhance the original CT image, (cid:2) 2 (cid:3) the segmentation of metal artifacts and metallic objects using a mutual information maximized segmentation algorithm (cid:2) MIMS (cid:3) , (cid:2) 3 (cid:3) a modiﬁed exemplar-based in-painting technique to restore the corrupted projection data in sinogram. The ﬁnal corrected image is then obtained by merging the segmented metallic object image with the ﬁltered back-projection (cid:2) FBP (cid:3) reconstructed image from the in-painted sinogram. Quantitative and qualitative experiments have been conducted on both a simulated phantom and clinical CT images and a comparative study has been led with Bal’s algorithm that proposed a similar segmentation-based method.


Introduction
Metallic objects like dental implants, surgical clips, or steel-hip prostheses lead to severe shadow and streak artifacts in CT images that superimpose the structures of interest and deteriorate image quality. The reason is that metallic objects have a very high density in the human body, which creates a barrier to the transmitted x-ray beam during CT examination. It results a lack of data in the projection data that lead to the production of streak artifacts in CT images 1, 2 . This photo deficiency caused by metallic object would become more severe under low dose scanning 2 . In the last decade, many approaches have been proposed to reduce these artifacts. These methods can be roughly classified into iterative and interpolation-based methods.
Iterative algorithms operate in a feed-back mode in both the image and projection data spaces 3-5 . However, three major difficulties can be pointed out for that methods: 1 the well-formatted original raw projection data are often unavailable because the leading manufacturers of CT imaging devices are often reluctant to provide it; 2 the involved high computational cost of iterative algorithms often requires an implementation on specialized processor units; 3 the iterative algorithm still need to be combined with sinogram correction method when the metal artifacts are rather severe.
Interpolation-based methods correct metal artifacts directly in sinogram space. Compared to iterative algorithms, these methods are less computationally expensive and can be implemented without the availability of the original raw projection data. They aim at identifying the corrupted segments in the sinogram and interpolating the data from noncorrupted neighboring projections. Some of these methods add other steps to improve the sinogram correction accuracy and design a four-stage process that consists of image enhancement, metallic object segmentation, image forward projection, and sinogram in-painting, final image reconstruction using FBP 6-12 . It is worth notable that the normalization operation suggested in 9, 10 has been considered to give a more accurate attenuation estimation.
We propose a new method to suppress the metal artifacts and improve the sinogram completeness that is based on the above described scheme. The major contribution states at the image enhancement, segmentation, and sinogram impainting levels with, respectively, the application of a large scale nonlocal means filter LS-NLM , a mutual information maximized segmentation MIMS , and a modified exemplar-based in-painting technique. The description of this method is given in Section 2. Comparative experiments with the method proposed by Bal and Spies in 12 are then provided in Section 3.T h i s method was chosen for comparison because it made use of a similar strategy to address the corrupted sinogram problem, that is, Tensor filtering, k-means clustering techniquebased segmentation, and linear interpolation-based sinogram in-painting. It will be referred thereafter in this paper by "Bal's algorithm." We will show that, compared to this algorithm, our method provides a better sinogram correction. Visual and quantitative analysis are also reported to highlight this superiority. A CUDA parallelization technique has been applied to accelerate the calculations of the patch distances involved in the image enhancement and sinogram inpainting steps, respectively.

Method
The proposed sinogram completeness algorithm is divided into four major stages.
Step 1 prefiltering . The original CT image including metal artifacts is enhanced with the edge-preserving LS-NLM filter.
Step 2 image segmentation . The metallic artifacts and objects are respectively extracted using the MIMS algorithm with a partitioning of the image into different regions.

Mathematical Problems in Engineering 3
Step 3 sinogram inpainting . Once the metallic objects and artifacts have been extracted, the segmented artifact image is forward projected to determine the projection data in the sinogram space which are affected by the artifacts. A subtraction is performed between the corrupted sinogram and the original one. The missing projection data in the subtracted sinogram are then restored using a modified exemplar-based in-painting technique.
Step 4 backward projection of the in-painted sinogram and image correction . The artifact compensated image is then reconstructed from the in-painted sinogram using the FBP algorithm. Afterward, the final corrected image is obtained by inserting the previously segmented metal component into the reconstructed image.
Of all the four steps, Step 2 and Step 3 are the two key steps in which the damaged sinogram data are estimated and corrected. The above stages are detailed in the following subsections with the flowchart displayed in Figure 1.

Image Enhancement
This first stage aims at applying an edge-preserving filtering operation to smooth and denoise the streak artifacts in the original CT images. The LS-NLM filter has been proven to be efficient for image denoising with edge preservation. It was, for instance, applied with success to suppress mottled noise in low-dose abdominal CT images 13 . The principle is to replace the value of a pixel by the weighted average of pixels located in a neighbourhood window of size N. Each weight expresses the similarity between the central pixel in the window and each neighboring pixel and is given by the pair-wise difference between patches surrounding each pair of considered pixels 14, 15 .
Let f O and f F , respectively, denote the images before and after filtering, the filter output is given by where f O i is the pixel located at the center of the neighborhood N and f O j the pixels located in the neighborhood of f O i . w ij denotes the weight between pixels f O i and f O j and is calculated as a similarity measure between the two patches s i and s j surrounding each pair of pixels i and j in the neighbourhood N i , respectively. The decay parameterh acts as a filtering parameter. In 2.2 , a Gaussian kernel of standard deviation α is used to take into account the distance between the central pixel and other pixels in the patch. The LS-NLM filter involves working with a large-size neighborhood N and a number of size patches s equal to the number of pixel in the neighborhood N, which implies high costs for calculating the distance between each patch pair in each neighborhood N. To accelerate the computation, a GPU parallelization using the CUDA framework was applied for the original pixel-wise processing based on 16 Figure 1: Flowchart of the proposed correction method.

Image Segmentation
The MIMS method is based on the maximum mutual information MMI and allows determining the class number n b a s e do nt h ed i fference of mutual information DMI 19 where P A , P B are the probability density distribution of A and B, respectively, and P A, B the joint probability density function. P A , P B ,a n dP A, B can be computed from the histogram quantization 20 . Based on 18 ,i n 2.5 ,M I f F ,S n is the normalized mutual entropy between the image f F and the current image S n segmented into n classes, and DMI n f F the normalized difference between the entropies MI f F ,S n and MI f F ,S n−1 .
When the class number increases, DMI n f F decreases while MI f F ,S n converges towards MI f F ,f F . This convergence is reached when DMI becomes smaller than a specified threshold ε. A local optimality can be obtained when DMI converges towards a local minimum while the mutual entropy synchronously reaches its maximum 19 . MIMS implementation only requires to set the maximum class numbers MCN MCN a and MCN m for the segmentations of artifacts and metals, resp. and the threshold ε. A description of the algorithm is given in Figure 2. The filtered image f F is classified into n classes by using an intensity threshold vector G k n k being the kth iteration and n referring to the threshold number . The thresholds in G k n and the class centers can be automatically computed within a simulated annealing-SA-based optimization process 21 .

Sinogram In-Painting
The output of the MIMS algorithm provides the segmented metal artifact and object images. We forward projected the segmented metal artifact image into the sinogram domain, subtracted then the original sinogram built from the original image that includes the metallic objects with the metal artifact sinogram to delete the corrupted projection data, and applied an in-painting technique on the subtracted sinogram to restore the missing projection data. The proposed in-painting technique refers to a modified exemplar-based in-painting method to find the best matched sinogram patch for the restoration of the missing projection data 22, 23 .
Let consider the following hold.
Output S Yes Yes

No
No n ≥MCN Figure 2: Flowchart of the MIMS algorithm: J corresponds to MI n f F ,S n and G k n is the intensity threshold vector for the kth iteration with current class number n. f F refers to the filtered image at the previous stage.
2 A patch P p centered on each pixel p located in the region Ω at the border of the frontier between the two regions Ω and θ so that the patch P p includes a certain number of pixels belonging to θ.
3 A set of patches P q centered on each pixel q of the region θ.
Then the patch P p is compared with each patch P q using the following similarity metric: where the distance d p, q is calculated between the corresponding pixels in patches P p and P q that belong to the region θ.P p and P q are of the same sizes, and · denotes the Euclidean distance between them. The patch P q that minimizes this distance is selected and its contents are copied into P p to restore the missing pixels of P p that are located in the region Ω.A n automatic region filling process is conducted based on 22 that introduces a filling order of the pixels located in the region Ω. It goes though the iteration of a three-stage process for each pixel.
1 The computation of priority for all the patches P p whose central pixels p are located in region Ω but just behind the border line of θ. The priority is computed with the following function Pr p for each of the pixels p: where D p and C p are, respectively, called data and fidelity term. The latter one is introduced to quantify the number of points, surrounding the target pixel p, that are known or have already been in-painted. This term tends to privilege those patches that have more pixels from the known region θ. D p is given by the scalar product between the vector normal to the front and the maximum gradient orientation at point p. Its objective is to encourage linear structures to be processed Mathematical Problems in Engineering 7 first against high curvature line. λ is a weighting coefficient 0 <λ<1 that controls the balance between C p and D p . Details of calculus of C p and D p can be found in 22 .
2 The choice of the pixel p with the highest priority calculated by 2.7 and the search for the best match between the patch P p and the set of patches P qdetermined by the similarity metric 2.6 . Then, the pixels of patch P p located in the region Ω are then restored with the corresponding known pixels of the patch P q .
3 The updates of the fidelity and data terms for pixels of P p that have been filled and are afterward located just behind the border line of θ.
In practice, to increase the accuracy in the in-painting process, the size of patch P p for inpainting is in fact set smaller than when applying the search for the best match in 2.6 .

Experiments
Experiments were conducted both on a simulated phantom and clinical CT images. All the images have a size of 512 × 512. Clinical CT datasets were acquired from a multidetector row CT unit with 16 detector rows Somatom Sensation 16; Siemens Medical Solutions . The scanning protocol was 100 mAs, 120 kVp, 5 mm slice thickness and the spatial resolution was 0.457 mm 2 . The images were reconstructed with a FBP algorithm using a convolution kernel "B40f." Here, convolution kernel is used to control the smoothing effect in CT images for Siemens Somatom Sensation 16 system, and B40f is the routine convolution kernel for brain CT reconstruction. Figures 3 a , 3 b ,a n d3 c display the three original clinical CT images that include metallic artifacts. The first dataset depicts a chest image in which metallic artifacts come from a metallic suture material. The two other datasets refer to a brain image where the metallic artifacts originate from golden earrings. A phantom image including metallic artifacts was also simulated to allow quantitative comparisons. It consists of a cylindrical metallic insert incorporated in a cylindrical water container Figure 4 a .T h e phantom was simulated from an artifact-free phantom CT image Figure 4 b applying the following intensity settings: air: 0; square container: 500; cylindrical water receptacle: 3000, water: 1000. Subsequently, we name the images in Figures 3 a , 3 b , 3 c ,a n d 4 a , ClinicalImage1, ClinicalImage2, ClinicalImage3 and PhantomImage, respectively. We applied a parallel geometry for the forward and backward projection operations involved in the metallic artifact reduction algorithm. This algorithm has been written in C language, MATLAB release R2006b , and NVIDIA CUDA libraries. It was then run on a PC with an Inter Core i5 processor, 2.68 × 4 GHZ, 6 G RAM, and GPU NVIDIA GTX465 . We compared our method with Bal's algorithm in 12 which made use of a similar strategy to reduce the metal artifacts.
1 A linear structure tensor LST filtering proposed in 24 is applied to reduce noise and smooth the artifacts in the original image with artifacts. Three parameters needed to be sets: the mask size υ, the scaling factor σ 0 , and the relationship between width and length of Gaussian filter. a b c Figure 3: Clinical CT images including metal artifacts. a ClinicalImage1 is a chest image with artifacts caused by metallic suture material; b ClinicalImage2 is a brain image with artifacts caused by one golden earring; c ClinicalImage3 is a brain image with artifacts caused by two golden earrings. Note that the image quality is severely degraded by metal artifacts. a b Figure 4: Simulated CT images including metal artifacts. a Phantom image including simulated metal artifact PhantomImage ; b the original artifact-free phantom image used to create PhantomImage. Note that the image quality is also severely degraded by metal artifacts.
3 An in-painting step in which the neighboring sinogram data from the projection of above K-means segmentation is used to complete the tagged metallic projection for the segmented metallic objects.
The proposed process involves also to specify a certain number of parameters: N, s, the decaying parameter h in the prefiltering step, the maximum class numbers MCN a ,MCN m to, respectively, extract the metallic artifacts and object, the threshold ε for the segmentation step, the size of the patches P p , and the value of λ for the in-painting step.
For both algorithms, the involved parameters were manually set to find the optimal parameter combination that led to the best qualitative results. This qualitative evaluation was carried out in collaboration with a radiologist Xindao Yin 10 years clinical experience . These optimal parameters are listed in Table 1. It is notable that we used the same parameter setting in Bal's method because we found similar results can be obtained when using the same parameter settings.

ClinicalImage1
ClinicalImage2 ClinicalImage3 PhantomImage Bal algorithm in 7 Prefilterig step υ 2, the scaling factor σ 0 is set to 10, the relationship between width and length of Gaussian filter is set to 2 σ ⊥ 2σ .  Figures 5 a and 5 b display for comparison each of the steps of each algorithm Bal's algorithm and our proposed method . ClinicalImage1 is used for validation. From Figures 5 a1 and 5 b1 , we can see that the LS-NLM filtering shows better properties of noise suppression and structure preservation than the ASF filtering used in Bal's algorithm. Also, only the metallic object was segmented with Bal's K-means algorithm Figure 5 a2 while our MIMS algorithm allowed both the artifacts and the metallic objects to be extracted Figure 5 b2 . Figures 5 a3 and 5 b3 point up the differences between the two resulting projected corrupted sinogram. In particular, the corrupted sinogram obtained from the segmented metallic artifact component in Figure 5 b2 is larger than the one issues from the segmented metallic component in Figure 5 a2 . Figures 5 a4 and 5 b4 provide the two resulting in-painted sinogram. In the one obtained with Bal's algorithm, we can see metallic shadows resulting from the metallic artifacts segmentation see the red arrow in the zoomed region in Figure 5 a4 . These shadows are absent in Figure 5 b4 ,t h a ti s ,t h e sinogram computed from our method. These metallic shadows create new artifacts in the reconstructed image in Figure 5 a5 . Finally, this example illustrates the best performance of our method on Bal's algorithm.

Quantitative Study
A quantitative study was then carried out on the simulated phantom of Figure 4 a with respect to the original artifact-free phantom image in Figure 4 b . We displayed the intensity profiles along a specified horizontal line in the original and corrected images to highlight the differences in the behavior of each method according to the crossed structure properties. We also computed the mean square error MSE and the standard deviation STD in two ROI located in two different homogeneous regions that is, inside region1 and outside region2 the phantom: The phantom images after correction by Bal's algorithm and our proposed method are, respectively, displayed in Figures 8 a and 8 b . The qualitative evaluation highlights the capacity of our algorithm to better suppress the metallic artifacts as to provide a superior consistency in the homogeneous region preservation. Figure 8 c plots the intensity profiles along the same given horizontal line in the original Figures 4 a and 4 b and corrected Figures 8 a and 8 b images, respectively. The profiles confirm that our method brings a better quality correction in the homogeneous region. Table 2 provides the MSE and STD measures for each method and in each ROI. The figures also confirm the supremacy of our approach with an MSE and an STD that are the lowest on the image set. Table 3 provides the total computation costs in CPU seconds for each step for each method. For our method, the computation cost is given without and with GPU acceleration. This method is rather expensive in computation time. The CUDA parallelization brings a substantial gain with an acceleration that can reach a rate of 10 to 30 depending on the complexity of the image. The parallelization makes then our method competitive in computation time with Bal's algorithm.

Discussion
The proposed strategy adopted for reducing the metallic artifacts in the reconstructed image relies on a four-stage process that consists of image enhancement, metallic object segmentation, image forward projection, and sinogram in-painting, final image reconstruction using FBP. The image enhancement makes use of an LS-NLM filter. This filter exploits a patch similarity measure to smooth the image while preserving the edges. Its response is not very sensitive to the size of the neighborhood N and patch s. However, the decaying parameter h, that quantifies the smoothing rate and how fast the weights decay with increasing dissimilarity of respective patches, is sensitive to the noise ratio in the image. Its value is set as a function of the noise variance. The MIMS-based segmentation involves only to set the maximum class number MCN MCN a and MCN m for the segmentations of artifacts and metals, resp. and the threshold ε that is applied for the convergence of the difference mutual information DMI . The choice of this threshold is not sensitive as we can see in the experiments. Its value is relatively stable on the set of the processed  images Table 3 . The MMI maximum mutual information and DMI are computed within a simulated annealing optimization process. This allows the thresholds in G k n and the class centers to be automatically computed. This guarantees an optimal choice of these parameters. Moreover, the metal artifacts and objects to be extracted have a very high density and if other highly contrasted structures are not located in the close neighborhood, the algorithm run quite well and provides satisfactory results. Considering now the sinogram inpainting stage, the proposed exemplar-based in-painting technique considers a global similarity measure and relies on redundant information present in the image. This modified exemplar-based in-painting method is reasonable for the CT sinogram completion because it is observed there are lots of repetitive structures in CT sinogram. When a large scale is selected, the proposed exemplar-based in-painting can give an effective restoration of repetitive structures in sinogram space. Parameters to be set relate to the sizes of the patches Table 1 and the weighting coefficient λ that is used to balance the fidelity and data term in 2.7 . We preferred to consider small patches 3 × 3 for the inpainting process in order not to lose subtle details.
For comparison, Bal's algorithm used a similar strategy: 1 a linear structure tensor LST filtering was first applied. Its response is highly dependent on the parametric streak strength and orientation quantification and several parameters were empirically tested to ensure a satisfying quantification and efficient filtering. 2 A K-means clustering technique was considered for the metallic artifact segmentation. It required choosing a suitable number of classes and the initial center of each cluster. The K-means algorithm is based on the intensity clustering of one single image and appears finally less efficient than the MIMS method to accurately segment the metallic artifacts. 3 An interpolation technique was then considered to recover missing or metallic data from neighboring projections from noncorrupted segments.
The global process is carried out in a small region surrounding the metallic object. Artifacts around the metal objects can be removed so as to remove projection data inconsistencies with data consistent with similar neighborhood. However, image details may be noticeably altered especially and some remaining artifacts still appear in the regions closest to the metal objects. The presence of pathology has not yet been considered and is something to be further evaluated in collaboration with our medical expert. The reason might be the metallic artefacts can in some situation completely superimpose the structures and the presence of a pathology may be hidden. Thus, the density change due to the artefact removing process may be difficult to evaluate.

Conclusion
This paper proposed a new strategy for reducing metallic artifacts in CT images. The proposed method outperforms Bal's algorithm in each of the three steps: image prefiltering, image segmentation, and sinogram inpainting. Visual and quantitative analyses on phantom and clinical data show that the proposed correction method provides a substantial reduction of the metallic artifacts in the corrected images. The pixel-wise operations in the pre-filtering and sinogram inpainting steps are greatly accelerated by a CUDA parallelization that makes the algorithm also competitive in computation time.
Although this algorithm demonstrated a good potential for the reduction of metallic artifacts in CT images, some improvements have to be further considered. First, segmentation accuracy might be further increased by applying more dedicate method such as the segmentation method in 25 . Second, the exemplar-based in-painting procedure is expensive in computation time due to the search for the patch priorities, which needs to be updated after the in-painting of each point within the corrupted sinogram. Third, some intensity inconsistencies can be still observed around the regions of the metallic objects in the corrected images, and the sinusoid property of sinogram has not be exploited in inpainting the missing sinogram data 26 . At last, the presence of pathology in the surrounding of the metal object has not been considered in the evaluation of the algorithm. Thus, further work will be devoted to solve the set of problems as to perform an extensive evaluation on clinical data.

CT:
Computed tomography LST: Linear structure tensor ASF: Adaptive Steering Filter NAST: Nonlinear anisotropic structure tensor LS-NLM: Large scale nonlocal means MIMS: Mutual information Maximized Segmentation MMI: Maximum mutual information DMI: Difference of mutual information.