Can Planning Images Reduce Scatter in Follow-Up Cone-Beam CT?

TL;DR

This paper evaluates different scatter estimation methods for cone-beam CT, demonstrating that incorporating accurate scatter estimates into statistical reconstruction improves image quality, with Monte-Carlo based estimates offering marginal benefits over simpler measurement-based methods.

Contribution

It introduces a method to match offline Monte-Carlo scatter estimates with new measurements and compares their effectiveness in statistical reconstruction for CBCT.

Findings

Statistical reconstruction with scatter estimates outperforms analytic and iterative pre-correction methods.

Monte-Carlo scatter estimates provide only slight improvements over measurement-based methods.

Incorporating scatter estimates significantly enhances image accuracy in CBCT reconstruction.

Abstract

Due to its wide field of view, cone-beam computed tomography (CBCT) is plagued by large amounts of scatter, where attenuated photons hit the detector, and corrupt the linear models used for reconstruction. Given that one can generate a good estimate of scatter however, then image accuracy can be retained. In the context of adaptive radiotherapy, one usually has a low-scatter planning CT image of the same patient at an earlier time. Correcting for scatter in the subsequent CBCT scan can either be self consistent with the new measurements or exploit the prior image, and there are several recent methods that report high accuracy with the latter. In this study, we will look at the accuracy of various scatter estimation methods, how they can be effectively incorporated into a statistical reconstruction algorithm, along with introducing a method for matching off-line Monte-Carlo (MC) prior estimates to the new measurements. Conclusions we draw from testing on a neck cancer patient are: statistical reconstruction that incorporates the scatter estimate significantly outperforms analytic and iterative methods with pre-correction; and although the most accurate scatter estimates can be made from the MC on planning image, they only offer a slight advantage over the measurement based scatter kernel superposition (SKS) in reconstruction error.

Figures

Tables

Table 1. Table 1 : Experimental results from various estimation and reconstructions

Experiment	FDK	NLL
oracle
none
SKS
offline-prior-MC

Equations

$$y_k\sim \mathrm{Poisson}\left\{b_k\exp (-{\left[\mathbf{\Phi }\mathbf{\mu }\right]}_i)+s_k\right\}\text{mbox}fork=1,\dots ,N,$$

$$\mathbf{s}\{i\}=\frac{\mathrm{SPR}}{N_{\mathrm{air}}(i)}\sum _{k\in {\mathcal{C}}_{\mathrm{air}}(i)}y(k)\text{mbox}fori=1,\dots ,P,$$

$$\mathbf{s}=\mathcal{F}(\mathbf{y}-\mathbf{b}\odot \exp (-\mathbf{\Phi }{\mathbf{\mu }}_{\mathrm{reg}})),$$

$$\tilde{\mathbf{s}}\sim \mathcal{MC}(\mathbf{b},\mathbf{\Phi },{\mathbf{\mu }}_{est.},N_{\mathcal{MC}}),$$

$$\mathbf{p}=\mathbf{\Phi }\mathbf{\mu }+\mathbf{n},$$

$$p_i=\log \left(\frac{b_i}{y_i-s_i}\right)\text{mbox}fori=1,...,N,$$

$$\hat{\mathbf{\mu }}=\underset{\mathbf{\mu }}{\mathrm{argmin}}{\left(\mathbf{\Phi }\mathbf{\mu }-\mathbf{p}\right)}^T\mathbf{W}\left(\mathbf{\Phi }\mathbf{\mu }-\mathbf{p}\right)+\lambda R(\mathbf{\mu }),$$

$$s_M\{i\}=\mathcal{I}\left(s\{i\};(u_M,v_M{)}_i\right)\text{mbox}fori=1,\dots ,P,$$

$$\left[\begin{array}[]{c}u^{\prime}\\ v^{\prime}\\ 1\end{array}\right]=\begin{bmatrix}w_{1,1}&w_{1,2}&w_{1,3}\\ w_{2,1}&w_{2,2}&w_{2,3}\\ 0&0&1\end{bmatrix}\left[\begin{array}[]{c}u\\ v\\ 1\end{array}\right],$$

