Superiorized method for metal artifact reduction

by   T. Humphries, et al.
University of Washington

Metal artifact reduction (MAR) is a challenging problem in computed tomography (CT) imaging. A popular class of MAR methods replace sinogram measurements that are corrupted by metal with artificial data. While these “projection completion” approaches are successful in eliminating severe artifacts, secondary artifacts may be introduced by the artificial data. In this paper, we propose an approach which uses projection completion to generate a prior image, which is then incorporated into an iterative reconstruction algorithm based on the superiorization framework. The prior image is reconstructed using normalized metal artifact reduction (NMAR), a popular projection completion approach. The iterative algorithm is a modified version of the simultaneous algebraic reconstruction technique (SART), which reduces artifacts by incorporating a polyenergetic forward model, least-squares weighting, and superiorization. The penalty function used for superiorization is a weighted average between a total variation (TV) term and a term promoting similarity with the prior image, similar to penalty functions used in prior image constrained compressive sensing. Because the prior is largely free of severe metal artifacts, these artifacts are discouraged from arising during iterative reconstruction; additionally, because the iterative approach uses the original projection data, it is able to recover information that is lost during the NMAR process. We perform numerical experiments modeling a simple geometric object, as well as several more realistic scenarios such as metal pins, bilateral hip implants, and dental fillings placed within an anatomical phantom. The proposed iterative algorithm is largely successful at eliminating severe metal artifacts as well as secondary artifacts introduced by the NMAR process, especially lost edges of bone structures in the neighborhood of the metal regions.



There are no comments yet.


page 11

page 27

page 28

page 29

page 30


Deep Sinogram Completion with Image Prior for Metal Artifact Reduction in CT Images

Computed tomography (CT) has been widely used for medical diagnosis, ass...

Metal Artifact Reduction in 2D CT Images with Self-supervised Cross-domain Learning

The presence of metallic implants often introduces severe metal artifact...

DuDoNet: Dual Domain Network for CT Metal Artifact Reduction

Computed tomography (CT) is an imaging modality widely used for medical ...

A Bayesian Approach to CT Reconstruction with Uncertain Geometry

Computed tomography is a method for synthesizing volumetric or cross-sec...

Mixed one-bit compressive sensing with applications to overexposure correction for CT reconstruction

When a measurement falls outside the quantization or measurable range, i...

Generative Mask Pyramid Network for CT/CBCT Metal Artifact Reduction with Joint Projection-Sinogram Correction

