Data consistency-driven scatter kernel optimization for x-ray cone-beam CT

The deconvolution method assumes that the scatter component can be estimated by a spatially invariant convolution of the primary signal with a scatter kernel (Love and Kruger 1987, Seibert and Boone 1988, Ohnesorge et al 1999, Li et al 2008, Ducote and Molloi 2010, Sun and Star-Lack 2010). The scatter kernel, which is based on the physical model, determines the magnitude and distribution of the scatter component in the projection data. The accuracy of the deconvolution method is largely a function of the shape of the scatter kernel. Various scatter kernel models have been proposed for estimating scatter components and for improving the accuracy of deconvolution method. In the present work, we used a symmetric scatter kernel proposed by Ohnesorge (Ohnesorge et al 1999) and modified by Sun (Sun and Star-Lack 2010), for convenience of implementation in the PSO framework. The symmetric scatter kernel or point spread function $\left(\text{PS}{{\text{F}}_{\text{sym}}}\right)$ follows the model consisting of an amplitude factor ${{A}_{f}}$ and a form function ${{h}_{s}}$ as shown below:

$$\begin{eqnarray}&&{{A}_{f}}=~A\cdot {{\left(\frac{{{I}_{p}}(x,y)}{{{I}_{0}}(x,y)}\right)}^{\alpha}}\cdot {{\left(\text{ln}\left(\frac{{{I}_{0}}(x,y)}{{{I}_{\text{p}}}(x,y)}\right)\right)}^{\beta}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{h}_{s}}=~\left[\text{exp}\left(\frac{-{{r}^{2}}}{2\sigma _{1}^{2}}\right)+B\,\text{exp}\left(\frac{-{{r}^{2}}}{2\sigma _{2}^{2}}\right)\right],\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\text{PS}{{\text{F}}_{\text{sym}}}~\left({{I}_{\text{p}}},{{I}_{0}},\ \text{r}\right)={{A}_{f}}\cdot {{h}_{s}}.\end{eqnarray} \tag{ 3 }$$

The spatial position r = xi + yj is defined on the detector, ${{I}_{0}}$ is the unattenuated intensity signal, and ${{I}_{\text{p}}}$ is the attenuated primary signal. The parameters $\alpha ,\beta ,A,B,{{\sigma}_{1}},$ and ${{\sigma}_{2}}$ are the fitting parameters that determine the scatter kernel shape. These parameters can be fit empirically from measurements or simulations of pencil beams directed through the water-equivalent phantoms, which can be cumbersome and suboptimal in various scanning tasks.

The scatter estimate at position $(x,y)$ can be acquired via the 2D convolution operation:

$$\begin{eqnarray}&&{{I}_{s}}(x,y)=\left(I_{p}^{\prime}(x,y){{A}_{f}}(x,y)\right)**{{h}_{s}}(x,y),\end{eqnarray} \tag{ 4 }$$

where $I_{p}^{\prime}(x,y)$ is the primary estimate, ${{A}_{f}}(x,y)$ is the amplitude factor defined in equation (1), and ${{h}_{s}}(x,y)$ is the form factor defined in equation (2).

From equation (4), one can estimate the scatter component of the projection data and can obtain scatter-corrected data by subtracting the estimated scatter from the measured data. In the present work, we used the measured intensity signal I as the initial primary estimate, after which we iteratively estimated the primary intensity. Typically, three iterations are known to be sufficient for the primary intensity recovery (Sun and Star-Lack 2010).

The key purpose of our work was to determine the fitting parameters of the symmetric scatter kernel using the DCC without any further empirical processes or hardware implementation. As the accuracy of scatter deconvolution is affected by the scatter kernel parameters, we optimized them iteratively using the DCC as the criterion.

Various versions of the DCC have been derived and utilized for diverse purposes (Patch 2002, Chen and Leng 2005, Yu et al 2006, Xu et al 2010, Tang et al 2011). The most well-known DCC are the Helgason–Ludwig (HL) conditions, which apply to the Radon transform in parallel projection geometry. The zeroth-order HL condition specifies that the total attenuation in a parallel-beam geometry such as that shown in figure 1 should remain constant from view to view. In other words, if we sum up the measured parallel projection data for each view, the sum should be a constant that is independent of the view-angle.

Zoom In Zoom Out Reset image size