$$\left[\begin{array}[]{c}r_{x}\\ r_{y}\\ r_{z}\end{array}\right]=\begin{bmatrix}\cos\varphi_{i}&-\sin\varphi_{i}&0\\ \sin\varphi_{i}&\cos\varphi_{i}&0\\ 0&0&1\end{bmatrix}\left[\begin{array}[]{c}t_{x}\\ t_{y}\\ t_{z}\end{array}\right].$$

$$\left[\begin{array}[]{c}u_{M}\\ v_{M}\\ 1\end{array}\right]=\begin{bmatrix}\frac{l_{\mathrm{SO}}}{l_{\mathrm{SO}}+r_{x}}&0&\frac{l_{\mathrm{SD}}}{l_{\mathrm{SO}}+r_{x}}r_{y}\\ 0&\frac{l_{\mathrm{SO}}}{l_{\mathrm{SO}}+r_{x}}&\frac{l_{\mathrm{SD}}}{l_{\mathrm{SO}}+r_{x}}r_{z}\\ 0&0&1\end{bmatrix}\left[\begin{array}[]{c}u\\ v\\ 1\end{array}\right],$$

$$\mathrm{NLL}(\mathbf{\mu };\mathbf{y},\mathbf{s})=\sum _{i=1}^N{b}_i\exp (-{\left[\mathbf{\Phi }\mathbf{\mu }\right]}_i)+s_i-y_i\log \left(b_i\exp (-{\left[\mathbf{\Phi }\mathbf{\mu }\right]}_i)+s_i\right).$$

$$\hat{\mathbf{\mu }}=\underset{\mathbf{\mu }\in \mathcal{C}}{\mathrm{argmin}}\mathrm{NLL}(\mathbf{\mu };\mathbf{y},\mathbf{s})+\lambda R(\mathbf{\mu }),$$

Full text

1 Introduction

Cone-beam computed tomography (CBCT) is an imaging modality that is seeing increased use for image guidance procedures, such as radiation therapy [1]. A key challenge of this geometry is the vast quantities of scattered photons that reach the detector [2], and contaminate other line-of-sight measurements. Usually in this context however, one has a planning scan of the same patient from a more accurate CT acquisition, such as a helical fan-beam, which has significantly lower scatter due to better collimation and narrower field-of-view. Typical approaches to CBCT scatter correction either form a self-consistent model based solely on the new measurements [3][4][5][6], or exploit the prior image as a basis for its estimation [7][8][9].

In this study, we will adopt the perspective that various scatter models allow one to make an estimate of its expectation, given a set of new CBCT projections, and potential access to a low-scatter prior CT image. We will then investigate two key aspects: how accurate can scatter be estimated in a moderate and low dose settings, and what is the robustness of various methods in this challenging scenario; and how should these estimates be incorporated into reconstruction.

In the first case, the scatter estimation methods we will look at fall into two distinct classes—methods that are blind to the planning image, and methods that exploit it. In our ‘prior blind’ category are a view dependent uniform estimation [5], the scatter kernel superposition (SKS) [4] and fast asymmetric SKS (fASKS) [6], along with simulating the scatter through a Monte-Carlo (MC) engine on a preliminary reconstruction with the fast Feldkamp–Davis–Kress (FDK) [10] algorithm. Conversely, based upon a rigid registration of the planning image, we will look at the effectiveness of taking scatter as a smooth projection difference [7][11], and from using the MC engine on this registered plan. Here, we will look at both calculating the planning MC estimate on-line, after registration [9], along with the notion of matching an off-line pre-calculated estimate to the measurements.

In the second case, in the subsequent reconstruction with each of the estimation methods, we make a distinction between ‘scatter correction’, where measurements a pre-processed to remove its effect, and ‘scatter-aware inference’, where the imaging operates based on the raw uncorrected measurements and knowledge of the scatter estimate. Most popular ‘analytic’ and ‘iterative’ techniques, such as FDK and PWLS [12], fall into the former category. Performing the inference, although more challenging due to its non-linearity and nonconvexity, represents a more accurate data model that may mitigate reconstruction artefacts and errors.

We begin this article with relevant background material, where we explain the system model in Section 2.1, an overview of scatter estimation methods in Section 2.2, and standard reconstruction based on pre-corrected measurements in Section 2.3. Next, we give details of matching an off-line planning MC scatter given a rigid translation of the specimen in Section 3.1, along with the model for statistical inference reconstruction in Section 3.2.

