Magnetic resonance-based computed tomography metal artifact reduction using Bayesian modelling

Medical x-ray computed tomography (CT) images of patients with metal implants often display major corruption from streak artifacts (Bal and Spies 2006, Jäkel and Reiss 2007), which affects both the visual quality of the images and the quantitative CT value accuracy. The latter is a potential hazard in radiotherapy (RT), where the CT values are used in treatment planning to provide electron density and relative stopping power estimates (Andersson et al 2014, Giantsoudi et al 2017). This is of particular concern in head-and-neck RT, where dental implants and fillings occur frequently and are simultaneously close to both the tumor site and critical organs. In this situation, metal artifact reduction (MAR) plays an important role in error management (Kim and Tomé 2007, Dietlicher et al 2014, Kwon et al 2015, Maerz et al 2015).

MAR is in general a difficult problem, as demonstrated by its approximately 40 year long history that has spawned numerous algorithmic approaches (Koehler et al 2011, Andersson et al 2014, Gjesteby et al 2016). A reason for this difficulty is the origin of the artifacts, which arise from multiple contributions that are amplified in the presence of metal (Beister et al 2012, Gjesteby et al 2016). Some of these contributions stem from incorrect assumptions in the CT reconstruction model that relates the image coefficients to the x-ray projection data (sinogram) (Goldman 2007, Buzug 2008). In particular, the mono-energetic approximation of the x-ray source spectrum in, e.g. the wide-spread filtered back projection (FBP) algorithm leads to incorrect modelling of the projections that are acquired through metal, and thus beam hardening artifacts (Stonestrom et al 1981, Hsieh et al 2000). Other contributions are more model-independent, such as the photon starvation in the metal projections that leads to noise artifacts (Nuyts et al 1998, Gjesteby et al 2016).

MAR algorithms may be categorized into three overall approaches. Image inpainting algorithms replace corrupted CT values with better estimates by post-processing already reconstructed images (Lell et al 2013). Sinogram inpainting algorithms, on the other hand, replace metal projections by estimates that fit the reconstruction model to a higher degree, which may be particularly effective in dealing with photon starvation (Kalender et al 1987, Meyer et al 2009, Zhang et al 2011, Lell et al 2012, Meyer et al 2012). Finally, model based iterative reconstruction (MBIR) algorithms change the CT reconstruction model itself to a more complex probabilistic forward model that better accounts for the artifact sources, at the cost of having to optimize a generally non-linear image functional in a slow, iterative algorithm (De Man et al 2000, 2001, Elbakri and Fessler 2002, Thibault et al 2007, Hamelin et al 2008, Stayman et al 2012, Van Slambrouck and Nuyts 2012, Nuyts et al 2013, Tilley et al 2016, Jin et al 2015, Stille et al 2016, Fu et al 2017).

An important part of many MAR algorithms is the inclusion of prior information about the image that is being estimated. While this is especially true for image inpainting algorithms, which impose such information directly in image space, prior information is also used in some of the most successful sinogram inpainting algorithms. An example is the 'normalized MAR' (nMAR) method (Meyer et al 2009, 2012), which replaces metal projections by (scaled) projection estimates simulated from a template image. MBIR models also include prior information as they consist of two parts: a sinogram data likelihood that addresses, e.g. noise and beam hardening artifacts by modelling the detector noise, the x-ray source spectrum and the implant material; and an image prior distribution that may be used to guide the reconstruction with statistical knowledge about the image being reconstructed (De Man et al 2001, Elbakri and Fessler 2002, Hamelin et al 2008, Van Slambrouck and Nuyts 2012, Nuyts et al 2013, Tilley et al 2016, Fu et al 2017).

In a general CT setting, the quality of the available prior information is rather limited. In sinogram inpainting, for example, the required template is typically generated by post-processing a CT image reconstruction that is heavily corrupted by metal artifacts (Meyer et al 2009, Lell et al 2012, Meyer et al 2012). Being generated from corrupted data, the resulting prior information may itself be compromised, thereby introducing new artifacts in complex, highly corrupted regions such as the head and neck near the teeth and oral cavity. Similarly, in MBIR the limited availability of accurate prior information often motivates the use of relatively simple functional priors that merely impose mathematical regularities in the reconstructed image (De Man et al 2000, Webster Stayman and Fessler 2001, Pan et al 2010, Makeev and Glick 2013).

In the specific context of RT, a potential source of prior information for improved MAR is the magnetic resonance (MR) image that is commonly co-acquired for tumor- and normal tissue delineation. Since metal artifacts are often more localized in MR compared to CT (Park et al 2015, Bell and Bashir 2018), the MR scan can provide superior anatomical prior information in regions that are heavily corrupted in CT (Jonsson et al 2010, Bell and Bashir 2018). Furthermore, since co-registration and acquisition of the MR and CT scans in the same patient fixation is already part of the tumor-delineation process, this can be done with minimal interruption to the existing clinical workflow.

Despite the obvious potential, to the best of our knowledge only a few prior attempts have been made to use MR for reducing CT metal artifacts, most notably the image inpainting algorithms described in Anderla et al (2013), Delso et al (2013), Park et al (2015). The principal difficulty faced by such MR-based approaches to CT MAR lies in the image contrast disparities between the two modalities, especially between bone and air, which are easily distinguishable in the CT but not in the MR scan unless dedicated sequences are used (Delso et al 2013). Additional difficulties arise from co-registration issues, mainly due to to inter-acquisitional motion, which limit the accuracy of MR-based predictions.

In order to overcome these difficulties, we propose a novel Bayesian algorithm to predict uncorrupted CT values from a combination of corrupted CT measurements and corresponding MR intensities. The contributions of this work can be summarized as follows:

• 1.  
  Due to the contrast disparity between MR and CT, single-voxel MR intensities are poor CT value predictors. In contrast to prior work, our model therefore relies on MR image patches (Andreasen et al 2015, 2016) instead, i.e. collections of MR intensities taken from local spatial neighborhoods that encode higher-order anatomical features.
• 2.  
  In contrast to existing methods that simply replace CT values that are deemed to be corrupt with entirely new values (Anderla et al 2013, Delso et al 2013, Park et al 2015), our method automatically blends MR-based predictions with the original, corrupted CT values. As we demonstrate experimentally, this way of retaining the information of the original CT helps to discern bone from air as well as to address co-registration issues.
• 3.  
  Rather than directly targetting only one particular MAR strategy such as image inpainting (Anderla et al 2013, Delso et al 2013, Park et al 2015), the predictive model developed here can be used flexibly within several widely-used MAR scenarios. In particular, we demonstrate experimentally the benefit of using the proposed model to directly calculate blended CT value replacements and thus perform image inpainting; to generate a template for sinogram inpainting; and to define an image reconstruction prior for MBIR.
• 4.  
  In order to facilitate clinical adoption, our model uses data from only the patient under consideration and automatically tunes its parameters to fit this specific patient. We demonstrate how the model thereby accounts for variations in patient-specific features, such as the severity of the artifact corruption as well as variations in the MR sequence parameters. It also removes the need for external training data of well-aligned MR and CT image pairs that is commonly required for MR-based CT prediction in the literature (Dowling et al 2012, Andreasen et al 2015, Torrado-Carvajal et al 2015, Han 2017, Torrado-Carvajal et al 2019).

