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

Authors: Yang Chen; Yinsheng Li; Hong Guo; Yining Hu; Limin Luo; Xindao Yin; Jianping Gu; Christine Toumoulin
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

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: (1) the application of a large-scale non local means filter (LS-NLM) to suppress the noise and enhance the original CT image, (2) the segmentation of metal artifacts and metallic objects using a mutual information maximized segmentation algorithm (MIMS), (3) a modified exemplar-based in-painting technique to restore the corrupted projection data in sinogram. The final corrected image is then obtained by merging the segmented metallic object image with the filtered back-projection (FBP) 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.


---

## Body

## 1. 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 Section2. Comparative experiments with the method proposed by Bal and Spies in [12] are then provided in Section 3. This 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 technique-based 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.


## 2. 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.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 Figure1.Figure 1
Flowchart of the proposed correction method.


### 2.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].LetfO and fF, respectively, denote the images before and after filtering, the filter output is given by
(2.1)fiF=∑j∈NiwijfjO∑j∈Niwij,(2.2)wij=exp⁡(-‖si-sj‖2,α2h),
where fiO is the pixel located at the center of the neighborhood N and fjO the pixels located in the neighborhood of fiO. wij denotes the weight between pixels fiOand fjO and is calculated as a similarity measure between the two patches si and sj surrounding each pair of pixels i and j in the neighbourhood Ni, 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–18].


### 2.2. Image Segmentation

The MIMS method is based on the maximum mutual information (MMI) and allows determining the class numbern based on the difference of mutual information (DMI) [19]. The mutual information MI(A,B) between image A and B is defined by
(2.3)MI(A,B)=H(A)+H(B)-H(A,B),
where H(A) and H(B), respectively, denote the entropies of the images A and B,and H(A,B) the mutual entropy. Based on [15, 16], the MI(A,B) can be rewritten as the joint probability density distribution of images A and B:
(2.4)MI(A,B)=∑P(A,B)log⁡P(A,B)P(A)P(B),(2.5)DMIn(fF)=MI(fF,Sn)-MI(fF,Sn-1)MI(fF,fF),
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), and P(A,B) can be computed from the histogram quantization [20]. Based on [18], in (2.5), MI(fF,Sn) is the normalized mutual entropy between the image fFand the current image Sn segmented into n classes, and DMIn(fF) the normalized difference between the entropies MI(fF,Sn) and MI(fF,Sn-1). When the class number increases, DMIn(fF) decreases while MI(fF,Sn) converges towards MI(fF,fF). 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 (MCNa and MCNm for the segmentations of artifacts and metals, resp.) and the threshold ε. A description of the algorithm is given in Figure 2. The filtered image fF is classified into n classes by using an intensity threshold vector Gnk (k being the kth iteration and n referring to the threshold number). The thresholds in Gnk and the class centers can be automatically computed within a simulated annealing- (SA-) based optimization process [21].Figure 2
Flowchart of the MIMS algorithm:J corresponds to MIn(fF,Sn) and Gnk is the intensity threshold vector for the kth iteration with current class number n. fF refers to the filtered image at the previous stage.


### 2.3. 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.(1)
The splitting of the subtracted sinogram into two regions: Ω the region to be filled (in which data are  missing) andθ the region where the information is complete.(2)
A patchPp centered on each pixel p located in the region Ω at the border of the frontier between the two regions Ω and θ so that the patch Pp includes a certain number of pixels belonging to θ.(3)
A set of patchesPq centered on each pixel q of the region θ.Then the patchPp is compared with each patch Pq using the following similarity metric:
(2.6)d(p,q)=‖Pp-Pq‖,
where the distance d(p,q) is calculated between the corresponding pixels in patches Pp and Pq that belong to the region θ. Pp and Pq are of the same sizes, and ∥·∥ denotes the Euclidean distance between them. The patch Pq that minimizes this distance is selected and its contents are copied into Pp to restore the missing pixels of Pp that are located in the region Ω. An 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 patchesPp 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:
(2.7)Pr(p)=λC(p)+(1-λ)D(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 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 pixelp with the highest priority (calculated by (2.7)) and the search for the best match between the patch Pp and the set of patches Pq - determined by the similarity metric (2.6). Then, the pixels of patch Pp located in the region Ω are then restored with the corresponding known pixels of the patch Pq.(3)
The updates of the fidelity and data terms for pixels ofPp 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 Pp for in-painting is in fact set smaller than when applying the search for the best match in (2.6).


## 2.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].LetfO and fF, respectively, denote the images before and after filtering, the filter output is given by
(2.1)fiF=∑j∈NiwijfjO∑j∈Niwij,(2.2)wij=exp⁡(-‖si-sj‖2,α2h),
where fiO is the pixel located at the center of the neighborhood N and fjO the pixels located in the neighborhood of fiO. wij denotes the weight between pixels fiOand fjO and is calculated as a similarity measure between the two patches si and sj surrounding each pair of pixels i and j in the neighbourhood Ni, 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–18].