From a dataset derived from repeat CT images of a neck caner patient, we then evaluate both the scatter estimation accuracy and reconstruction error with the range the methods under test. The results are then presented in Section 4.3.

2 Background

2.1 System Description

The system we will directly study in this work is a circular scan CBCT. This consists of a point source and flat panel detector, which rotate throughout $360^{\circ}$ around the specimen where photons are emitted and measured after a fixed angular increment. Assuming that recorded x-rays are drawn from independent Poisson distributions given as [13][14], then we can write the distribution using monoenergetic Beer-Lambert law with additive scatter component as

[TABLE]

where $\hat{\boldsymbol{y}}\in\mathbb{R}^{N}$ is a column vector of measurements, $\boldsymbol{b}\in\mathbb{R}^{N}$ is a vector of input source fluxes, $\boldsymbol{\mu}\in\mathbb{R}^{M}$ is a vector of linear attenuation coefficients, and $\boldsymbol{\Phi}\in\mathbb{R}^{N\times M}$ is a system matrix describing the path of each ray through the specimen and onto the detector.

2.2 Estimating Scatter

In this section, we give brief overviews of the various scatter estimation techniques that we evaluate in this study.

2.2.1 Uniform:

A simple method that calculates a constant scatter at each projection angle, we denote as ‘uniform’ [5][15]. Here, using the assumption that a scatter-to-primary ratio (SPR) is known a-priori, along with a distinction between air and object containing projections, we can write

[TABLE]

where $\mathcal{C}_{\mathrm{air}}(i)$ is the set of air containing measurements at the $i^{\mathrm{th}}$ angle, where the set satisfies $y(k)\geq t_{\mathrm{air}}\,\forall k\in\mathcal{C}_{\mathrm{air}}$ with some scalar threshold $t_{\mathrm{air}}$. This essentially calculates the mean scatter given a constant ratio. In practice, this SPR can be found by observing the magnitude of signal in the air region with and without a specimen present, and assume the difference is scatter. To ensure the scatter is less than the minimum value in a given projection, a non-negativity constraint can be added [15].

2.2.2 SKS/ASKS:

The scatter kernel superposition (SKS) [4] and asymmetric SKS (ASKS) [6], perform estimation as a convolution of the scatter free incident beam with an appropriate kernel. Since the incident beam is itself unknown, the methods iteratively estimate this as the difference between raw measurements and updated scatter estimate from the previous iteration. Due to the ability calculate convolution rapidly through the FFT, both of these methods are relatively fast, especially when the projections are sub-sampled. Although they implicitly model the scatter media as homogeneous, the estimates are accurate in practice.

2.2.3 Diff. filt.:

A simple concept for predicting the scatter contribution based upon a planning image is though the smooth difference—‘diff. filt.’—between CBCT measurements and projections of a registered plan [7]. This model can be expressed as

[TABLE]

where $\boldsymbol{\mathcal{F}}$ is a projection-wise filter, or set of filters—subsequent median and Gaussian filters are used in [7]—and $\boldsymbol{\mu}_{\mathrm{reg}}$ is the registered planning CT onto a preliminary reconstruction of the CBCT measurements.

2.2.4 Monte-Carlo (MC):

Monte-Carlo scatter estimation techniques essentially draw a number of samples from an accurate probabilistic model of physical interactions. Given that the model is a faithful representation of reality, then the true expectation of scatter can be found with infinite samples.

We denote this process as

[TABLE]

where $\tilde{\boldsymbol{s}}$ is the estimate after $N_{\mathcal{MC}}$ photon simulations distributed throughout several projection angles and $\boldsymbol{\mu}_{\mathrm{est.}}$ is the image onto which the estimation is based. We will test the ability of estimating MC scatter onto both a preliminary FDK with appropriate prior blind estimation, such as SKS, and onto the planning image.

If the MC simulation is made after the measurements are taken, then it may be appropriate to sub-sample both the image $\tilde{\boldsymbol{s}}$ and the number of photons, in order to complete the calculation on-line [9]. In the off-line setting, there is no immediate limitation on computational time, as it can be performed days or weeks ahead of the follow-up CBCT. Eventual matching of this off-line estimate is detailed in Section 3.1.

2.3 Reconstruction from Scatter Correction

