# A Re-weighted Joint Spatial-Radon Domain CT Image Reconstruction Model for Metal Artifact Reduction

###### Abstract

High density implants such as metals often lead to serious artifacts in the reconstructed CT images which hampers the accuracy of image based diagnosis and treatment planning. In this paper, we propose a novel wavelet frame based CT image reconstruction model to reduce metal artifacts. This model is built on a joint spatial and Radon (projection) domain (JSR) image reconstruction framework with a built-in weighting and re-weighting mechanism in Radon domain to repair degraded projection data. The new weighting strategy used in the proposed model makes the regularization in Radon domain by wavelet frame transform more effective. The proposed model, which will be referred to as the re-weighted JSR model, combines the ideas of the recently proposed wavelet frame based JSR model [20] and the normalized metal artifact reduction model [38], and manages to achieve noticeably better CT reconstruction quality than both methods. To solve the proposed re-weighted JSR model, an efficient alternative iteration algorithm is proposed with guaranteed convergence. Numerical experiments on both simulated and real CT image data demonstrate the effectiveness of the re-weighted JSR model and its advantage over some of the state-of-the-art methods.

Keywords: Computerized tomography, metal artifact reduction, tight wavelet frame, joint spatial-Radon domain reconstruction.

## 1 Introduction

X-ray based computerized tomography (CT) is one of the most important imaging modalities in diagnosis and treatment planning of various diseases such as cancer. Although magnetic resonance imaging (MRI) is a safer alternative to CT, the acquisition time of MRI is normally much longer than CT, and it is prohibited for patients with metal implants. However, metal implants often lead to serious degradations in the reconstructed CT images. This is mainly due to the commonly adopted assumption that the CT imaging process is linear with monochromatic X-ray source. This assumption works well when the imaged region contains similar types of materials such as soft tissue and muscles which have similar attenuation coefficient [30]. In practice, however, the X-ray source available is multi-chromatic and the detected photons are only the mean energy after attenuation. When the X-ray beam passes through high density materials such as bones and metals, the low energy components of the X-ray are attenuated more than high energy components. As a result, the reconstructed CT image by a simple model such as the widely used filtered backprojection will suffer serious degradations [5]. Therefore, more careful design of the CT image reconstruction model is necessary to suppress metal artifacts. In this paper, we focus on metal artifact reduction in 2-dimensional (2D) CT image reconstruction. Generalization to 3D is straightforward in theory, though requires further numerical studies.

Based on the Lambert-Beer Law, the measured photon intensity from multi-chromatic energy X-ray can be described by

(1.1) |

where denotes the imaging process, with being the incident photons of X-ray at energy level , and is the linear attenuation coefficient at energy level . In this case, the imaging process is nonlinear if the photon intensity is not constant and is varied at each energy level . In practice, the exact at each energy level is not readily available and difficult to measure, and solving a huge nonlinear system is overwhelmingly challenging. Therefore, monochromatic energy is a common assumption adopted in the literature. In the monochromatic energy X-ray imaging, model (1.1) is simplified to

(1.2) |

where is the measured projection data with obtained from (1.1). However, since (1.2) is merely an approximation of (1.1), when high density components such as metal implants are present, simple models solving the linear inverse problem (1.2) may lead to severe metal artifacts in the reconstructed images. The regions corresponding to metals in the Radon domain, which are referred to as the metal trace, are nonlinearly degraded. Most existing metal artifact reduction methods introduce different techniques to repair the degraded projection data or/and to suppress metal artifacts in spatial domain using regularization.

Classical CT image reconstruction algorithms include filtered backprojection(FBP) [33], FDK [25] and Simultaneous Algebraic Reconstruction Technique(SART) [1, 49]. Since the problem (1.2) is an ill-posed linear inverse problem, when metals are present in the reconstructed images, images reconstructed by these classical approaches may be seriously degraded. In the past three decades, different strategies for example, hardware filtering, dual energy CT (DECT) [17, 3] and statistical reconstruction methods are combined with the aforementioned classical reconstruction algorithms to reduce the metal artifacts. Nonetheless, suppressing metal artifact in CT images remains a great challenge.

In image processing, regularization based approach has been playing an important role, since most of the image processing tasks, such as super-resolution, deblurring, inpanting, etc., are ill-posed inverse problems. Total variation (TV) [44, 2, 11] and wavelet frame based approach [9, 10, 6, 21] are two successful examples among others. In low-dose CT imaging, the iterative reconstruction model with TV regularization [47, 41, 32] and wavelet frame based regularization [48, 36, 31, 54, 58, 20, 52, 13] were proposed and achieved better visual quality of the reconstructed images than traditional methods. However, these methods are effective in suppressing artifacts caused by noise and lack of projection data, while they are less effective in suppressing metal artifacts.

Recently in [38], the authors proposed a normalized metal artifact reduction (NMAR) model which outperforms classical reconstruction methods such as FBP. It first normalizes the measured projection data using the projection of a segmented image which is obtained from thresholding of a roughly reconstructed image. Then, interpolation is conducted on the normalized projection data to repair the degradations caused by metals. Finally, it reconstructs the CT image from the repaired projection data using FBP. However, the NMAR model only works well when noise level of the projection data is low. In [37], a nonconvex CT image metal artifact reduction model with sparsity based interpolation in wavelet domain was shown to outperform the NMAR model. But the convergence of this model is not guaranteed since the -norm was used to promote sparsity. A TV regularization based metal artifact reduction model was proposed in [53] and it outperforms the NMAR method on simulated and real data set at the presence of noise. Other TV based metal artifact reduction models can be found in [56, 55] and wavelet frame based model in [57]. Note that, the TV based model works well for piecewise constant images, whereas textures in the image may be smeared out. The prior image constrained iterative compressed sensing(PICCS) reconstruction technique is able to reduce metal artifacts [4] by estimating the removed metal degraded projection data with prior image constrains, which also can be used to reduce the limited angle artifacts [12].