## 2.2. Image Segmentation

The MIMS method is based on the maximum mutual information (MMI) and allows determining the class numbern based on the difference of mutual information (DMI) [19]. The mutual information MI(A,B) between image A and B is defined by
(2.3)MI(A,B)=H(A)+H(B)-H(A,B),
where H(A) and H(B), respectively, denote the entropies of the images A and B,and H(A,B) the mutual entropy. Based on [15, 16], the MI(A,B) can be rewritten as the joint probability density distribution of images A and B:
(2.4)MI(A,B)=∑P(A,B)log⁡P(A,B)P(A)P(B),(2.5)DMIn(fF)=MI(fF,Sn)-MI(fF,Sn-1)MI(fF,fF),
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), and P(A,B) can be computed from the histogram quantization [20]. Based on [18], in (2.5), MI(fF,Sn) is the normalized mutual entropy between the image fFand the current image Sn segmented into n classes, and DMIn(fF) the normalized difference between the entropies MI(fF,Sn) and MI(fF,Sn-1). When the class number increases, DMIn(fF) decreases while MI(fF,Sn) converges towards MI(fF,fF). 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 (MCNa and MCNm for the segmentations of artifacts and metals, resp.) and the threshold ε. A description of the algorithm is given in Figure 2. The filtered image fF is classified into n classes by using an intensity threshold vector Gnk (k being the kth iteration and n referring to the threshold number). The thresholds in Gnk and the class centers can be automatically computed within a simulated annealing- (SA-) based optimization process [21].Figure 2
Flowchart of the MIMS algorithm:J corresponds to MIn(fF,Sn) and Gnk is the intensity threshold vector for the kth iteration with current class number n. fF refers to the filtered image at the previous stage.


## 2.3. 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.(1)
The splitting of the subtracted sinogram into two regions: Ω the region to be filled (in which data are  missing) andθ the region where the information is complete.(2)
A patchPp centered on each pixel p located in the region Ω at the border of the frontier between the two regions Ω and θ so that the patch Pp includes a certain number of pixels belonging to θ.(3)
A set of patchesPq centered on each pixel q of the region θ.Then the patchPp is compared with each patch Pq using the following similarity metric:
(2.6)d(p,q)=‖Pp-Pq‖,
where the distance d(p,q) is calculated between the corresponding pixels in patches Pp and Pq that belong to the region θ. Pp and Pq are of the same sizes, and ∥·∥ denotes the Euclidean distance between them. The patch Pq that minimizes this distance is selected and its contents are copied into Pp to restore the missing pixels of Pp that are located in the region Ω. An 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 patchesPp 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:
(2.7)Pr(p)=λC(p)+(1-λ)D(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 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 pixelp with the highest priority (calculated by (2.7)) and the search for the best match between the patch Pp and the set of patches Pq - determined by the similarity metric (2.6). Then, the pixels of patch Pp located in the region Ω are then restored with the corresponding known pixels of the patch Pq.(3)
The updates of the fidelity and data terms for pixels ofPp 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 Pp for in-painting is in fact set smaller than when applying the search for the best match in (2.6).


## 3. Experiments

Experiments were conducted both on a simulated phantom and clinical CT images. All the images have a size of512×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 mm2. 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), and 3(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)). The 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), and 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).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)
(c)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.
(a)
(b)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.(2)
A cluster-based K-means method to segment the metallic objects. This algorithm requires a suitable setting of the class numbers and the initial center of each cluster.(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 MCNa, MCNm to, respectively, extract the metallic artifacts and object, the threshold ε for the segmentation step, the size of the patches Pp, 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 Table1. 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.Table 1
Parameter setting for Bal’s algorithm and our proposed method. For our proposed correction,MCNa and MCNm depict the maximum class number (MCN) in the MIMS segmentation of artifacts and metal components, respectively.
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σ∥).
Segmentation step
K-means segmentation using 5 classes with CT values: −950 (air),
0 (soft tissue), 200 (normal tissue), 750 (bone), 5000 (metal).
Inpainting step
Linear interpolation using 5 points in each symmetric side
Our proposed method
Pre-filterig step
41
×
41
N
,
11
×
11
s
,
h
=
500
41
×
41
N
,
11
×
11
s
,
h
=
450
41
×
41
N
,
11
×
11
s
,
h
=
450
41
×
41
N
,
11
×
11
s
,
h
=
350
Segmentation step
MCN
a
=
16
,
MCN
m
=
22,ε=600
MCN
a
=
12
,
MCN
m
=
18,ε=600
MCN
a
=
12
,
MCN
m
=
18,ε=600
MCN
a
=
5
,
MCN
m
=
7,ε=600
In-painting step
11
×
11
P
p
,
3
×
3
P
i
q
,
λ
=
0.65
11
×
11
P
p
,
3
×
3
P
i
q
,
λ
=
0.3
11
×
11
P
p
,
3
×
3
P
i
q
,
λ
=
0.3
11
×
11
P
p
,
3
×
3
P
i
q
,
λ
=
0.8


