Image reconstruction for interrupted-beam x-ray CT on diagnostic clinical scanners

Our system model relies on that of previous work (Thibault et al 2007). The p th data point, b_p can be modeled as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:origmodel} b_p &= f_p + \epsilon_p, \nonumber \\ f_p &\sim {\rm Poisson}(I_p e^{-\bar{y}_p}), \nonumber \\ \epsilon_p &\sim \mathcal{N}(0,\sigma^2), \nonumber \\ \bar{y}_p &= {\bf a}_p^T{\bf x},\nonumber \end{align} \tag{ 1 }$$

where ${\bf x} \in \mathbb{R}^N$ is the object being imaged with N voxels and I_p is the source intensity (Thibault et al 2007). The line integral projections, $\bar{y}_p$, are attenuated by the object, ${\bf x}$, as described by the ray path in ${\bf a}_p$. The Poisson variable f _p with mean and variance $\bar{y}_p$ is called the quantum noise, while the Gaussian variable $\newcommand{\e}{{\rm e}} \epsilon_p$ with variance $\sigma^2$ is called the electronic noise.

For simplification, we work with the postlog data, y _p, instead of working with (1) directly. We use the following linear model (Thibault et al 2007):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:realmodel} {\bf y} = {\bf Ax} + {\bf n}, \nonumber \end{align} \tag{ 2 }$$

where ${\bf y} = \left [y_1, ..., y_P \right ]^T$ is the preprocessed data and ${\bf A} = \left [ {\bf a}_1, ..., {\bf a}_P \right ]^T$ is the system matrix of line integrals. Section 2.2 describes the construction of ${\bf A}$ that takes into account partial obstruction of the x-ray source. ${\bf n} \sim \mathcal{N}(0, {\bf W}^{-1})$ is Gaussian noise designed to approximate the combined effects of both noise sources and the preprocessing steps. For clarity, we note that $\bar{y}_p$ is the ideal, noiseless version of y _p, whereas y _p is the estimate of $\bar{y}_p$ that arises from log transform preprocessing. The covariance matrix, ${\bf W}^{-1}$, is assumed to be diagonal, with its inverse being

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf W} = \left[\begin{array}{@{}ccc@{}} w_1 & & \nonumber \\ & \ddots & \nonumber \\ & & w_P \end{array}\right], \nonumber \end{align} \tag{ 3 }$$

where w_p is the variance of the p th data point. Based on previous results (Thibault et al 2007), we use $w_p=b_p$, which leads to high-intensity pre-log data being more important in the reconstruction. Calculating y _p from b_p requires a logarithm and some spectral corrections that we describe in section 3.2.2.

Given the signal model in (2), one could obtain the penalized weighted least squares (PWLS) estimate of ${\bf x}$ by solving the following optimization problem:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\defeq}{:=} \newcommand{\mattnorm}[1]{\left \| {#1} \right \|} \newcommand{\argmin}[1]{\underset{#1}{{\rm argmin}}} \displaystyle \label{eq:costfun} \hat{{\bf x}} = \argmin{{\bf x} \succeq {\bf 0}} \left \{\Psi({\bf x}) \defeq \frac{1}{2} \mattnorm{{\bf y}-{\bf Ax}}_{{\bf W}}^2 + R({\bf x}) \right \}, \nonumber \end{align} \tag{ 4 }$$

where the above estimate also enforces a non-negativity constraint on the attenuation coefficients in ${\bf x}$. To compensate for aliasing due to MSC-based undersampling, we incorporate $R({\bf x})$, a regularization function that promotes sparsity. We define $R({\bf x})$ to have the following general form:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R({\bf x}) = \sum_{j=1}^J \psi_j \left ( \left [{\bf Cx} \right ]_j \right), \nonumber \end{align} \tag{ 5 }$$

where ${\bf C} \in \mathbb{R}^{J \times N}$ is a finite difference operator and $\left [ {\bf Cx} \right ]_j$ is the j th output of ${\bf Cx}$. The function $\psi_j(t)$ is a potential function. The classic compressed sensing choice for $\psi_j(t)$ would be the absolute value function, thus yielding the $\newcommand{\e}{{\rm e}} \ell_1$ regularizer. However, in our experiments we found that this could lead to salt-and-pepper noise artifacts as demonstrated in previous results (Thibault et al 2007). Instead, we use the hyperbola function (Long et al 2010),

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:psieq} \psi_j(t) = \beta_j \delta^2 \left [ \sqrt{1 + (t/\delta)^2} - 1 \right ]. \nonumber \end{align} \tag{ 6 }$$

With small values of $\delta$, (6) approximates the classic Total Variation function used for compressed sensing. $\beta_j$ is a space-varying regularization weight. We allow the regularization parameter to vary since the data fidelity weighting can induce resolution variation on the reconstructed image (Fessler and Rogers 1996). In particular, low weighting of projections through the center of the image volume can lead to over-regularization in this area. We compensate for this in our selection of the $\beta_j$ weights. First, we calculate

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \kappa_n = \sqrt{\left. \left ( \sum_{p=1}^P a_{pn}w_p \right) \middle/ \left ( \sum_{p=1}^P a_{pn} \right)\right.}, \nonumber \end{align} \tag{ 7 }$$

where a_pn is the p th, nth value of ${\bf A}$. We then apply

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:varbeta} \beta_j = \beta \prod_{\forall n, c_{jn} \neq 0} {\rm max} \left (\kappa_n, 0.01\kappa_{{\rm max}} \right), \nonumber \end{align} \tag{ 8 }$$