In most cases, reconstruction is performed by inferring the attenuation coefficient given the model in (1), which follows an effective ‘correction’ of scatter. In the crudest form, this involves simply subtracting the scatter estimate from the measurements. An advantage of pre-correcting for scatter in this manner, allows a linear system to be exposed and solved, of the form

[TABLE]

where $\boldsymbol{n}\in\mathbb{R}^{N}$ is noise, $\boldsymbol{p}\in\mathbb{R}^{N}$ is the linearised projection, calculated by

[TABLE]

where it can be solved by analytic filtered back-projection methods, such as FDK for CBCT [10], or with iterative methods, that approximate the noise model in (1) and incorporate regularisation, such as penalised weighted-least-squares (PWLS). Reconstruction through PWLS involves solving the problem

[TABLE]

where $\boldsymbol{W}\in\mathbb{R}^{N\times N}$ is a diagonal weighting matrix with entries $w_{ii}=(y_{i}-s_{i})^{2}/y_{i}$, $R(\boldsymbol{\mu})$ is a regularisation function to promote desirable structure in $\boldsymbol{\mu}$, and $\lambda$ is usually a scalar constant trade-off between data fit and regularisation.

3 Method

3.1 Off-line Scatter Matching

We propose that the expectation of scatter may be calculated off-line to a high accuracy based upon a prior image, then matched to the measurements during replanning. Conceptually, this is very similar to the notion of SKS/ASKS [6], where the scatter point spread function of a scanner are measured through blocks of material, and combined with convolution. Instead here, the entire global scatter profile is estimated, and simply shifted to fit the current pose of the patient. Our framework for this off-line scheme is illustrated in Figure 1.

For the match, we seek a transformation in the coordinates of detecting elements at each projection angle, for which we adopt the notation

[TABLE]

where $P$ is the number of projections, $\mathcal{I}(\cdot)$ is the 2D linear interpolation of the image $s\{i\}$ corresponding to the $i^{\mathrm{th}}$ projection angle of scatter estimate, $s_{M}\{i\}$ is the matched estimate, and $(u_{M},v_{M})_{i}$ are the transformed 2D coordinates according to

[TABLE]

where $u$ and $v$ are the original vector of coordinates for the detector, and $w$ are parameters we wish to calculate through the matching process.

With reference to Figure 2, each scatter projection image with coordinates $u,v$ is updated based on a rigid translation of an object at the centre of rotation by $[t_{x},t_{y},t_{z}]^{T}$ onto a FDK reconstruction of the measurements. We define a new set of coordinates $[r_{x},r_{y},r_{z}]^{T}$, due to rotation of the source around by the $i^{\mathrm{th}}$ angle $\varphi_{i}$ by

[TABLE]

We note that due to the projection, a translation along $r_{x}$ leads to a change in scaling, and translations in $r_{y},r_{z}$ lead to shifting. If we define the distances $l_{\mathrm{SO}},l_{\mathrm{SD}}$ as the lengths from source to origin and the detector respectively, the transformation to adjust the projection is

[TABLE]

3.2 Reconstruction with Scatter Estimate

Instead using the approximate linearisation of the model in (1), one could use it exactly. This is done for general additive Poisson noise in [16], and we will repeat it here explicitly for the incorporation of a scatter expectation in CBCT. In this case, reconstruction is taken as the maximum likelihood given (1), with a regularisation function to impose desirable properties in the image, as in PWLS (6).

As in [16], we pursue finding the maximum likelihood, by minimising the negative log-likelihood of (1), which is denoted

[TABLE]

It should be noted that for $s_{i}>0$, (9) is nonconvex, so it may not be minimised with the same ease of the PWLS. Nevertheless, it is continuously differentiable with respect to $\boldsymbol{\mu}$, and can therefore be treated with an appropriate first order method. We note that reconstruction is then solution of

[TABLE]

where $\mathcal{C}$ is a set of box constraints on $\boldsymbol{\mu}$ so that $0\leq\mu_{i}\leq\zeta\mbox{ for }i=1,...,N$, where $\zeta$ is the maximum allowable attenuation coefficient.

Although some may consider the difference between our notions of ‘correction’ with PWLS and ‘estimation’ to be trivial, there is a compelling distinction. Whilst in the corrective case, one must carefully design the process to well approximate the model used, in estimation, the expectation of scatter can be used directly, and reconstruction may be considered as the direct inference from the raw measurements. How this translates into practical reconstruction accuracy will be studied in the experimental section.