Another strategy to reduce the metal artifact is based on the analysis of the geometric property of the metal trace in Radon domain, which does not require any interpolation. Metal artifacts in the reconstructed CT image is characterized mathematically by the wavefront set through microlocal analysis of the projection data [39]. Based on such characterization, a beam hardening corrector is proposed in [40] to reduce the metal artifacts. However, this model requires the knowledge of the X-ray energy spectrum which is generally hard to acquire in practice. One may refer to [27] for a systematic review of the metal artifact reduction techniques in CT image reconstruction.

In this work, we propose a new model for metal artifact reduction in CT image reconstruction. This model, which will be referred to as the re-weighted JSR model, has a weighting and re-weighting mechanism naturally embedded in a framework of joint Radon domain inpainting and spatial domain regularization. The re-weighted JSR model is composed of different terms specialized in different tasks, and yet highly collaborative with each other. In one of the terms, weights are applied to the projection data, which are calculated from the projection of a segmented CT image. After weighting the projection data, the inpainting in Radon domain becomes easier since the weighted projection data has sparser representation in the wavelet frame domain than the original projection data. The effect of weighting is cancelled out in another term of the re-weighted JSR model to ensure a correct reconstruction in spatial domain. Such weighting strategy also ensures that the corrected projection data and the unknown CT image obey the linear model (1.2) better than the original projection data. Regularization by wavelet frame transforms is applied in both spatial and Radon domain to suppress noise and preserve features, and to facilitate a stable image reconstruction.

The rest of the paper is organized as follows. In Section 2, we briefly review the concept of wavelet frames and wavelet frame transforms, followed by a detailed description of how metal trace and weights that will be needed in the re-weighted JSR model are estimated. Section 3 introduces the re-weighted JSR model and an efficient optimization algorithm solving the proposed model. Numerical experiments are presented in Section 4 for phantoms and in Section 4.2 for a real scan data. Finally, we summarize our contributions and findings in Section 5.

## 2 Preliminaries

### 2.1 Tight wavelet frames

In this subsection, we briefly recall some of the basics of tight wavelet frames that will be used in later sections. Tight wavelet frame systems are redundant systems. Their redundancy not only grants robustness to the system, but also grants vast flexibility in designing frame systems satisfying properties that are desirable in various applications (e.g. short support, symmetry and high order of vanishing moments). Together with their multiscale structure, tight wavelet frame systems can robustly decompose piecewise smooth functions (such as images) into smooth and sparse components, which is the key to their success in image restoration. The interested readers should consult [42, 43, 15, 16] for theories of wavelet frames, [46, 22] for a short survey on the theory and applications of wavelet frames, and [21] for a more detailed survey.

Given a set of compactly supported functions , with , the quasi-affine system generated by is defined by the collection of dilations and shifts of as

(2.1) |

where is defined by

(2.2) |

Then, is called a tight wavelet frame of if

In the tight wavelet frame system , each generator function is called a framelet.

The constructions of compactly supported and desirably (anti)symmetric framelets are usually based on the multiresolution analysis (MRA) generated by some refinable function with refinement mask satisfying

The idea of an MRA-based construction of framelets is to find masks , which are finite sequences (or filters), such that

(2.3) |

The sequences are called wavelet frame masks, or the high pass filters associated to the tight wavelet frame system, and is also known as the low pass filter.

The unitary extension principle (UEP) [42] provides a rather general characterization of MRA-based tight wavelet frames. Roughly speaking, as long as are finitely supported and their Fourier series satisfy

(2.4) |

for all and , the quasi-affine system with defined by (2.3) forms a tight frame of . Note that the quasi-affine system is shift-invariant and it corresponds to the undecimated wavelet (frame) transform [42, 14], which performs better than the traditional wavelet system in image restoration. Furthermore, undecimated wavelet frame transform has a more natural link to differential operators in the framework of variational and PDEs models [6, 7, 19, 23].

###### Example 1.

###### Example 2.

In the discrete setting, the -level wavelet frame transform/decomposition with filter banks is defined by

The wavelet frame coefficients of image are computed by where denotes the convolution operator with a certain boundary condition, e.g., periodic boundary condition, and is defined by

(2.7) |

Denote the inverse wavelet frame transform (or wavelet frame reconstruction operator) as , which is the adjoint operator of , then the UEP leads to a perfect reconstruction formula as

For the 2D image processing, fast decomposition/reconstruction algorithms are available for the quasi-affine system , which is constructed by the tensor product from the univariate wavelet frame. We briefly recall how a set of 2D filters can be constructed from a given set of 1D filters by tensor product. Given a set of univariate masks , define the 2D masks as

It is known that if the set of 1D filters satisfies the UEP conditions, so does the corresponding set of 2D filters constructed by tensor product [21].