where c_jn is the kth, nth value of ${\bf C}$ and $\beta$ is a baseline value that governs the regularization level in the entire image volume. This choice of $\beta_j$ values promotes uniform resolution in the reconstructed image (Fessler and Rogers 1996, Kim et al 2015).

We construct the data operator, ${\bf A}$, based on the classic Siddon line integral method (Siddon 1985). We choose to do this since the Siddon method introduces minimal parameterization to the system geometry, allowing adaptation to correct for the effects of partial source obstruction, a distortion the MSC introduces to the standard system geometry. Partial source obstruction results in each row of the detector only seeing a portion of the source. This shifts the effective location of the ideal point source, resulting in row-variant cone beam magnification effects. To compensate for partial source obstruction effects, we shift the location of the source for each row of the detector. We determine the size of the shift from offline simulation experiments on the scanner geometry for each MSC position and MSC design. Then, in our Siddon ray-tracing procedure (Siddon 1985), we apply this shift for each row of the detector and MSC position.

Figure 2 shows a schematic of the partial source obstruction effect. To estimate what portion of the source was observable, we approximated the source as a summation of Gaussian basis functions based on experimental data. Then, we calculated the centroid of the observable source and used this information to perturb the point source location for each row of the detector in the projector.

There are several approaches for minimizing the cost function in (4): our chosen method was to use the ordered subsets/momentum algorithm (OS-momentum) (Kim et al 2015) due to its speed and ease of implementation. The OS-momentum method develops an optimization transfer algorithm for (4) and accelerates the optimization transfer algorithm using ordered subsets and Nesterov momentum (Nesterov 1983). We review this procedure for (4).

The first step is the creation of the surrogate, which can be done by upper-bounding the Hessian of (4), ${\bf H}$ with a diagonal matrix, ${\bf D}$. The Hessian can be decomposed into two parts given by the data-fitting term and the regularizer, i.e. ${\bf H}={\bf H}_f + {\bf H}_R$. Then, the following is a diagonal upper bound:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf D}_f + {\bf D}_R = {\bf D} \succeq {\bf H} = {\bf H}_f + {\bf H}_R, \nonumber \end{align} \tag{ 9 }$$

provided that ${\bf D}_f \succeq {\bf H}_f$ and ${\bf D}_R \succeq {\bf H}_R$ where ${\bf D} \succeq {\bf H}$ implies that ${\bf D} - {\bf H}$ is positive semidefinite. Thus, we can consider the data-fitting and regularizer terms separately in two stages. A diagonal upper bound for the data-fitting part of the Hessian, ${\bf A}^T{\bf W}{\bf A}$, can be found using the method of De Pierro (1995):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf D}_f = {\rm diag} \left \{\left | {\bf A}^T \right | {\bf W} \left | {\bf A} \right | {\bf 1} \right \} \succeq {\bf A}^T {\bf W}{\bf A}, \nonumber \end{align} \tag{ 10 }$$

where $|\cdot|$ is the element-wise absolute value operation. This bound is commonly used for quadratic surrogates (SQS) methods (Erdoğan and Fessler 1999). The method of De Pierro also gives a bound for the regularizer (De Pierro 1995),

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf D}_R = {\rm diag} \left \{\left | {\bf C}^T \right | {\rm diag} \left \{\ddot{\psi}_j(0) \right \} \left | {\bf C} \right | {\bf 1} \right \}. \nonumber \end{align} \tag{ 11 }$$