4 Experimentation

4.1 Data

The data set we are using is derived from repeat CT scans of a neck cancer patient from the Cancer Image Archive [17][18]. With these, we will use the first CT scan as the planning image, then synthesise CBCT measurements on the follow up after 5 months—these are shown in Figures 3(a) and 3(b) respectively. A strong advantage of using this approach is that one has access to a ground truth, against which one can perform valid quantitative assessments.

To generate the CBCT data, we used the Monte-Carlo simulation tool Gate [19] with a 60 keV monoenergetic source on the oracle image, where we did runs with $5\times 10^{10}$ and $1\times 10^{10}$ photons over 160 projection angles to represent two levels of dose.

4.2 Methods Under Test

4.2.1 Scatter Method Implementation

• •

Oracle: this is using the true scatter signal from the measurement synthesis, to represent the ultimate conceivable estimate, and ground-truth for assessment.

• •

None: scatter signal is not estimated at all.

• •

Uniform: calculated using (2) [5] with $\mathrm{SPR}=0.04$ and $t_{\mathrm{air}}=4000,800$ for moderate and low doses.

• •

SKS/fASKS: implemented with same parameters as [6] for ‘full-fan’ acquisition and 20 iterations each.

• •

Diff. filt.: using (3) based upon the registered plan in Figure 3(c).

• •

FDK-MC: a sub-sampled MC estimate based upon preliminary FDK with SKS.

• •

Online-prior-MC: the work-flow in [9], with sub-sampled MC applied to the registered planning image in Figure 3(c).

• •

Offline-prior-MC: using a detailed MC of the unregistered planning image in Figure 3(b), and matching as in Section 3.1 from same registration parameters as other planning methods.

For illustrative purposes, select scatter estimates, along with the ground-truth are shown in Figure 4.

4.2.2 Reconstruction Implementation

The reconstruction methods under test were FDK, PWLS according to (6), and inference (NLL) according to (10). All iterative methods were run for 200 iterations, which was deemed ample for convergence, and all $\lambda$ in the case of PWLS and NLL was set to $2\times 10^{5}$, which was numerically tuned for good performance in both cases.

4.3 Results

Results are summarised in Tables 3 and 2, and selected reconstruction images are shown in Table 1.

The first observation that can be made from the scatter accuracy in high and low dose cases, is that the prior-MC methods are the best. This is interestingly opposed to the FDK-MC, especially since this is based upon the SKS FDK, over which it has a significantly worse scatter estimate and only a slight decrease of analytic reconstruction accuracy. Another significant result is the very poor performance of the diff. filt. method, giving the worst estimates of scatter, which is likely due to large errors from mismatches between the registered plan. Perhaps this would decrease with a non-rigid registration as in [7], though this will inevitably be increasingly difficult and unstable in the lower dose settings.

In general, the relationship between relative errors in SKS/fASKS and the prior MC methods is enlightening. Although fASKS is the best performer in FDK, this does propagate through the rest. Apart from this however, the relative performance of these methods is very similar within the iterative results, all of which become rather close to the oracle scatter reconstruction in NLL. Of these, SKS may be the most appealing due to its fast computation and no reliance to planning registration.

One global trend in both the moderate and low dose results in Tables 3 and 2 is that NLL is more accurate than PWLS on every count. This is unsurprising, since PWLS may be considered an approximation to NLL, but is motivating for pursuing inference methods rather than pre-corrected in the high scatter setting of CBCT. Furthermore, as the low reduces, the difference between these reconstruction models grows.

5 Conclusions

In this study, we have provided valid evidence for differences between various scatter estimation strategies, and how these may best be incorporated into reconstruction. The most conclusive message is that opting for the NLL is more accurate than pre-correcting for scatter and using the PWLS, and in a lower dose setting this difference becomes significant. In terms of scatter estimation, several of the methods aided rather accurate reconstruction: SKS, fASKS, and both on-line/off-line planning MC estimating. Which approach specifically to chose will be dependent on the application, and further evaluations may be valuable for other imaging scenarios and acquisitions.

6 Acknowledgements

The authors would like to thank the Maxwell Advanced Technology Fund, EPSRC DTP studentship funds and ERC project: C-SENSE (ERC-ADG-2015-694888) for supporting this work.