Among many different choices of framelets, the ones constructed from B-splines are popular in image processing. This is mainly due to their short supports and symmetry which are desirable in many applications. A tight frame system constructed from a low order B-spline has less filters and each filter has shorter support than the tight frame systems constructed from high order B-splines. Thus, low order B-spline framelets are more computationally efficient to use than high order B-spline framelets, whereas the latter can capture richer image features than the former. Furthermore, since high order B-spline framelets have larger supports, they may introduce more numerical viscosity which often lead to smoother reconstructions in image restoration tasks such as denoising, deblurring, inpainting, etc. Therefore, the choice of B-spline framelets really depends on the task and the computation load one can afford. In all the models we use in the following sections, we choose the Haar framelets for spatial domain regularization and the piecewise cubic B-pline framelets for projection domain regularization. The reason for such choice is simply because most CT images can be well approximated by piecewise constant functions and their corresponding projection images have higher regularity. Throughout the rest of the paper, we fix the level of decomposition to . We finally note that, choices of (e.g. choice of framelets and levels of decomposition) indeed affect the reconstruction. For example, using a data-driven tight frame can generate better reconstructions than B-spline framelets in sparse view [52] and limited view [13] CT reconstruction. However, we forgo a more detailed discussion on the choices of in order not to dilute the main focus of this paper.

### 2.2 Initialization

The proposed re-weighted JSR model requires a pre-estimation of the metal trace and weights in projection domain. They can be obtained fairly easily from a roughly reconstructed CT image using a simple reconstruction model. In this paper, we use the tight wavelet frame based analysis model [8]. This subsection describes the details on how metal trace and weights are computed using the NURBS-based cardiac-torso(NCAT) phantom [45]. Two metal components are implanted in the NCAT phantom as shown in Figure 1(a) and the simulated projection data is obtained from a multi-chromatic X-ray source. Details on the settings of the imaging system are postponed to Section 4.1.1.

#### 2.2.1 Metal trace estimation

Given the measured projection data from the multi-chromatic X-ray source imaging system (1.2) and the projection operator , the unknown CT image can be approximately reconstructed by the following tight wavelet frame based analysis model

(2.8) |

where is the tight wavelet frame transform as reviewed in the previous subsection, and is the norm. The second term of (2.8) is defined by

The definition of the norm was first introduced by [6], whose corresponding proximal operator is the isotropic soft shrinkage (3.5).

The optimization problem (2.8) can be solved by the split Bregman algorithm [29, 8] efficiently, which is also equivalent to the alternating direction method of multipliers (ADMM) [24, 26, 28]. The reconstructed phantom image by model (2.8), denoted by , is shown in Figure 1(b). Metal location in image domain can be robustly estimated by the summation of the high frequency wavelet frame coefficients (Figure 1(c)) followed by a simple thresholding. Then, the index of the metal trace in Radon domain, denoted by , can be identified by the projection of the indicator function associated to the metal location (Figure 1(d)).

Note that one may estimate the metal location by simply thresholding the initially reconstructed image. However, the metal artifacts may have a significant influence on the estimation if the threshold is not properly chosen. Thanks to the multiscale structure of the wavelet frame transform, we are able to robustly detect features from poorly reconstructed images based on the summation of high frequency tight framelet coefficients. This has already been observed in the past [18, 7]. In Table 1, we show that the quality of the reconstructed image using the proposed re-weighted JSR model is not very sensitive to the choice of the threshold (denoted by ) on the summation of high frequency tight framelet coefficients. Furthermore, the proposed approach is better than directly thresholding on the initially reconstructed image.

Thresholding wavelet | Threshold | SSIM | RelErr |
---|---|---|---|

frame coefficients | 0.9826 | 0.0835 | |

(proposed) | 0.9835 | 0.0842 | |

0.9835 | 0.0845 | ||

0.9839 | 0.0860 | ||

0.9830 | 0.0877 | ||

0.9821 | 0.0887 | ||

Thresholding initially | 0.9758 | 0.0873 | |

reconstructed image | 0.9816 | 0.0883 | |

0.9824 | 0.0892 | ||

0.9823 | 0.0894 | ||

0.9824 | 0.0890 | ||

0.9825 | 0.0893 |

#### 2.2.2 Weights estimation

To suppress the effect of data inconsistency in the Radon domain near the metal trace , we will weight the projection data in the proposed re-weighted JSR model. This subsection describes how the weights are estimated.

Given the metal trace obtained from the previous subsection, we first solve the following tight wavelet frame based analysis model

(2.9) |

where is the restriction operator with respect to the domain . The restriction operator is defined as: for ; and for . The optimization problem (2.9) can be solved similar to (2.8) by the split Bregman algorithm. We denote the solution of (2.9) as (Figure 2(a)).

The weights will be computed by projecting a segmented image that approximates tissue classification of the unknown CT image. One may obtain such approximated tissue classification by segmenting image obtained from (2.8) or obtained from (2.9). However, the approximation from either image will be rather inaccurate since has severe artifacts in between metal locations (Figure 1(b)), while the metal components are missing from (Figure 2(a)) though there are less artifacts in between metal locations. Therefore, we propose to segment a combined image defined by

(2.10) |

with a tuning parameter. In this paper, the segmentation of is obtained by the algorithm proposed by [34, 35]. The segmented image, denoted as (Figure 2(b)), contains three components: the air, the low density components such as soft tissues, and the high density components such as bones and metals. The intensity values of the segmented image from prior image are assumed to be constant for each segmented component. We propose to use the mean values of in the segmented regions as the constants. After obtaining , the weight that will be used in our re-weighted JSR model is defined by . Finally, we note that the reconstructed image using the re-weighted JSR model is relatively insensitive to the choice of the parameter in (2.10), which is demonstrated in Table 2 using NCAT phantom.

