Metal Artifact Reduction in X-ray CT via Ray Profile Correction

Abstract

In computed tomography (CT), metal implants increase the inconsistencies between the measured data and the linear assumption of the Radon transform made by the analytic CT reconstruction algorithm. The inconsistencies appear in the form of dark and bright bands and streaks in the reconstructed image, collectively called metal artifacts. The standard method for metal artifact reduction (MAR) replaces the inconsistent data with interpolated data. However, sinogram interpolation not only introduces new artifacts but it also suffers from the loss of detail near the implanted metals. With the help of a prior image that is usually estimated from the metal artifact-degraded image via computer vision techniques, improvements are feasible but still no MAR method exists that is widely accepted and utilized. We propose a technique that utilizes a prior image from a CT scan taken of the patient before implanting the metal objects. Hence, there is a sufficient amount of structural similarity to cover the loss of detail around the metal implants. Using the prior scan and a segmentation or model of the metal implant, our method then replaces sinogram interpolation with ray profile matching and estimation, which yields much more reliable data estimates for the affected sinogram regions. Experiments with clinical dataset obtained using surgical imaging CT scanner show very promising results.

1. Introduction

X-ray computed tomography (CT) is a leading cross-sectional imaging technique lauded for its high image resolution and rapid speed of acquisition. It is based on the differing photoelectric X-ray absorption properties of tissue and assumes that the X-rays are attenuated exponentially according to the Beer–Lambert law:

$$I_E(L)\approx I_0(E)exp\left\{-{\int }_0^L\mu (\mathit{s}+\lambda ·\mathit{r};E)d\mathit{\lambda }\right\}.$$

(1)

In this equation, at some energy level E, $I_0(E)$ is the number (or intensity) of incident photons and $\mu (\mathit{x};E)$ is the linear attenuation coefficient at location $\mathit{x}$. Then, $I_E(L)$ is the number (or intensity) of detected photons along a linear ray described as a line equation, $\mathit{s}+\lambda ·\mathit{r}$, and its length is L. In this equation, $\mathit{s}$ is the position of the X-ray source, $\mathit{r}$ is the directional vector of the ray, and $\lambda$ is a step size. From this, the Radon transform [1] is derived by taking the logarithm of Equation (1) such that

$$M(L;E)=-ln\frac{I_E(L)}{I_0(E)}\approx {\int }_0^L\mu (\mathit{s}+\lambda ·\mathit{r};E)d\lambda ,$$

(2)

This expression forms the underlying principle of most existing CT reconstruction algorithms. It assumes that the attenuation responses along a ray are commutative and accumulate in an additive fashion.

The approximation as formulated in Equation (2) only holds when the X-ray source is monochromatic. When it comes to (more widely employed) polychromatic X-ray sources, the polychromatic projection data can be experimentally calculated [2] as

$$P(L;E)=-ln\frac{\int I_0(E)exp\left\{-{\int }_0^L\mu (\mathit{s}+\lambda ·\mathit{r};E)d\mathit{\lambda }\right\}dE}{\int I_0(E)dE}$$

(3)

The integration over the X-ray energy spectrum in Equation (3) breaks the linearity between attenuation and penetration length during the reconstruction procedure assumed in Equation (2). Furthermore, low-energy photons are attenuated more easily than high-energy ones (see, for example, [3]), and, as a result, the X-ray beam is hardened as it passes through the material, shifting the energy spectrum toward higher energy.

This becomes a particular problem when CT scanning is conducted for patients bearing metallic implants, which have dramatically increased attenuation properties for lower energies. The selective photon absorption not only increases the amount of dose absorbed by the patient, it also amplifies the X-ray beam’s hardness. A direct effect is that more photons are detected than expected. A failure to consider the non-linear behavior of the beam hardening effect in Equation (2) results in dark streaks along the lines of greatest attenuation [4,5]. The high pass filter used in a Filtered-back projection (FBP) reconstruction algorithm [6] further exaggerates the differences between adjacent detector elements where one has received a hardened beam and the other one has not. This unintended contrast produces bright streaks in other directions. As a consequence, by ways of these adverse mechanisms, metal artifacts obscure information about anatomical structures, making it difficult for radiologists to correctly interpret the affected CT images.

There have been extensive efforts in developing metal artifact reduction (MAR) algorithms to compensate the approximation errors in Equation (2) caused by implanted metals or high density objects. These efforts can be largely divided into two types of approaches: iterative reconstruction and sinogram correction.

