I Introduction
Artifacts caused by metal objects, such as prostheses, rods, screws, and fillings, are a wellknown 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. Dualenergy 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 socalled 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 “denormalized” 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 interpolationbased projection completion is performed. Ref. WK04 uses a distancedependent spatial weighting to combine the LI image with one processed using multidimensional 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 metalcorrupted data are ignored rather than replaced) with the inclusion of regularization terms to mitigate the illposedness 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 ; ZWX11and used optimizationbased reconstruction with machinelearned 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 twodimensional attenuation distribution, , which depends on position, , and the energy of the incident Xray beam, . If the Xray beam is monoenergetic with energy , the idealized measurement along a line is modeled as
(1) 
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 logtransforming the data,
(2) 
Discretizing as an image with pixels, and taking measurements along lines, then yields a system of linear equations
(3) 
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 Xray beams, with a spectrum of energies typically ranging from zero up to roughly 150 keV. In this instance, the measurement model is expressed as
(4) 
where is the beam spectrum. Logtransforming the data and applying the mean value theorem for integrals then gives
(5) 
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 socalled 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 Xray interactions with matter, the measurement along line
is typically modeled as a Poisson random variable
, with mean equal to. The measurement therefore has a signaltonoise 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.
Water correction essentially assumes that the object consists only of water or waterlike 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 socalled secondorder 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 Xrays DNDMS01 ; EF02 ; VVDBS11 ; LS14 .
As discussed in the Introduction, a standard method of metal artifact reduction (MAR) is to replace metalcontaminated 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 nonmetal 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 reinserted if desired.
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 “denormalize” the sinogram. The denormalization 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.
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 blockiterative 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 (BISART) algorithm is then given by
(6)  
(7)  
The subscript indicates that only rows of and corresponding to the measurements in are used, including when forming the matrices and . The operator ensures nonnegativity of the solution after every iteration.
The BISART algorithm accounts for neither the statistical uncertainty in the measured data, nor the polyenergetic nature of the Xray measurements. Additionally, images reconstructed by BISART may be of poor quality when the data are noisy. In Ref. HG17 , we introduced three enhancements to the BISART 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
(8) 
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
(9) 
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
(11) 
We can then apply BISART 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
(12) 
where ^{1}^{1}1Matrix 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
(13) 
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 http://math.haifa.ac.il/YAIR/bibsuperiorizationcensor.html.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 sparseview and limitedangle 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:
(14) 
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 wPSARTTV (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.
ii.4 Prior image constrained superiorized algorithm
Our prior image constrained superiorized algorithm (Algorithm 5), which we denote as wPSARTPICS, 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:
(15) 
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.
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 NMARbased 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 wPSARTPICS 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 builtin 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 SheppLogan phantom) in the center of the object. The lowconstrast 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 parallelbeam 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 wPSARTTV 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 lowcontrast features. Pixel values within the ROI were rescaled to the interval using the formula
(16) 
where and were set to of the background soft tissue value. The peak signaltonoise ratio (PSNR) of the region was then calculated as
(17)  
(18) 
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 wPSARTTV 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 lowcontrast features. Images E–H, reconstructed using pSARTPICS with , demonstrate better definition about the features.
The PSNR values given in Table 1 corroborate these observations. At all four noise levels, the best PSNR values are obtained using wPSARTPICS 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 
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 Archive^{2}^{2}2https://www.cancerimagingarchive.net/. 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 pSARTPICS 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 wPSARTTV 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 wPSARTPICS (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.
It is apparent that the wPSARTPICS 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 Xrays 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 NMARbased 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 (wPSARTTV) HG17 . Numerical experiments modeling several different anatomical scenarios were perfomed, and indicate that the proposed algorithm (wPSARTPICS) is able improve on both NMAR and wPSARTTV. 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 wPSARTTV.
While the wPSARTPICS 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 Xrays. Additionally, we note that the simulated sinogram data for the pelvis phantom also included measurements along which no photons were detected; nonetheless, the wPSARTPICS 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 wPSARTPICS 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.
Acknowledgments
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.
References
 [1] Bruno De Man, Johan Nuyts, Patrick Dupont, Guy Marchal, and Paul Suetens. Metal streak artifacts in Xray 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 interpolationbased 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 tissueclass 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 projectionbased metal artifact reduction method for flatdetector 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 xray 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 Xray 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. Xray CT highdensity 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. Xray 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. Xray 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 xray 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 Xray 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 maximumlikelihood 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 blockiterative 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 variationbased 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 xray 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. Priorguided metal artifact reduction for iterative xray 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 Xray 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 Xray 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. Blockiterative 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 sparseview and limitedangle polyenergetic data. Phys. Med. Biol., 62(16):6762, 2017.
 [47] P. Sukovic and N.H. Clinthorne. Penalized weighted leastsquares image reconstruction for dual energy Xray 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. Timeresolved cardiac interventional conebeam 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. http://web.eecs.umich.edu/fessler/code/.
Comments
There are no comments yet.