SSIM | RelErr | |
---|---|---|

0.9611 | 0.1127 | |

0.9678 | 0.1020 | |

0.9764 | 0.0918 | |

0.9764 | 0.0918 | |

0.9810 | 0.0938 | |

0.9794 | 0.1003 |

## 3 A re-weighted joint spatial-Radon domain reconstruction for metal artifact reduction

In this section, we propose a new re-weighted joint spatial-Radon domain reconstruction (re-weighted JSR) model to reconstruct CT images from the multi-chromatic X-ray imaging system with reduced metal artifacts. Efficient algorithm is adopted to solve the proposed model. Intuition of the weighting strategy is discussed at the end of this section.

### 3.1 The re-weighted JSR model

To reduce metal artifacts and reconstruct high quality CT images, we propose the following re-weighted JSR model

(3.1) |

Here, is the measured projection data from the multi-chromatic X-ray source imaging system which is contaminated by Poisson noise, is the Haar framelet transform (multilevel with ), is the piecewise cubic B-spline framelet transform (multilevel with ), is the restriction operator that extracts the projection data in the complement of metal trace, and . The division is defined point-wisely.

The novelty of the proposed re-weighted JSR model (3.1) is how the weights are introduced to the joint spatial-Radon reconstruction framework. In (3.1), the projection data is weighted by so that the inpainting mechanism induced by the last two terms of (3.1) is more effective. The weighting is cancelled out during the reconstruction of from the first two terms of (3.1) to ensure a correct reconstruction. An empirical explanation of such weighting strategy will be presented in Section 3.2.

A special case of the re-weighed JSR model (3.1) is

(3.2) |

where the weight is removed. The model (3.2) was proposed in [20] for sparse angle CT image reconstruction. Although in principle, model (3.2) can also be applied to metal artifact reduction, it introduces new artifacts in the reconstructed images. For convenience, we refer to the model (3.2) as the unweighted JSR model. In Section 4, we shall compare the unweighted JSR model with the re-weighted JSR model to demonstrate the effectiveness of the proposed weighting strategy.

To solve the re-weighted JSR model, we first rewrite it into a more succinct form. For convenience of notation, we temporarily let and . Denote Then, model (3.1) can be rewritten as

(3.3) |

where the norm is defined block-wisely. Model (3.3) is a standard analysis based model and can be efficiently solved by the split Bregman algorithm [29, 8] with guaranteed convergence. The split Bregman algorithm solving (3.3) takes the following form

(3.4) |

Each subproblem has a closed form solution.

Details of the algorithm is presented in Algorithm 1. Since the sub-variables are coupled together in the -subproblem, they can be updated alternatively. The linear system in the step solving for is solved using the conjugate gradient (CG) algorithm. The thresholding operator , , is the isotropic soft shrinkage operator [6]. Now, we recall the definition of the isotropic soft shrinkage operator. Given the wavelet frame coefficients and thresholds with , the shrinkage operator used in Algorithm 1 is the isotropic shrinkage operator [6] defined by

(3.5) |

where . As convention, we choose .

For a better presentation of the proposed CT image reconstruction method with reduced metal artifacts, we summarize the entire procedure in the flow chart shown in Figure 3.

### 3.2 Why does the weighting strategy work ?

In this subsection, we give some empirical observations on the weighting strategy used in the proposed model (3.1). For simplicity, all the numerical results in this subsection are obtained using the NCAT phantom. The re-weighted JSR model can be viewed as a combination of the following two models: the spatial domain reconstruction model

(3.6) |

and the Radon domain inpainting model

(3.7) |

Obviously, (3.6) is an analysis model with tight wavelet frame regularization. The Radon domain inpainting model (3.7) is another analysis based model aims at repairing the degraded projection data around the metal trace.

We observe that the weighted projection data admits a sparse representation than the measured projection data in wavelet frame domain. The distributions of wavelet frame coefficients (using piecewise linear tight wavelet frame system) are plotted in Figure 4. To have a fair comparison, and are normalized to the same scale, and the coefficients near the metal trace are removed. From the distributions we can see that the weighted projection data has a sparser representation than the original projection data . Therefore, the sparsity prior in the inpainting model (3.7) is more effective than directly using as input.

By combining (3.7) with (3.6), the proposed re-weighted JSR model is able to repair the degraded projection so that the repaired data is closer to , where is the Radon domain solution from (3.1) and is the unknown ground-truth CT image. Note that it is hard to obtain the ground truth image for polychromatic energy CT with varying attenuation coefficients with respect to the energy level, we choose the NCAT phantom (Figure 1(a)) as the approximated ground truth CT image . Empirically, the repaired projection data fits the linear inverse problem (1.2) better than the repaired projection data and computed from the NMAR and the unweighted JSR model (3.2) respectively. This is the key to the success of the re-weighted JSR model since the linear model (1.2) is what we commonly assume for CT imaging. Such linear assumption is not correct (though reasonable) for a multi-chromatic imaging system.

To support such claim, we present comparisons of with , , and in Figure 5 using the NCAT phantom. We observe that is a better approximation to the projection of the reference image than the repaired projection data from the JSR and NMAR model (see Figure 5(c),(d)). The NMAR model also generates a better repaired projection than the unweighted JSR model due to its re-weighting strategy. However, the unweighted JSR model is still able to reduce the majority of the metal artifacts in the reconstructed image due to its sparsity based joint regularization. These observations, together with the reconstruction results in Section 4, show that the re-weighted JSR model combines the merits of the NMAR model’s weighting strategy and the sparsity based joint regularization of the unweighted JSR model.