Combining these into ${\bf D}={\bf D}_f+{\bf D}_R$ gives the following separable quadratic surrogate (Erdoğan and Fessler 1999) at iteration k for the cost function in (4):

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\mattnorm}[1]{\left \| {#1} \right \|} \displaystyle \label{eq:surrogate} \phi_k({\bf x}) = \frac{1}{2} \mattnorm{{\bf x} - \left ( {\bf k}^{(k)} - {\bf D}^{-1} \nabla \Psi \left ({\bf x}^{(k)} \right) \right) }_{\bf D}^2, \nonumber \end{align} \tag{ 12 }$$

where $\nabla \Psi \left ({\bf x}^{(k)} \right)$ is the gradient computed at the point, ${\bf x}^{(k)}$. Iteratively applying the update in k will descend the cost function in (4).

Such an algorithm is typically slow, so we add ordered subsets and momentum accelerations (Kim et al 2015). The ordered subsets acceleration assumes that the overall cost function can be decomposed into sums of functions that have the same shape. If this is satisfied, then we apply a similar update to that in (12), but instead using the ordered subset cost function:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\defeq}{:=} \newcommand{\mattnorm}[1]{\left \| {#1} \right \|} \displaystyle \Psi_m({\bf x}) \defeq \frac{1}{2} \mattnorm{{\bf y}_m - {\bf A}_m {\bf x}}_{{\bf W}_m}^2 + \frac{1}{M} R({\bf x}), \nonumber \end{align} \tag{ 13 }$$

for $m=1, ..., M$. Further, we assume that $\nabla \Psi \left ({\bf x}^{(k)} \right) \approx M \nabla \Psi_1 \left ({\bf x}^{(k)} \right) \approx \hdots \approx M \nabla \Psi_M \left ({\bf x}^{(k)} \right)$ and replace $\nabla \Psi \left ({\bf x}^{(k)} \right)$ in (12) with $\nabla \Psi_m \left ({\bf x}^{(k)} \right)$. This can also be combined with Nesterov momentum by introducing auxiliary variables ${\bf v}$ and ${\bf z}$ (Kim et al 2015), giving the overall OS-Momentum Algorithm. All subsets are updated in the FFT bit-reversal order to ensure incoherence between consecutive updates.

Our simulation experiments used data from the Low Dose CT Grand Challenge (McCollough et al 2017). The data consisted of a central $736 \times 64 \times 4006$ (detector channels $\times$ detector rows $\times$ projections) section of the original helical sinogram of a liver examination from patient 67. This section of data was cropped from the original data set from z  =  125 mm to z  =  165 mm, yielding 4,006 projections. From this, we reconstructed a $768 \times 768 \times 32$ image volume with 0.66 mm $\times$ 0.66 mm $\times$ 3 mm resolution. The images displayed in subsequent figures are from slice 15 of the image volume.

Although the tungsten MSC has a very high attenuation factor and completely blocks the beam, the size of the focal spot leads to a non-binary sampling pattern corrupted by penumbra effects (Chen et al 2017, 2019). The signal in the penumbra region may experience different physical effects than the main signal region—this can be caused by a spectral shift due to the source location on the anode or a partially attenuated path through a corner of the MSC. For these reasons, we decided to exclude elements from the penumbra region that had more than an 80% attenuation due to the MSC. We corrupted the data using a noise model that incorporated penumbra attenuation effects (Chen et al 2017) based on the W4S16 and W2S16 sampling patterns (4-fold and 8-fold reduction factors). Then, we reconstructed by minimizing (4) using the OS-Momentum algorithm (Kim et al 2015). We reconstructed images across a large array of values for the regularization parameter. We used a $\delta=5$ HU based on a previous study (Long et al 2009). These results were compared to quarter and eighth dose tube-current reduction simulations based on the Siemens ReconCT software (Yu et al 2012), also reconstructed by minimizing (4) with the OS-Momentum algorithm. A reference standard was reconstructed with uncorrupted data. Lastly, our experiments used a sinogram oversampling factor of 3 to reduce ringing from the Siddon ray-tracing procedure.

We did not simulate the effects of partial source obstruction in our simulation experiments. Our main task was to isolate the sampling and noise properties of the MSC to ascertain feasibility of the approach on in vivo data. However, these effects are fully-realized in the real-world experiments.

3.2.1. Data Collection

Our real-world prototype consisted of a Siemens SOMATOM Force scanner equipped with a W4S16 collimator in place of the standard adaptive collimator. Figure 3 shows the W4S16 collimator secured in the adaptive collimator jaws of the Siemens scanner.

Zoom In Zoom Out Reset image size

In our prototype experiments, we used the American College of Radiology CT (ACR) phantom. The phantom has inserts at known locations that can be used to determine the resolution of the system. Four data sets were acquired with the prototype scanner: a no-MSC, 1/4 tube current air scan, a no-MSC, 1/4 tube current ACR scan, an MSC air scan, and an MSC ACR scan. The air scans were acquired for log-domain preprocessing (described in section 3.2.2). From these data, we retrospectively simulated a scan in which the MSC has a piecewise linear motion from the first to the sixteenth position over one rotation. This generated a $920 \times 96 \times 2100$ data set. Of this data set, elements that had greater than 80% attenuation due to the MSC were excluded. From this, we reconstructed a $768 \times 768 \times 32$ image volume at 0.66 mm $\times 0.66$ mm $\times 3$ mm resolution. The reconstruction process used the custom preprocessing pipeline described in section 3.2.2 and the partial source obstruction model shown in section 2.2.

3.2.2. Prototype preprocessing

The prototype experiments used a custom preprocessing pipeline to convert the prelog data points b_p to postlog data points y _p. Preprocessing pipelines consist of several corrections for non-ideal system factors. These factors include alterations to the spectral distribution of the x-ray source due to beam hardening, negative values that can arise in b_p due to the electronic noise, geometric distortions, and anatomy-specific factors. Corrections for these factors can vary between vendors, and the details are often proprietary. Our goal is to build upon past literature on preprocessing while making modifications necessary for SparseCT. These include air calibration, log transformation, and beam-hardening correction. We apply these modifications to the data in our prototype experiments.

We first correct the negative values due to the electronic noise in b_p. Current systems have signal-dependent filters (SDFs) for handling small values that occur behind highly-attenuating structures such as bone. However, the MSC itself is an abnormal high-attenuating structure, so we do not expect these filters to work in our setting. Instead, we apply clipping:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle s_p = {\rm max}(b_p, b_{p,{\rm min}}), \nonumber \end{align} \tag{ 14 }$$

where $b_{p,{\rm min}}=0.000\,01I_p$. This limits the minimum value to be five orders of magnitude below the incident value. Clipping can create some bias in low-signal regions, but we design the reconstruction weights in (4) such that low-signal regions have relatively small effects on the final reconstruction.

We next apply the logarithm. This requires knowledge of the fluence, I_p, which is measured with an air scan. We collected air scans for each MSC position with our prototype. With these scans, we apply the logarithm as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u_p = - {\log}\left (s_p/I_p \right). \nonumber \end{align} \tag{ 15 }$$

The final step is beam-hardening correction. For this, we relied on a propriety 2-compartment model from Siemens (related, Raupach et al (2001)). Since we apply the method in the post-log domain after the air scan calibration, we assume that the beam-hardening method is performing in normal operating conditions and that no further modifications are necessary due to the MSC. This gives the final postlog data for reconstruction:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y_p = {\rm BHC}(u_p), \nonumber \end{align} \tag{ 16 }$$

where ${\rm BHC}(\cdot)$ is the beam-hardening correction function. We did not perform scatter correction due to the anti-scatter grid on the detector.

In our simulation experiments we compared to the full-dose reference scan shown in figure 4(a). We reconstructed this image with mild regularization to reduce the influence of the noise present in the original data. After calculating the reference image, we simulated SparseCT acquisitions with the W4S16 and the W2S16 collimator design. For comparison, we also simulated 1/4th dose and 1/8th dose tube-current acquisitions and reconstructed them with the same algorithm as the full-dose scan. Then, we calculated the normalized root-mean-square difference (NRMSD) and structural similarity measure (SSIM) across six slices between the low-dose image volumes and the reference image volume. The SSIM calculations used 400 HU for the dynamic range since this was the viewing window size. All displayed images in figures 4 and 5 were reconstructed with the NRMSD-optimal $\beta$ for that method. The results of these reconstructions are shown in figures 4 and 5. Both tube-current reduction and SparseCT acquisitions show image degradation as the dose is reduced in the simulated images in figures 4 and 5. At high dose reduction levels, the SparseCT acquisitions showed less image degradation, as illustrated in figure 5. To more fully characterize the methods, we also plotted SSIM and NRMSD versus the ground truth image as a function of regularization parameter. These are shown in figure 6. The SparseCT methods perform better for the metrics at lower regularization levels. Among our experiments the minimum W4S16 NRMSD was 2.56% at $\beta=1270$ with an SSIM of 0.8752. The minimum 1/4 tube current NRSMD was 2.69% at $\beta=3002$ with an SSIM of 0.8810. The minimum W2S16 NRMSD was 3.19% at $\beta=1270$ with an SSIM of 0.8336. The minimum 1/8 tube current MSE was 4.19% at $\beta=3999$ with an SSIM was 0.8015.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2.1. Partial source obstruction

Figure 7 shows images illustrating the effects of partial source obstruction on one view of the sinogram. The original, open-collimator scan is shown in figure 7(a). Figure 7(b) shows how the operating table becomes slanted in the presence of the MSC—this is due to variations in the cone beam magnification effect caused by changes in the effective source location. Figure 7(c) shows one view from ${\bf Ax}_{{\rm no~msc}}$, where ${\bf x}_{{\rm no~msc}}$ is the reconstructed image volume from no-collimator sinogram (i.e. that shown in figure 7(a)). Figure 7(c) shows the effect of the simulated partial source obstruction model. The model reproduces the slant of the scanner table. The model does not smoothly interpolate into regions where the MSC incurs full obstruction, but the effects of this region are negligible due to sinogram-domain masking and the weights matrix, ${\bf W}$.

Zoom In Zoom Out Reset image size

4.2.2. Resolution

Figure 8 compares results with a quarter-dose open collimator versus the MSC with linear motion by zooming in on the 7 line pairs per cm (lp cm⁻¹) insert of the ACR phantom. For these reconstructions, a $\beta$ of 128 was used for all methods. The open collimator results were obtained with the same dose reduction achieved for the SparseCT method and with the same reconstruction algorithm. Both methods resolve the 7 lp cm⁻¹ insert. Although both methods were able to resolve the 7 lp cm⁻¹ insert in figure 8, we found inserts finer than 7 lp cm⁻¹ were more difficult.

Zoom In Zoom Out Reset image size

Figures 4 and 5 suggest that SparseCT and tube current have different error characteristics. Although SSIM and NRMSD are similar for both methods, the errors for the tube current reduction method are heavily concentrated on the interior of the abdomen, whereas the errors for SparseCT are more diffuse. This effect becomes greater at higher dose reduction factors. Model-based reconstruction methods for low-dose CT often rely on detailed system models and priors to compensate for this increased noise. Figure 6 suggests that with these models, the extent of regularization necessary to maintain high image quality is reduced with an interrupted-beam design. Nonetheless, the trends suggested by figures 4–6 should be considered a proof-of-concept of the potential of an interrupted-beam design. Further validation is necessary in future work.

Our reconstruction pipeline shares a number of aspects with other low-dose reconstruction methods. Although data fidelity weighting and regularization for uniform resolution were originally developed for tube-current reduction schemes (Fessler and Rogers 1996, Thibault et al 2007), we found they readily extend to the SparseCT setting to consider penumbra effects. We also found that much of the standard preprocessing pipeline can be preserved after air calibration.

Still, SparseCT introduces some new considerations that are not present in other sparse sampling schemes such as view-based undersampling. In particular, we discovered that partial source obstruction can affect the data. Our correction for this required simulations of the source distribution and an alteration to the projection operator. Many implementations of modern projection operators (such as separable footprints (Long et al 2010)) are highly parameterized, but the jittering necessary for partial source obstruction precludes use of detailed parameterizations. This was the primary motivation for our use of the very general Siddon method (Siddon 1985). In the future, we will examine whether more parameterized projection operators can be adapted for SparseCT.

Here, we simulated a linear motion of the MSC by retrospectively combining scans at different MSC positions. These results readily extend to a prospective linear motion with a few modifications under development. First, the air calibration scan must be registered to each location of the MSC in the sinogram. This requires a process for detecting the location of the MSC, then interpolating between the two closest corresponding air scan locations in order to perform the appropriate air calibration. Another modification will register the partial source obstruction effects to the detected location of the MSC by first fitting a polynomial for partial source obstruction corrections based on simulations, then applying this polynomial to the appropriate position determined from the MSC registration procedure. After these modifications are applied, the rest of the pipeline described in this paper can be used for reconstruction.

Another possibility for extension is that of dual-energy CT. The MSC naturally segments the sinogram into high and low-energy areas. In our case we entirely blocked some segments of the sinogram in order to achieve dose reduction, but an alternative approach would be to only partially attenuate these areas to achieve dual-energy CT imaging, an application that has been explored for CBCT systems (Lee et al 2017). The partial source obstruction models developed in this paper would continue to be relevant to reconstruction in this setting. Joint reconstruction techniques (Knoll et al 2014, Rigie and La Rivière 2015, Lee et al 2017) could be used to mitigate the image degradation from the low-energy filter by combining information from both segments of the sinogram.