### 3.1. Qualitative Study

Figures5(a) and 5(b) and 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), that is, the 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.Artifact correction on ClinicalImage1 (Figure3(a)) using: (a), Bal’s algorithm [7]; (b), our proposed method.
(a)
(b)Figures6 and 7 illustrate these results on the other two clinical CT images ClinicalImage2 and ClinicalImage3, respectively. A region of interest (ROI) delineated by red lines is zoomed to emphasize the behavior of each algorithm. Severe streak artifacts can be observed in the original image that spread from the ears at the metallic component location. The resulting image analysis (Figures 6(c), 6(d), 7(c) and 7(d)) makes appear that each method provides a substantial reduction of the metallic artifacts. However, results depicted in (Figures 6(d) and 7(d)) appear smoother with less accentuated artifacts. Structures are better preserved while in (Figures 6(c) and 7(c)) some artifacts remain that deteriorate the quality of the image. These artifacts are more pronounced in the close metallic object surrounding.Artifact correction on ClinicalImage2. (a) Original CT image with artifacts; (b) zoomed ROI in the original image; (c) zoomed ROI in the result image obtained from Bal’s algorithm; (d) zoomed ROI in the result image obtained from the proposed method.
(a)
(b)
(c)
(d)Artifact correction on ClinicalImage3. (a) Original CT image with artifacts; (b) zoomed ROI in the original image; (c) zoomed ROI in the result image obtained from Bal’ algorithm; (d) zoomed ROI in the result image obtained from the proposed method.
(a)
(b)
(c)
(d)We can also see in Figures6(d) and 7(d) that some new streak artifacts (pointed by black arrows in Figures 6(d) and 7(d)), though not strong, were introduced into the processed images. The introduced new artifacts might come from the errors of the segmentation and inpainting in the proposed approach.


### 3.2. Quantitative Study