The iterative reconstruction methods adapt existing X-ray CT systems by incorporating one or more of the following types of a-priori knowledge: (1) low-level information (such as tissue classified image or metal region classified sinogram) of the images to be reconstructed [7,8,9,10]; (2) the X-ray spectrum of the source [11]; (3) the attenuation functions of the base materials [12,13,14]; and (4) the composition of the metal components [15]. More recently proposed iterative algorithms attempt to reduce beam hardening effects without the need of any prior knowledge by decomposing the image to be reconstructed into low and high density components [16,17]. Nevertheless, iterative reconstruction techniques have the downside that they are computationally demanding and so using them in clinical practice remains a challenge.

On the other hand, the sinogram correction methods aim to directly correct the metal shadow regions in the projection data in which the corresponding rays have interacted with metal objects. One early approach replaces the corrupted data with their neighbors using linear [18] or high-order interpolation schemes [19,20,21]. However, the interpolation-based MARs often suffer from loss of detail around the metal objects, and they also have high propensity to introduce new streak artifacts [22]. To address the lack of structural information, Prell et al. [23] and Meyer et al. [24] attempted to build prior CT images by roughly segmenting the uncorrected or pre-corrected CT image into soft-tissue, air, and bone equivalent materials. This has been a promising idea and further efforts have emerged that seek to produce a better prior image with the help of advanced computer vision techniques [25,26].

In this paper, we present a new MAR method that also uses the general approach of correcting a contaminated sinogram by substituting corrupted data by cleaner data available in prior images. The synthesis process we propose is not unlike the one often used in image-guided surgery (IGS) [27]. These methods perform a real-time correlation of the operative field with a preoperative imaging dataset to show the precise location of a selected surgical instrument in the surrounding anatomic structures. To realize this, before the surgery, the patient undergoes a series of CT scans that reveal the soft tissues and bony structures. In our scenario, these preoperative CT images serve as prior images to help remove the metal artifacts that appear in intra- or post-surgery CT scans due to the implanted metal objects. Since such prior images have been acquired from the same patient, they will likely contain very similar internal structures, especially around metal implants. Furthermore, as these regions are often at least partially surrounded by bone, it is unlikely that they are markedly deformed during the surgery. Thus, to find surrogate values to replace unreliable data in the sinogram, we first search ray paths in the prior images that have very similar density profiles along the ray passing through the metal objects. Then, the best matched prior ray profiles are used to correct the ray paths profiles that are corrupted by metal artifacts. Finally, the unreliable data are replaced with the re-projections of corrected ray profiles. We explored this general idea in [28] but this preliminary work was limited to 2D fan-beam CT geometry. In this paper, we generalize our method to 3D cone-beam CT geometry and also present a significantly refined and mature framework.

In the following, Section 2 describes the proposed method and its technical details. Then, in Section 3, we show metal artifact reduction results. Section 4 concludes the paper with a discussion on future research directions for the proposed method.

2. Methods

In the following, we use spine surgery as an example but our method applies in any scenario where a prior patient scan is available. In fact, it even applies in single-scan scenarios when metal is only present in a restricted region of the object.

Our method requires two CT datasets (or one that has a sufficient portion free of metal). One dataset is an artifact-free prior CT scan taken before the (spine) surgery. The other is obtained during the surgery, containing metal artifacts due to the implanted pedicle screws. Although the two CT scans are taken from the same patient, they might be obtained in different conditions (e.g., patient’s pose, X-ray dosage amount, field-of-view, etc.). Therefore, it is necessary to register one volume to the other before applying the proposed ray profile correction method. The volume registration is the process of aligning two or more volumes of the same scene. This process involves designating one image as the reference volume, also called the fixed volume, and applying geometric transformations or local displacements to the other volumes so that they align with the reference. This registration step is described in Section 2.1. Ray profile correction is also required to know which parts of a ray path belong to metal objects, and which ones do not. For this, we extract the implanted pedicle screws from the uncorrected CT volume (the volume suffering from metal artifacts). This metal localization step is described in Section 2.2. In the following, we call a set of sample points along a ray the ray profile. A line integral is then computed as the weighted sum of all sample points of a ray profile. In addition, we define the regions where the corresponding rays pass through metal objects as metal shadow. The projection values under the metal shadow are unreliable because of beam-hardening, photon starvation, and so on, and they will result in metal artifacts. Our goal is to compute surrogate values in the metal shadow regions by correcting the corresponding ray profiles using the aligned prior CT volume and geometric information of the implanted metal (here the pedicle screws). This new correction scheme is explained in Section 2.3. Finally, the corrected metal shadow is combined smoothly with the original CT projection data, as described in Section 2.4. The overall process is illustrated in Figure 1.