To quantitatively measure the difference between and the repaired projection data from different models and the measured projection data , we calculate the -norms , , and . Since the region of the metal trace has major contribution to these quantities, we also compute the -norms excluding the regions of the metal trace. Results are shown in Table 3. Obviously, the repaired projection data from the re-weighted JSR model is closer to than that from the NMAR and JSR model. However, although is closer to in regions outside of , is overall closer to than due to the inaccurate recovery of the projection data inside the metal trace by the re-weighted JSR model (see Figure 5(c)). This is probably why the re-weighted JSR model still cannot fully remove metal artifacts, though it improves over the unweighted JSR and the NMAR model.

NCAT | Models | Including | Excluding |

Measured Data | 1547.6 | 713.8 | |

NMAR | 2802.1 | 713.4 | |

JSR | 3048.2 | 709.7 | |

Re-weighted JSR | 2777.5 | 688.3 |

## 4 Numerical experiments

In this section, we validate the effectiveness of the re-weighted JSR model using both simulated phantoms and real data. All numerical experiments are implemented in MATLAB running on a platform with 16 GB RAM and Intel(R) Core(TM) i7-6700T CPU at 2.8-GHz with 4 cores.

The regularization parameter appeared in any of the models presented in the previous sections takes the form where is a tuning parameter. The stopping criterion for Algorithm 1 is

The unweighted JSR model and the TV-FADM model are solved using the algorithms proposed in [20] and [53] respectively with the same stopping criterion as above. We set the maximum allowable number of iteration of all the iterative algorithms to 700. We use the relative error (RelErr) and structural similarity (SSIM) index [51] to quantitative evaluate the reconstructed images from different methods.

### 4.1 Experiments on image phantoms

#### 4.1.1 Experimental design

The fan-beam CT imaging system similar to [50] is chosen for all the simulations in this subsection. The source to detector distance is 949.075 mm, the distance from the source to the iso-center is 541 mm and the strip width is 1.024 mm. The source trajectory covers the full circular orbit of with 984 angular views and the number of bins per view is 888. The tube potential is 140kvp with 2.5mm aluminum and 0.5mm copper filters.

The NCAT phantom (shown in Figure 1(a)) and the cerebral phantom^{1}^{1}1http://see.xidian.edu.cn/vipsl/database_CTMR.html (shown in Figure 8(a)) are chosen as image phantoms. For the NCAT phantom, it has pixels. Two metal components (Titanium) are implanted in the image, which is shown in Figure 1(a) with red curves labeling the locations of the metals. For the cerebral phantom, it has pixels and three metal components (Titanium) are implanted. Both of the phantoms contain three major components, i.e. soft tissue, bone and metal components, and their linear attenuation coefficients can be found in [30].

The projection data is obtained from a multi-chromatic X-ray imaging system definition by (1.1) with the energy spectrum shown in Figure 6. The measured projection data contaminated by Poisson noise is generated in the following way

(4.1) |

where is used to add the Poisson noise, is the incident photons’ number and the term is used to replace the pixel value after adding Poisson noise. In this subsection, we select in (4.1).

#### 4.1.2 Reconstruction results

Table 4 is the parameter setting of re-weighted JSR model (3.1) and unweighed JSR model (3.2) for both phantoms. The analysis model (2.8), the inpainting model (2.9), the unweighted JSR model (3.2) and the re-weighted JSR model (3.1) are all initialized with . The number of CG iteration in all the algorithms involved is set to . We compare the re-weighted JSR model (3.1) with the classical FBP algorithm [33], the unweighted JSR model (3.2), NMAR [38] and TV-FADM [53]. The codes of TV-FADM are provided by the authors of [53]. Comparisons of the aforementioned methods using the NCAT and cerebral phantoms are shown in Figure 7 and Figure 8, respectively.

Figure 7(a) shows that the reconstructed NCAT phantom from FBP has severe metal artifacts and is noisy. The reconstructed image from the unweighted JSR model (3.2) shown in Figure 7(b) has a better visual effect with noticeably less noise and metal artifacts. Sharp edges are also well preserved except for the blurry effects in the region surrounding the metals. The reconstructed image from NMAR shown in Figure 7(c) also has most of the metal artifacts suppressed and the regions surrounding the metals are much less blurry than the unweighted JSR. However, the unweighted JSR does a better job than NMAR in suppressing noise and preserving sharp image features away from the metals. TV-FADM is able to reconstruct images with minimum metal artifacts and noise, as shown in Figure 7(d). However, the metal components are fused with nearby structures which is highlighted by the red arrow. The reconstructed image from proposed re-weighted JSR model has the best overall quality with rather minor metal artifacts.