The paper is structured as follows. We first derive and discuss the proposed predictive model, cover the automatic parameter estimation and present three MAR algorithms that are based on this model (section 2 'Methods'). Next, (in section 3 'Experiments') we describe the practical implementation of the three MAR algorithms, and benchmark their performance in comparison to conventional alternatives; the results are afterwards presented in sequence for the three MAR algorithms (section 4 'Results'). We end with a final summary of our findings and main conclusions, followed by a discussion of, in particular, the potential of clinical implementation of our methods (section 5 'Conclusion and Discussion').

An early version of our approach, with only qualitative (visual) results on a small set of patients, was previously presented (Nielsen et al 2018). The current manuscript provides substantially more technical derivations and analysis, in particular putting more focus on the self-tuning of the model parameters; proposes a novel sinogram inpainting variant as well as an improved MBIR algorithm; and contains detailed quantitative results on a larger set of patients.

Let a set of voxels covering a patient volume be assigned indices i from the index set $\newcommand{\cl}[1]{\mathcal{#1}} \cl{T}$ ($\newcommand{\cl}[1]{\mathcal{#1}} i \in \cl{T}$), and denote the voxel center locations as $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\cl}[1]{\mathcal{#1}} \{\vc{x}_i\}_{i \in \cl{T}}$. To avoid cluttered notation, we will suppress the index set and use $\newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\e}{{\rm e}} \{\cdot \} \equiv \{\cdot \}_{i \in \cl{T}}$ in the remaining. Consider now a CT volume reconstructed using FBP (Buzug 2008), with per-voxel FBP values {t_i}. Assuming the availability of an MR image that is co-registered with the FBP reconstruction, we extract a set of small cuboidal MR volumes ('patches') centered at $\newcommand{\vc}[1]{{\bf #1}} \{\vc{x}_i \}$, with intensities stacked in the vectors $\newcommand{\vc}[1]{{\bf #1}} \{\vc{m}_i \}$. The values in {t_i} are potentially corrupted by metal artifacts, and our task is to estimate the underlying true CT values, {y _i}, given the observation of $\newcommand{\vc}[1]{{\bf #1}} \{t_i, \vc{m}_i\}$.

To accomplish this, we infer the posterior predictive distribution $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(\{y_i \}|\{t_i, \vc{m}_i \}, \vcs{\sigma})$, where $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \vcs{\sigma}$ denotes the set of model parameters, from a generative, probabilistic model of the corrupted CT values contingent on the observed MR scan. We specifically model the two-step process by which the uncorrupted CT values are first sampled from a prior distribution, assumedly in an independent fashion for each voxel:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle p(\{y_i\}|\{\vc{m}_i\}, \vcs{\sigma}) = \prod_{i \in \cl{T}} p(y_i|\vc{m}_i, \vcs{\sigma}), \nonumber \end{align*}$$

and then corrupted by sampling from a noise model $p(\{t_i\}|\{y_i\})$, again independently for each voxel:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle p(\{t_i\}|\{y_i\}, \vcs{\sigma}) = \prod_{i \in \cl{T}} p(t_i|y_i, \vcs{\sigma}). \nonumber \end{align*}$$

Under this generative model, applying Bayes' rule yields:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle p(\{y_i\}|\{t_i, \vc{m}_i \}, \vcs{\sigma}) = \prod_{i \in \mathcal{T}} p(y_i|t_i, \vc{m}_i, \vcs{\sigma}), \label{eq:posterior} \nonumber \end{align} \tag{ 1 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle {{\mathsf w\mathsf h\mathsf e\mathsf r\mathsf e}} \quad p(y_i|t_i, \vc{m}_i, \vcs{\sigma}) = Z^{-1} p(t_i|y_i, \vcs{\sigma})\,p(y_i|\vc{m}_i,\vcs{\sigma}) \nonumber \end{align} \tag{ 2 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle { {\mathsf w\mathsf i\mathsf t\mathsf h}} \quad Z = \int_{-\infty}^{\infty} p(t_i|y_i, \vcs{\sigma})\,p(y_i|\vc{m}_i, \vcs{\sigma}) dy_i \label{eq:partition}. \nonumber \end{align} \tag{ 3 }$$

The predictive distribution of equation (1) depends on an FBP-dependent likelihood function $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(t_i |y_i, \vcs{\sigma})$ and an MR-determined prior distribution $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(y_i|\vc{m}_i,\vcs{\sigma})$, which we detail below.

2.1.1. Likelihood

The likelihood $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(t_i|y_i, \vcs{\sigma})$ models the distribution of the corrupted CT value given the underlying true CT value. Introducing the maximal noise variance ${\sigma_t^*}^{2}$, we model the artifact corruption as additive Gaussian noise with variance $\sigma_t^{2}$ smoothly decreasing with the distance to the implants, such that it is equal to ${\sigma_t^{*}}^{2}$ within the implants but 0 far away where the FBP is free of artifacts. In particular:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\normd}[3]{\mathcal{N}(#1 | #2, #3)} \newcommand{\norm}[3]{\mathcal{N}(#1 | #2, #3^{-1})} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle p(t_i|y_i, \vcs{\sigma}) = \normd{t_i}{y_i}{\sigma_t^{2}} \quad \mathsf{with} \quad \sigma_t^{2} = {\sigma_t^{*}}^{2} f(\vc{x}_i)\label{eq:likelihood} \nonumber \end{align} \tag{ 4 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle \mathsf{and} \quad f(\vc{x}_i) = 1 + \tanh(-\frac{D^2_{\perp}(\vc{x}_i)}{\kappa})\label{eq:sigmoid}. \nonumber \end{align} \tag{ 5 }$$

Here, $\newcommand{\vc}[1]{{\bf #1}} D_{\perp}(\vc{x}_i)$ is the perpendicular spatial distance to a segmentation of the metal implants; $\newcommand{\norm}[3]{\mathcal{N}(#1 | #2, #3^{-1})} \newcommand{\normd}[3]{\mathcal{N}(#1 | #2, #3)} \normd{\cdot}{\psi}{\nu^2}$ denotes a Gaussian with mean $\psi$ and variance $\nu^2$; and $\kappa$ sets the decrease rate of the sigmoidal function $\newcommand{\vc}[1]{{\bf #1}} f(\vc{x}_i)$. As an illustration, consider the two voxels at different distances to the metal in figure 1(a), for which the noise models (figure 1(b)) have different widths due to the modulation with $\newcommand{\vc}[1]{{\bf #1}} f(\vc{x}_i)$ (figure 1(c)).

Zoom In Zoom Out Reset image size

The metal segmentation and $\kappa$ parameter must be defined by the user; section 3.2 will present the definitions used for our experiments. The maximal artifact noise variance ${\sigma_t^*}^2$ will however be chosen automatically given the observed data, as described in section 2.3.

2.1.2. Prior

Letting $\newcommand{\cl}[1]{\mathcal{#1}} \cl{T}_u \subseteq \cl{T}$ denote an assumed uncorrupted part of the patient volume (see figure 2(a)), we learn the prior distribution $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(y_i|\vc{m}_i, \vcs{\sigma})$ from all matched CT/MR pairs $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\cl}[1]{\mathcal{#1}} \{y_n, \vc{m}_n \}_{n \in \cl{T}_u}$ in $\newcommand{\cl}[1]{\mathcal{#1}} \cl{T}_u$ using kernel regression (Bishop 2006). In particular, we estimate the joint distribution $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(y_i, \vc{m}_i|\vcs{\sigma})$ with kernel density estimation (Bishop 2006) using Gaussian kernels with diagonal covariance matrices:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\normd}[3]{\mathcal{N}(#1 | #2, #3)} \newcommand{\norm}[3]{\mathcal{N}(#1 | #2, #3^{-1})} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle & p(y_i, \vc{m}_i|\vcs{\sigma}) \nonumber \\ & =|\cl{T}_u|^{-1} \sum_{n \in \cl{T}_u} \normd{y_i}{y_n}{\sigma_y^{2}} \normd{\vc{m}_i}{\vc{m}_n}{\sigma_m^{2} \vc{I}_{M}}\label{eq:kde} , \nonumber \end{align} \tag{ 6 }$$

and obtain the prior distribution $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(y_i|\vc{m}_i, \vcs{\sigma})$ accordingly:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle p(y_i|\vc{m}_i, \vcs{\sigma}) = \frac{p(y_i, \vc{m}_i | \vcs{\sigma})}{p(\vc{m}_i | \vcs{\sigma})} = \sum_{n \in \cl{T}_u}w_n^i\mathcal{N}(y_i | y_n, \sigma_y^{2})\label{eq:priorcalced} \nonumber \end{align} \tag{ 7 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\normd}[3]{\mathcal{N}(#1 | #2, #3)} \newcommand{\norm}[3]{\mathcal{N}(#1 | #2, #3^{-1})} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle \mathsf{with} \quad w_n^i = \frac{\normd{\vc{m}_i}{\vc{m}_{n}}{\sigma_m^{2} I_M}}{\sum_{n' \in \cl{T}_u}\normd{\vc{m}_i}{\vc{m}_{n'}}{\sigma_m^{2} I_M}}.\label{eq:mriweights} \nonumber \end{align} \tag{ 8 }$$

Zoom In Zoom Out Reset image size

Here, $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\norm}[3]{\mathcal{N}(#1 | #2, #3^{-1})} \newcommand{\normd}[3]{\mathcal{N}(#1 | #2, #3)} \normd{\cdot}{\vcs{\psi}}{\vcs{\Sigma}}$ denotes a multivariate Gaussian with mean $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \vcs{\psi}$ and covariance matrix $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \vcs{\Sigma}$, $\newcommand{\vc}[1]{{\bf #1}} \vc{I}_M$ is the identity matrix of dimension M (the number of voxels in a patch) and $|\cdot|$ denotes set cardinality. $\sigma_y^2$ and $\sigma_m^2$ are the kernel variances, which, together with the artifact noise variance ${\sigma_t^*}^{2}$ (see section 2.1), will be automatically estimated given the observed data as described in section 2.3. We obtain the required uncorrupted part of the patient volume $\newcommand{\cl}[1]{\mathcal{#1}} \cl{T}_u$ (as well as its complement $\newcommand{\cl}[1]{\mathcal{#1}} \cl{T}_c$) by thresholding $\newcommand{\vc}[1]{{\bf #1}} f(\vc{x}_i)$ as shown in figure 1(c), leading to the following sets:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle \cl{T}_u \equiv \{i \in \cl{T}|f(\vc{x}_i) \leqslant 0.5 \} \quad \mathsf{and} \quad \cl{T}_c \equiv \cl{T} \setminus \cl{T}_u \label{eq:setdef} , \nonumber \end{align} \tag{ 9 }$$

which are illustrated in figure 2(c).

The kernel regression process is illustrated in figure 2(b) for the special case where 1x1x1 patches are used. After estimating $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(y_i, \vc{m}_i | \vcs{\sigma})$ using uncorrupted CT/MR samples $\newcommand{\vc}[1]{{\bf #1}} \{y_n, \vc{m}_n \}$, the prior distribution $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(y_i|\vc{m}_i, \vcs{\sigma})$ corresponds to drawing a trace at the observed $\newcommand{\vc}[1]{{\bf #1}} \vc{m}_i$ and normalizing, leading to a one-dimensional Gaussian mixture model that displays several peaks along the CT-axis. The mixture contains a large number of Gaussians, i.e. one for every sample in $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\cl}[1]{\mathcal{#1}} \{y_n,\vc{m}_n \}_{n \in \cl{T}_u}$. Since the samples in practice are clustered corresponding to tissue types, as illustrated in figure 2(b), the effective number of modes is however automatically reduced. In this way, kernel regression automatically performs an implicit tissue classification and reflects it in the prior model.

The expectation value of $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(y_i|\vc{m}_i, \vcs{\sigma})$ presents a way to calculate a (Bayesian) prior estimate of the uncorrupted CT:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle {\hat{y}^{\mathsf{pCT}}}_i \equiv \int_{-\infty}^{\infty}y_i p(y_i|\vc{m}_i, \vcs{\sigma}) dy_i= \sum_{n \in \cl{T}_u} w_n^i y_n, \quad \forall i \in \cl{T}. \label{eq:pct} \nonumber \end{align} \tag{ 10 }$$

This straighforward pseudo-CT (Andreasen et al 2015, Hofmann et al 2008) (pCT) estimate, which is dependent only on the MR scan, will serve as a baseline for our experiments.

Using the prior and likelihood that we have now defined, the posterior in equation (1) may be calculated. Defining the parameter set $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\e}{{\rm e}} \vcs{\sigma} \equiv \{\sigma_t^*, \sigma_y, \sigma_m \}$ and plugging in equations (7) and (4):

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\normd}[3]{\mathcal{N}(#1 | #2, #3)} \newcommand{\norm}[3]{\mathcal{N}(#1 | #2, #3^{-1})} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle p(y_i|t_i, \vc{m}_n, \vcs{\sigma}) = Z^{-1} \sum_{n \in \cl{T}_u} w_n^i \normd{y_i}{y_n}{\sigma_y^{2}} \normd{t_i}{y_i}{\sigma_t^{2}}. \label{eq:temp} \nonumber \end{align} \tag{ 11 }$$

Using that:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\normd}[3]{\mathcal{N}(#1 | #2, #3)} \newcommand{\norm}[3]{\mathcal{N}(#1 | #2, #3^{-1})} \displaystyle &\normd{t_i}{y_i}{\sigma_t^{2}} \normd{y_i}{y_n}{\sigma_y^{2}} \nonumber \\ &= \normd{y_i}{\mu_n^i}{(\sigma_t^{-2}+\sigma_y^{-2})^{-1}} \normd{t_i}{y_n}{\sigma_t^{2} + \sigma_y^{2}}\label{eq:temp2} \nonumber \end{align} \tag{ 12 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu_n^i \equiv \frac{1}{1+\frac{\sigma_t^{2}}{\sigma_y^{2}}} t_i + \frac{1}{1+\frac{\sigma_y^2}{\sigma_t^2}} y_n \label{eq:lincomb}, \nonumber \end{align} \tag{ 13 }$$

the normalizing factor (equation (3)) becomes:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\normd}[3]{\mathcal{N}(#1 | #2, #3)} \newcommand{\norm}[3]{\mathcal{N}(#1 | #2, #3^{-1})} \newcommand{\cl}[1]{\mathcal{#1}} \displaystyle Z = \sum_{n' \in \cl{T}_u} w_n^i \normd{t_i}{y_n}{\sigma_t^{2}+\sigma_y^{2}} \nonumber \end{align*}$$

and we get:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle p(y_i|t_i, \vc{m}_i, \vcs{\sigma}) = \sum_{n \in \cl{T}_u}v_n^i\mathcal{N}(y_i | \mu_n^i, (\sigma_y^{-2} + \sigma_t^{-2})^{-1}) \label{eq:posterior_calculated}, \nonumber \end{align} \tag{ 14 }$$

with weights:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \displaystyle v_n^i\equiv \frac{w_n^i \mathcal{N}(t_i | y_n, \sigma_t^{2} + \sigma_y^{2})}{\sum_{n' \in \cl{T}_u} w_{n'}^i \mathcal{N}(t_i | y_{n'}, \sigma_t^{2} + \sigma_y^{2})}\label{eq:weights}. \nonumber \end{align} \tag{ 15 }$$

Equation (14) is our main result. We will later use it to calculate a Bayesian estimate of {y _i} via its expectation:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle \hat{y}_i \equiv \int_{-\infty}^{\infty}y_i p(y_i|t_i,\vc{m}_i, \vcs{\sigma}) dy_i = \sum_{n \in \cl{T}_u} v_n^i \mu_n^i, \quad \forall i \in \cl{T}, \label{eq:kermar_estimate} \nonumber \end{align} \tag{ 16 }$$

which will be used in the proposed MAR methods.

2.2.1. Interpreting the posterior

The means and weights $\newcommand{\cl}[1]{\mathcal{#1}} \{\mu_n^i, v_n^i\}_{n \in \cl{T}_u}$ in the Gaussian mixture of equation (14) are similar to the ones in the prior model of equation (7), but are now determined by both the potentially corrupted FBP value t_i and the regression points $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\cl}[1]{\mathcal{#1}} \{y_n, \vc{m}_n \}_{n \in \cl{T}_u}$. How the two are blended is determined by the exact values of the parameters $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \vcs{\sigma}$, which are therefore critical for the shape of the posterior.

Kernel regression parameters

The influence of the kernel regression variances $\sigma_y^2$ and $\sigma_m^2$ is isolated to the prior $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(y_i|\vc{m}_i, \vcs{\sigma})$, and is as illustrated in figure 3(a).

Zoom In Zoom Out Reset image size

The $\sigma_y^2$ parameter controls the width of the Gaussians in the mixture of equation (7). This affects the modelled CT value variance within the implicit tissue classes defined by the kernel regression, which is reflected by the width of the peaks in the prior.

The $\sigma_m^2$ parameter scales the squared MR image patch differences in the Gaussian weights $w_n^i$, such that decreasing it translates to more weight on a smaller set of Gaussians with smaller $\newcommand{\vc}[1]{{\bf #1}} \|\vc{m}_i - \vc{m}_n\|$. Compared to larger settings (figure 3(a), right), a small value (figure 3(a), left) means that only patches in uncorrupted areas that have very similar intensity profiles contribute to the prediction, which shows as an amplification of the corresponding peaks in the prior.

Artifact noise variance

The influence of ${\sigma_t^{*}}^{2}$ is isolated to the likelihood $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(t_i|y_i, \vcs{\sigma})$ and governs the expected level of corruption near the metal voxels (near the metal, $\sigma_t^2 \approx {\sigma_t^*}^{2}$). For a given voxel, $\sigma_t^2$ directly determines the width of the likelihood and thus how the prior is modified to yield the posterior. Typical values of $\sigma_t^2$ lead to posteriors such as those in figure 3(b), where peaks near the observed FBP value t_i are amplified.

In the special case where $\sigma_t^2$ is very large, however, the likelihood becomes effectively a constant, so that the prior dominates the estimate and $\mu_n^i \approx y_n^i$, $v_n^i \approx w_n^i$, and the estimate in equation (16) therefore approaches the pCT in equation (10) ($\hat{y}_i \approx \hat{y}_i^{\mathsf{pCT}}$). The purely MR-based pCT thus arises as a special case of equation (16).

In the other extreme where $\sigma_t^2 \to 0$, the likelihood becomes a sharp function, while $\mu_n^i \to t_i$ and $\hat{y}_i \to t_i$. In that scenario the estimate therefore simply copies the FBP without reference to the MR scan. The FBP is thus another special case of equation (16).

2.3.1. Parameter tuning

As we have seen, the $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \vcs{\sigma}$ parameters critically alter the shape of the posterior as well as the predictions derived from it, and so it is important that they are chosen appropriately. The observed data themselves may be used as a guide for this process; for instance, a highly corrupted FBP would suggest a large ${\sigma_t^*}^{2}$ while a certain degree of CT value variation between voxels containing the same tissue type would suggest a certain setting for $\sigma_y^2$. We will now use a Bayesian method to automatically pick the parameters that best explain the observed data given our probabilistic model.

Equations (6) and (4) together define a joint distribution over the observed data $\newcommand{\vc}[1]{{\bf #1}} \{t_i, \vc{m}_i \}$ together with the unobserved variables {y _i}:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle &p(\{t_i, {\bf m}_i, y_i\}| \vcs{\sigma}) = \prod_{i \in \cl{T}}p(t_i, \vc{m}_i, y_i|\vcs{\sigma}) \nonumber \\ &\mathsf{with} \quad p(t_i, \vc{m}_i, y_i|\vcs{\sigma}) = p(y_i, \vc{m}_i|\vcs{\sigma})\,p(t_i|y_i,\vcs{\sigma}) \label{eq:generative}. \nonumber \end{align} \tag{ 17 }$$

Marginalizing equation (17) over the unobserved variables {y _i}, using equations (6) and (4) for the components, yields

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle &p(\{t_i, {\bf m}_i\}| \vcs{\sigma})= \prod_{i \in \cl{T}} p(t_i, {\bf m}_i| \vcs{\sigma}) \label{eq:marginal} \quad \mathsf{with}\nonumber \\ &p(t_i, \vc{m}_i|\vcs{\sigma}) = \int_{-\infty}^{\infty} p(t_i, \vc{m}_i, y_i|\vcs{\sigma}) d y_i\nonumber \\ &= |\cl{T}_u|^{-1} \sum_{n \in \cl{T}_u}\mathcal{N}({\bf m}_i | {\bf m}_n, \sigma_m^{2} {\bf I}_M)\nonumber \\ &\left[ \int_{-\infty}^{\infty}\mathcal{N}(t_i | y_i, f(\vc{x}_i) {\sigma_t^*}^{2})\mathcal{N}(y_i | y_n, \sigma_y^{2})dy_i \right]\nonumber \\ &= |\cl{T}_u|^{-1}\sum_{n \in \cl{T}_u} \mathcal{N}(t_i | y_n, f(\vc{x}_i) {\sigma_t^*}^{2} + \sigma_y^{2})\mathcal{N}({\bf m}_i | {\bf m}_n, \sigma_m^{2} {\bf I}_M) \nonumber . \nonumber \end{align} \tag{ 18 }$$

In the last step, we used equation (12) to carry out the integral. Given the observed variables $\newcommand{\vc}[1]{{\bf #1}} \{t_i, \vc{m}_i \}$, appropriate parameter values $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \vcs{\sigma}$ may now be obtained by maximizing equation (18), effectively fitting the model to the available data.

We perform the optimization using an iterative expectation-maximization (EM) algorithm (Dempster et al 1977) that only requires closed-form updates. Starting from an initial setting of $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \vcs{\sigma}$, EM alternately applies an expectation (E) and a maximization (M) step. In the E-step, a lower bound to the objective function $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \log p(\{t_i, {\bf m}_i\}| \vcs{\sigma})$ is constructed given the current estimate of $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \vcs{\sigma}$, derived directly from Jensen's inequality (Borman 2009):

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle \log p(\{t_i, {\bf m}_i\}| \vcs{\sigma}) & = \sum_{i \in \cl{T}} \log \sum_{n \in \cl{T}_u} \phi_n^i \nonumber \\ & = \sum_{i \in \cl{T}} \log \sum_{n \in \cl{T}_u} v_n^i \left(\frac{\phi_n^i}{v_n^i} \right) \nonumber \\ & \geqslant \underbrace{\sum_{i \in \cl{T}} \sum_{n \in \cl{T}_u} v_n^i \log \left(\frac{\phi_n^i}{v_n^i} \right)}_{(*)}, \nonumber \end{align*}$$

where $v_n^i$ are the weights in equation (15) and we defined:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle \phi_n^i = |\cl{T}_u|^{-1} \mathcal{N}(t_i | y_n, f(\vc{x}_i) {\sigma_t^*}^{2} + \sigma_y^{2})\mathcal{N}({\bf m}_i | {\bf m}_n, \sigma_m^{2} {\bf I}_M). \nonumber \end{align*}$$

The lower bound ($*$) is, by design, equal to $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \log p(\{t_i, {\bf m}_i\}| \vcs{\sigma})$ at the current parameter estimate, so when it is maximized during the M-step, the objective is guaranteed to increase (Dempster et al 1977). In general, the non-analytic $\newcommand{\vc}[1]{{\bf #1}} f(\vc{x}_i)$ prevents a closed-form formulation of this maximization, and so we simplify the objective by thresholding at f   =  0.5, as illustrated in figure 1⁵. This removes $\newcommand{\vc}[1]{{\bf #1}} f(\vc{x}_i)$ from the equations while splitting $\newcommand{\cl}[1]{\mathcal{#1}} \cl{T}$ into the corrupted and uncorrupted sets $\newcommand{\cl}[1]{\mathcal{#1}} \cl{T}_c$ and $\newcommand{\cl}[1]{\mathcal{#1}} \cl{T}_u$, leading to a closed-form M-step and algorithm 1:

Algorithm 1. Automatic parameter estimation.

1: Choose an initial estimate of the parameters

($\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \vcs{\sigma} \gets \{10^6, 10^6, 10^6 \}$), and set $\delta \leftarrow 1$.

2: while $\delta > 10^{-3}$ do

3:    $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \vcs{\sigma}_0 \leftarrow \vcs{\sigma}$

4:    E-step: Calculate ${v_n^i}$, $\newcommand{\cl}[1]{\mathcal{#1}} \forall i \in \cl{T}$ and $\newcommand{\cl}[1]{\mathcal{#1}} \forall n \in \cl{T}_u$,

using equation (15).

5:    M-step: Update the parameter estimates:

6:

$\newcommand{\cl}[1]{\mathcal{#1}} [{\sigma_y}]^{2} \leftarrow \frac{1}{|\cl{T}_u|} \sum_{i \in \cl{T}_u}\sum_{n \in \cl{T}_u} {v_n^i} (t_i-y_n){}^2$

$\newcommand{\cl}[1]{\mathcal{#1}} [{\sigma_t}^*]^{2} \leftarrow \frac{1}{|\cl{T}_c|} \sum_{i \in \cl{T}_c}\sum_{n \in \cl{T}_u} {v_n^i} (t_i - y_n){}^2 - [\sigma_y]^{2}$

$\newcommand{\vc}[1]{{\bf #1}} \newcommand{\cl}[1]{\mathcal{#1}} [\sigma_m]^{2} \leftarrow \frac{1}{|\cl{T}|} \sum_{i \in \cl{T}}\sum_{n \in \cl{T}_u} {v_n^i} \frac{\| (\vc{m}_i-\vc{m}_n)\|^2}{M}$.

7:  $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \delta \leftarrow \| \vcs{\sigma} - \vcs{\sigma}_0\| / 3$

8: end while

The calculation of $v_n^i$ in the E-step may be viewed as a probabilistic assignment of the data points to the tissue classes that were implicitly defined during kernel regression. Given this classification, the update equations in the M-step estimate the within-class variance, over different parts of the patient volume:

• $\sigma_y^{2}$ is calculated over the uncorrupted set $\newcommand{\cl}[1]{\mathcal{#1}} \cl{T}_u$, which reflects the observed noise in the CT values within a tissue class in the absence of the metal artifacts.
• ${\sigma_t^*}^{2}$ updates to the additional variance in the metal artifact corrupted volume $\newcommand{\cl}[1]{\mathcal{#1}} \cl{T}_c$, reflecting the level of artifact corruption.
• $\sigma_m^{2}$ updates to the MR image patch variance over the entire volume.

During the iterations, the update of the parameters in the M-step improves the soft classification in the E-step as the model increasingly fits the data, which in turn improves the parameter estimates. This continues until the objective is maximized and the parameters stop changing.

We propose to use the posterior predictive distribution $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(y_i|t_i, \vc{m}_i, \vcs{\sigma})$, specified in equation (14) and with parameters $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \vcs{\sigma}$ automatically tuned using algorithm 1, as the basis for three different MAR algorithms. Specifically, we define an image inpainting method that simply uses the mean of the predictive distribution; a sinogram inpainting method that uses the image inpainting result as a template to estimate the metal projections; and an MBIR algorithm that uses the full distribution $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(y_i|t_i, \vc{m}_i, \vcs{\sigma})$ as a reconstruction prior together with a Poisson likelihood model of the x-ray intensity measurements.

2.4.1. Image inpainting

2.4.2. Sinogram inpainting

2.4.3. MBIR

MBIR maximizes the posterior distribution of a CT reconstruction given the acquired x-ray intensity data and a reconstruction prior, for which we propose to use our predictive distribution:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle p(\{y_i\}|\{n_j\}) \propto p(\{n_j\}|\{y_i\}) \prod_{i \in \mathcal{T}} p(y_i|t_i, \vc{m}_i, \vcs{\sigma}) \label{eq:reconstruction_posterior}, \nonumber \end{align} \tag{ 19 }$$

where $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(y_i|t_i, \vc{m}_i, \vcs{\sigma})$ is given by equation (14). Here {n_j} is the set of x-ray intensity measurements for difference paths through the patient volume. Following (Nuyts et al 1998), we use a Poisson likelihood:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \displaystyle p(\{n_j\}|\{y_i\}) = & \prod_{j} p(n_j|\{y_i\}), \nonumber \\ p(n_j|\{y_i\}) = \mathsf{Poisson}\{\lambda_j\} & , \quad \lambda_j = \Gamma_j \mathsf{e}^{-\sum_{i \in \cl{T}} l_{j, i} y_i} \nonumber , \nonumber \end{align*}$$

where $\Gamma_j$ is the emitted x-ray count toward detector j , and $\newcommand{\vc}[1]{{\bf #1}} \vc{L}$ is the system matrix whose entry l_j ,i defines the intersection between the x-ray path j  and voxel i. To maximize the reconstruction posterior in equation (19) we use the iterative MLTR algorithm (Nuyts et al 1998, Van Slambrouck and Nuyts 2012) in its extended form De Man et al (2000), De Man (2001), which allows for a general reconstruction prior. Starting from an initial reconstruction estimate, the algorithm maximizes the log-posterior by iteratively applying an additive term. In addition to the x-ray data likelihood, this term depends on the first and second derivatives of the log-prior (De Man et al 2000, De Man 2001), which in our case become

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle \ln(\,p(y_i|t_i,\vc{m}_i, \vcs{\sigma}))^{\prime} = (\sigma_t^{-2} + \sigma_y^{-2})[\sum_{n \in \cl{T}_u} \tilde{v}_n^i \mu_n^i - y_i] \quad \mathsf{and} \quad \label{eq:mltr_first} \nonumber \end{align} \tag{ 20 }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\cl}[1]{\mathcal{#1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \newcommand{\vc}[1]{{\bf #1}} \displaystyle &\ln(\,p(y_i|t_i,\vc{m}_i, \vcs{\sigma}))^{\prime \prime} \nonumber \\ &=(\sigma_t^{-2}+\sigma_y^{-2})^2\bigg(\sum_{n \in \cl{T}_u}\tilde{v}_n^i {\mu_n^{i}}^2 - [\sum_{n \in \cl{T}_u} \tilde{v}_n^i \mu_n^i]^2\bigg) - (\sigma_t^{-2}+\sigma_y^{-2}) \label{eq:mltr_second} \nonumber \end{align} \tag{ 21 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\normd}[3]{\mathcal{N}(#1 | #2, #3)} \newcommand{\norm}[3]{\mathcal{N}(#1 | #2, #3^{-1})} \newcommand{\cl}[1]{\mathcal{#1}} \displaystyle \mathsf{where} \quad \tilde{v}_n^i = \frac{v_n^i \normd{y_i}{\mu_n^i}{(\sigma_t^{-2}+\sigma_y^{-2})^{-1}}}{\sum_{n' \in \cl{T}_u} v_{n'}^i \normd{y_i}{\mu_{n'}^i}{(\sigma_t^{-2}+\sigma_y^{-2})^{-1}} }\label{eq:mltr_weights}. \nonumber \end{align} \tag{ 22 }$$

We will refer to the resulting reconstruction algorithm as MLTR-k.

We evaluated the three proposed algorithms on an anonymized retrospective data set of nine head-and-neck radiotherapy patients containing dental implants and/or fillings. The image sets had a resolution of $1.2 \times 1.2 \times 2.0$ mm (CT) and $0.5 \times 0.5 \times 5.5$ mm (MR). The MR images were resampled to the CT resolution after rigid, multi-modal co-registration by mutual information (Wells et al 1996, Maes et al 1997) using the MatLab (version R2019a) image processing toolbox (Mathworks 2019). The CT scanner was a Philips Brilliance Big Bore with kVp120, the MR scanner a Philips Panorama 1.0T HFO. The patients were MR scanned with a T1-weighting 2D spin-echo sequence at TE  =10 ms and TR  =520.2–572.2 ms (varying between scans). The patients were positioned in the same fixation for both MR and CT scans. For MBIR and sinogram inpainting, we exported the sinograms from the CT scanner.

A generic procedure for implementing the three MAR algorithms is described in algorithm 2:

Algorithm 2. Summary of implementation.

1: Calculate the FBP.

2: Segment the metal implants.

3: Calculate $\newcommand{\vc}[1]{{\bf #1}} f(\vc{x}_i)$ according to equation (5).

4: Obtain the uncorrupted set $\newcommand{\cl}[1]{\mathcal{#1}} \cl{T}_u$ by thresholding

$\newcommand{\vc}[1]{{\bf #1}} f(\vc{x}_i)$ at 0.5 (see equation (9)).

5: Estimate $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} \vcs{\sigma}$ using algorithm 1, then calculate the

weights $\newcommand{\cl}[1]{\mathcal{#1}} \{w_n^i\}_{n \in \cl{T}_u}$ using equation (8) and $\newcommand{\cl}[1]{\mathcal{#1}} \{\mu_n^i, v_n^i \}_{n \in \cl{T}_u}$

using equations (13) and (15), $\newcommand{\cl}[1]{\mathcal{#1}} \forall i \in \cl{T}$.

6: Apply the three MAR algorithms:

1. For image inpainting, use equation (16) and the metal segmentation of step 2.

2. For sinogram inpainting, use the image inpainting result as a template and proceed as outlined in section 2.4.

3. For MBIR, use the derivatives in equations (20) and (21) along with equation (22) in an MLTR implementation.

For our experiments, we implemented the various steps as follows:

• Step 1:  
  We acquired the FBPs as reconstructed by the vendor-provided scanner software.
• Step 2:  
  We performed the metal segmentation automatically using Otsu's method (Otsu 1979).
• Step 3:  
  We found the exact value of $\kappa$ in the expression for $\newcommand{\vc}[1]{{\bf #1}} f(\vc{x}_i)$ to be non-critical, and used $\kappa = (10~{\rm mm}){}^2$ for all our head-and-neck patients.
• Step 5:  
  The vast majority of the weights $\newcommand{\cl}[1]{\mathcal{#1}} \{w_n^i\}_{n \in \cl{T}_u}$ will in practice attain very small values and therefore not contribute to the model in a meaningful way. In order to speed up computations, we therefore used a fast patch matching algorithm (Ta et al 2014) to identify 200 regression points for each voxel i with particularly small patch differences $\newcommand{\vc}[1]{{\bf #1}} \| \vc{m}_n - \vc{m}_i \|$ and therefore large weights $w_n^i$, effectively clamping the weights of all other regression points to zero. We used $5 \times 5 \times 5$ patches, i.e. $6 \times 6 \times 10$ mm.
• Step 6:  
  For sinogram inpainting, we used 3D spiral forward projection to detect the metal projections and calculate the prior projections. For the interpolation of the metal projections, we used 2D barycentric linear interpolation over the cylindrical detector array on a triangular grid. For image reconstruction, since the MBIR implementation in particularrequired numerous slow but highly parallellizable forward and back projection operations, speed was a priority both for potential clinical implementation and to facilitate the experiments. We therefore used the GPU-accelerated primitives in the ASTRA (van Aarle et al 2015, 2016) package, as well as the bundled FBP implementation (with a Shepp–Logan filter). We performed the reconstructions in 2D, after rebinning and interpolating the 3D spiral sinograms from the scanner to sets of 2D sinograms with a linear detector geometry.We initialized the iterative reconstruction process in MBIR with uniform images with an attenuation coefficient of 10⁻⁶ mm⁻¹. We stopped iterating once the voxel-averaged change between iterations fell below 10⁻⁶ mm⁻¹.

Our entire computational framework used Python. The MAR algorithms, parameter estimation algorithm and sinogram rebinning methods in particular used NumPy and SciPy, while the reconstruction algorithms used ASTRA in its Python-wrapped form.

For each MAR algorithm, we evaluated the quantitative accuracy of the artifact-reduced CT images as follows: first, we acquired manual delineations of the oral cavity and the teeth for each patient, drawn in the FBP CT images guided by the MR scans. We split these delineations into a corrupted and uncorrupted part using our definitions of $\newcommand{\cl}[1]{\mathcal{#1}} \cl{T}_c$ and $\newcommand{\cl}[1]{\mathcal{#1}} \cl{T}_u$, and calculated reference mean CT values for each structure and each patient using the CT values in the uncorrupted parts. Around these mean values, we calculated the CT value standard deviations (STDs) over the corrupted parts; this provided an image quality metric for each structure, with each MAR, and for each patient. We performed the calculations in Hounsfield Units (HU) (Buzug 2008). We finally contrasted the STD observations for the different MAR algorithms using a repeat measurements Student's t-test, and report the patient-averaged STD results for all MAR algorithms and structures, as well as the p-values of the tests.

The use of the STD metric as an image quality metric relies on the assumption that different partitions of each ROI do not display large differences in mean HU values in the absence of artifacts as a result of large anatomical variations; such differences would lead to a systematic error in the reference mean, directly affecting the STD results. To estimate the magnitude of this error, we used the uncorrupted parts of our nine-patient cohort to emulate the artifact-free tissue distributions in the oral cavity and teeth. Using that the corrupted and uncorrupted sets in our experiments were, on average, not far from equal-sized, we then created 10⁶ random bipartions of the extracted set of CT values, for each calculating the difference in mean between the partitions. These differences between means turned out to be near-Gaussian distributed; for the oral cavity the average was  ∼0 HU and the standard deviation 1.1 HU, while for the teeth the average was  ∼0 HU and standard deviation 6.3 HU. Even after multiplying the standard deviations by a factor 10 for a more conservative estimate, the error on the reference means in our experiments should thus be in the tens of HU, which is an order of magnitude smaller than the variations between MAR algorithms that we present in the following section. Consequently, the STD metric appears to be a valid indicator of artifact suppression.

Figure 4 shows visual results of the MR-based image inpainting algorithm kerMAR for representative head-and-neck patients, at a narrow, soft tissue enhanced window level (a) and a wider level (b). The T1w MR images are shown in (c). The green arrows in (a) indicate regions where kerMAR provided notable artifact reduction by suppressing both high and low intensity streaks, while at the same time preserving complex structures where the CT values are difficult to predict from the MR scan, such as the teeth.

Zoom In Zoom Out Reset image size

The region where kerMAR was least successful was the windpipe, indicated by the rings. Here, the sometimes severe misalignments of the CT and MR images, due to inter-scan patient motion, led to a highly inaccurate MR-based prior and thus anatomical errors. However, in similar but less severe cases, such as with smaller misalignments, mechanisms in the kerMAR algorithm allowed it to avoid such anatomical errors. In particular, as we saw in section 2.2.1, the kerMAR image is an intermediate between the extreme special cases of pCT (${\sigma_t^*}^{2} \to \infty$) and FBP (${\sigma_t^*}^{2} \to 0$), and thus corresponds to a non-trivial blending of the FBP and MR-based prediction. Figure 5 shows the kerMAR image along with these two special cases for a patient where the MR and CT scans were misaligned in the teeth. In this less severe case of misalignment, while we see poor pCT performance and thus a poor prior-based prediction, the blending with the FBP led to a much improved kerMAR. Similar results occurred in the spinal cord (red square), where similar inaccuracies led to an apparent introduction of MR features in the pCT, which were successfully suppressed in the kerMAR.

Zoom In Zoom Out Reset image size

To summarize the results of the visual inspection, we saw (1) clear visual artifact reduction with both kerMAR and pCT and (2) a clear performance difference in bone and air regions between the two MAR algorithms. For the quantitative analysis of the oral cavity and teeth, these observations are reflected in the STDs shown in figure 6, which measure the CT value standard deviation in the corrupted parts of the structures around the mean intensity in the uncorrupted parts. In the oral cavity, kerMAR and pCT were equally effective, with significant STD reductions relative to the FBP of  ∼150 HU. In the teeth, however, the pCT did not improve on the FBP, consistent with our observations of poor pCT performance in such thick, bony regions. Here, only kerMAR showed a benefit, providing a significant improvement over the pCT of  ∼100 HU.

Zoom In Zoom Out Reset image size

4.1.1. Influence of automatic parameter tuning

We finally considered the influence of the automatic patient-specific parameter tuning (algorithm 1) by intentionally swapping the estimated parameters between patients. Figure 7 shows the results for two representative patient-pairs, where we swapped the parameters for patient 1 and 3 with those of patient 2 and 4, respectively. Patients 1 and 3 displayed relatively smaller settings of $\sigma_y^2$ and ${\sigma_t^*}^{2}$, and swapping the parameters with those of patient 2 and 4 therefore led to increased $\sigma_y^2$ and ${\sigma_t^*}^2$. This substantially altered the predictive model; as illustrated in figure 3, a ${\sigma_y}^2$ increase widens the prior peaks while a ${\sigma_t^*}^{2}$ increase widens the likelihood, leading to a prediction that was both less precise and less accurate in areas where the MR-based prediction was compromised. These changes in the predictive model caused errors especially in the teeth (rings and arrow), as well as blur for patient 1, as compared to the automatically tuned parameters.

Zoom In Zoom Out Reset image size

Patients 2 and 4 correspondingly displayed larger settings of $\sigma_y^2$ and ${\sigma_t^*}^{2}$. Upon swapping, the ensuing decrease in $\sigma_y^2$ and ${\sigma_t^*}^2$, and thus narrowing of the posterior, led to image noise and less effective artifact reduction; this is clearly visible for patient 2. Along with the larger ${\sigma_t^*}^{2}$ and ${\sigma_y}^2$, patients 2 and 4 further displayed smaller values of $\sigma_m^2$. The swap therefore also led to an increase in $\sigma_m^2$ and thus a prior model with a less sensitive MR patch-based distinction between its regression points. Compared to using the automatically tuned parameters, this led to errors for patient 2 in the thick molars indicated by the circles, where the relevant MR patch differences are subtle.

Figure 8 shows visual results of the proposed MR-based sinogram inpainting method nMAR-k, alongside those produced with a standard nMAR implementation based on a conventional, FBP-based template image. As seen at the narrow, soft tissue enhancing window level in figure 8(a), nMAR and nMAR-k performed similarly in terms of artifact reduction; neither algorithm reconstructed an artifact-free oral cavity, while both reduced artifacts elsewhere to similar degrees (see regions indicated by arrows).

Zoom In Zoom Out Reset image size

At the wider window level in figure 8(b), we see a clearer difference between nMAR and nMAR-k: The sometimes lower quality of the FBP-based template in standard nMAR led to the visible artifacts (arrows), which are absent using the proposed MR-based method. A clear benefit of using the higher quality MR-based template thus appears to be to avoid certain artifacts introduced by the CT-based one. The influence of such improvements on the quantitative STD metric was comparatively minor, as figure 9 shows significant improvement of both methods over the FBP, but little quantitative difference between them.

Zoom In Zoom Out Reset image size

In order to evaluate the benefit of including the predictive distribution $\newcommand{\vc}[1]{{\bf #1}} \newcommand{\vcs}[1]{\boldsymbol{#1}} p(\{y_i\}|\{t_i, \vc{m}_i \}, \vcs{\sigma})$ as an image reconstruction prior in the MLTR algorithm, we ran the algorithm both with (the proposed MLTR-k algorithm) and without this prior (referred to in the remainder as simply 'MLTR'). We ran each algorithm until the voxel-averaged difference accross two subsequent iterations of the reconstructed image decreased below a predefined threshold (10⁻⁶). Figure 10 displays the obtained reconstructions at the specified convergence threshold for both methods in five representative patients. Results are shown at a soft tissue enhancing window level (a) and at a wider level (b). At the narrow window level, the main detectable benefit of MLTR-k compared to MLTR was less blurry reconstructions, which indicate that the algorithm was closer to full convergence; in terms of artifact reduction, MLTR-k more successfully eliminated low intensity streaks but did not wholly eliminate high intensity ones. At the wider window level, MLTR-k more clearly improved the MAR over MLTR (arrows in (b)), although these benefits came, in a few cases, at the expense of some newly introduced artifacts in the teeth (see black rings), that however were not a general trend.

Zoom In Zoom Out Reset image size

Figure 11(a) further shows the logarithm of the voxel-averaged difference used as convergence criterion for both algorithms and our nine patients across iterations. As we can see, the rate of convergence is universally larger for MLTR-k with the MR-based prior, as compared to the purely likelihood-based MLTR. In particular, the convergence speed was essentially doubled with MLTR-k, decreasing the number of iterations from  ∼400–600 to  ∼200–300.

Zoom In Zoom Out Reset image size

Visually, MLTR-k with its MR-based prior thus provided superior MAR to the prior-free MLTR, as well as an image closer to full convergence at the specified convergence threshold. These results are supported by the quantitative results in figure 11(b) that show significant standard deviation improvements compared to MLTR, in both the oral cavity and teeth, of respectively  ∼200 HU and  ∼150 HU.

The nMAR sinogram inpainting algorithm (with the CT-based template) may be used as a gold standard for benchmarking and comparing the three proposed MR-based algorithms between themselves. To more easily perform this comparison, we summarize the results of our quantitative evaluations in figure 12 along with p-values of the comparison to nMAR; the results for the MR-based results are highlighted. We see no case of significantly worse performance compared to the gold standard of the three MR-based algorithms, and further see a quantitative improvement in the oral cavity using the image space approach (kerMAR). Additionally, the quantitative results should be taken together with the visual comparison, and when we compared the sinogram inpainting method (nMAR-k) to nMAR in section 4.2, we found the negligible quantitative results to be accompanied by visual improvements. The comparison is more difficult with the MBIR approach (MLTR-k), where we again see insignificant STD results but where, as we saw in section 4.3, the visual quality may have been compromised by additional artifacts.

Zoom In Zoom Out Reset image size

In summary, the image and sinogram inpainting algorithms at least show benefits over the gold standard nMAR algorithm, while the MBIR results are more inconclusive.

While our MR-based method is applicable whenever an MR scan is co-acquired with a CT, it should be noted that such cases may be of limited frequency outside of head-and-neck radiotherapy, as MR at the time of this study is only used infrequently in the case of metal artifacts at other sites, such as hip implants. Due to the continually increasing importance of MR in radiotherapy, this may however change in the future (Lagendijk et al 2014).

The inclusion of the FBP for the CT value prediction in our image inpainting algorithm kerMAR helped to mitigate co-registration and alignment issues between the MR and CT images, but for clinical adoption the associated errors in especially the windpipe should probably further reduced (see figures 4(a) and (b)). A straightforward way to achieve this is to improve the co-registration, e.g. by using a deformable registration methods, which are routinely used in modern radiotherapy clinics. Further accuracy gains may be achieved by improving the noise model, which in its current form simply assumes that the noise introduced by the artifacts decreases sigmoidally with the distance to the metal implants. This is a rough assumption as especially streak artifacts may persist throughout the entire CT image (see, for instance, figure 4(a)). Our current model may also not translate well to other types of implants, such as metallic dual hip prostheses, where the severity of the artifacts also depends heavily on whether or not they are located between the prostheses. Future experiments may therefore investigate different formulations of the function $\newcommand{\vc}[1]{{\bf #1}} f(\vc{x}_i)$ that better capture the spatial dependency of the artifact noise.

In order to successfully translate the proposed techniques into clinical routine, the required computation time will need to be further reduced. The most time-consuming part of our model is the calculation of patch distances over regression points, which is currently implemented on a subset of 200 points $\newcommand{\cl}[1]{\mathcal{#1}} \forall i \in \cl{T}$, found using a patch matching algorithm (Ta et al 2014). In our current Python implementation on a single CPU core (Intel Core i7-4712HQ $@$ 2.30 GHz), this process takes between 10–30 min. The algorithm is however parallelizable, and on a similar-sized dataset, Ta et al report results on the order of  ∼1 min on a multi-CPU cluster (Ta et al 2014). In the future, we therefore intend to speed up our algorithm in a similar manner.

Our MBIR experiments used a relatively simple Poisson likelihood. While this model helps address the artifacts stemming from photon starvation of the metal projections (Nuyts et al 1998, Buzug 2008), it does not account for e.g. the important effect of beam hardening. This may have been the source of the artifacts that were sometimes introduced with MLTR-k (see e.g. the black rings in figure 10) (Nuyts et al 1998, De Man et al 2000, 2001, Buzug 2008, Van Slambrouck and Nuyts 2012). Future work may therefore consider using a more accurate likelihood model for the MBIR technique.