2.1. Rigid Volume Registration

For the ray profile correction, we need to find matched prior profiles from the prior CT volume generated in the pre-operative CT scan. This prior CT volume is usually significantly misaligned with the CT volume obtained during or after the surgery. The patient may be in different pose or the CT scan may cover a different range of the spine region (or have a different field-of-view). All these two volumes might have in common is the surgical region itself. One naïve approach for finding matched ray profiles would be to exhaustively search the prior CT volume with the ray profiles extracted from the uncorrected CT volume. This approach would be computationally very demanding as the search space is almost infinite.

Instead, we align the two volumes and then find the set of matched prior profiles. The challenge in aligning two CT volumes taken at different times is that there can be large discrepancies in the soft internal structures (e.g., tissues). To resolve this problem, we first extract bone structures which are quite robust to deformation and also less contaminated by metal artifacts, in contrast to soft structures (see Figure 2 for a visualization). For the bone structure extraction, assuming there are only low and high density materials, we employ the balanced histogram thresholding (BHT) method [29]. Figure 3 shows two CT volumes obtained before and after surgery along with bone structures extracted using the BHT method. The prior volume is then rigidly registered to the uncorrected one by solving the following minimization problem:

$$\widehat{\mathit{\theta }}=\underset{\mathit{\theta }}{\mathrm{argmin}}\sum _{i=1}^N{w}_i·|f_i^{unc}-T_{\mathit{\theta }}{(f^{pri})}_i|.$$

(4)

Here, f is the bone-only CT volume and its superscripts $unc$ and $pri$ indicate the uncorrected and the prior CT volumes, respectively, which consist of N voxels in total. For brevity, f is a flattened 1D array of a volume in row-major order and i is an index at location (x, y, z) in the volume. The $T_{\mathit{\theta }}(·)$ is a rigid volume transformation operator with parameter vector, $\theta$, which includes three translations and three rotations in a 3D Cartesian coordinate system. Using Equation (4), we find the optimal parameter vector, $\widehat{\mathit{\theta }}$, which has a minimum weighted sum of the voxel-wise absolute difference between the two volumes. The weight term, $w_i$, is there to further penalize the data mismatch term if two voxels came from different anatomical structures (materials). It is formalized as follows:

$$w_i=\left\{\begin{array}{cc}1\hfill & \mathrm{if}{f}_i^{unc}\mathrm{and}{T}_{\mathit{\theta }}{(f^{pri})}_i\in \mathrm{same}\mathrm{material}\hfill \\ c\hfill & \mathrm{if}\mathrm{otherwise}\hfill \end{array}\right.$$

(5)

Here, there are two material types, low and high, and they are pre-computed using the BHT segmentation method. We give more penalty, $c(\ge 1)$, if the two voxels are not in the same material. This minimizes the contribution of mis-categorized voxels in the uncorrected volume, such as voxels in the bright band or in implanted metals that are regarded as bone after applying the BHT segmentation.

We use a GPU-accelerated Hybrid-PSO (particle swarm optimization) algorithm to solve the minimization problem in Equation (4), which avoids a convergence to a local minima [30,31]. More specifically, in each generation, we randomly choose half of the particles and randomly adjust either a translation or rotation parameter with uniform probability. In every third generation, we pick half of the worst particles. The first half of these are replaced with completely new random values. Among the remainder, three-fifth of the particles are randomized as we do in each generation and the crossover is applied to the others. These types of randomization strategies have proven effective in finding the global solution in different optimization tasks [31,32]. A new generation is spawned until either there are no change in the best solution compared to previous one or it reaches the maximum number of generations pre-defined. Figure 3 shows an example of the rigid registration result.

2.2. Localization of Implanted Metal Objects

Before we can correct the ray profiles (see Section 2.3), we need to know whether a sample point in a given ray profile originates from metal (or not). For this purpose, we segment the screws from the uncorrected CT volume in the following two steps. In the first step, we use a simple threshold-based segmentation approach to get coarsely metal-only segmented result. After that, we apply the DBSCAN (Density-based spatial clustering of applications with noise) algorithm [33] to the segmented structures. DBSCAN is a popular clustering algorithm which classifies points that are not well connected to a cluster as outliers. We found that DBSCAN did very well to remove any remaining noise and obtain an accurate segmentation of the metal objects (in our case, the screws). The clean segmentation also allows us to precisely determine how many screws were implanted. Optionally, we might also include prior geometric knowledge to accelerate the process and to further improve the clustering accuracy. Figure 4 shows a screw extracted from an uncorrected CT volume. Figure 4c is what we refer to as the metal-only CT volume.

2.3. Ray Profile Correction

To compute the profile’s surrogate values (the values subject to replacement), we use the observation that metal artifacts usually appear around implanted metals and that the degree of corruption tends to decrease with distance from the metal region. Using this observation, the noisy ray profiles are corrected by dividing it into two regions: metal and non-metal regions. A metal region is the part of a profile that traverses a metal-only CT volume. For these regions, since we usually know the material of the implanted metals and their linear attenuation coefficients, the surrogate values are replaced with the linearly interpolated values of the two profiles extracted from the prior and the metal-only CT volume. We use linear interpolation to take into account the partial volume effect around the metal boundaries. For the non-metal regions, the surrogate values are computed by linear interpolating between the noisy and prior profiles extracted from the uncorrected and the prior CT volume, respectively. The interpolation weight is given by the distance from the nearby metal boundaries along the ray path. It takes into account that, when a point in a ray profile is close to a metal, it is more likely deteriorated by metal artifacts and thus we put more emphasis on prior information. Conversely, when a profile point is sufficiently far away from metal, we can safely rely on its value in the uncorrected volume. As such, our method smoothly blends prior image information into the currently acquired imagery but only at locations where the current image information is likely unreliable due to metal artifacts.

Our ray profile correction scheme is formulated as follows:

$$p_i^{new}=\left\{\begin{array}{cc}\mathrm{lerp}(\rho ,p_i^{clean},p_i^{metal}/\rho )\hfill & \mathrm{if}i\in \mathrm{metal}\mathrm{region}\hfill \\ \mathrm{lerp}(p_i^{clean},p_i^{noisy},k·exp(-p_i^{dt}/h))\hfill & \mathrm{if}i\in \mathrm{non-metal}\mathrm{region}\hfill \end{array}\right.$$

(6)

where lerp($\alpha ,\beta ,\omega$) is the linear interpolation operator such that $\omega ·\alpha +(1-\omega )·\beta$. In this equation, $p_i$ represents the sampled value of a ray profile at position i while its superscripts $new$, $metal$, $noisy$, and $clean$ indicate the corrected profile and the profiles extracted from the metal-only, uncorrected, and aligned prior volumes, respectively. The superscript $dt$, denotes the distance transform of $p^{metal}$ and henceforth, $p_i^{dt}$ is the distance from the position i to the closest metal boundary along the ray path [34]. The value $\rho$ is the linear attenuation coefficient of the implanted metals while h is a scalar that controls the smoothness of the weight factor and k is a constant value representing how similar the aligned prior volume is to the uncorrected one. In our study, k is experimentally determined as $0.7$ and it will vary depending on the quality of the volume registration, as discussed in Section 2.1. Automatic determination of the values, h and k, will be our future research topic.

Figure 5 shows an example of the proposed ray profile correction scheme that is applied along a ray depicted as a red line in a zoomed-in image at the upper-left corner. In the non-metal region, as the location along the ray is far from the metal region, the two profiles, $p^{noisy}$ and $p^{clean}$, become similar to each other; on the other hand, as the location is close to the metal region, $p^{noisy}$ shows much higher intensities due to the metal artifacts. Thus, $p^{new}$ in the non-metal region is calculated by relying more on the $p^{clean}$ as the location is close to the metal region. The weight function in the non-metal region is shown in Figure 5c. The corrected values in the metal region are calculated by interpolating between $\rho$ and $p^{clean}$; if a probability of being the implanted metal is high, the corrected values are more biased on $\rho$ and vice versa. The probability (or weight function) is derived from $p^{metal}$ and shown in Figure 5b. It indicates the degree of confidence that a location is filled with the implanted metal.

2.4. Seamless In-Painting

Recall that sample profiles are only computed for rays that traverse the segmented metal in the metal-only CT volume. We integrate these rays and store them in a corrected sinogram. Figure 6a shows a portion of an uncorrected sinogram while Figure 6b shows the same region with the corrected profiles only. The final task is to replace the metal shadow regions of the uncorrected sinogram with these corrected regions. However, a direct replacement of the data can lead to undesired discontinuities around the boundary of the metal shadow, resulting in the generation of new artifacts [22]. Therefore, it is important to seamlessly combine the new data with the existing ones at the boundary while internally keeping the relative contrast and the details of the data. To achieve such a seamless in-painting, we solve the following minimization problem [35]:

$${\widehat{P}}^{new}=\underset{P^{new}}{\min }\sum _{i\in R}\left(\sum _{j\in N_i\cap R}{(\nabla P_j^{new}-\nabla P_j^{corr})}^2+\sum _{j\in N_i\cap \sim R}{(P_j^{new}-P_j^{orig})}^2\right),$$

(7)

where P is projection data and its superscripts $new$, $corr$, and $orig$ represent the in-painted, corrected, and original projection data, respectively. R denotes the metal shadow region and $N_i$ is the eight-connected neighborhood of a pixel, i. In this equation, the first term aims to preserve the gradients of the original (uncorrected) projection data in the metal shadow regions while lowering the intensities to non-metal values. Preserving the original gradients ensures that the detail and contrast in the projection data are maintained for the subsequent reconstruction. The second term in the equation affects the metal shadow region’s boundary only and ensures a smooth transition to the outside regions. Figure 6c shows a result of this seamless in-painting process using the same region as in panel Figure 6a,b with Gauss–Seidel iteration and successive over-relaxation.

3. Results and Discussion

To test the proposed MAR framework, we used clinical CT projection data obtained during an image-guided surgery procedure on a cervical spine region using a Medtronic O-Arm surgical imaging CT scanner. The scanner has a source to axis distance of $647.7$ mm and a source to detector distance of $1147.7$ mm. It is equipped with a flat X-ray detector with $1024\times 384$ bins and an active area of $393.432\times 290.224$ mm${}^2$. During the scan, 360 projections were collected uniformly distributed over 360${}^{\circ }$. The 3D reconstruction used a filtered back-projection algorithm [1] and produced a $512\times 512\times 192$ volume with a voxel size of $0.415\times 0.415\times 0.83$ mm${}^3$.

Figure 7 and Figure 8 show some results we obtained using our metal artifact reduction algorithm. The spine has two pedicle screws implanted. Figure 7 shows one of them in transverse, sagittal, and coronal views. Note that the sagittal and coronal views are horizontal and vertical cut slices passing through the screws, respectively. Figure 8 shows the other implanted screws in the same manner. Overall, the proposed method effectively removes the metal artifacts (dark/bright bands and streaks) and reveals clear outlines of the implanted pedicle screws, which are suitable for evaluating their placements after the surgery. For example, in Figure 7 (top and bottom rows), the yellow arrow indicates a pedicle screw where only the corrected image (Column b) can reveal that is has been correctly inserted into the bone without extending into the tissue. Likewise, the yellow arrow in the middle row in Figure 7 shows a volume feature that was previously hidden by the beam hardening artifacts (Column a) but is now readily visible.

One side effect of the proposed method is the tendency of blurring the anatomical structures near metal objects. The difference images (Column c) between the uncorrected and corrected images show: (1) the removed artifacts; (2) a bright version of the metal pieces (as mentioned, our method lowers their projection values in the sinogram); and (3) some incorrectly removed details. The latter causes the blurring effects (annotated by the red arrow in Figure 7 and Figure 8). We think this is primarily because of the distance-based artifact region prediction model in Equation (6) where the model estimates the artifact regions based on the distance ($p_i^{dt}$) from a point to the nearest metal boundary along a ray path regardless of whether the point is corrupted by metal artifacts or not. One way to mitigate the blurring effect is by adding an additional stage at the end of our MAR framework to exploit the information hidden in low- and high-pass filtered sinogram [36] or images [37].

4. Conclusions

We present a new method for metal artifact reduction (MAR). It assumes that a prior CT scan taken before implanting the metal objects into the patient is available. Using this prior scan and a segmentation or a model of the metal implant, we utilize a novel ray profile correction scheme that computes an accurate estimate of the rays traversing the projection regions affected by the metal artifacts. Our experiments with clinical CT data indicate that the proposed method can clearly reveal the placements of implanted pedicle screws that were ambiguous before due to the significant metal artifacts.

Future research will extend this work along three directions. Firstly, we will further investigate the behavior of the phenomenon of beam hardening for better metal artifact region prediction. Secondly, using a CAD model available for commercial pedicle screws in Figure 9, we believe a better localization of the implanted metal can be accomplished, leading to a better ray profile correction. Lastly, we also believe that utilizing existing algorithms that can reveal features hidden by metal artifacts in conjunction with a MAR method will be able to help resolve the occasional blurring effects around metals. Finally, we also plan to conduct more clinical case studies which will contain a wide selection of different implanted metals and different amounts of metal artifacts as well as through comparison studies with other existing methods in terms of restoration quality and time performance.