A conventional approach to computed tomography (CT) or cone beam CT (CBC...

Generative Mask Pyramid Network forCT/CBCT Metal Artifact Reduction with Joint Projection-Sinogram Correction

A conventional approach to computed tomography (CT) or cone beam CT (CBC...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

Artifacts caused by metal objects, such as prostheses, rods, screws, and fillings, are a well-known source of image artifacts in computed tomography (CT). These artifacts are caused by numerous physical factors, including beam hardening, increased measurement noise, photon starvation, partial volume and exponential edge gradient effects, and scatter DNDMS98 . They typically appear as dark streaks or bands between metal objects, as well as thin, alternating dark and light streaks emanating from these objects. These artifacts tend to be dramatic and may obscure important features of the image.

Despite the long history of proposed metal artifact reduction (MAR) techniques (see Ref. gjest16 for a recent review), it remains a challenging problem. Dual-energy systems are able to significantly reduce metal artifacts BDNRBJ11 , but require special equipment. On conventional systems, a common approach, often referred to as projection completion, regards measurements that have passed through a metal object as unreliable, and replaces them with artificially generated data. Typically, a preliminary, artifact corrupted image is first reconstructed from the raw sinogram to identify metal regions in the image space, followed by a simulated forward projection to identify the so-called metal trace in the sinogram domain. Once data inside the metal trace have been replaced, an image is reconstructed from the modified data using a standard technique such as filtered backprojection (FBP).

Early projection completion methods GP81 ; KHE87 used linear interpolation (LI) to replace the metal trace within each column of the sinogram. While the LI approach is effective in removing the most severe metal artifacts, it also produces secondary artifacts, due mainly to the loss of edge information about other structures (e.g. bone or air pockets) where their traces intersect with the metal trace MB09 . A popular way to address this issue is to use the preliminary image to not only identify the metal trace, but also generate a prior image, which omits metal but retains information about the other structures such as bone. A number of authors BS06 ; PKBK09 ; KCWM12 ; LFN09 ; VS12 ; ZYJYJ13 ; THZ13 ; ZY18 propose replacing the metal trace with the corresponding measurements obtained from simulated forward projection of the prior, rather than LI data. Alternatively, in normalized metal artifact reduction (NMAR) methods, the prior is used to normalize the projections MB09 ; MRLSK10 ; MRLSK12 before a sinogram correction using LI; the sinogram is subsequently “de-normalized” to reintroduce lost edge information. The prior image is typically generated using tissue classification of an image reconstructed using FBP BS06 ; PKBK09 ; MRLSK10 ; KCWM12 or by iterative reconstruction with regularization to suppress artifacts LFN09 ; VS12 ; ZYJYJ13 ; THZ13 ; one recent work ZY18

has also proposed generating the prior using a convolutional neural network.

Some approaches based on projection completion employ additional heuristics to reduce secondary artifacts. Ref. WCSLX04 proposes segmenting bone from the preliminary reconstruction (replacing bone pixels with smoothed soft tissue values) and then reintroducing bone to the image after an interpolation-based projection completion is performed. Ref. WK04 uses a distance-dependent spatial weighting to combine the LI image with one processed using multi-dimensional adaptive filtering to reduce noise. Ref. BF11 begins with an LI reconstruction and then performs several iterations in which the reconstruction is forward projected, averaged with experimental data using spatial weighting, then reconstructed again. Other methods retain edge information by performing projection completion in the Fourier KWB12 or wavelet ZRWWB00 ; MARZ13 ; BC17 domains, rather than on the original sinogram. Another recent approach GSYXC18 uses a convolutional neural network to combine information from an NMAR image and the uncorrected image. Finally, projection completion can also be performed directly in the sinogram domain without reconstruction of a preliminary image VJVG10 , avoiding issues of mismatch between the simulated forward projection and true system geometry.

An alternative to projection completion is to use the original data and reconstruct the image using an iterative algorithm. Iterative algorithms, while more computationally expensive than FBP, have the advantage of being able to accurately model physical effects contributing to metal artifacts, such as beam hardening and noise. They also may be able to solve the exterior problem (where the metal-corrupted data are ignored rather than replaced) with the inclusion of regularization terms to mitigate the ill-posedness of the problem. Ref. WSOV96

applied both the algebraic reconstruction technique (ART) and expectation maximization (EM), finding that both outperformed FBP reconstruction. Later work has incorporated polyenergetic modeling into the EM framework 

DNDMS01 ; VN12 , applied EM to projection completed data OB07 , used ART in conjunction with total variation (TV) penalties RBFK11 ; ZWX11

and used optimization-based reconstruction with machine-learned regularization 

BC17 . In Ref. CYSTS19

, the authors simultanously estimate the image and the mismatch between polyenergetic and (idealized) monoenergetic data, using a regularized least squares approach incorporating a prior image and several different regularizers.

Projection completion approaches and iterative methods each have advantages and disadvantages. Projection completion is especially effective in removing severe streaking artifacts from images, but the use of artificially generated data carries an inherent risk of creating new artifacts. Iterative methods can make use of the original projection data, but the poor quality of the metal contaminated measurements make obtaining a high quality image difficult. In this paper we combine projection completion with an iterative method to attempt to mitigate these issues. We use the NMAR method to construct a prior image which is free of severe artifacts, but may contain some secondary artifacts. This image is then used both as the initial estimate and as a prior in an iterative method which reconstructs an image from the original projection data. The iterative method incorporates the superiorization methodology HGDC12 with a secondary objective guided by the prior image as well as a total variation (TV) minimization term. Numerical experiments indicate that the proposed algorithm is successful in eliminating both severe artifacts due to streaking, as well as secondary artifacts introduced during the NMAR process.

Ii Methodology

ii.1 Mathematical model

We consider a two-dimensional attenuation distribution, , which depends on position, , and the energy of the incident X-ray beam, . If the X-ray beam is monoenergetic with energy , the idealized measurement along a line is modeled as


where is the initial intensity of the beam and is the (idealized) intensity measured by the detector, in counts per second. The line integral, , which is a sample from the Radon transform of , can be obtained by log-transforming the data,


Discretizing as an image with pixels, and taking measurements along lines, then yields a system of linear equations


where is the discretized image of , is the line integral data, and is the system matrix whose th element is the length of intersection of the th line with the th pixel of .

Clinical CT systems, however, generate polyenergetic X-ray beams, with a spectrum of energies typically ranging from zero up to roughly 150 keV. In this instance, the measurement model is expressed as


where is the beam spectrum. Log-transforming the data and applying the mean value theorem for integrals then gives


where , and is an unknown energy in the range of the spectrum, depending on the path . The inconsistency with the monoenergetic model (2) leads to so-called beam hardening artifacts BD76

, including cupping and streaking. From the physical perspective, as the polyenergetic beam passes through the object, photons with lower energy are attenuated at a higher rate than photons with greater energy, causing the spectrum of the beam to “harden” as it becomes skewed towards higher energies. Metal objects induce particularly severe beam hardening artifacts due to their high attenuation coefficients, particularly at the low end of the energy spectrum.

Another contributing factor to metal artifacts is increased image noise. Due to the stochastic nature of X-ray interactions with matter, the measurement along line

is typically modeled as a Poisson random variable

, with mean equal to

. The measurement therefore has a signal-to-noise ratio (SNR) of

; i.e., the SNR deterioriates as becomes small. This situation occurs if the initial beam intensity is small, or if the line integral through is large. Measurements passing through metal objects thus tend to be very noisy, degrading image quality. In the extreme case may be so small that (no counts are detected along the line), resulting in an effective attenuation measurement of infinity.

ii.2 Beam hardening and metal artifact reduction

A standard approach to reduce beam hardening artifacts is to perform water correction (sometimes known as soft tissue correction) on the measured data JS78 . Water correction (Algorithm 1) simulates a monoenergetic measurement from the corresponding polyenergetic measurement via a two step process. First, the effective length, , of water through which the polyenergetic beam would have to pass to generate is estimated. This can be accomplished by solving a nonlinear equation (Line 3 of Algorithm 1, with denoting the attenuation coefficient of water) or, more commonly, by interpolating from a table of measured values based on known thicknesses of water. Once is known, the corresponding measurement at some reference energy can be obtained straightforwardly.

1 Given: Polyenergetic measurements , . for  to  do
2         Solve for ;
3         Set ;
5      end for
return .
Algorithm 1 Water correction.

Water correction essentially assumes that the object consists only of water or water-like materials. This assumption breaks down if the object contains dense material like bone, contrast agents, or metal, whose attenuation curves do not behave like scaled versions of the attenuation curve of water. Water correction is therefore effective in reducing cupping artifacts, but not the so-called second-order artifacts caused by those materials. Correcting for these artifacts typically requires more advanced data corrections JS78 ; JR97 ; KMPK10 or iterative methods which directly model polyenergetic X-rays DNDMS01 ; EF02 ; VVDBS11 ; LS14 .

As discussed in the Introduction, a standard method of metal artifact reduction (MAR) is to replace metal-contaminated measurements with linearly interpolated data GP81 ; KHE87 . Pseudocode for this method, which we denote as MAR, is provided in Algorithm 2. A preliminary image is reconstructed from the water corrected sinogram, and the image is segmented into non-metal and metal parts using thresholding. The metal index is forward projected to obtain the metal trace in the sinogram domain. The sinogram is then processed columnwise (i.e. with respect to each projection angle) to replace measurements corresponding to the metal trace with their interpolated values. An image can then be reconstructed from the corrected sinogram, and the metal objects can be re-inserted if desired.

1 Given: Water-corrected sinogram, , and system matrix, . Reconstruct image from ;
2 Segment to obtain metal index ;
3 Set (forward projection);
4 Set wherever ;
5 for every column of  do
6      Fill in zero values by linear interpolation to obtain corresponding column of ;
8      end for
Algorithm 2 MAR method

The MAR method reduces severe artifacts caused by metal, but produces secondary artifacts, as the interpolation results in the loss of edge information about other structures. The normalized metal artifact reduction (NMAR) method MRLSK10 (Algorithm 3) addresses this issue through the use of a prior image. Following the preliminary reconstruction, is segmented into air, soft tissue, bone, and metal. The prior image is created by setting regions containing air and soft tissue to their respective attenuation values at the reference energy. Regions containing bone are assigned the corresponding attenuation coefficients from , to accurately capture the varying attenuation of bone. The values assigned to the metal regions are unimportant MRLSK10 ; we use soft tissue for simplicity. The prior image is then forward projected, and the resulting sinogram is divided elementwise () into the original sinogram to normalize it. The MAR algorithm is then applied to the normalized sinogram, followed by an elementwise multiplication () with to “de-normalize” the sinogram. The de-normalization step effectively recovers edge information about bone and other structures that may be lost during the MAR step MB09 ; MRLSK10 .

Fig. 1 highlights the differences between MAR and NMAR. In the sinogram profiles (top right) we can observe that the interpolated profile used by MAR (red line) does not accurately capture the traces of the bone and air pocket, while the NMAR method does. While the MAR correction (bottom center) removes the large streak between the two metal objects that appears in the water corrected (bottom left) image, secondary streaking artifacts are created between the metal objects and the bone and air pockets. NMAR is able to remove virtually all artifacts in this simple experiment, though a mild artifact persists surrounding the bone region.

1 Given: Water-corrected sinogram, , system matrix, , reference energy . Reconstruct image from ;
2 Segment to obtain indices and (air, soft tissue, bone, metal);
3 Set wherever ;
4 Set wherever ;
5 Set wherever ;
6 ;
7 ;
8 ;
9 ;
Algorithm 3 NMAR method
Figure 1: Illustration of MAR and NMAR methods. Top row: Left: sinogram of a simple mathematical phantom including metal objects (small circular regions) as well as bone (oblong region in bottom right) and air pocket (top right). Right: Profile along white dashed line in sinogram, showing water corrected sinogram (black line), MAR correction (red line) and NMAR correction (blue line). Profiles have been vertically offset to highlight differences. Bottom: Images reconstructed from water corrected (left), MAR corrected (center) and NMAR corrected (right) sinograms using filtered backprojection. Phantom specifications are given in Section III.

ii.3 Iterative reconstruction

Standard iterative reconstruction methods are based on solving the linear system as described in (3). Our approach is based on a block-iterative variant of the simulataneous algebraic reconstruction technique (SART) AK84 , as described in Ref. CE02 . We first define diagonal matrices and as:

To implement block iteration, the columns of the sinogram are evenly partitioned into subsets, indexed by ; for example, if , then would correspond to the first, sixth, eleventh, etc. columns, to the second, seventh, twelfth, and so on.This accelerates the convergence of the algorithm. One iteration of the block iterative SART (BI-SART) algorithm is then given by


The subscript indicates that only rows of and corresponding to the measurements in are used, including when forming the matrices and . The operator ensures non-negativity of the solution after every iteration.

The BI-SART algorithm accounts for neither the statistical uncertainty in the measured data, nor the polyenergetic nature of the X-ray measurements. Additionally, images reconstructed by BI-SART may be of poor quality when the data are noisy. In Ref. HG17 , we introduced three enhancements to the BI-SART algorithm to address these issues, which are especially important in the context of metal artifact reduction. We summarize these enhancements below:

ii.3.1 Polyenergetic forward model

To account for the polyenergetic measurements, we adopt the polyenergetic SART (pSART) method LS14b . pSART defines a polyenergetic forward projection operation using linear interpolation between tabulated attenuation curves for known basis materials such as soft tissue, bone, and metal. Specifically, letting denote the attenuation image at the reference energy , we define


where and are the linear attenuation coefficient (LAC) curves of the two basis materials whose LACs at bracket the pixel value . For example, if has an LAC halfway between that of soft tissue and bone at , its LAC is assumed to be halfway between that of soft tissue and bone at all other energies as well. The polyenergetic forward projection operation is then defined as


where is a discretization of the continuous beam spectrum into energy levels, and is the th row of (cf. (4) and (5)). The pSART iteration is achieved by simply replacing the monoenergetic forward projection in SART with in (7).:


While pSART has not been proven to converge in general, it has been shown to be effective in reducing beam hardening artifacts in numerical experiments LS14b ; H15 ; HWF17 .

ii.3.2 Weighted least squares

To model measurement uncertainty, a weighted least squares approach (WLS) SC00 can be employed. We define the diagonal weighting matrix as


We can then apply BI-SART to the system , which assigns proportionally higher weight to the less noisy measurements. Incorporating the polyenergetic forward projection as well, the modified iteration is given by


where 111Matrix does not need to be modified, as the factor of is cancelled during multiplication with ..

ii.3.3 Superiorization

Superiorization HGDC12 is an optimization heuristic in which solutions generated by an iterative algorithm are perturbed in every iteration, to improve the quality of the solution in some respect. For example, if (6) is taken to be the basic iterative algorithm with defined in (12), then the superiorized version is of the form


where , , and the

are a sequence of bounded perturbation vectors. The key result underlying superiorization is that if the basic algorithm is

perturbation resilient, the superiorized version is able to find a solution that is equally satisfactory with respect to satisfying the constraints of the problem ( in our case), though typically this requires more iterations. The perturbation vectors are usually chosen to be nonascending directions of some penalty function at the current iterate, for example, . Thus, the solution found by the superiorized algorithm can be expected to be superior with respect to , while equally compatible with the constraints. The methodology has been featured in a special edition of Inverse Problems CHJ17 , and a continously updated online bibliography is maintained at

All three of the proposed enhancements (polyenergetic forward projection, weighted least squares, and superiorization) can be implemented independently of one another. In HWF17 we found that polyenergetic forward projection and TV superiorization could be used to eliminate beam hardening artifacts as well as artifacts due to sparse-view and limited-angle data. In HG17 , we found that including all three enhancements in the reconstruction algorithm provided the best results when metal artifacts were present. Pseudocode for this algorithm is presented in Algorithm 4. In the implementation of the algorithm, we allow for a total of perturbations to be applied in every iteration, to further improve the solution with respect to the penalty function. (This is of course equivalent to a single perturbation, albeit one that cannot be determined a priori). The step sizes decrease geometrically with rate to ensure that the perturbations are summable.

In HG17 , total variation (TV) was used as the penalty function employed by superiorization:


where denotes the pixel in the th row and th column of the image . The small parameter is introduced to avoid singularity of at points where the image is piecewise constant. While the algorithm, which we denoted as wPSART-TV (weighted polyenergetic SART with TV superiorization) was effective in reducing metal artifacts, we found that it was difficult to remove strong streaking artifacts between metal objects; for example, in the case of a bilateral hip implant. This was especially true at low count rates, where photon starvation was a significant factor. This motivates the development of the hybrid algorithm, discussed in the next section.

1 Given: Log-transformed polyenergetic sinogram , system matrix , number of subsets , reference energy , spectrum , tabulated attenuation curves , initial image , stopping criteria and/or , superiorization parameters , . ;
2 ;
3 while  and  do
4      ;
5      ;
6      while  do
7           ;
8           while true do
9                ;
10                ;
11                ;
12                if  and  then
13                     ;
14                     break;
16                     end if
18                     end while
19                    ;
21                     end while
22                    , where ;
23                     ;
25                     end while
Algorithm 4 Block-iterative, weighted, polyenergetic SART with superiorization

ii.4 Prior image constrained superiorized algorithm

Our prior image constrained superiorized algorithm (Algorithm 5), which we denote as wPSART-PICS, consists of several steps. First, the water corrected sinogram is corrected using NMAR to obtain a metal corrected sinogram. A modified version of Algorithm 4 is then applied to this sinogram to reconstruct a prior image, . In the modified version of the algorithm, monoenergetic forward projection (i.e. multiplication by ) is used instead of , and the weighting matrix is omitted, since the data have been water corrected and the noisy metal trace removed by NMAR. Total varation (14) is used as the penalty function for superiorization, in order to generate a smooth prior image.

The prior image is then incorporated into a second penalty function:


where controls the weighting of the two terms. This type of penalty is used in the prior image constrained compressed sensing (PICCS) approach described in CTL08 ; LTC12 ; we use the term “superiorization” instead of “compressed sensing” as it more accurately describes our method. Algorithm 4 is then applied with the original polyenergetic projection data and as the penalty function. The prior image is also used as the initial estimate, as it is expected to be a good approximation to the true image.

1 Given: Water corrected sinogram , all inputs needed for Algorithm 4. Apply NMAR to to obtain ;
2 Use Algorithm 4 (modified) with as sinogram to reconstruct prior image ;
3 Set ;
4 Use Algorithm 4 with as sinogram, , and as penalty function to reconstruct ;
Algorithm 5 wPSART-PICS algorithm

By penalizing the TV of the difference between the reconstructed image and the prior, effectively promotes similar edge structure between the two images. Our motivation for using in this application is that if the NMAR-based prior is largely free of severe streaking artifacts, these artifacts will also be penalized in the image being reconstructed. Since the final step of Algorithm 5 uses the original polyenergetic data, however, it should be possible to recover image details that are lost due to the interpolation performed by NMAR.

Iii Results

We validated our approach using several numerical experiments. The wPSART-PICS algorithm was implemented in Matlab using the Michigan Image Reconstruction Toolbox (MIRT) MIRT to simulate the CT spectrum, material attenuation curves, and generate the system matrix for iterative reconstruction. We implemented our own version of NMAR in Matlab as well, using the built-in Image Processing Toolbox methods to segment the CT image.

iii.1 Simple phantom experiments

We first generated the simple mathematical phantom shown in Fig. 1 to fine tune the algorithm parameters, before studying more realistic data. The phantom is pixels with a pixel size of 0.75 mm. It consists of a background region of soft tissue, two small circular regions containing titanium, and a larger oblong region containing bone. The linear attenuation coefficients of soft tissue, bone and titanium at the reference energy of 70 keV are , , and 2.44 cm, respectively. The phantom also consists of an air pocket and several low contrast features (modeled after the Shepp-Logan phantom) in the center of the object. The low-constrast features have attenuation values of relative to the background, and are heavily obscured by metal streaking artifacts if no metal artifact correction is performed (see Fig. 1).

We analytically computed 720 parallel-beam views over a 180 arc around the phantom. A 130 kVp spectrum was simulated to generate polyenergetic data. Poisson noise was subsequently added to the data, proportionally to initial count rates of , , , and . To reconstruct the data using Algorithm 5, we ran 24 iterations of the modified algorithm (Line 3) to generate . Algorithm 4 was then run for 32 iterations to produce a final image. In both instances, subsets of projection data were used. The parameters controlling the gradient descent steps within the superiorized algorithm were set to , . To test the sensitivity of the resulting image to the choice of in (15), we performed reconstructions with , and ; the first value corresponding to the case of ordinary TV regularization. For comparison, we also reconstructed an image using the wPSART-TV approach of  Ref. HG17 ; this is equivalent to applying Algorithm 4 with and .

For quantitative comparisons, we defined a pixel region of interest (ROI) in the center of the object, surrounding the six small low-contrast features. Pixel values within the ROI were rescaled to the interval using the formula


where and were set to of the background soft tissue value. The peak signal-to-noise ratio (PSNR) of the region was then calculated as


The rescaling was performed so that the PSNR values were better able to capture differences between the reconstructed values within the ROI.

Fig. 2 shows the images reconstructed using all approaches at the noise level of . While the prior image reconstructed from the NMAR data (B) is free of metal streaking artifacts, there is a small artifact around the oblong bone region in the bottom right due to inaccuracy in the segmentation applied during NMAR. This artifact is absent from images C–H, which are reconstructed from the original data. Image C, reconstructed using the wPSART-TV algorithm, retains a mild streaking artifact between the two metal objects. This is likely due to the fact that the image is reconstructed without any prior knowledge from the NMAR correction, which removes the streak entirely. Image D begins from , but uses only TV as the penalty function during superiorized reconstruction (since ), resulting in blurred edges around the low-contrast features. Images E–H, reconstructed using pSART-PICS with , demonstrate better definition about the features.

Figure 2: Images reconstructed at noise level of . A: True image. B: Prior image reconstructed from NMAR data. C: Image reconstructed using wPSART-TV. D – H: Images reconstructed using wPSART-PICS with and , respectively. Top two rows: full image; bottom two rows: zoomed in on ROI highlighted in red box. Grayscale window is the value of the soft tissue background.

The PSNR values given in Table 1 corroborate these observations. At all four noise levels, the best PSNR values are obtained using wPSART-PICS with an value of 0.0 (column H). This indicates that (for this experiment, in any event), there is no benefit to including the TV of the image in the penalty function (15). This is perhaps not surprising, since the prior image itself is reconstructed using TV superiorization and is therefore quite smooth itself. That being said, the difference in PSNR values between columns E–H is only on the order of 1–2%, indicating that the algorithm is not particularly sensitive to the choice of , as long as . Comparatively, the difference in PSNR values between columns D and E ( versus ) ranges from 2 – 10%.

Counts B C D E F G H
1.0e5 21.36 21.41 21.24 21.73 21.92 21.97 21.97
2.0e5 21.96 22.00 21.30 22.35 22.49 22.58 22.62
5.0e5 23.28 22.24 21.37 23.56 23.82 23.93 23.97
1.0e6 23.95 20.46 21.67 24.20 24.52 24.62 24.66
Table 1: PSNR values of reconstructed images within the ROI shown in Fig. 2. Column headings refer to images shown in Fig. 2.

iii.2 Anatomical phantom experiments

To test the performance of the algorithm on more realistic data, slices of clinical CT images were downloaded from the Cancer Imaging Archive222 Four pixel slices were selected from the dataset, representing the abdominal, shoulder, pelvic and head regions. Metal objects representing pins, dental fillings, and a bilateral hip implant were subsequently manually inserted into the image slices; the pins and hip implants were modeled as titanium, while the fillings were modeled as gold. The objects were modeled after those appearing in Ref. ZY18 . Fan beam data corresponding to 900 views over 360 were simulated using the MIRT toolbox, using the same 130 kVp spectrum and reference energy of 70 keV as for the simple phantom experiment. Soft thresholding ZY18 with base materials of fat, soft tissue, and bone were used to generate attenuation coefficients for nonmetal regions at every energy level when simulating the polyenergetic data. Noise was then added to the data proportional to count rate of . The four true images and their reconstructions using FBP are shown in the first two rows of Fig. 3, illustrating the severity of the metal artifacts when no correction is performed.

Images were reconstructed using the same algorithms as in the previous section. In light of the results from the simple phantom experiment, we only used values of , and for the pSART-PICS algorithm. All superiorized algorithms were run with and ; while the choice of , gave good results for the simple phantom used in the previous section, we found that it tended to oversmooth the images based on real CT data. By reducing both and , we smooth the images less because we perform fewer gradient descent iterations in Algorithm 4, while also reducing the step size more quickly.

Reconstructed images of the four phantoms are shown in Fig. 3 and Fig. 4. Row C of the figures shows the prior images reconstructed from the NMAR data. It is apparent in Fig. 4 that some information has been lost due to the interpolation performed by NMAR; for example, the area around the spine and backmost ribs in the abdomen image, the scapula on the right side of the shoulder image, and the edges of the hip bones in the pelvis image. While the wPSART-TV algorithm (Row D) does a better job of recovering these features, streaking artifacts caused by metal persist in every image, most noticeably between the two implants in the hip image. The three images reconstructed using wPSART-PICS (rows E–G) are able to eliminate streaks caused by metal while also recovering features which were occluded in the prior image. As in the simple phantom experiment, there is little apparent difference in the images reconstructed using (rows F and G), while the image using (i.e., TV as the penalty function) is noticeably smoother. In the abdomen and shoulder images, including the prior in the penalty function has actually introduced some mild streaking into the image along the horizontal direction, since these streaks are also present in the prior image.

Figure 3: Phantoms and reconstructions based on clinical CT image slices. Row A: True phantoms at reference energy of 70 keV with metal inserts indicated in pink. Row B: reconstructions generated using FBP and water corrected data. Row C: Prior images reconstructed from NMAR data. Row D: Images reconstructed using wPSART-TV. Rows E–G: Images reconstructed using wPSART-PICS with and , respectively. Grayscale window is cm. Images were reconstructed at 512512 pixels, but have been cropped to eliminate black space. Green boxes indicate zoomed-in regions shown in Fig. 4.
Figure 4: Zoomed-in regions indicated in Fig. 3.

It is apparent that the wPSART-PICS algorithm does not give an acceptable reconstructed image in the case of the head phantom (last column of Fig. 3). The reason for this is that the metal simulated in the head phantom experiment (gold) has such a high value that the three objects modeling dental fillings block all X-rays passing through them at the simulated count rate of counts per ray. This results in a measured intensity of , translating to infinite attenuation along the line. Due to the projection weighting employed by the algorithm (11), these measurements are ignored during reconstruction; however, the lack of information about pixel intensities along those lines produce artifacts in the final image. The artifacts are actually more severe in the images which begin with the prior as an initial estimate (Rows E–G); this is because there is a slight mismatch between the metal trace in the forward projection of the prior image and that of the measured sinogram, which creates severe streaking in the reconstructed image after just one iteration of the algorithm.

Iv Conclusions

In this paper we present an iterative reconstruction algorithm for CT imaging which uses the superiorization methodology to perform metal artifact reduction (MAR). The algorithm uses a prior image reconstructed using the normalized metal artifact reduction (NMAR) method of Meyer et al. MRLSK10 , which eliminates most severe streaking artifacts caused by metal, but may introduce secondary artifacts in the vicinity of the metal regions. The NMAR-based prior is then incorporated into a penalty function akin to that used in prior image constrained compressive sensing (PICCS) algorithms CTL08 , which is used within a superiorized iterative algorithm based on our previous work (wPSART-TV) HG17 . Numerical experiments modeling several different anatomical scenarios were perfomed, and indicate that the proposed algorithm (wPSART-PICS) is able improve on both NMAR and wPSART-TV. In particular, it is able to recover details that are lost during the NMAR process (particularly with respect to the structure of bone around the metal regions) while also removing streaking artifacts that persist in images reconstructed using wPSART-TV.

While the wPSART-PICS algorithm was effective in removing metal artifacts in most cases, the numerical experiments do highlight some limitations. In an experiment simulating gold dental fillings (fourth column of Fig. 3), the algorithm failed to produce acceptable results, due to the total photon starvation induced by the dense metal objects. In a separate experiment (not shown), we simulated titanium fillings rather than gold; this eliminated the artifacts as the titanium fillings no longer blocked all X-rays. Additionally, we note that the simulated sinogram data for the pelvis phantom also included measurements along which no photons were detected; nonetheless, the wPSART-PICS algorithm was successfully able to reconstruct this image. As indicated in Fig. 5, the key difference appears to be that not all measurements through the metal objects in the pelvis phantom are totally attenuated; only those passing through both objects. Thus, the algorithm is able to reconstruct the object accurately based on the remaining information. We conclude that while the wPSART-PICS algorithm is able to perform well in challenging scenarios including some photon starvation, the NMAR algorithm may be preferable in cases where extreme photon starvation occurs. Additional limitations include the use of spectral knowledge in the forward projection step ((9)), which may not always be available, and the need to tune parameters such as , , and/or in Algorithm 4 to obtain good results.

Figure 5: Fan-beam sinograms corresponding to pelvis and head phantoms with initial count rate of . Measurements corresponding to infinite attenuation are shown in pink.


The authors would like to thank Aviv Gibali (ORT Braude College, Karmiel, Israel) for his initial suggestion to apply superiorization to metal artifact reduction, and contributions to earlier work.


  • [1] Bruno De Man, Johan Nuyts, Patrick Dupont, Guy Marchal, and Paul Suetens. Metal streak artifacts in X-ray computed tomography: a simulation study. In Nuclear Science Symposium, 1998. Conference Record. 1998 IEEE, volume 3, pages 1860–1865. IEEE, 1998.
  • [2] Lars Gjesteby, Bruno De Man, Yannan Jin, Harald Paganetti, Joost Verburg, Drosoula Giantsoudi, and Ge Wang. Metal artifact reduction in CT: where are we after four decades? IEEE Access, 4:5826–5849, 2016.
  • [3] Fabian Bamberg, Alexander Dierks, Konstantin Nikolaou, Maximilian F Reiser, Christoph R Becker, and Thorsten RC Johnson. Metal artifact reduction by dual energy computed tomography using monoenergetic extrapolation. European radiology, 21(7):1424–1429, 2011.
  • [4] Gary H Glover and Norbert J Pelc. An algorithm for the reduction of metal clip artifacts in CT reconstructions. Medical physics, 8(6):799–807, 1981.
  • [5] Willi A Kalender, Robert Hebel, and Johannes Ebersberger. Reduction of CT artifacts caused by metallic implants. Radiology, 164(2):576–577, 1987.
  • [6] Jan Müller and Thorsten M Buzug. Spurious structures created by interpolation-based CT metal artifact reduction. In Medical Imaging 2009: Physics of Medical Imaging, volume 7258, page 72581Y. International Society for Optics and Photonics, 2009.
  • [7] Matthieu Bal and Lothar Spies. Metal artifact reduction in CT using tissue-class modeling and adaptive prefiltering. Medical physics, 33(8):2852–2859, 2006.
  • [8] Daniel Prell, Yiannis Kyriakou, Marcel Beister, and Willi A Kalender. A novel forward projection-based metal artifact reduction method for flat-detector computed tomography. Physics in Medicine & Biology, 54(21):6575, 2009.
  • [9] Seemeen Karimi, Pamela Cosman, Christoph Wald, and Harry Martz. Segmentation of artifacts and anatomy in CT metal artifact reduction. Medical physics, 39(10):5857–5868, 2012.
  • [10] Catherine Lemmens, David Faul, and Johan Nuyts. Suppression of metal artifacts in CT using a reconstruction procedure that combines map and projection completion. IEEE transactions on medical imaging, 28(2):250–260, 2009.
  • [11] Joost M Verburg and Joao Seco. CT metal artifact reduction method correcting for beam hardening and missing projections. Physics in Medicine & Biology, 57(9):2803, 2012.
  • [12] Yanbo Zhang, Hao Yan, Xun Jia, Jian Yang, Steve B Jiang, and Xuanqin Mou. A hybrid metal artifact reduction algorithm for x-ray CT. Medical physics, 40(4), 2013.
  • [13] Zhiwei Tang, Guangshu Hu, and Hui Zhang. Efficient metal artifact reduction method based on improved total variation regularization. J Med Biol Eng, 34(3):261–268, 2013.
  • [14] Yanbo Zhang and Hengyong Yu. Convolutional neural network based metal artifact reduction in X-ray computed tomography. IEEE transactions on medical imaging, 37(6):1370–1381, 2018.
  • [15] Esther Meyer, Rainer Raupach, Michael Lell, Bernhard Schmidt, and Marc Kachelrieß. Normalized metal artifact reduction (NMAR) in computed tomography. Medical physics, 37(10):5482–5493, 2010.
  • [16] Esther Meyer, Rainer Raupach, Michael Lell, Bernhard Schmidt, and Marc Kachelrieß. Frequency split metal artifact reduction (FSMAR) in computed tomography. Medical physics, 39(4):1904–1916, 2012.
  • [17] Jikun Wei, Laigao Chen, George A Sandison, Yun Liang, and Lisa X Xu. X-ray CT high-density artefact suppression in the presence of bones. Physics in Medicine & Biology, 49(24):5407, 2004.
  • [18] Oliver Watzke and Willi A Kalender. A pragmatic approach to metal artifact reduction in CT: merging of metal artifact reduced images. European radiology, 14(5):849–856, 2004.
  • [19] F Edward Boas and Dominik Fleischmann. Evaluation of two iterative techniques for reducing metal artifacts in computed tomography. Radiology, 259(3):894–902, 2011.
  • [20] Bärbel Kratz, Imke Weyers, and Thorsten M Buzug. A fully 3d approach for metal artifact reduction in computed tomography. Medical physics, 39(11):7042–7054, 2012.
  • [21] Shiying Zhao, DD Robertson, Ge Wang, Bruce Whiting, and Kyongtae T Bae. X-ray CT metal artifact reduction using wavelets: an application for imaging total hip prostheses. IEEE transactions on medical imaging, 19(12):1238–1247, 2000.
  • [22] Abolfazl Mehranian, Mohammad Reza Ay, Arman Rahmim, and Habib Zaidi. X-ray CT metal artifact reduction using wavelet domain sparse regularization. IEEE transactions on medical imaging, 32(9):1707–1722, 2013.
  • [23] Parisa Babaheidarian and David Castañón. A randomized approach to reduce metal artifacts in x-ray computed tomography. Electronic Imaging, 2017(17):24–29, 2017.
  • [24] Lars Gjesteby, Hongming Shan, Qingsong Yang, Yan Xi, Bernhard Claus, Yannan Jin, Bruno De Man, and Ge Wang.

    Deep neural network for ct metal artifact reduction with a perceptual loss function.

    In Fifth international conference on image formation in X-ray computed tomography, 2018.
  • [25] Wouter JH Veldkamp, Raoul MS Joemai, Aart J van der Molen, and Jacob Geleijns. Development and validation of segmentation and interpolation techniques in sinograms for metal artifact suppression in CT. Medical physics, 37(2):620–628, 2010.
  • [26] Ge Wang, Donald L Snyder, Joseph A O’Sullivan, and Michael W Vannier. Iterative deblurring for CT metal artifact reduction. IEEE transactions on medical imaging, 15(5):657–664, 1996.
  • [27] Bruno De Man, Johan Nuyts, Patrick Dupont, Guy Marchal, and Paul Suetens. An iterative maximum-likelihood polychromatic algorithm for CT. IEEE transactions on medical imaging, 20(10):999–1008, 2001.
  • [28] Katrien Van Slambrouck and Johan Nuyts. Metal artifact reduction in computed tomography using local models in an image block-iterative scheme. Medical physics, 39(11):7080–7093, 2012.
  • [29] M Oehler and TM Buzug. Statistical image reconstruction for inconsistent CT projection data. Methods of information in medicine, 46(03):261–269, 2007.
  • [30] Ludwig Ritschl, Frank Bergner, Christof Fleischmann, and Marc Kachelrieß. Improved total variation-based CT image reconstruction applied to clinical data. Physics in Medicine & Biology, 56(6):1545, 2011.
  • [31] Xiaomeng Zhang, Jing Wang, and Lei Xing. Metal artifact reduction in x-ray computed tomography (CT) by constrained optimization. Medical physics, 38(2):701–711, 2011.
  • [32] Z. Chang, D. H. Ye, S. Srivastava, J. Thibault, K. Sauer, and C. Bouman. Prior-guided metal artifact reduction for iterative x-ray computed tomography. IEEE Transactions on Medical Imaging, 38(6):1532–1542, June 2019.
  • [33] Gabor T Herman, Edgar Garduño, Ran Davidi, and Yair Censor. Superiorization: An optimization heuristic for medical physics. Medical Physics, 39(9):5532–5546, 2012.
  • [34] R. A. Brooks and G. Di Chiro. Beam hardening in X-ray reconstructive tomography. Phys. Med. Biol., 21(3):390–398, 1976.
  • [35] P. M. Joseph and R. D. Spital. A method for correcting bone induced artifacts in computed tomography. Journal of Computer Assisted Tomography, 2(1):100–108, 1978.
  • [36] P.M. Joseph and C. Ruth. A method for simultaneous correction of spectrum hardening in CT images containing both bone and iodine. Med. Phys., 24(10), 1997.
  • [37] Y. Kyriakou, E. Meyer, D. Prell, and M. Kachelrieß. Empirical beam hardening correction (EBHC) for CT. Med. Phys., 37(10):5179–5187, 2010.
  • [38] I. A. Elbakri and J. A. Fessler. Statistical image reconstruction for polyenergetic X-ray computed tomography. IEEE Trans. Med. Imag., 21(2):89–99, 2002.
  • [39] G. Van Gompel, K. Van Slambrouck, M. Defrise, K.J. Batenburg, J. de Mey, J. Sijbers, and J. Nuyts. Iterative correction of beam hardening artifacts in CT. Med. Phys., 38(7):S36–S49, 2011.
  • [40] Yuan Lin and Ehsan Samei. An efficient polyenergetic sart (psart) reconstruction algorithm for quantitative myocardial CT perfusion. Medical physics, 41(2), 2014.
  • [41] A. H. Andersen and A. C. Kak. Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm. Ultrasonic Imaging, 6(1):81–94, 1984.
  • [42] Y. Censor and T. Elfving. Block-iterative algorithms with diagonally scaled oblique projections for the linear feasibility problem. SIAM Journal on Matrix Analysis and Applications, 24(1):40–58, 2002.
  • [43] T. Humphries and A. Gibali. Superiorized polyenergetic reconstruction algorithm for reduction of metal artifacts in CT images. In 2017 IEEE Nuclear Science Symposium Conference Record, 2017.
  • [44] Y. Lin and E. Samei. An efficient polyenergetic SART (pSART) reconstruction algorithm for quantitative myocardial CT perfusion. Med. Phys., 41(2):021911–1 – 021911–14, 2014.
  • [45] T. Humphries. Technical note: Convergence analysis of a polyenergetic SART algorithm. Medical Physics, 42(7):1407–1404, 2015.
  • [46] T. Humphries, J. Winn, and A. Faridani. Superiorized algorithm for reconstruction of CT images from sparse-view and limited-angle polyenergetic data. Phys. Med. Biol., 62(16):6762, 2017.
  • [47] P. Sukovic and N.H. Clinthorne. Penalized weighted least-squares image reconstruction for dual energy X-ray transmission tomography. IEEE Trans. Med. Imag., 19(11):1075–1081, Nov 2000.
  • [48] Y. Censor, G. T. Herman, and M. Jiang. Special issue on superiorization: theory and applications. Inverse Problems, 33(4), 2017.
  • [49] G.-H. Chen, J. Tang, and S. Leng. Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets. Med. Phys., 35(2):660–663, 2008.
  • [50] P. T. Lauzier, J. Tang, and G.-H. Chen. Time-resolved cardiac interventional cone-beam CT reconstruction from fully truncated projections using the prior image constrained compressed sensing (PICCS) algorithm. Phys. Med. Biol., 57(9):2461, 2012.
  • [51] J. Fessler. Michigan image reconstruction toolbox.