In order to apply this condition to scatter kernel optimization, we used the mid-plane of cone-beam projection data after each deconvolution. We converted the fan-beam projection data to parallel-beam data by use of fan-parallel rebinning (Kak et al 1988), because the zeroth-order HL condition is valid only in parallel-beam geometry.

Both fan- and parallel-beam geometries are shown in figure 2. β represents the fan-beam source angle, $\gamma$ represents the fan-angle, $\theta$ represents the corresponding parallel-beam angle, and $t$ represents the radial parallel coordinate. $s$ represents the virtual fan-beam detector coordinate and $D$ is the distance of the source point from the origin.

Zoom In Zoom Out Reset image size

The relationship between ($\theta ,t$) and ($\beta ,\gamma$) is given by equations (5) and (6), below.

$$\begin{eqnarray}&&\left\{\begin{array}{*{35}{l}} \theta \,\,=\beta +\gamma \\ t\,\,\,=s\,\text{cos}\,\gamma \\ \end{array}\right.\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\left\{\begin{array}{l} \theta =~\beta +\text{arctan}\left(\frac{s}{D}\right) \\ t=~\frac{sD}{\sqrt{\left({{s}^{2}}+{{D}^{2}}\right)}}.~ \\ \end{array}\right.\end{eqnarray} \tag{ 6 }$$

The key step in this work is to find, using the DCC as the criterion, the optimum scatter kernel parameters. While the parameters values are iteratively changed, the optimum that best satisfies the DCC is sought. However, this strategy can be challenging with respect to its computational burden if all of the possible parameter values in a certain range need to be checked in a combinatorial fashion.

For efficient implementation of the optimization process, we incorporated the PSO algorithm (Kennedy and Eberhart 1995, Yuhui and Eberhart 1998, Parsopoulos and Vrahatis 2002, Trelea 2003). The PSO algorithm is an optimization algorithm that finds the best element from sets of available alternatives. PSO simulates the foraging process of certain animals; it is a meta-heuristic, in that it makes few assumptions for problem optimization and can search large spaces of candidate solutions. The PSO algorithm, in comparison with other algorithms, offers the advantages of simplicity in implementation and relatively fast convergence.

In PSO, an individual solution is expressed as a 'particle', and a solution group is expressed as a 'swarm'. There are two conditions that need to be met in the PSO algorithm. Particles, which stand for individuals in a swarm, should share information with each other during iterations. Also, they should move based on that shared information. When better positions of agents are found, the movement of the swarm will be updated and guided accordingly. This process is repeated until a satisfactory solution is reached within a reasonable number of iterations. In our study, we had six parameters in the scatter kernel; therefore, we set the number of swarms as 6, and for each swarm we set the number of particles as 20. A schematic illustration of a particle and a swarm is provided in figure 3, and the movement of a particle in an update step is indicated in figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 3, $i$ represents the index of a parameter, and $k$ represents the index of a single candidate of a parameter. Each particle or candidate of a parameter has an initial position vector and velocity vector. The position vector $S$ represents the value of a parameter, and the velocity vector $V$ depicts the increment of the parameter value. The superscript 'now' and 'next' represent the generations of the iteration. Each particle is evaluated by the cost function, which is the satisfaction level of data consistency as will be explained in the following section. The movements of particles are determined by equations (7) and (8).

$$\begin{eqnarray}&&V_{k}^{\text{next}}=~{{a}_{1}}*V_{k}^{\text{now}}+{{a}_{2}}*\text{rand}*\left(\text{pbes}{{\text{t}}_{k}}-S_{k}^{\text{now}}\right)+{{a}_{3}}*\text{rand}*\left(\text{gbest}-S_{k}^{\text{now}}\right),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&S_{k}^{\text{next}}=S_{k}^{\text{now}}+V_{k}^{\text{next}},\end{eqnarray} \tag{ 8 }$$

where $\text{pbes}{{\text{t}}_{k}}$ represents the position vector that has the best value of the $k\text{th}$ particle, and $\text{gbest}$ represents the position vector that has the best value in the swarm. $\text{rand}$ represents the random number uniformly distributed in the range [0, 1]. ${{a}_{1}}$ is a parameter called the inertia weight parameter, which is employed to accommodate the impact of the previous velocity on the current one. We initialized ${{a}_{1}}$ by 1 and decreased it toward 0 as the number of iterations increased. The inertia weight is known to help global searching in the early iterations as well as fine searching in the last search toward the global optimum. ${{a}_{2}}$ and ${{a}_{3}}$ are positive constants that are used to limit the velocity vectors; in this study, we set both ${{a}_{2}}$ and $~{{a}_{3}}$ to 0.5, which was the suggested value (Parsopoulos and Vrahatis 2002). As meta-heuristic algorithms, including the PSO algorithm, do not guarantee the global optimum in general, we incorporated the inertia weight and the positive constants in the algorithm implementation in this work to gear toward the global optimum. With these enforcing schemes, the PSO is known to have a better chance of finding the global optimum within a reasonable number of iterations, as discussed in Parsopoulos and Vrahatis (2002). In this way, it is likely that the global optimum can be reached empirically for the six parameters of the scatter kernel. We randomized the position and velocity vectors in the initialization step within a reasonable range, referencing Sun's work (Sun and Star-Lack 2010) to guide the ranges of the positions and velocities.

A conceptual workflow of the proposed scatter kernel optimization method is summarized in figure 5.

Zoom In Zoom Out Reset image size

To validate the proposed method, we performed both simulation and experimental studies. In the simulation study, the thoracic part of the XCAT numerical phantom (Segars et al 2010) was used. We acquired scatter-free projection data for 360 views, and added scatter noise by convolving the scatter-free data with a specific scatter kernel. For the scatter calculation, a spatially variant kernel that considers the object thickness was employed, referencing Sun (Sun and Star-Lack 2010) again to determine the parameters of multiple thickness groups. Although scatter kernel information is thus available, for the purposes of scatter correction, we assumed it to be unknown. The recovery of the scatter kernel parameters in the simulation study appeared to be very useful for at-a-glance assessment of optimization algorithm performance. However, since the scatter kernel model used for correction in this work is composed of a multiplicative part to and a convolving part with the primary signals, one-shot deconvolution in the Fourier domain would not work in this approach; instead an iterative deconvolution, as described above, is usually used. Since the iterative deconvolution method applies consecutive convolution and subtraction to the measured data, which contains both primary and scatter signals, the convolution kernel is likely to be different from the kernel that was used for scatter generation in the simulation study. It is particularly so because the kernel used in the scatter generation process employs a thickness-dependent model, whereas the kernel used in the deconvolution process assumes a spatially invariant model in this work. The x-ray tube voltage in the simulation study was set to 140 kVp. A system geometry for the simulation study is indicated in table 1.

In the experimental studies, we used a bench-top CBCT system, the scanning parameters of which are summarized in table 2. The system consists of an x-ray tube (NDI-451, Varian Co., Salt Lake City, Utah, USA) and a flat panel detector (1642AP, Perkin Elmer, Santa Clara, California, USA). An aluminum wedge filter was used in order to make a uniform beam quality. Neither a bowtie nor an anti-scatter grid was employed in the experiments. A PH-3 angiographic CT head phantom ACS (Kyoto Kagaku Co., Ltd, Kyoto, Japan) and the pelvic part of the RT-humanoid phantom (the Rando phantom, Humanoid Systems, Carson, California, USA) were utilized in the experiments. Using these phantoms, we acquired not only cone-beam projection data but also fan-beam data by collimating the x-ray beam to secure, by inherent scatter suppression, almost scatter-free reference data. The x-ray tube voltage was 140 kVp and the number of projection views was 360, as in the simulation study.

In the scatter-free projection data, the DCC was well satisfied as shown in figure 6(a), though there was a slight variation, possibly due to the beam-hardening and the partial volume effects. Contrastingly, in the scatter-contaminated case, there was a large variation in the sum of line integrals from view to view, and in fact, lowered values can be seen in figure 6(b). This variation and reduction in the sum of line integrals was due mostly to the added scatter. Figure 6(c) plots the sum of line integrals of the scatter-corrected data by use of the proposed method; as is apparent, the DCC was fairly well recovered. As a quantitative metric that can be used for a cost function, we introduced the inconsistency level defined in equation (9). The inconsistency level represents the percentage of standard deviation with respect to the average. A low inconsistency level implies that the DCC is highly satisfied.

$$\begin{eqnarray}&&\text{Inconsistency}\ \text{level}=~\frac{\text{std(value)}}{\text{average(value)}}\times 100~(\%).\end{eqnarray} \tag{ 9 }$$

Zoom In Zoom Out Reset image size

Table 3 shows the inconsistency levels for the respective cases. The difference in the inconsistency level of the scatter-contaminated data with respect to the scatter-free data was calculated as $\left(2.781-0.529\right)/0.529\times 100=426~\%$, and similarly, the corresponding value of the scatter-corrected data was $\left(0.637-0.529\right)/0.529\times 100=20~\%$. Table 4 shows the selected parameters of the scatter kernel for the optimum DCC condition in the simulation study.

Table 4. Optimally selected parameters of scatter kernel in simulation study.

$A~\,\,\left(\text{c}{{\text{m}}^{-2}}\right)$	$\alpha$	$\beta$	$B$	${{\sigma}_{1}}~\,\,\left(\text{cm}\right)$	${{\sigma}_{2}}~\,\,\left(\text{cm}\right)$
$6.62\times {{10}^{-4}}$	$-0.14$	$1.33$	$2.38$	$22.20$	$5.36$

Figure 7 shows the converging behavior of the inconsistency level as the number of iterations increases by use of the PSO algorithm. We used $0.5$ for $\epsilon$ to demonstrate the converging nature, but a higher value can be used in practical applications so as to stop at earlier iterations.

Zoom In Zoom Out Reset image size

The reconstructed images on the mid-plane and off-mid-plane from the scatter-free, scatter-contaminated, and scatter-corrected projection data and their line profiles along the solid lines are shown in figures 8 and 9, respectively. Both on the mid-plane and off-mid-plane, image artifacts caused by the scatter component are clearly seen in figure 8(b) and in the dotted line profile in figure 9; such artifacts were substantially mitigated by use of the proposed method.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For a quantitative assessment of image quality, we used structure similarity (SSIM), which is a metric that evaluates the similarity between two images. For comparison of the scatter effects between cases, contrast also was calculated. Region of interests (ROIs) that we used for the assessments are indicated in figure 8. The results of the quantitative assessment of each case are summarized in tables 5 and 6. On the mid-plane images, the optimally selected scatter kernel improved the contrast up to 99.4% and the SSIM up to 99.8% from 75.1% and 70.4% in the scatter-contaminated images, respectively. Similarly on the off-mid-plane images, the scatter correction improved the contrast up to 99.5% and the SSIM up to 99.8% from 82.3% and 66.0% in the scatter-contaminated images.

For direct evaluation of scatter, the scatter-to-primary ratio (SPR) in the thoracic region of the projection data also was calculated. The results of this direct assessment are summarized in table 7; as indicated, the scatter correction results reduced the SPR value in the projection data from 0.282 to 0.011. These quantitative assessment results show that the image accuracy had been much improved by application of the proposed scatter correction method across the entire field of view.

Figure 10 shows the converging behavior of the inconsistency level in the experimental study. As in the simulation, we used $0.5$ for $\epsilon$, and a similar convergence was observed. In this study, for the purposes of comparison, we also acquired scatter-corrected data with a suboptimal kernel at the 15th iteration.

Zoom In Zoom Out Reset image size

In the fan-beam projection data, the DCC was satisfactorily met, as shown in figure 11(a), though there existed a weak variation of data sum, which was thought to have been due to beam-hardening and residual fan-beam-mode scatter. As is shown in figure 11(b), the DCC was much degraded due to the cone-beam scatter contribution. Figures 11(c) and (d) show the data sum profiles of the scatter-corrected data by use of the suboptimal kernel and optimized kernel, respectively. The results clearly demonstrate that the DCC can provide an effective criterion for scatter correction adjustment in CBCT. The inconsistency level measurements, meanwhile, are summarized in table 8. In figures 11(a) and (b) and table 8, it is worth noting that the beam-hardening in CBCT was, compared with the scatter, relatively minor in its effects on the data consistency. The difference in inconsistency level of the scatter-contaminated data with respect to the scatter-free data was calculated to be $\left(2.407-0.414\right)/0.414\times 100=481~\%$, and similarly, the corresponding value of the optimally scatter-corrected data was 30%. The selected parameters of the scatter kernel for the optimum DCC condition in the experimental head phantom study are summarized in table 9.

Zoom In Zoom Out Reset image size

Table 9. Optimally selected parameters of scatter kernel in experimental head phantom study.

$A~\,\,\left(\text{c}{{\text{m}}^{-2}}\right)$	$\alpha$	$\beta$	$B$	${{\sigma}_{1}}~\,\,\left(\text{cm}\right)$	${{\sigma}_{2}}~\,\,\left(\text{cm}\right)$
$4.35\times {{10}^{-4}}$	$-0.04$	$0.33$	$1.47$	$11.11$	$2.20$

In figure 12, the reconstructed images on the mid-plane and off-mid-plane from scatter-free, scatter-contaminated, and scatter-corrected projection data are shown, respectively. The line profiles of the reconstructed images are also shown in figure 13; the cupping artifact from the scatter component was substantially mitigated. The sagittal and coronal views without and with scatter correction are shown in figure 14. It is apparent that the scatter artifacts were successfully removed by the proposed method.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 12 indicates the ROIs used in our quantitative assessments of image quality, while tables 10 and 11 summarize the results for SSIM and contrast, respectively. On the mid-plane images, the optimally selected scatter kernel improved the SSIM up to 92.3% and the contrast up to 91.8% from 60.3% and 56.1%, respectively, in the scatter-contaminated images. Meanwhile, on the off-mid-plane images, the optimally selected scatter kernel improved the SSIM up to 92.6% and the contrast up to 96.9% from 51.5% and 51.8%, respectively. Additionally, table 12 shows that the optimally selected scatter kernel decreased the SPR in the projection data from 0.403 to 0.099.

The overall results confirm that the image accuracy was much improved after application of the proposed scatter correction method. Particularly, the results of suboptimal scatter correction compared with those of optimum correction highlight the importance of the optimization of scatter kernel parameters.

For further validation of the proposed method, we also performed an experimental study using the pelvic part of the Rando phantom which is less isotropic in its shape compared to the head phantom. Figure 15 shows the converging behavior of the inconsistency level. Unlike previous studies, we used 2 for $\epsilon$, as the inconsistency level of the pelvis phantom turned out to be higher than that of the head phantom. This discrepancy was thought to be due predominantly to the asymmetric shape of the pelvis phantom, which will be discussed to detail in the discussion. We also acquired the scatter correction data with a suboptimal kernel for comparison whose parameters were determined at the 5th iteration.

Zoom In Zoom Out Reset image size

In the fan-beam projection data, the DCC is relatively poorly met, as shown in figure 16(a). It was conjectured that the asymmetric shape of the phantom leads to noticeably varying degrees of beam-hardening contribution to the data consistency. Figure 16(b), however, shows that the DCC is much more challenged by the cone-beam scatter contribution. Note that the sum of the projection data drops by a substantially larger amount than the fluctuation seen in the fan-beam case. Figures 16(c) and (d) reveal that as the kernel reaches the optimum, the sum profiles of the scatter-corrected data tend to recover that of the fan-beam case. The results clearly demonstrated that even in the pelvic scans, the DCC can provide an effective criterion for scatter correction in CBCT. The inconsistency level measurements are summarized in table 13. The difference of the inconsistency level of the scatter-contaminated data from that of the scatter-free data was calculated as $\left(5.080-1.913\right)/1.913\times 100=165~\%$, and similarly, the corresponding value of the optimally scatter-corrected data was 28%. Table 14 shows the selected scatter kernel parameters for the optimum DCC condition in the experimental pelvis phantom study.

Zoom In Zoom Out Reset image size

Table 14. Optimally selected parameters of scatter kernel in experimental pelvis phantom study.

$A~\,\,\left(\text{c}{{\text{m}}^{-2}}\right)$	$\alpha$	$\beta$	$B$	${{\sigma}_{1}}~\,\,\left(\text{cm}\right)$	${{\sigma}_{2}}~\,\,\left(\text{cm}\right)$
$8.50\times {{10}^{-4}}$	$-0.04$	$1.68$	$2.09$	$16.72$	$3.61$

In figures 17 and 18, the reconstructed images on the mid-plane and on the off-mid-plane from scatter-free, scatter-contaminated, and scatter-corrected projection data and their line profiles are shown, respectively. The sagittal and coronal views without and with scatter correction are shown in figure 19. They serve to demonstrate that the scatter artifacts were successfully mitigated by means of the proposed method.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 17 marks the ROIs that we utilized for quantitative assessment, while tables 15 and 16 summarize the results for SSIM and contrast, respectively. On the mid-plane images, the SSIM and contrast were improved up to 91.7% and 83.8%, respectively, by use of the optimally selected scatter kernel. On the off-mid-plane images meanwhile, the SSIM and contrast were improved up to 89.9% and 85.0%, respectively. According to the SPR measurement, use of the optimally selected scatter kernel reduced the value of the SPR from 0.508 to 0.115 (table 17). Indeed, all of the assessments showed that the image accuracy had been much improved by use of the proposed scatter correction method.