Figure 8 shows the reconstructed cerebral phantom from different methods. We highlight some regions with more distinct differences with red contours. Since the cerebral phantom contains more textures, it is more challenging than the NCAT phantom. The pros and cons of these methods are mostly the same as the previous example. However, we note that the reconstructed image from TV-FADM shown in Figure 8(e) has severe artifact, which is due to the well-known staircase artifact of TV regularization. We found that TV-FADM is relatively sensitive to the choice of its parameters. It is not easy to balance between sharpness of image features and metal artifacts reduction. The soft tissue around metal components is also not well preserved by the NMAR method as indicated by the blue arrow in Figure 8(d). Furthermore, the circled areas in Figure 8(d) show that there are still some artifacts around the metal. Same as the NCAT phantom, the proposed re-weighted JSR model has the best overall performance. Notice that the intensity of metals in Figure 8(d) and 8(f) seems lower than the rest of the reconstructed images. This is because we set the intensity of the metal components in the segmentation with the same mean value as that of bones. Increasing the value of metal components of can increase the intensity of metals in the reconstructed images, whereas it also introduces more artifacts around the metals.

Quantitative assessments of the reconstruction quality of these methods are given in Table 5, where the relative errors, SSIM values and computation time are presented. Although Algorithm 1 that solves the re-weighted JSR model is relatively time consuming, we are able to gain on quality of the reconstructed images. Finally, to numerically demonstrate the convergence of the Algorithm 1, we present the decay of , and the cost function of the re-weighted JSR model in Figure 9(a), 9(b) and 9(c) respectively.

### 4.2 Numerical experiments: real data

We perform a CT scan of a chicken leg placed in a disposable cup (Figure 10(a)). We first scan the chicken leg without metals (Figure 10(b)) to create a reference image using FBP algorithm. Then, we place two steel thread nails on each side of the chicken leg and scan the subject again using the same scanning protocol (Figure 10(c)). The projection data is acquired from a MicroCT scanner equipped at the Division of Nuclear Technology and Applications, Institute of High Energy Physics, Chinese Academy of Sciences. The X-ray source is with 90 kV and 70 mA energy and the flat plane detector contains pixels. The scanning trajectory is a full circle with equally spaced views at per view. The physical size of each detector unit is . The distance from the X-ray source to the detector is . In order to conduct a 2D experiment, we choose the 512th row of the detector array.

Figure 11 shows the images reconstructed using FBP, the analysis model (2.8), the inpainting model (2.9) and the segmented image from the image obtained by (2.10). The reference image without metal implants are shown in Figure 11(a). All the images in this subsection are displayed within the grayscale interval . The segmented image shown in Figure 11(e) is used to estimate the weights needed in NMAR and the re-weighted JSR model.

Figure 12 shows a comparison between the reconstructed image from NMAR and the unweighted JSR model. Figure 13 shows a comparison between the reconstructed images from TV-FADM and the proposed re-weighted JSR model. Zoom-in views are provided in both Figure 12 and Figure 13 for a better visual assessment. As one can see that the reconstructed images from the unweighted JSR model and TV-FADM are less noisy than NMAR as indicated by the blue ellipse curve, whereas NMAR does a better job in preserving image features and suppressing metal artifacts. However, there are also new artifacts around the metal on the right as shown in Figure 12(d). The proposed re-weighted JSR model has best overall performance in terms of feature preservation, noise and metal artifact reduction.

## 5 Conclusion

In this paper, we proposed a new model for metal artifact reduction in multi-chromatic X-ray CT imaging. The proposed model had a weighting and re-weighting mechanism naturally embedded in a framework of joint Radon domain inpainting and spatial domain regularization. Regularization by wavelet frame transforms was applied in both spatial and Radon domain to facilitate a high quality and stable image reconstruction. The proposed model was then rewritten into a more compact form so that the split Bregman algorithm could be directly applied to solve the model efficiently. Our numerical experiments using both image phantoms and real data showed that the proposed model was more effective than the classical FBP method [33], the popular NMAR method [38], and the recently proposed TV-FADM method [53]. Comparisons with the unweighted JSR model [20] further illustrated the effectiveness of the weighting and re-weighting mechanism.

## Acknowledgements

We would like to thank the anonymous reviewers for their constructive suggestions and comments that helped tremendously in improving the presentation of this paper. We would also like to thank Professor Qibin Fan from the School of Mathematics and Statistics, Wuhan University for his support on this project.

## References