A quantitative study was then carried out on the simulated phantom of Figure4(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:
(3.1)MSE=1NΩ∑j∈Ω(fjc-fjo)2,(3.2)STD=1NΩ∑j∈Ω(fjc-f-Ωc)2,
where NΩ denotes the pixel number in the chosen region Ω (Region1 or Region2). fjc and fjo are the intensity of the pixel within the region Ω in the corrected image fc and the original reference image fo, respectively. f-Ωc characterizes the mean intensity in the ROI, in the corrected images.The phantom images after correction by Bal’s algorithm and our proposed method are, respectively, displayed in Figures8(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 2
MSE and STD on the original corrupted and corrected phantom images. We can see the calculated MSE and STD are lower for the images corrected with our method.
Original uncorrected image (Figure 4(a))
Image corrected by the Bal algorithm (Figure 8(a))
Image corrected by the proposed method (Figure 8(a))
MSE of the whole images
1.37
1.31
1.26
STD
Region1
302.98
64.21
42.20
Region2
293.52
99.04
56.81
MSE
Region1
15.15
3.27
2.16
Region2
6.22
2.66
1.13Line profile comparison for different artifact corrections. (a) corrected image using Bal’s algorithm; (b) corrected image using our method; (c) line profiles along a given line in the original image and each corrected image. The location of the line is drawn in red in the medallion image that appears at the upper left part of the plot.
(a)
(b)
(c)Table3 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.Table 3
Computation cost (in CPU seconds) for each steps in the Bal algorithm and the proposed correction method.
ClinicalImage1
ClinicalImage2
ClinicalImage3
PhantomImage
Bal’s Method in [7]
Pre-filtering step
0.98
0.97
0.96
0.97
Segmentation step
67.23
23.44
24.21
14.13
Inpainting step
1.32
1.57
1.22
0.85
Total
69.53
25.98
26.39
15.95
Proposed correction (unparallelized)
Pre-filtering step
237.56
236.92
237.17
236.32
Segmentation step (metal)
34.45
26.26
28.34
12.39
Segmentation step (artifacts)
26.21
16.96
16.92
4.93
Inpainting step
44.02
49.09
64.78
52.81
Total
342.24
329.23
347.21
306.45
Proposed correction (CUDA parallelized)
Pre-filtering step
6.23
5.97
6.39
5.98
Segmentation step (metal)
34.45
26.26
28.34
12.39
Segmentation step (artifacts)
26.21
16.96
16.92
4.93
Inpainting step
1.42
1.58
2.06
1.74
Total
33.86
24.51
25.37
12.65


## 3.1. Qualitative Study

Figures5(a) and 5(b) and 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), that is, the 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.Artifact correction on ClinicalImage1 (Figure3(a)) using: (a), Bal’s algorithm [7]; (b), our proposed method.
(a)
(b)Figures6 and 7 illustrate these results on the other two clinical CT images ClinicalImage2 and ClinicalImage3, respectively. A region of interest (ROI) delineated by red lines is zoomed to emphasize the behavior of each algorithm. Severe streak artifacts can be observed in the original image that spread from the ears at the metallic component location. The resulting image analysis (Figures 6(c), 6(d), 7(c) and 7(d)) makes appear that each method provides a substantial reduction of the metallic artifacts. However, results depicted in (Figures 6(d) and 7(d)) appear smoother with less accentuated artifacts. Structures are better preserved while in (Figures 6(c) and 7(c)) some artifacts remain that deteriorate the quality of the image. These artifacts are more pronounced in the close metallic object surrounding.Artifact correction on ClinicalImage2. (a) Original CT image with artifacts; (b) zoomed ROI in the original image; (c) zoomed ROI in the result image obtained from Bal’s algorithm; (d) zoomed ROI in the result image obtained from the proposed method.
(a)
(b)
(c)
(d)Artifact correction on ClinicalImage3. (a) Original CT image with artifacts; (b) zoomed ROI in the original image; (c) zoomed ROI in the result image obtained from Bal’ algorithm; (d) zoomed ROI in the result image obtained from the proposed method.
(a)
(b)
(c)
(d)We can also see in Figures6(d) and 7(d) that some new streak artifacts (pointed by black arrows in Figures 6(d) and 7(d)), though not strong, were introduced into the processed images. The introduced new artifacts might come from the errors of the segmentation and inpainting in the proposed approach.


## 3.2. Quantitative Study

A quantitative study was then carried out on the simulated phantom of Figure4(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:
(3.1)MSE=1NΩ∑j∈Ω(fjc-fjo)2,(3.2)STD=1NΩ∑j∈Ω(fjc-f-Ωc)2,
where NΩ denotes the pixel number in the chosen region Ω (Region1 or Region2). fjc and fjo are the intensity of the pixel within the region Ω in the corrected image fc and the original reference image fo, respectively. f-Ωc characterizes the mean intensity in the ROI, in the corrected images.The phantom images after correction by Bal’s algorithm and our proposed method are, respectively, displayed in Figures8(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 2
MSE and STD on the original corrupted and corrected phantom images. We can see the calculated MSE and STD are lower for the images corrected with our method.
Original uncorrected image (Figure 4(a))
Image corrected by the Bal algorithm (Figure 8(a))
Image corrected by the proposed method (Figure 8(a))
MSE of the whole images
1.37
1.31
1.26
STD
Region1
302.98
64.21
42.20
Region2
293.52
99.04
56.81
MSE
Region1
15.15
3.27
2.16
Region2
6.22
2.66
1.13Line profile comparison for different artifact corrections. (a) corrected image using Bal’s algorithm; (b) corrected image using our method; (c) line profiles along a given line in the original image and each corrected image. The location of the line is drawn in red in the medallion image that appears at the upper left part of the plot.
(a)
(b)
(c)Table3 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.Table 3
Computation cost (in CPU seconds) for each steps in the Bal algorithm and the proposed correction method.
ClinicalImage1
ClinicalImage2
ClinicalImage3
PhantomImage
Bal’s Method in [7]
Pre-filtering step
0.98
0.97
0.96
0.97
Segmentation step
67.23
23.44
24.21
14.13
Inpainting step
1.32
1.57
1.22
0.85
Total
69.53
25.98
26.39
15.95
Proposed correction (unparallelized)
Pre-filtering step
237.56
236.92
237.17
236.32
Segmentation step (metal)
34.45
26.26
28.34
12.39
Segmentation step (artifacts)
26.21
16.96
16.92
4.93
Inpainting step
44.02
49.09
64.78
52.81
Total
342.24
329.23
347.21
306.45
Proposed correction (CUDA parallelized)
Pre-filtering step
6.23
5.97
6.39
5.98
Segmentation step (metal)
34.45
26.26
28.34
12.39
Segmentation step (artifacts)
26.21
16.96
16.92
4.93
Inpainting step
1.42
1.58
2.06
1.74
Total
33.86
24.51
25.37
12.65


## 4. 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 neighborhoodN 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 (MCNa and MCNm 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 Gnk 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 non-corrupted 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.


## 5. 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.



---

*Source: 786281-2012-08-14.xml*