- [1] A. H. Andersen and A. C. Kak, Simultaneous algebraic reconstruction technique (sart): a superior implementation of the art algorithm, Ultrasonic imaging, 6 (1984), pp. 81–94.
- [2] G. Aubert and P. Kornprobst, Mathematical problems in image processing: partial differential equations and the calculus of variations, Springer, 2006.
- [3] F. Bamberg, A. Dierks, K. Nikolaou, M. F. Reiser, C. R. Becker, and T. R. Johnson, Metal artifact reduction by dual energy computed tomography using monoenergetic extrapolation, European radiology, 21 (2011), pp. 1424–1429.
- [4] P. Bannas, Y. Li, U. Motosugi, K. Li, M. Lubner, G.-H. Chen, and P. J. Pickhardt, Prior image constrained compressed sensing metal artifact reduction (piccs-mar): 2d and 3d image quality improvement with hip prostheses at ct colonography, European radiology, 26 (2016), pp. 2039–2046.
- [5] T. M. Buzug, Computed tomography: from photon statistics to modern cone-beam CT, Springer Science & Business Media, 2008.
- [6] J.-F. Cai, B. Dong, S. Osher, and Z. Shen, Image restoration: total variation, wavelet frames, and beyond, Journal of the American Mathematical Society, 25 (2012), pp. 1033–1089.
- [7] J.-F. Cai, B. Dong, and Z. Shen, Image restoration: a wavelet frame based model for piecewise smooth functions and beyond, Applied and Computational Harmonic Analysis, 41 (2016), pp. 94–138.
- [8] J.-F. Cai, S. Osher, and Z. Shen, Split bregman methods and frame based image restoration, Multiscale modeling and simulation, 8 (2009), pp. 337–369.
- [9] R. H. Chan, T. F. Chan, L. Shen, and Z. Shen, Wavelet algorithms for high-resolution image reconstruction, SIAM Journal on Scientific Computing, 24 (2004), pp. 1408–1432.
- [10] R. H. Chan, S. D. Riemenschneider, L. Shen, and Z. Shen, Tight frame: an efficient way for high-resolution image reconstruction, Applied and Computational Harmonic Analysis, 17 (2004), pp. 91–115.
- [11] T. Chan and J. Shen, Image processing and analysis: variational, PDE, wavelet, and stochastic methods, Society for Industrial Mathematics, 2005.
- [12] G.-H. Chen and Y. Li, Synchronized multiartifact reduction with tomographic reconstruction (smart-recon): A statistical model based iterative image reconstruction method to eliminate limited-view artifacts and to mitigate the temporal-average artifacts in time-resolved ct, Medical physics, 42 (2015), pp. 4698–4707.
- [13] J. K. Choi, B. Dong, and X. Zhang, Limited tomography reconstruction via tight frame and simultaneous sinogram extrapolation, Journal of Computational Mathematics, 34 (2016), pp. 575–589.
- [14] R. Coifman and D. Donoho, Translation-invariant de-noising, Wavelets and statistics, 103 (1995), p. 125.
- [15] I. Daubechies, Ten lectures on wavelets, SIAM, 1992.
- [16] I. Daubechies, B. Han, A. Ron, and Z. Shen, Framelets: Mra-based constructions of wavelet frames, Applied and Computational Harmonic Analysis, 14 (2003), pp. 1–46.
- [17] F. Dilmanian, X. Wu, E. Parsons, B. Ren, J. Kress, T. Button, L. Chapman, J. Coderre, F. Giron, D. Greenberg, et al., Single-and dual-energy ct with monochromatic synchrotron x-rays, Physics in medicine and biology, 42 (1997), p. 371.
- [18] B. Dong, A. Chien, and Z. Shen, Frame based segmentation for medical images, Communications in Mathematical Sciences, 9 (2010), pp. 551–559.
- [19] B. Dong, Q. Jiang, and Z. Shen, Image restoration: wavelet frame shrinkage, nonlinear evolution pdes, and beyond, Multiscale Modeling and Simulation, 15 (2017), pp. 606–660.
- [20] B. Dong, J. Li, and Z. Shen, X-ray ct image reconstruction via wavelet frame based regularization and radon domain inpainting, Journal of Scientific Computing, 54 (2013), pp. 333–349.
- [21] B. Dong and Z. Shen, MRA-Based Wavelet Frames and Applications, IAS Lecture Notes Series, Summer Program on “The Mathematics of Image Processing”, Park City Mathematics Institute, (2010).
- [22] B. Dong and Z. Shen, Image restoration: a data-driven perspective, in Proceedings of the International Congress of Industrial and Applied Mathematics (ICIAM), 2015, pp. 65–108.
- [23] B. Dong, Z. Shen, and P. Xie, Image restoration: a general wavelet frame based model and its asymptotic analysis, SIAM Journal on Mathematical Analysis, 49 (2017), pp. 421–445.
- [24] E. Esser, Applications of lagrangian-based alternating direction methods and connections to split bregman, CAM Report 09-31, Department of Mathematics, University of California Los Angeles, (2009).
- [25] L. Feldkamp, L. Davis, and J. Kress, Practical cone-beam algorithm, JOSA A, 1 (1984), pp. 612–619.
- [26] D. Gabay and B. Mercier, A dual algorithm for the solution of nonlinear variational problems via finite element approximation, Computers and Mathematics with Applications, 2 (1976), pp. 17–40.
- [27] L. Gjesteby, B. De Man, Y. Jin, H. Paganetti, J. Verburg, D. Giantsoudi, and G. Wang, Metal artifact reduction in ct: Where are we after four decades?, IEEE Access, (2016).
- [28] R. Glowinski and A. Marrocco, Sur l’approximation par éléments finis d’ordre un, et la résolution par pénalisation-dualité, d’une classe de problèmes de dirichlet non linéares, Journal of Equine Veterinary Science, 2 (1975), pp. 41–76.
- [29] T. Goldstein and S. Osher, The split bregman method for -regularized problems, SIAM journal on imaging sciences, 2 (2009), pp. 323–343.
- [30] J. H. Hubbell and S. M. Seltzer, Tables of x-ray mass attenuation coefficients and mass energy-absorption coefficients 1 kev to 20 mev for elements z= 1 to 92 and 48 additional substances of dosimetric interest, tech. report, National Inst. of Standards and Technology-PL, Gaithersburg, MD (United States). Ionizing Radiation Div., 1995.
- [31] X. Jia, B. Dong, Y. Lou, and S. B. Jiang, Gpu-based iterative cone-beam ct reconstruction using tight frame regularization., Physics in Medicine and Biology, 56 (2011), pp. 3787–807.
- [32] X. Jia, Y. Lou, R. Li, W. Y. Song, and S. B. Jiang, Gpu-based fast cone beam ct reconstruction from undersampled and noisy projection data via total variation, Medical physics, 37 (2010), pp. 1757–1760.
- [33] A. Katsevich, Theoretically exact filtered backprojection-type inversion algorithm for spiral ct, SIAM Journal on Applied Mathematics, 62 (2002), pp. 2012–2026.
- [34] C. Li, R. Huang, Z. Ding, J. C. Gatenby, D. N. Metaxas, and J. C. Gore, A level set method for image segmentation in the presence of intensity inhomogeneities with application to mri, IEEE Transactions on Image Processing, 20 (2011), pp. 2007–2016.
- [35] C. Li, C. Xu, C. Gui, and M. D. Fox, Distance regularized level set evolution and its application to image segmentation, IEEE transactions on image processing, 19 (2010), pp. 3243–3254.
- [36] Y. Liu, J. Ma, Y. Fan, and Z. Liang, Adaptive-weighted total variation minimization for sparse data toward low-dose x-ray computed tomography image reconstruction, Physics in medicine and biology, 57 (2012), p. 7923.
- [37] A. Mehranian, M. R. Ay, A. Rahmim, and H. Zaidi, X-ray ct metal artifact reduction using wavelet domain sparse regularization, IEEE transactions on medical imaging, 32 (2013), pp. 1707–1722.
- [38] E. Meyer, R. Raupach, M. Lell, B. Schmidt, and M. Kachelrieß, Normalized metal artifact reduction (nmar) in computed tomography, Medical physics, 37 (2010), pp. 5482–5493.
- [39] H. S. Park, J. K. Choi, and J. K. Seo, Characterization of metal artifacts in x-ray computed tomography, Communications on Pure and Applied Mathematics, (2017). doi:10.1002/cpa.21680.
- [40] H. S. Park, D. Hwang, and J. K. Seo, Metal artifact reduction for polychromatic x-ray ct based on a beam-hardening corrector, IEEE transactions on medical imaging, 35 (2016), pp. 480–487.
- [41] L. Ritschl, F. Bergner, C. Fleischmann, and M. Kachelrieß, Improved total variation-based ct image reconstruction applied to clinical data, Physics in medicine and biology, 56 (2011), p. 1545.
- [42] A. Ron and Z. Shen, Affine systems in : The analysis of the analysis operator, Journal of Functional Analysis, 148 (1997), pp. 408–447.
- [43] A. Ron and Z. Shen, Affine systems in ii: Dual systems, Journal of Fourier Analysis and applications, 3 (1997), pp. 617–637.
- [44] L. I. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear Phenomena, 60 (1992), pp. 259–268.
- [45] W. P. Segars, M. Mahesh, T. J. Beck, E. C. Frey, and B. M. Tsui, Realistic ct simulation using the 4d xcat phantom, Medical physics, 35 (2008), pp. 3800–3808.
- [46] Z. Shen, Wavelet frames and image restorations, in Proceedings of the International congress of Mathematicians, vol. 4, Hindustan Book Agency, 2010, pp. 2834–2863.
- [47] E. Y. Sidky and X. Pan, Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization, Physics in medicine and biology, 53 (2008), p. 4777.
- [48] Z. Tian, X. Jia, K. Yuan, T. Pan, and S. B. Jiang, Low-dose ct reconstruction via edge-preserving total variation regularization, Physics in medicine and biology, 56 (2011), p. 5949.
- [49] G. Wang, D. L. Snyder, J. A. O’Sullivan, and M. W. Vannier, Iterative deblurring for ct metal artifact reduction, IEEE transactions on medical imaging, 15 (1996), pp. 657–664.
- [50] J. Wang, T. Li, H. Lu, and Z. Liang, Penalized weighted least-squares approach to sinogram noise reduction and image reconstruction for low-dose x-ray computed tomography, IEEE transactions on medical imaging, 25 (2006), pp. 1272–1283.
- [51] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, Image quality assessment: from error visibility to structural similarity, IEEE transactions on image processing, 13 (2004), pp. 600–612.
- [52] R. Zhan and B. Dong, Ct image reconstruction by spatial-radon domain data-driven tight frame regularization, SIAM Journal on Imaging Sciences, 9 (2016), pp. 1063–1083, https://doi.org/10.1137/16M105928X.
- [53] H. Zhang, L. Wang, L. Li, A. Cai, G. Hu, and B. Yan, Iterative metal artifact reduction for x-ray computed tomography using unmatched projector/backprojector pairs, Medical Physics, 43 (2016), pp. 3019–3033.
- [54] Y. Zhang, B. Dong, and Z. Lu, minimization for wavelet frame based image restoration, Mathematics of Computation, 82 (2013), pp. 995–1015.
- [55] Y. Zhang, Y.-F. Pu, J.-R. Hu, Y. Liu, and J.-L. Zhou, A new ct metal artifacts reduction algorithm based on fractional-order sinogram inpainting, Journal of X-ray science and technology, 19 (2011), pp. 373–384.
- [56] Y. Zhang, H. Yan, X. Jia, J. Yang, S. B. Jiang, and X. Mou, A hybrid metal artifact reduction algorithm for x-ray ct, Medical physics, 40 (2013), p. 041910.
- [57] S. Zhao, D. Robeltson, G. Wang, B. Whiting, and K. T. Bae, X-ray ct metal artifact reduction using wavelets: an application for imaging total hip prostheses, IEEE transactions on medical imaging, 19 (2000), pp. 1238–1247.
- [58] W. Zhou, J.-F. Cai, and H. Gao, Adaptive tight frame based medical image reconstruction: a proof-of-concept study for computed tomography, Inverse problems, 29 (2013), p. 125006.